User:Eml4507.s13.team3.steiner/Team Negative Damping (3): Repart 6

6.4
= Problem 1: Pb. 3.1, p. 3-8 Section 3 Notes using Matlab and CALFEM = On my honor, I have neither given nor received unauthorized aid in doing this assignment.

Given
Plane truss as shown. A force, F, of 20000 N is applied as shown. The length of each element is 1 meter. Element 2 is vertical and the elements are spaced evenly around node 1. Each element has the same physical properties: Young's modulus of 206 GPa and cross sectional area of .0001 m. Additionally, node 2 experiences an initial displacement of 2 cm in the x direction and -1 cm in the y direction; node 3 experiences an initial displacement of -3 cm in the x direction and 5 cm in the y direction.

Find
Find the unknown displacements and forces at each global node as well as the member forces of each element using Matlab and verify using CALFEM. Plot the deformed and undeformed shape on the same plot.

Using Matlab
Enter the number of global nodes. >> G=10;

Enter the x, y, and z coordinates. >> x=[-1.5 1.5 -1.5 1.5 1.5 -1.5 -4 4 4 -4]; >> y=[0 0 1.5 1.5 -1.5 -1.5 4 4 -4 -4]; >> z=[8 8 4 4 4 4 0 0 0 0];

Enter the input force matrix. The first element of each row corresponds to the global degree of freedom; the second element corresponds to the force associated with that degree of freedom. Only the known forces are entered into this matrix. >> F_in=[1 0; 2 60000; 3 0; 4 0; 5 60000; 6 0; 7 0; 8 0; 9 0; 10 0; 11 0; 12 0; 13 0; 14 0; 15 0; 16 0; 17 0; 18 0];

Enter the input displacement matrix. The first element of each row corresponds to the global degree of freedom; the second element corresponds to the displacement associated with that degree of freedom. >> Q_in=[19 0; 20 0; 21 0; 22 0; 23 0; 24 0; 25 0; 26 0; 27 0; 28 0; 29 0; 30 0];

Enter the number of elements. >> E=25;

Enter the global node corresponding to the first local node of each element. The column number corresponds to the element number of the given node. >> l1=[1 1 2 1 2 2 2 1 1 3 4 3 5 3 6 4 5 4 3 5 6 6 3 4 5];

Enter the global node corresponding to the second local node of each element. >> l2=[2 4 3 5 6 4 5 3 6 6 5 4 6 10 7 9 8 7 8 10 9 10 7 8 9];

Enter the Young's Modulus associated with each element. The column number corresponds to the element number of the given node. >> Modulus=3e7*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];

Enter the Area associated with each element. The column number corresponds to the element number of the given node. >> Area=(pi/4)*2^2*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];

Enter the Mass Density associated with each element. The column number corresponds to the element number of the given node. >> Ro=7.3e-4*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; Run the truss program. >> [F_bar] = truss(G, x, y, z, F_in, Q_in, E, l1, l2, Modulus, Area, Ro)

The internal spring force matrix F_bar F_bar = 1.0e+04 * Columns 1 through 9 0.0000   3.5997    3.5997   -3.5997   -3.5997    5.5981   -5.5981    5.5981   -5.5981      -0.0000   -3.5997   -3.5997    3.5997    3.5997   -5.5981    5.5981   -5.5981    5.5981     Columns 10 through 18 -0.0000  -0.0000   -0.9053    0.9053    1.8114   -1.8114    1.8114   -1.8114    3.4748       0.0000    0.0000    0.9053   -0.9053   -1.8114    1.8114   -1.8114    1.8114   -3.4748     Columns 19 through 25 3.4748  -3.4748   -3.4748   -6.7822    6.7822    6.7822   -6.7822      -3.4748    3.4748    3.4748    6.7822   -6.7822   -6.7822    6.7822

Row 1 gives the member forces with respect to the Global nodes and row 2 gives the member forces with respect to the members.

The axial stresses are given by dividing the axial forces by the area. stress = 1.0e+04 * Columns 1 through 9 -0.0000  -1.1458   -1.1458    1.1458    1.1458   -1.7819    1.7819   -1.7819    1.7819     Columns 10 through 18 0.0000   0.0000    0.2882   -0.2882   -0.5766    0.5766   -0.5766    0.5766   -1.1061     Columns 19 through 25 -1.1061   1.1061    1.1061    2.1588   -2.1588   -2.1588    2.1588



Deformed (green) and undeformed (blue) plot of the Truss. The truss doesn't deform much due to the given forces.

Matlab Code
This code is equivalent to that used in R2.3 except that it uses the given spring constants rather than than the modulus and area to determine the stiffness matrix.

CALFEM Verification
First begin by constructing the Edof  matrix. Column 1 corresponds to the element number; columns 2 and 3 correspond to the global node that corresponds to the first and second local nodes respectively. >> Edof=[1 1 2 5 6 2 1 2 3 4           3 1 2 7 8];

