User:Eml4507.s13.team3.steiner

/Team Negative Damping (3): Report 2/ /Team Negative Damping (3): Report 3/ /Team Negative Damping (3): Report 4/ /Team Negative Damping (3): Report 5/ /Team Negative Damping (3): Report 6/ /Team Negative Damping (3): Report 7/

= Problem 2: Equation of Motion of a Spring-Mass-Dashpot System =

Problem Statement
Derive the equation of motions for the two spring-mass-dashpot systems in series.



Solution
In order to derive an equation of motion, all of the forces must be analyzed. In the spring-mass-dashpot model there are forces from the spring, driving force, and a damping force.



The force of a spring oriented in the y-direction is given by Hooke's Law:
 * {| style="width:100%" border="0"

$$\displaystyle F_{k}= ky$$
 * style="width:95%" |
 * style="width:95%" |
 * }

The force of a damper oriented in the y-direction is:
 * {| style="width:100%" border="0"

$$\displaystyle F_{c}=cy'$$
 * style="width:95%" |
 * style="width:95%" |
 * }

Summing the forces with respect to their directions in the Free Body Diagram:
 * {| style="width:100%" border="0"

$$\displaystyle -F_{k}- F_{c}- my''+ F(t)=0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Plugging in spring and damper values:
 * {| style="width:100%" border="0"

$$\displaystyle -ky -cy'- my''+ F(t)=0 $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Reorganizing  the equation we get:$$\displaystyle $$
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |
 * }

Solution
In order to derive an equation of motion, all of the forces must be analyzed. Once again, in the spring-mass-dashpot model there are forces from the spring, driving force, and damping force.



The forces directions with respect to the Free Body Diagram:
 * {| style="width:100%" border="0"

$$\displaystyle my'' + F_{1} = F(t)$$
 * style="width:95%" |
 * style="width:95%" |


 * }

To find the force from the spring and damper, we set them equal to eachother:


 * {| style="width:100%" border="0"

$$\displaystyle F_{k}=F_{c} $$
 * style="width:95%" |
 * style="width:95%" |
 * }


 * {| style="width:100%" border="0"

$$\displaystyle ky_{k}=Cy'_{c} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Solving for dampers velocity: Summing the forces with respect to their directions in the Free Body Diagram:
 * {| style="width:100%" border="0"

$$\displaystyle y''_{c}= y'_{k}\frac{k}{c}$$
 * style="width:95%" |
 * style="width:95%" |
 * }

Plugging into acceleration equations after taking the derivative:
 * {| style="width:100%" border="0"

$$\displaystyle y = y_{c} + y''_{k} $$
 * style="width:95%" |
 * style="width:95%" |
 * }
 * {| style="width:100%" border="0"

$$\displaystyle y = y_{k} + y'_{k}\frac{k}{c} $$
 * style="width:95%" |
 * style="width:95%" |
 * }

Plugging into the original:$$\displaystyle $$
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * style="width:95%" |
 * }

Problem Statement
Part 1: Provide the numerical values for the element stiffness matrices $$ \mathbf{k^{(1)}} $$ and $$ \mathbf{k^{(2)}} $$.

Part 2: Find the coordinates of the global nodes 1 and 3 when the global node 2 is at the origin of the global coordinate system.

The following are physical parameters:

$$ E^{(1)} = 3, E^{(2)} = 5, $$

$$ A^{(1)} = 1, A^{(2)} = 2, $$

$$ L^{(1)} = 4, L^{(2)} = 2 $$

$$ \theta^{(1)}=30^{o}, \theta^{(2)}=-45^{o} $$

Part 1
$$ \mathbf{k^{(e)}} = k^{(e)} \begin{bmatrix} (l^{(e)})^{2} & l^{(e)}m^{(e)} & -(l^{(e)})^{2} & -l^{(e)}m^{(e)}\\ l^{(e)}m^{(e)} & (m^{(e)})^{2} & -l^{(e)}m^{(e)} & -(m^{(e)})^{2}\\ -(l^{(e)})^{2}& -l^{(e)}m^{(e)} & (l^{(e)})^{2} & l^{(e)}m^{(e)}\\ -l^{(e)}m^{(e)} & -(m^{(e)})^{2} & l^{(e)}m^{(e)} & (m^{(e)})^{2} \end{bmatrix} $$

where

$$ l^{(e)} = cos(\theta^{(e)}), m^{(e)} = sin(\theta^{(e)}) $$

and

$$ k^{(e)}={\frac{E^{(e)}A^{(e)}}{L^{(e)}}} $$

As given in the problem statement,

$$ E^{(1)} = 3, E^{(2)} = 5, $$

$$ A^{(1)} = 1, A^{(2)} = 2, $$

$$ L^{(1)} = 4, L^{(2)} = 2 $$

$$ \theta^{(1)}=30^{o}, \theta^{(2)}=-45^{o} $$

Therefore,

$$ k^{(1)}={\frac{3*1}{4}}={\frac{3}{4}}, k^{(2)}={\frac{5*2}{2}}=5, $$

$$ (l^{(1)})^{2}=.75, (m^{(1)})^{2}=.25, l^{(1)}m^{(1)}=.433 $$

and

$$ (l^{(2)})^{2}=.5, (m^{(2)})^{2}=.5, l^{(2)}m^{(2)}=-.5 $$

Substituting these values back into the original stiffness matrix, $$ \mathbf{k^{(e)}} $$ yields:

$$ \mathbf{k^{(1)}} = .75 \begin{bmatrix} .75 & .433 & -.75 & -.433\\ .433 & .25 & -.433 & -.25\\ -.75 & -.433 & .75 & .433\\ -.433 & -.25 & .433 & .25 \end{bmatrix} = \begin{bmatrix} .5625 & .3248 & -.5625 & -.3248\\ .3248 & .1875 & -.3248 & -.1875\\ -.5625 & -.3248 & .5625 & .3248\\ -.3248 & -.1875 & .3248 & .1875 \end{bmatrix} $$

and

$$ \mathbf{k^{(2)}} = 5 \begin{bmatrix} .5 & -.5 & -.5 & .5\\ -.5 & .5 & .5 & -.5\\ -.5 & .5 & .5 & -.5\\ .5 & -.5 & -.5 & .5 \end{bmatrix} = \begin{bmatrix} 2.5 & -2.5 & -2.5 & 2.5\\ -2.5 & 2.5 & 2.5 & -2.5\\ -2.5 & 2.5 & 2.5 & -2.5\\ 2.5 & -2.5 & -2.5 & 2.5 \end{bmatrix} $$

Part 2
Global Nodes 1 and 2 can be located in our coordinate system using the angle and length of each respective element.

However, each angle must be measured from the right horizontal on global node 2. $$ \theta^{(2)} $$ remains the same, but $$ \theta^{(1)}=210^{o} $$.

Trigonometry shows that the coordinates of global nodes 1 and 3,

$$ (x_{1},y_{1}) $$ and $$ (x_{2},y_{2}) $$

can be written as

$$ (L^{(1)}cos(\theta^{(1)}),L^{(1)}sin(\theta^{(1)})) $$ and $$ (L^{(2)}cos(\theta^{(2)}),L^{(2)}sin(\theta^{(2)})) $$ respectively.

Therefore, global nodes 1 and 3 have the following coordinates:

$$ (-2\sqrt{3},-2) $$ and $$ (2\sqrt{2},-2\sqrt{2}) $$

= R2.3 Exercise 21, p. 103 Kim and Sankar: Matlab, Statics, and CALFEM =

Given
Determine the global stiffness matrix, the unknown displacement degrees of freedom, the reaction forces, and the member forces for the given truss system using Matlab. Verify the reaction and member forces with those computed using a Statics approach. Verify the global stiffness matrix, the unknown displacement degrees of freedom, the reaction forces, and the member forces using CALFEM.

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

Enter the x, y, and z coordinates (this program was designed solve 1D, 2D, and 3D problems, so the z coordinates must be entered. As long as all nodes have the same z coordinates, the solution will match a 2D problem). >> x=[ 0  12  12  24   0 -12]; >> y=[ 9   0   9   0   0   0]; >> z=[ 0   0   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    0             4     0             5 -1200             7   400             8     0            10     0            13     0            14     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. Only the known displacements are entered into this matrix. If the force for a particular degree of freedom is unknown, the its displacement is known. Because there is no force in the z direction, it is assumed that all displacements in the z direction are equal to 0. >> Q_in=[ 3 0 6 0            9 0            11 0            12 0            15 0            16 0            17 0            18 0];

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

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=[2 3 4 1 5 5 1 6 3];

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

Enter the spring constant associated with each element. The column number corresponds to the element number of the given node. >> Modulus=[10000 10000 10000 10000 10000 10000 10000 10000 10000];

Enter the spring constant associated with each element. The column number corresponds to the element number of the given node. >> Area=[3.14159 3.14159 3.14159 3.14159 3.14159 3.14159 3.14159 3.14159 3.14159];

Run the displacement program. >> [k Q F Q_bar F_bar]=displacement(G, x, y, z, F_in, Q_in, E, l1, l2, spring_const);

The stiffness matrix k: k = 1.0e+03 * 5.2988        0   -1.3404    1.0053   -2.6180         0         0         0         0         0   -1.3404   -1.0053            0    4.9986    1.0053   -0.7540         0         0         0         0         0   -3.4907   -1.0053   -0.7540      -1.3404    1.0053    6.5764   -1.0053         0         0   -2.6180         0   -2.6180         0         0         0       1.0053   -0.7540   -1.0053    4.2446         0   -3.4907         0         0         0         0         0         0      -2.6180         0         0         0    3.9584   -1.0053   -1.3404    1.0053         0         0         0         0            0         0         0   -3.4907   -1.0053    4.2446    1.0053   -0.7540         0         0         0         0            0         0   -2.6180         0   -1.3404    1.0053    3.9584   -1.0053         0         0         0         0            0         0         0         0    1.0053   -0.7540   -1.0053    0.7540         0         0         0         0            0         0   -2.6180         0         0         0         0         0    5.2360         0   -2.6180         0            0   -3.4907         0         0         0         0         0         0         0    3.4907         0         0      -1.3404   -1.0053         0         0         0         0         0         0   -2.6180         0    3.9584    1.0053      -1.0053   -0.7540         0         0         0         0         0         0         0         0    1.0053    0.7540

The displacement matrix Q: Q = 1.0000   0.8260       2.0000   -1.4992       3.0000         0       4.0000    0.6112       5.0000   -2.1836       6.0000         0       7.0000    0.5204       8.0000   -1.9258       9.0000         0      10.0000    1.0695      11.0000         0      12.0000         0      13.0000    0.3056      14.0000   -1.4992      15.0000         0      16.0000         0      17.0000         0      18.0000         0

Every 3 rows correspond to a separate global node. (ie. global node 1 experiences a positive displacement in the x direction, a negative displacement in the y direction, and no displacement in the z direction) It can be seen that global node 4 (rows 10-12) experiences only a positive displacement in the x direction.

The Force matrix, F F = 1.0e+03 * -0.0000      0.0000            0       0.0000      -1.2000            0       0.4000       0.0000            0      -0.0000       0.9000            0       0.0000            0            0      -0.4000       0.3000            0

This shows that at global node 4, a positive vertical force of 900 lb is exerted on the truss. Additionally, at global node 6 a negative horizontal force of 400 lb and a positive vertical force of 300 lb is exerted on the truss. This is in accord with a force balance of the entire truss.

The local displacement matric Q_bar

Q_bar = -1.7991   1.5718   -1.0695    1.4992    0.3056   -0.3056    0.8260         0    1.9258      -1.5604    0.8556   -0.6112    1.4992    0.6112         0    0.5204   -0.2387    2.1836

The internal spring force matrix F_bar F_bar = 1.0e+03 * -0.5000   1.5000   -1.2000         0   -0.8000   -0.8000    0.8000    0.5000   -0.9000       0.5000   -1.5000    1.2000         0    0.8000    0.8000   -0.8000   -0.5000    0.9000

Matlab Code

This declares the function and a few variables.

function [k, Q, F, Q_bar, F_bar] = displacement(G, x, y, z, F_in, Q_in, E, l1, l2, spring_const) %Element Computation L = zeros(1,E); l = zeros(1,E); m = zeros(1,E); n = zeros(1,E); k_elem = zeros(3,3,E); This for loop computes the length, direction cosines, and the top left quarter of the global k matrix for each element.

for i=1:E L(i)=sqrt((x(l2(i))-x(l1(i)))^2+(y(l2(i))-y(l1(i)))^2+(z(l2(i))-z(l1(i)))^2); l(i)=(x(l2(i))-x(l1(i)))/L(i); m(i)=(y(l2(i))-y(l1(i)))/L(i); n(i)=(z(l2(i))-z(l1(i)))/L(i); k_elem(:,:,i) = (Modulus(i)*Area(i))/L(i)*[l(i)*l(i), l(i)*m(i), l(i)*n(i); m(i)*l(i), m(i)*m(i), m(i)*n(i); n(i)*l(i), n(i)*m(i), n(i)*n(i)]; end

This nested for loop constructs the stiffness matrix.

%Stiffness Matrix Construction k = zeros(3*G,3*G); for i=1:E for j=1:3 for h=1:3 k(3*(l1(i)-1)+h,3*(l1(i)-1)+j) = k(3*(l1(i)-1)+h,3*(l1(i)-1)+j) + k_elem(h,j,i); k(3*(l1(i)-1)+h,3*(l2(i)-1)+j) = k(3*(l1(i)-1)+h,3*(l2(i)-1)+j) - k_elem(h,j,i); k(3*(l2(i)-1)+h,3*(l1(i)-1)+j) = k(3*(l2(i)-1)+h,3*(l1(i)-1)+j) - k_elem(h,j,i); k(3*(l2(i)-1)+h,3*(l2(i)-1)+j) = k(3*(l2(i)-1)+h,3*(l2(i)-1)+j) + k_elem(h,j,i); end end end

This for loop removes the rows and columns associated with known displacements. It results in the reduced stiffness matrix.

%remove rows and columns corresponding to known displacements removed = 0; k_red = k;  for i=1:length(Q_in(:,1)) k_red(Q_in(i,1)-removed,:)=[]; k_red(:,Q_in(i,1)-removed)=[]; removed = removed + 1; end This computes the unknown displacements from the reduced k matrix.

%Get unknown displacements Q_out=[F_in(:,1) inv(k_red)*F_in(:,2)]; This concatenates the original known displacements with those newly computed and then sorts them.

%Construct (and order) full displacement matrix Q=[Q_out; Q_in]; sort_done=0; temp = zeros(1); while(sort_done==0) sort_done=1; for i=1:length(Q(:,1))-1 if(Q(i,1)>Q(i+1,1)) temp=Q(i,:); Q(i,:)=Q(i+1,:); Q(i+1,:)=temp; sort_done=0; end end end This computes the unknown forces using the global stiffness matrix and the now-known displacement matrix.

%Get Force Matrix F = k*Q(:,2); This for loop translates the global displacements back to local coordinates and then computes the internal forces for each element.

%Get Member forces %Get local displacements Q_bar=zeros(2); for i=1:E T = [l(i) m(i) n(i) 0   0    0 0   0    0    l(i) m(i) n(i)]; Q_bar(:,i)=T*transpose([Q(3*l1(i)-2, 2) Q(3*l1(i)-1, 2) Q(3*l1(i), 2) Q(3*l2(i)-2, 2) Q(3*l2(i)-1, 2) Q(3*l2(i), 2)]); %Get member forces k_bar= (Modulus(i)*Area(i))/L(i)*[1 -1 -1 1];       F_bar(:,i)=k_bar*Q_bar(:,i); end This set of loops removes rows and columns associated with extra degrees of freedom.

%remove rows of zeros if problem is 1D or 2D DOF = 3; %assume 3 degrees of freedom if(x==0) for i=1:G k(DOF*i-(1+i),:)=[]; k(:,DOF*i-(1+i))=[]; end DOF = DOF - 1; end if(y==0) for i=1:G k(DOF*i-i,:)=[]; k(:,DOF*i-i)=[]; end DOF = DOF - 1; end if(z==0) for i=1:G k(DOF*i-(i-1),:)=[]; k(:,DOF*i-(i-1))=[]; end DOF = DOF - 1; end

Statics (Hand Calculation)
Begin by solving for the reaction force at global node 4 using a moment balance about global node 6.

$$\Sigma M_{6}=-24R_{2y}-9R_{3x}+36R_{4y}$$

$$R_{4y}=(24*1200+9*400)/36=900$$

Sequentially move from node to node solving the force balance equations at each node.

The following table shows the equivalence between the reaction forces computed by Matlab and those computed from the Statics approach.

The following table shows the equivalence between the member forces computed by Matlab and those computed from the Statics approach.

= Problem 2.6 =

Exercise 2, p. 98 Kim and Sankar using Matlab and CALFEM
Spring-rigid-body arrangement as shown. Global nodes 1 and 5 are supports. A force of 1000 N is applied to global node 3 in the negative x direction. find the displacements of global nodes 2, 3, and 4, the forces in the springs, and the reactions at the walls. Motion is constrained to the horizontal direction. Spring constants for springs 1-6 are given as 500, 400, 600, 200, 400, and 300 N/mm respectively.



Solution
Enter the number of global nodes. >> G=5;

Enter the x, y, and z coordinates (this program was designed solve 1D, 2D, and 3D problems, so the y and z coordinates must be entered. As long as all nodes have the same y and z coordinates, the solution will match a 1D problem). The exact x coordinates for this problem are not important. As long as the relative location of each nodes remains the same, the solution will be correct. >> x=[0 1 2 3 4]; >> y=[0 0 0 0 0]; >> z=[0 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=[4 0; 7 -1000; 10 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. Only the known displacements are entered into this matrix. If the force for a particular degree of freedom is unknown, the its displacement is known. >> Q_in=[1 0; 2 0; 3 0; 5 0; 6 0; 8 0; 9 0; 11 0; 12 0; 13 0; 14 0; 15 0];

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

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 2 2 1 3 4];

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

Enter the spring constant associated with each element. The column number corresponds to the element number of the given node. >> spring_const=[500 400 600 200 400 300];

Run the displacement program. >> [k Q F Q_bar F_bar]=displacement(G, x, y, z, F_in, Q_in, E, l1, l2, spring_const);

The stiffness matrix k: k = 700       -500        -200           0           0        -500        1500        -600        -400           0        -200        -600        1200        -400           0           0        -400        -400        1100        -300           0           0           0        -300         300

The displacement matrix Q: Q = 1.0000        0    2.0000         0    3.0000         0    4.0000   -0.8542    5.0000         0    6.0000         0    7.0000   -1.5521    8.0000         0    9.0000         0   10.0000   -0.8750   11.0000         0   12.0000         0   13.0000         0   14.0000         0   15.0000         0

This shows that the first degree of freedom (x direction) for global node 2 experiences a negative displacement of .8542 mm, the first degree of freedom (x direction) for global node 3 experiences a negative displacement of 1.5521 mm, and the first degree of freedom (x direction) of global node 4 experiences a negative displacement of .875 mm.

The Force matrix, F F = 1.0e+03 * 0.7375        0         0   -0.0000         0         0   -1.0000         0         0         0         0         0    0.2625         0         0

This shows that the first degree of freedom (x direction) for global node 1 exerts a positive force of 737.5 N and the first degree of freedom (x direction) for global node 5 exerts a positive force of 262.5 N.

The internal spring force matrix F_bar F_bar = 427.0833   8.3333  418.7500  310.4167 -270.8333 -262.5000 -427.0833   -8.3333 -418.7500 -310.4167  270.8333  262.5000

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.

This declares the function and a few variables.

function [k, Q, F, Q_bar, F_bar] = displacement(G, x, y, z, F_in, Q_in, E, l1, l2, spring_const) %Element Computation L = zeros(1,E); l = zeros(1,E); m = zeros(1,E); n = zeros(1,E); k_elem = zeros(3,3,E); This for loop computes the length, direction cosines, and the top left quarter of the global k matrix for each element.

for i=1:E L(i)=sqrt((x(l2(i))-x(l1(i)))^2+(y(l2(i))-y(l1(i)))^2+(z(l2(i))-z(l1(i)))^2); l(i)=(x(l2(i))-x(l1(i)))/L(i); m(i)=(y(l2(i))-y(l1(i)))/L(i); n(i)=(z(l2(i))-z(l1(i)))/L(i); k_elem(:,:,i) = spring_const(i)*[l(i)*l(i), l(i)*m(i), l(i)*n(i); m(i)*l(i), m(i)*m(i), m(i)*n(i); n(i)*l(i), n(i)*m(i), n(i)*n(i)]; end This nested for loop constructs the stiffness matrix.

%Stiffness Matrix Construction k = zeros(3*G,3*G); for i=1:E for j=1:3 for h=1:3 k(3*(l1(i)-1)+h,3*(l1(i)-1)+j) = k(3*(l1(i)-1)+h,3*(l1(i)-1)+j) + k_elem(h,j,i); k(3*(l1(i)-1)+h,3*(l2(i)-1)+j) = k(3*(l1(i)-1)+h,3*(l2(i)-1)+j) - k_elem(h,j,i); k(3*(l2(i)-1)+h,3*(l1(i)-1)+j) = k(3*(l2(i)-1)+h,3*(l1(i)-1)+j) - k_elem(h,j,i); k(3*(l2(i)-1)+h,3*(l2(i)-1)+j) = k(3*(l2(i)-1)+h,3*(l2(i)-1)+j) + k_elem(h,j,i); end end end This for loop removes the rows and columns associated with known displacements. It results in the reduced stiffness matrix.

%remove rows and columns corresponding to known displacements removed = 0; k_red = k;  for i=1:length(Q_in(:,1)) k_red(Q_in(i,1)-removed,:)=[]; k_red(:,Q_in(i,1)-removed)=[]; removed = removed + 1; end This computes the unknown displacements from the reduced k matrix.

%Get unknown displacements Q_out=[F_in(:,1) inv(k_red)*F_in(:,2)]; This concatenates the original known displacements with those newly computed and then sorts them. %Construct (and order) full displacement matrix Q=[Q_out; Q_in]; sort_done=0; temp = zeros(1); while(sort_done==0) sort_done=1; for i=1:length(Q(:,1))-1 if(Q(i,1)>Q(i+1,1)) temp=Q(i,:); Q(i,:)=Q(i+1,:); Q(i+1,:)=temp; sort_done=0; end end end This computes the unknown forces using the global stiffness matrix and the now-known displacement matrix.

%Get Force Matrix F = k*Q(:,2); This for loop translates the global displacements back to local coordinates and then computes the internal forces for each element.

%Get Member forces %Get local displacements Q_bar=zeros(2); for i=1:E T = [l(i) m(i) n(i) 0   0    0 0   0    0    l(i) m(i) n(i)]; Q_bar(:,i)=T*transpose([Q(3*l1(i)-2, 2) Q(3*l1(i)-1, 2) Q(3*l1(i), 2) Q(3*l2(i)-2, 2) Q(3*l2(i)-1, 2) Q(3*l2(i), 2)]); %Get member forces k_bar= spring_const(i)*[1 -1 -1 1];       F_bar(:,i)=k_bar*Q_bar(:,i); end This set of loops removes rows and columns associated with extra degrees of freedom.

%remove rows of zeros if problem is 1D or 2D DOF = 3; %assume 3 degrees of freedom if(x==0) for i=1:G k(DOF*i-(1+i),:)=[]; k(:,DOF*i-(1+i))=[]; end DOF = DOF - 1; end if(y==0) for i=1:G k(DOF*i-i,:)=[]; k(:,DOF*i-i)=[]; end DOF = DOF - 1; end if(z==0) for i=1:G k(DOF*i-(i-1),:)=[]; k(:,DOF*i-(i-1))=[]; end DOF = DOF - 1; end

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 2 2 4           3 2 3            4 1 3            5 3 4            6 4 5];

Construct an empty stiffness matrix K that has side lengths equal to the number of dimensions times the number of global nodes. For this problem, motion is constrained to 1 dimension and there are 5 global nodes. >> K=zeros(5,5);

Construct the force matrix F. The forces corresponding to global nodes 2, 3, and 4 are known; the forces corresponding to global nodes 1 and 5 are unknown. >> F=zeros(5,1); F(3)=-1000;

Construct the displacement matrix Q. Global nodes 1 and 5 have known displacements of 0. >> Q=[1 0 5 0];

The function spring1e generates the element stiffness matrices. Its parameter is the spring stiffness for each respective element. The stiffness are given in the problem statement.

The commands: >> Ke2=spring1e(200); >> Ke3=spring1e(300); >> Ke4=spring1e(400); >> Ke5=spring1e(500); >> Ke6=spring1e(600);

Yield the following stiffness matrices: >> Ke2=[200 -200 -200 200];   >> Ke3=[300 -300 -300 300];   >> Ke4=[400 -400 -400 400];   >> Ke5=[500 -500 -500 500];   >> Ke6=[600 -600 -600 600];

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

This results in the following stiffness matrix. >> K = 700       -500        -200           0           0        -500        1500        -600        -400           0        -200        -600        1200        -400           0           0        -400        -400        1100        -300           0           0           0        -300         300

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     -0.8542      -1.5521      -0.8750            0   r = 737.5000     -0.0000       0.0000       0.0000     262.5000

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

The function extract is used to computer the displacement of each spring, which is needed to determine the force in each spring. >> 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);

This results in the following displacements. The number following 'ed' corresponds to the element number. The first value represents the displacement of the left side of each spring. The second column corresponds to the displacement of the right side of each spring. ed1 =      0   -0.8542 ed2 = -0.8542  -0.8750 ed3 = -0.8542  -1.5521 ed4 =      0   -1.5521 ed5 = -1.5521  -0.8750 ed6 = -0.8750        0

The function spring1s can be used to determine the force in each spring. Its input parameters are the spring stiffness and the spring displacements. >> es1=spring1s(500,ed1) >> es2=spring1s(400,ed2) >> es3=spring1s(600,ed3) >> es4=spring1s(200,ed4) >> es5=spring1s(400,ed5) >> es6=spring1s(300,ed6)

This results in the following spring forces. es1 = -427.0833 es2 =  -8.3333 es3 = -418.7500 es4 = -310.4167 es5 = 270.8333 es6 = 262.5000

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