University of Florida/Eml4507/Fead-s13-team4-R2
Problem R2.1 Model 2-D Truss With 2 Inclined Elements and Verify Assembly
editGiven 2-D Truss Model and Constants
editConsider the 2D truss system
Use the following constants to very calculations
Find:
edit1) Explain the Assembly Process of the Above 2D Truss System
edit2) Verify With CALFEM
editSolution: Assembly Process for 2d Truss System
editConstruct local stiffness matrices
editThe stiffness matrices will be constructed for individual elements 1 and 2. The rows and columns will be numbered to represent the force and degree of freedom. For element 1
For element 2
The variable K can be determined by
If the variables are plugged in, we get the following matrices:
For element 1
For element 2
It is clear that if we wish to combine the matrices, it will place coordinate 3,3 of element 2 on the location of 3,3 of element 1 due to them referencing the same point in the same direction. When combining the matrices, addition will be used since the matrices represent directional forces and forces always combine in the same direction by addition.
Combining local stiffness matrices into a global matrix
editThe local matrices from element 1 and 2 will be combined as described above to give
Solution: Verify With CALFEM
editMatLab Code
editEDU>> Edof = [1 1 2 3 4
2 3 4 5 6];
EDU>> K = zeros(6);
EDU>> f = zeros(6,1);
EDU>> f(4) = -7;
EDU>> E1 = 3; E2 = 5;
EDU>> L1 = 4;L2 = 2;
EDU>> A1 = 1;A2 = 2;
EDU>> ep1 = [E1 A1]
ep1 =
3 1
EDU>> ep2 = [E2 A2]
ep2 =
5 2
EDU>> ex1 = [-3.464 0]; ex2 = [0 1.414];
EDU>> ey1 = [-2 0]; ey2 = [0 -1.414];
EDU>> Ke1 = bar2e(ex1,ey1,ep1)
Ke1 =
0.5625 0.3248 -0.5625 -0.3248 0.3248 0.1875 -0.3248 -0.1875 -0.5625 -0.3248 0.5625 0.3248 -0.3248 -0.1875 0.3248 0.1875
EDU>> Ke2 = bar2e(ex2,ey2,ep2)
Ke2 =
2.5004 -2.5004 -2.5004 2.5004 -2.5004 2.5004 2.5004 -2.5004 -2.5004 2.5004 2.5004 -2.5004 2.5004 -2.5004 -2.5004 2.5004
EDU>> K = assem(Edof(1,:),K,Ke1);
EDU>> K = assem(Edof(2,:),K,Ke2)
K =
0.5625 0.3248 -0.5625 -0.3248 0 0 0.3248 0.1875 -0.3248 -0.1875 0 0 -0.5625 -0.3248 3.0629 -2.1756 -2.5004 2.5004 -0.3248 -0.1875 -2.1756 2.6879 2.5004 -2.5004 0 0 -2.5004 2.5004 2.5004 -2.5004 0 0 2.5004 -2.5004 -2.5004 2.5004
EDU>> bc = [1 0;2 0;5 0;6 0];
EDU>> [a,r]=solveq(K,f,bc)
a =
0 0 -4.3519 -6.1268 0 0
r =
4.4378 2.5622 0 -0.0000 -4.4378 4.4378
Conclusion
editBy comparing the value from the hand calculation and the value from the CALFEM program, we can see that the two matrices are identical. This verifies the assembly process of using addition when combining local matrices by addition.
Problem R2.2
editDisccription:
editWe are asked to consider the L2-ODE-CC system referenced in (1) p.53-2 to have the following complex conjugate roots:
We are asked to find the homogeneous solution in standard form and to find the damping ratio,
We are then to consider the initial conditions
We are to plot this solution and then use that plot to find the log decrement and compare it to the log decrement that would be obtained using (6) p.53-5b.
Finally, we're to find the quality factor, Q, and the loss factor,
Solution:
editFrom R1.5, we have already obtained that for the above conditions,
with coefficients to the characteristic equation of
and
From (1) and (2) p.53-5, we have that
Knowing both a and b we get
We then have from equation (1) p.53-5d
Using the above graph of the homogeneous solution, accurate estimates were able to be made as to the various values of y at the peaks
Using equation (3) p.53-5c we are able to find
Conclusion
editWe then arrive at
For comparison, the equation (6) p.53-5b gives
We now find the quality factor, Q. From (2) p.53-5d, we have
From (5) p.53-5d, we know that the loss factor is just the inverse of the quality factor so
Problem R2.3
editProblem Statement
editKS 2008 p.103 sec.2.7 pb.21 write a matlab program to assemble the global stiffness matrix, compute the unknown disp dofs, the reaction forces, the member forces; compare the results with exact hand calculation Statics, Mechanics of Materials).
run CALFEM to verify the results, but now do this problem in a different way (as in commercial FE codes).
first, establish 2 arrays: (1) the global node coordinate array “coord”, and (2) the elem connectivity array “conn”. next, write a matlab code to loop over the number of element to call the function “bar2e” of CALFEM to compute the element stiffness matrices. with this method, you let your matlab code construct the arrays ex1, ey1, ex2, ey2, etc. for you, using the info from the arrays “coord” and “conn”.
Solution
editMATLAB Code
editThe following is a MATLAB code that assembles the global stiffness matrix, computes the unknown displacement dofs, reaction forces, and member forces
function FEA_R2p3 E = 10000; % Young's modulus A = pi*2^2/4; % Area % Lengths of each element L1 = 12; L2 = 12; L3 = 12; L4 = sqrt((12^2)+(9^2)); L5 = 9; L6 = sqrt((12^2)+(9^2)); L7 = 12; L8 = 9; L9 = sqrt((12^2)+(9^2)); % topology matrix % [element# node1x node1y node2x node2y] Edof = [1 1 2 3 4;2 3 4 5 6;3 5 6 7 8; 4 1 2 9 10;5 9 10 3 4;6 9 10 5 6; 7 9 10 11 12;8 11 12 5 6;9 11 12 7 8]; % element stiffness matrix in local coordinates for elements 1 through 9 esmloc1 = (E*A/L1)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc2 = (E*A/L2)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc3 = (E*A/L3)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc4 = (E*A/L4)*[(L1/L4)^2 (L1/L4)*(L5/L4) -(L1/L4)^2 -(L1/L4)*(L5/L4); (L1/L4)*(L5/L4) (L5/L4)^2 -(L1/L4)*(L5/L4) -(L5/L4)^2; -(L1/L4)^2 -(L1/L4)*(L5/L4) (L1/L4)^2 (L1/L4)*(L5/L4); -(L1/L4)*(L5/L4) -(L5/L4)^2 (L1/L4)*(L5/L4) (L5/L4)^2]; % [cos^2 cossin -cos^2 -cossin;cossin sin^2 -cossin -sin^2; % -cos^2 -cossin cos^2 cossin;-cossin -sin^2 cossin sin^2] esmloc5 = (E*A/L5)*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1]; esmloc6 = (E*A/L6)*[(L1/L4)^2 -(L1/L4)*(L5/L4) -(L1/L4)^2 (L1/L4)*(L5/L4); -(L1/L4)*(L5/L4) (L5/L4)^2 (L1/L4)*(L5/L4) -(L5/L4)^2; -(L1/L4)^2 (L1/L4)*(L5/L4) (L1/L4)^2 -(L1/L4)*(L5/L4); (L1/L4)*(L5/L4) -(L5/L4)^2 -(L1/L4)*(L5/L4) (L5/L4)^2]; % [cos^2 -cossin -cos^2 cossin;-cossin sin^2 cossin -sin^2; % -cos^2 cossin cos^2 -cossin;cossin -sin^2 -cossin sin^2] esmloc7 = (E*A/L7)*[1 0 -1 0;0 0 0 0;-1 0 1 0;0 0 0 0]; esmloc8 = (E*A/L8)*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1]; esmloc9 = (E*A/L9)*[(L1/L4)^2 -(L1/L4)*(L5/L4) -(L1/L4)^2 (L1/L4)*(L5/L4); -(L1/L4)*(L5/L4) (L5/L4)^2 (L1/L4)*(L5/L4) -(L5/L4)^2; -(L1/L4)^2 (L1/L4)*(L5/L4) (L1/L4)^2 -(L1/L4)*(L5/L4); (L1/L4)*(L5/L4) -(L5/L4)^2 -(L1/L4)*(L5/L4) (L5/L4)^2]; % [cos^2 -cossin -cos^2 cossin;-cossin sin^2 cossin -sin^2; % -cos^2 cossin cos^2 -cossin;cossin -sin^2 -cossin sin^2] % element nodal coordinates % ex1 and ey1 are the x and y coordinates to the nodes of element 1 ex1 = [0 12]; ey1 = [0 0]; ex2 = [12 24]; ey2 = [0 0]; ex3 = [24 36]; ey3 = [0 0]; ex4 = [0 12]; ey4 = [0 9]; ex5 = [12 12]; ey5 = [9 0]; ex6 = [12 24]; ey6 = [9 0]; ex7 = [12 24]; ey7 = [9 9]; ex8 = [24 24]; ey8 = [9 0]; ex9 = [24 36]; ey9 = [9 0]; K = zeros(12); F = zeros(12,1); F(6)=-1200; F(11)=400; ep = [E A]; % assembly of the global stiffness matrix K = assem(Edof(1,:),K,esmloc1); K = assem(Edof(2,:),K,esmloc2); K = assem(Edof(3,:),K,esmloc3); K = assem(Edof(4,:),K,esmloc4); K = assem(Edof(5,:),K,esmloc5); K = assem(Edof(6,:),K,esmloc6); K = assem(Edof(7,:),K,esmloc7); K = assem(Edof(8,:),K,esmloc8); K = assem(Edof(9,:),K,esmloc9); % boundary conditions bc = [1 0;2 0;8 0]; % reaction force vector Q = solveq(K,F,bc); % finding element displacements ed1 = extract(Edof(1,:),Q); ed2 = extract(Edof(2,:),Q); ed3 = extract(Edof(3,:),Q); ed4 = extract(Edof(4,:),Q); ed5 = extract(Edof(5,:),Q); ed6 = extract(Edof(6,:),Q); ed7 = extract(Edof(7,:),Q); ed8 = extract(Edof(8,:),Q); ed9 = extract(Edof(9,:),Q); disp = [ed1;ed2;ed3;ed4;ed5;ed6;ed7;ed8;ed9]; % reaction forces at each element; member forces (stress) N1 = bar2s(ex1,ey1,ep,ed1); MF1 = N1/A; N2 = bar2s(ex2,ey2,ep,ed2); MF2 = N2/A; N3 = bar2s(ex3,ey3,ep,ed3); MF3 = N3/A; N4 = bar2s(ex4,ey4,ep,ed4); MF4 = N4/A; N5 = bar2s(ex5,ey5,ep,ed5); MF5 = N5/A; N6 = bar2s(ex6,ey6,ep,ed6); MF6 = N6/A; N7 = bar2s(ex7,ey7,ep,ed7); MF7 = N7/A; N8 = bar2s(ex8,ey8,ep,ed8); MF8 = N8/A; N9 = bar2s(ex9,ey9,ep,ed9); MF9 = N9/A; RF = [N1;N2;N3;N4;N5;N6;N7;N8;N9]; MF = [MF1;MF2;MF3;MF4;MF5;MF6;MF7;MF8;MF9];
With this, the output is as follows:
GlobalStiffnessMatrix = K dispDofs = disp reactionForces = RF memberForces = MF
GlobalStiffnessMatrix = 1.0e+003 * 3.9584 1.0053 -2.6180 0 0 0 0 0 -1.3404 -1.0053 0 0 1.0053 0.7540 0 0 0 0 0 0 -1.0053 -0.7540 0 0 -2.6180 0 5.2360 0 -2.6180 0 0 0 0 0 0 0 0 0 0 3.4907 0 0 0 0 0 -3.4907 0 0 0 0 -2.6180 0 6.5764 -1.0053 -2.6180 0 -1.3404 1.0053 0 0 0 0 0 0 -1.0053 4.2446 0 0 1.0053 -0.7540 0 -3.4907 0 0 0 0 -2.6180 0 3.9584 -1.0053 0 0 -1.3404 1.0053 0 0 0 0 0 0 -1.0053 0.7540 0 0 1.0053 -0.7540 -1.3404 -1.0053 0 0 -1.3404 1.0053 0 0 5.2988 0 -2.6180 0 -1.0053 -0.7540 0 -3.4907 1.0053 -0.7540 0 0 0 4.9986 0 0 0 0 0 0 0 0 -1.3404 1.0053 -2.6180 0 3.9584 -1.0053 0 0 0 0 0 -3.4907 1.0053 -0.7540 0 0 -1.0053 4.2446
dispDofs = 0 0 0.3056 -1.4992 0.3056 -1.4992 0.6112 -2.1836 0.6112 -2.1836 1.0695 0 0 0 0.8260 -1.4992 0.8260 -1.4992 0.3056 -1.4992 0.8260 -1.4992 0.6112 -2.1836 0.8260 -1.4992 0.5204 -1.9258 0.5204 -1.9258 0.6112 -2.1836 0.5204 -1.9258 1.0695 0
reactionForces = 1.0e+003 * 0.8000 0.8000 1.2000 -0.5000 0 0.5000 -0.8000 0.9000 -1.5000
memberForces = 254.6479 254.6479 381.9719 -159.1549 0 159.1549 254.6479 286.4789 -477.4648
It can be seen here that the problem done by hand corresponds to the force output of the MATLAB program.
CALFEM MATLAB
edit% CALFEM % global coordinate matrix Coord = [0 0;12 0;24 0;36 0;12 9;24 9]; % element connectivity array conn = Edof; % global nodal dof matrix Dof = [1 2;3 4;5 6;7 8;9 10;11 12]; % contructing element arrays [ex,ey]=coordxtr(conn,Coord,Dof,2); K=zeros(12); EA=[E A]; % Young's modulus, Area % compute element stiffness matrix Ke and global stiffness matrix for i=1:9 Ke = bar2e(ex(i,:),ey(i,:),EA); % element stiffness matrix K = assem(conn(i,:),K,Ke); % global stiffness matrix end % element displacement vector ed=extract(conn,Q); % reaction forces for i=1:9 N(i)=bar2s(ex(i,:),ey(i,:),ep,ed(i,:)); end % member forces MF1C = N(1)/A; % stress in AC MF2C = N(2)/A; % stress in CE MF3C = N(3)/A; % stress in EF MF4C = N(4)/A; % stress in AB MF5C = 0; % stress in BC MF6C = N(6)/A; % stress in BE MF7C = N(7)/A; % stress in BD MF8C = N(8)/A; % stress in DE MF9C = N(9)/A; % stress in DF MFC = [MF1C;MF2C;MF3C;MF4C;MF5C;MF6C;MF7C;MF8C;MF9C];
The output to the CALFEM code is as follows:
CGlobalStiffnessMatrix = K CdispDofs = ed CreactionForces = N' CmemberForces = MFC
CGlobalStiffnessMatrix = 1.0e+003 * 3.9584 1.0053 -2.6180 0 0 0 0 0 -1.3404 -1.0053 0 0 1.0053 0.7540 0 0 0 0 0 0 -1.0053 -0.7540 0 0 -2.6180 0 5.2360 0 -2.6180 0 0 0 0 0 0 0 0 0 0 3.4907 0 0 0 0 0 -3.4907 0 0 0 0 -2.6180 0 6.5764 -1.0053 -2.6180 0 -1.3404 1.0053 0 0 0 0 0 0 -1.0053 4.2446 0 0 1.0053 -0.7540 0 -3.4907 0 0 0 0 -2.6180 0 3.9584 -1.0053 0 0 -1.3404 1.0053 0 0 0 0 0 0 -1.0053 0.7540 0 0 1.0053 -0.7540 -1.3404 -1.0053 0 0 -1.3404 1.0053 0 0 5.2988 0 -2.6180 0 -1.0053 -0.7540 0 -3.4907 1.0053 -0.7540 0 0 0 4.9986 0 0 0 0 0 0 0 0 -1.3404 1.0053 -2.6180 0 3.9584 -1.0053 0 0 0 0 0 -3.4907 1.0053 -0.7540 0 0 -1.0053 4.2446
CdispDofs = 0 0 0.3056 -1.4992 0.3056 -1.4992 0.6112 -2.1836 0.6112 -2.1836 1.0695 0 0 0 0.8260 -1.4992 0.8260 -1.4992 0.3056 -1.4992 0.8260 -1.4992 0.6112 -2.1836 0.8260 -1.4992 0.5204 -1.9258 0.5204 -1.9258 0.6112 -2.1836 0.5204 -1.9258 1.0695 0
CreactionForces = 1.0e+003 * 0.8000 0.8000 1.2000 -0.5000 0 0.5000 -0.8000 0.9000 -1.5000
CmemberForces = 254.6479 254.6479 381.9719 -159.1549 0 159.1549 254.6479 286.4789 -477.4648
The results of all methods verifies their accuracy.
Problem R2.4
editDescription:
editWe are to find an expression for the excitation referenced in (1) p.53-6 and then find the actual solution and the amplification factor, A for the following conditions:
We are then to find the complete solution.
Solution:
edit
Let us guess the solution is:
then
We then write
It follows that
and that
We solve for A and B to be
Once given the conditions
we can find
Plugging these values in yields
The particular solution is then
Equations (3) p. 53-8 and (2) p. 53-8 give us
The following is the graph of the homogeneous solution:
The following is the graph of the particular solution:
The following is the graph of the total solution, :
A graph of all three plots is produced below:
Problem R2.5 CALFEM 3.4 Manual exs1 on a Spring System and Verification
editGiven a Linearly Elastic Spring System
editConsider a system of 3 linearly elastic springs
Use the following diagram for the node locations
Find:
edit1.Find the results when running exs1 from CALFEM
edit2.Verify the results of CALFEM to manual calculations
editSolution: Results of CALFEM
editConstruct the topology matrix
editThis matrix contains the node numbers and degrees of freedom.
Construct the global stiffness matrix and load vector f
edit
The load vector is in position 2 with F=100
Generating the element stiffness matrices
editThe element stiffness matrices are generated using the CALFEM function spring1e. The matrices are generated using k and 2k where k=1500
Assembling the element stiffness matrices into the global stiffness matrix
editThe assembly is done using the assem fucntion from CALFEM.
Assembling the second element at the first row:
Assembling the first element at the second row:
Assembling the second element at the third row:
Solve the system of equations using boundary conditions
editThe CALFEM fucntion solveq is used to solve the system of equations subject to certain boundary conditions.
Evaluate the element forces from the element displacements
edit
Evaluate the spring forces from function spring1s
edit
Solution: Results of verification
editThis section verifies the CALFEM code by manual calculations
Nodal equilibrium equations
edit
Setting up the equation to solve
edit
Implementing the boundary conditions
edit
Inverting the global stiffness matrix to find displacements
edit
Displacements 1 and 3 are zero due to boundary conditions.
Calculate forces in springs from equation
editThe forces in the springs are found by the following equations. Negative forces display compressions and positive tensions.
Discuss if the solution is verified
editThe solution by CALFEM is verified since both the displacements and the forces are the same by both methods.
Problem R2.6
editProblem Statement
editFind the displacements of the nodes as well as the forces in the springs(tension/compression). Also determine the reaction forces at the walls.
Nodal Equilibrium Equations
edit
Using Matlab
Solutions
Displacements
These are the displacements at , , and respectively. The negative sign shows that the nodes are moving in the direction of the applied force which is against the convention of positive displacement left to right.
Now it is possible to find the forces in the spring. Negative forces denote compression and positive forces denote tension.
Forces in Springs
editNow,
Reaction Forces
edit
Verification with Calfem
edita values are the nodal displacements 1 through 5.
r values are the reaction forces at the walls.
es1-es6 are the forces in the springs. Negative values denote compression and positive values denote tension.
These answers match those previously calculated.
===Code===
Edof=[1 1 2;2 2 4;3 2 3;4 1 3;5 3 4;6 4 5];
ep1=500;
ep2=400;
ep3=600;
ep4=200;
ep5=400;
ep6=300;
Ke1=spring1e(ep1);
Ke2=spring1e(ep2);
Ke3=spring1e(ep3);
Ke4=spring1e(ep4);
Ke5=spring1e(ep5);
Ke6=spring1e(ep6);
K=zeros(5,5);
K=assem(Edof(1,:),K,Ke1);
K=assem(Edof(2,:),K,Ke2);
K=assem(Edof(3,:),K,Ke3);
K=assem(Edof(4,:),K,Ke4);
K=assem(Edof(5,:),K,Ke5);
K=assem(Edof(6,:),K,Ke6);
bc= [1 0; 5 0];
f=zeros(5,1);
f(3)=-1000;
[a,r]=solveq(K,f,bc)
ed1=extract(Edof(1,:),a);
ed2=extract(Edof(2,:),a);
ed3=extract(Edof(3,:),a);
ed4=extract(Edof(4,:),a);
ed5=extract(Edof(5,:),a);
ed6=extract(Edof(6,:),a);
es1=spring1s(ep1,ed1)
es2=spring1s(ep2,ed2)
es3=spring1s(ep3,ed3)
es4=spring1s(ep4,ed4)
es5=spring1s(ep5,ed5)
es6=spring1s(ep6,ed6)
'''Solution'''
a =
0
-0.8542
-1.5521
-0.8750
0
r =
737.5000
-0.0000
0.0000
0.0000
262.5000
es1 =
-427.0833
es2 =
-8.3333
es3 =
-418.7500
es4 =
-310.4167
es5 =
270.8333
es6 =
262.5000
--EML4507.s13.team4ever.Wulterkens (discuss • contribs) 07:21, 6 February 2013 (UTC)