Construct an empty stiffness matrix K that has side lengths equal to the number of degrees of freedom times the number of global nodes. For this problem, there are 4 global nodes and motion is constrained to 2 degrees of freedom. >> K=zeros(8,8);

Construct the force matrix F. The forces corresponding to global node 1 are known; the forces corresponding to global nodes 2, 3, and 4 are unknown. >> F=zeros(8,1); >> F(1)=20000/sqrt(2); >> F(2)=20000/sqrt(2);

Construct the displacement matrix Q. Global nodes 2, 3, and 4 and known displacements as shown. >> Q=[3 .02 4 -.01        5 -.03         6 .05         7 0         8 0];

The function bar2e generates the element stiffness matrices for a bar with 2 degrees of freedom. Its parameters are the x and y coordinates of the element and a vector containing the Young's modulus and cross sectional area. This information is given in the problem statement.

The Vector containing the Young's modulus and cross sectional area. >> E=206e9; >> A=1e-4; >> ep=[E A];

The x and y coordinates for each element. >> ex1=[0 sqrt(3)/2]; >> ey1=[0 -.5]; >> ex2=[0 0]; >> ey2=[0 1]; >> ex3=[0 -sqrt(3)/2]; >> ey3=[0 -.5];

The commands: >> Ke1=bar2e(ex1, ey1, ep); >> Ke2=bar2e(ex2, ey2, ep); >> Ke3=bar2e(ex3, ey3, ep);

Yield the following stiffness matrices: >> Ke1 = 1.0e+07 * 1.5450  -0.8920   -1.5450    0.8920          -0.8920    0.5150    0.8920   -0.5150          -1.5450    0.8920    1.5450   -0.8920           0.8920   -0.5150   -0.8920    0.5150      Ke2 = 0          0           0           0                0    20600000           0   -20600000                0           0           0           0                0   -20600000           0    20600000   >> Ke3 = 1.0e+07 * 1.5450   0.8920   -1.5450   -0.8920          0.8920    0.5150   -0.8920   -0.5150         -1.5450   -0.8920    1.5450    0.8920         -0.8920   -0.5150    0.8920    0.5150

The global stiffness matrix can be assembled with the element stiffness matrices using the following command. >> K=assem(Edof(1,:),K,Ke1); >> K=assem(Edof(2,:),K,Ke2); >> K=assem(Edof(3,:),K,Ke3);

This results in the following stiffness matrix. >> K = 1.0e+07 * 3.0900        0         0         0   -1.5450    0.8920   -1.5450   -0.8920               0    3.0900         0   -2.0600    0.8920   -0.5150   -0.8920   -0.5150               0         0         0         0         0         0         0         0               0   -2.0600         0    2.0600         0         0         0         0         -1.5450    0.8920         0         0    1.5450   -0.8920         0         0          0.8920   -0.5150         0         0   -0.8920    0.5150         0         0         -1.5450   -0.8920         0         0         0         0    1.5450    0.8920         -0.8920   -0.5150         0         0         0         0    0.8920    0.5150

This stiffness matrix is in accord with that computed in Matlab in the previous section.

The function solveq takes the stiffness matrix, known force matrix, and displacement matrix and returns the complete force and displacement matrices. >> [q r]=solveq(K,F,Q) q = -0.0290      0.0108       0.0200      -0.0100      -0.0300       0.0500            0            0   r = 1.0e+05 * 0.0000      0.0000            0      -4.2816      -3.6562       2.1109       3.5148       2.0293

The displacement and force matrices are in accord with those computed in Matlab in the previous section.

The function extract is used to evaluate the displacements of the local nodes for each beam, which is needed to determine the force in each beam. >> ed1=extract(Edof(1,:),q); >> ed2=extract(Edof(2,:),q); >> ed3=extract(Edof(3,:),q);

This results in the following displacements. The number following 'ed' corresponds to the element number. Each column corresponds to the respective degree of freedom for the beam. ed1 = -0.0290   0.0108   -0.0300    0.0500 ed2 = -0.0290   0.0108    0.0200   -0.0100 ed3 = -0.0290   0.0108         0         0

The function bar2s can be used to determine the force in each beam. Its input parameters are the x and y coordinates for the beam (ex1 and ey1), the physical properties of the beam (Young's modulus and cross sectional area - ep) and the displacement (ed). >> es1=bar2s(ex1, ey1,ep,ed1) >> es1=bar2s(ex2, ey2,ep,ed2) >> es1=bar2s(ex3, ey3,ep,ed3)

This results in the following spring forces. es1 = -4.2219e+05 es2 = -4.2816e+05 es3 = -4.0586e+05

This is in accord with the forces computed in Matlab in the previous section.