View source for University of Florida/Eml4507/s13 Team 3 Report 1

{{Big3|Report 1}}


==Problem 1: Completing a Matlab Primer==
[[File:Matlab_Primer.pdf|Matlab Primer]]

==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.

[[File:FEAHW1.JPG|framed|Spring-mass-dashpot system<ref>http://upload.wikimedia.org/wikiversity/en/b/b1/Fead.s13.sec53a.djvu FEAD: sec.53a, p.53-2b, Pb-53.1]</ref>]]
  
[[File:FEAHW2.JPG|framed|Spring-dashpot-mass system<ref>http://upload.wikimedia.org/wikiversity/en/d/d5/Iea.s12.sec1.djvu From IEA S12:p.1-4]</ref>]]

===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.

[[File:FEAHW1FBD.jpg|FEAHW1FBD]]

The force of a spring oriented in the y-direction is given by Hooke's Law:
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle F_{k}= ky</math></div>
|}

The force of a damper oriented in the y-direction is:
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle F_{c}=cy'</math></div>
|}



Summing the forces with respect to their directions in the Free Body Diagram:
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle -F_{k}- F_{c}- my''+ F(t)=0 </math></div>
|}



Plugging in spring and damper values: 
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle -ky -cy'- my''+ F(t)=0 </math></div>
|}


Reorganizing   the equation we get:<math>\displaystyle </math>
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math style="padding: 5px; border: 2px solid gray;">\displaystyle my''+cy'+ky= F(t)</math></div>
|}

===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.

[[File:FEAHW2fbd.jpg|FEAHW2fbd]]

The forces directions with respect to the Free Body Diagram: 
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle my'' + F_{1} = F(t)</math></div>

|}


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

:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle F_{k}=F_{c} </math></div>
|}


:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle ky_{k}=Cy'_{c} </math></div>
|}


Solving for dampers velocity:
Summing the forces with respect to their directions in the Free Body Diagram:
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle y''_{c}= y_{k}\frac{k}{c}</math></div>
|}



Plugging into acceleration equations after taking the derivative: 
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle y'' = y''_{c} + y''_{k} </math></div>
|}
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math>\displaystyle y'' = y''_{k} + y'_{k}\frac{k}{c} </math></div>
|}


Plugging into the original:<math>\displaystyle </math>
:{| style="width:100%" border="0"
|-
| style="width:95%"  |
<div style="text-align: left;"><math style="padding: 5px; border: 2px solid gray;">\displaystyle m(y''_{k} + y_{k}\frac{k}{c}) + ky_{k} = F(t)</math></div>
|}

== Problem 3: Finding the numerical values for the element stiff matrices k<sup>(1)</sup> and k<sup>(2)</sup> ==
[[File:R1.3 FEAD s13.PNG|framed|Plane truss with two members]]
=== Problem Statement ===

Part 1: Provide the numerical values for the element stiffness matrices <math> \mathbf{k^{(1)}} </math> and <math> \mathbf{k^{(2)}} </math>.

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:

<math>
E^{(1)} = 3, E^{(2)} = 5,
</math>

<math>
A^{(1)} = 1, A^{(2)} = 2,
</math>

<math>
L^{(1)} = 4, L^{(2)} = 2
</math>

<math>
\theta^{(1)}=30^{o}, \theta^{(2)}=-45^{o}
</math>

=== Problem Solution ===
==== Part 1 ====
<math>
\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}
</math>

where 

<math>
l^{(e)} = cos(\theta^{(e)}), m^{(e)} = sin(\theta^{(e)})
</math>

and

<math>
k^{(e)}={\frac{E^{(e)}A^{(e)}}{L^{(e)}}}
</math>

As given in the problem statement,

<math>
E^{(1)} = 3, E^{(2)} = 5,
</math>

<math>
A^{(1)} = 1, A^{(2)} = 2,
</math>

<math>
L^{(1)} = 4, L^{(2)} = 2
</math>

<math>
\theta^{(1)}=30^{o}, \theta^{(2)}=-45^{o}
</math>

Therefore,

<math>
k^{(1)}={\frac{3*1}{4}}={\frac{3}{4}}, k^{(2)}={\frac{5*2}{2}}=5,
</math>

<math>
(l^{(1)})^{2}=.75, (m^{(1)})^{2}=.25, l^{(1)}m^{(1)}=.433
</math>

and

<math>
(l^{(2)})^{2}=.5, (m^{(2)})^{2}=.5, l^{(2)}m^{(2)}=-.5
</math>

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

<math>
\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}
</math>

and

<math>
\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}
</math>

==== 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. <math> \theta^{(2)} </math> remains the same, but <math> \theta^{(1)}=210^{o} </math>.

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

<math> (x_{1},y_{1}) </math> and <math> (x_{2},y_{2}) </math> 

can be written as 

<math> (L^{(1)}cos(\theta^{(1)}),L^{(1)}sin(\theta^{(1)})) </math> and <math> (L^{(2)}cos(\theta^{(2)}),L^{(2)}sin(\theta^{(2)})) </math> respectively.

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

<math> (-2\sqrt{3},-2) </math> and <math> (2\sqrt{2},-2\sqrt{2}) </math>

==Problem 4: Analysis of a truss with three members in CALFEM==
   On my honor, I have neither given nor received unauthorized aid in doing this assignment. 
[[File:R14.PNG|framed|Plane truss with three members]]

===Problem Statement===
Consider a plane truss consisting of tree bars with the properties E = 200 GPa, A<sub>1</sub> = 6.0 · 10<sup>-4</sup> m<sup>2</sup>, A<sub>2</sub> = 3.0 · 10<sup>-4</sup> m<sup>2</sup> and A<sub>3</sub> = 10.0 · 10<sup>-4</sup> m<sup>2</sup>, and loaded by a single force P = 80 kN. The corresponding finite element model consists of three elements and eight degrees of freedom. 
<ref>[http://www.solid.lth.se/fileadmin/hallfasthetslara/utbildning/kurser/FHL064_FEM/calfem34.pdf Analysis of a plane truss, CALFEM: A Finite Element Toolbox, Version 3.4, p.9.3-9 (pdf p.234)]</ref>

===Solution===
The free body diagram is defined by the matrix 
   Edof=[1 1 2 5 6; 
         2 5 6 7 8; 
         3 3 4 5 6];
Where the first number in each row is the element, and the four other numbers are the degrees of freedom at each node of the element.

Next, create the n x n stiffness matrix '''K''', and the n x 1 load vector '''f'''
   >> K=zeros(8);
      f=zeros(8,1); f(6)=-80e3;

The axial rigidity of the elements are defined by the element property vectors
   >> E=2.0e11;
   >> A1=6.0e-4;    A2=3.0e-4;   A3=10.0e-4;
   >> ep1=[E A1];   ep2=[E A2];   ep3=[E A3];

and the element coordinate vectors '''ex1''', '''ex2''', '''ex3''', '''ey1''', '''ey2''', '''ey3'''
   >> ex1=[0 1.6];   ex2=[1.6 1.6];   ex3=[0 1.6];
   >> ey1=[0 0];   ey2[0 1.2];   ey3=[1.2 0];

Using '''bar2e''', compute '''Ke1''', '''Ke2''', and '''Ke3'''
   >> Ke1=bar2e(ex1,ey1,ep1)
   Ke1 = 
   1.0e+007 * 
   7.5000    0   -7.5000    0
        0    0         0    0
  -7.5000    0    7.5000    0
        0    0         0    0
   
   >> Ke2=bar2e(ex2,ey2,ep2)
   Ke2 = 
   1.0e+007 * 
   0         0    0         0
   0    5.0000    0   -5.0000
   0   -5.0000    0    5.0000
   
   >> Ke3=bar2e(ex3,ex3,ep3)
   Ke3 = 
   1.0e+007 *
   6.4000   -4.8000   -6.4000   4.8000
  -4.8000    3.6000    4.8000  -3.6000
  -6.4000    4.8000    6.4000  -4.8000
   4.8000   -3.6000   -4.8000   3.6000

Assemble the element stiffness matrices to create the global stiffness matrix

   >> K=assem(Edof(1,:),K,Ke1);
   >> K=assem(Edof(2,:),K,Ke2);
   >> K=assem(Edof(3,:),K,Ke3)
   
   K = 
   1.0e+008 *
   0.7500    0        0        0   -0.7500         0    0         0
        0    0        0        0         0         0    0         0
        0    0   0.6400   -4.800   -0.6400    0.4800    0         0
        0    0        0        0         0         0    0         0
        0    0  -0.6400    4.800    1.3900   -0.4800    0         0
        0    0   0.4800   -3.600   -0.4800    0.8600    0   -0.5000
        0    0        0        0         0         0    0         0
        0    0        0        0         0   -0.5000    0    0.5000

Create a prescribed displacement matrix '''bc'''
   >> bc=[1 0;2 0;3 0;4 0;7 0;8 0];
Then solve the system of equations using '''solveq''', which results in displacements '''a''' and support forces '''r'''
   >> [a,r]=solveq(K,f,bc)
   a =
   1.0e-002 *
         0
         0
         0
         0
   -0.0398
   -0.1152
         0
         0
   
   r =
   1.0e+004 *
    2.9845
         0
   -2.9845
    2.2383
    0.0000
    0.0000
         0
    5.7617

From the displacement matrix, it can be seen that the vertical displacement at the point of loading is 1.15 mm. Obtain the element displacements '''ed1''', '''ed2''', and '''ed3''' using '''extract'''. 
   >> ed1=extract(Edof(1,:),a);
   >> ed2=extract(Edof(2,:),a);
   >> ed3=extract(Edof(3,:),a);
Then calculate the normal forces '''N1''', '''N2''', and '''N3''' using '''bar2s'''.
   >> N1=bar2s(ex1,ey1,ep1,ed1)
   N1 =
   -2.9845e+004
   
   >> N2=bar2s(ex2,ey2,ep2,ed2)
   N2 =
   5.7617e+004
      
   >> N3=bar2s(ex3,ey3,ep3,ed3)
   N3 =
   3.7306e+004
The normal forces are '''N<sub>1</sub> = −29.84 kN''', '''N<sub>2</sub> = 57.62 kN''' and '''N<sub>3</sub> = 37.31 kN'''.

==Problem 5: Derivation of standard L2-ODE-CC==

===Problem Statement===
<math>\lambda_1 = -0.5 \lambda_2 = -1.5</math><br>
<math>\lambda_1 = \lambda_2 = -0.5</math><br>
<math>\lambda_1,_2 = -0.5+i2</math><br>

For each case find the standard L2-ODE-CC and the homogeneous solution. Plot each individual solution, and plot all solutions together
for comparison<br>

===Case 1===
<math> (\lambda-\lambda_1)(\lambda-\lambda_2)=0</math><br>
<math> (\lambda+.5)(\lambda+1.5)=0</math><br>
<math> \lambda^2 + 2\lambda + 0.75 = 0</math><br>
<br>
The standard L2-ODE-CC is then<br>
<math>y'' + 2y' +.75 = 0 </math> <br>
<br>
Homogeneous solution<br>
<math>y_h = c_1e^{-.5x} + c_2e^{-1.5x}</math><br>
<math>y'_h = -.5c_1e^{-.5x} - 1.5c_2e^{-1.5x}</math><br>
Plug in initial values<br>
<math> y(0)= 1 = c_1+c_2</math><br>
<math> y'(0) = 0 = -.5c_1 - 1.5c_2</math><br>
<br>
<math> c_1 = 1-c_2 </math><br>
<math> 0 = -.5(1-c_2)-1.5c_2 </math><br>
<math> .5 = -1c_2 </math><br>
<math> c_2 = -.5 </math><br>
<math> c_1 = 1.5 </math><br>
<br>
The solution is then<br>
<math>y_h = 1.5e^{-.5x} - .5e^{-1.5x}</math>

===Case 2===
<math> (\lambda-\lambda_1)(\lambda-\lambda_2)</math><br>
<math> (\lambda+.5)(\lambda+.5)</math><br>
<math> \lambda^2 + 1\lambda + 0.25 = 0</math><br>
<br>
The standard L2-ODE-CC is then<br>
<math>y'' + y' +.25 = 0 </math> <br>

Homogeneous solution<br>
<math>y_h = c_1e^{-.5x} + c_2xe^{-.5x}</math><br>
<math>y'_h = -.5c_1e^{-.5x} + c_2(e^{-.5x} - .5xe^{-.5x}</math><br>
Plug in initial values<br>
<math> y(0)= 1 = c_1</math><br>
<math> y'(0) = 0 = -.5c_1 + c_2</math><br>
<br>
<math> c_1 = 1 </math><br>
<math> 0 = -.5c_1 + c_2 </math><br>
<math> 0 = -.5 + c_2 </math><br>
<math> c_2 = .5 </math><br>
<br>
The solution is then<br>
<math>y_h = e^{-.5x} - .5xe^{-.5x}</math>

===Case 3===
<math> (\lambda-\lambda_1)(\lambda-\lambda_2)</math><br>
<math> (\lambda+.5+i2)(\lambda+.5-i2)</math><br>
<math> \lambda^2 + .5\lambda - i2\lambda + .5\lambda + .25 + .5i2 +i2\lambda - .5i2 + 4</math><br>
<math> \lambda^2 + 1\lambda + 4.25 = 0</math><br>
<br>
The standard L2-ODE-CC is then<br>
<math>y'' + y' + 4.25y = 0 </math> <br>

Homogeneous solution<br>
<math>y_h = e^{-.5x}(c_1cos(2x) + c_2sin(2x))</math><br>
<math>y'_h = -.5e^{-.5x}(c_1cos(2x) + c_2sin(2x)) + e^{-.5x}(-2c_1sin(2x) + 2c_2cos(2x)) </math><br>
Plug in initial values<br>
<math> y(0)= 1 = 1*(c_1 + 0)</math><br>
<math> y'(0) = 0 = -.5(c_1 + 0) + 1*(0 + 2c_2)</math><br>
<br>
<math> c_1 = 1 </math><br>
<math> 0 = -.5c_1 + 2c_2 </math><br>
<math> .5 = 2c_2 </math><br>
<math> c_2 = .25 </math><br>
<br>
The solution is then<br>
<math>y_h = e^{-.5x}(cos(2x) + .25sin(2x)) </math>

==Problem 6: Computing the element stiffness matrices for elements in a truss system==

=== Matlab Code ===

=== Define Input/Output Variables ===
The program begins with the definition of the input and output variables and then clears any existing variables.
 function [K] = stiffness( G )
 
 clear


=== Global Node Information ===
Then the number of global nodes in the problem is input into the variable "G". An array "g" is then created of size "G" to stor the numbers of each node. Then three arrays "x", "y" and "z" are created of size "G" to store the value of each of the corresponding coordinates of each global node. The user is prompted to enter this information using a for loop that loops "G" times.
 %Global Node Information
 G = input('Give number of Global Nodes: ');
 
 g = [1:G];
 x = zeros(1,G);
 y = zeros(1,G);
 z = zeros(1,G);
 
 for i=1:G
     fprintf('\nFor Global Node %g \n',i) 
     x(i)=input('x coordinate: ');
     y(i)=input('y coordinate: ');
     z(i)=input('z coordinate: ');
 end


=== Element Information ===
Here the user is prompted to enter the number of elements in the system into variable "E." Then the user is asked if the Young's Modulus and Area is the same for each element in the problem. Afterwards five arrays are constructed; "e" stores the element number information, "l1" stores the global node corresponding to local node one for each element, "l2" stores the global node corresponding to local node two for each element, "Modulus" stores the Young's Modulus for each element and "Area" stores the Area for each element.
 %Element Information
 E = input('\nGive number of Elements: ');
 same_modulus = input('Does each element have the same Modulus (1-yes, 0-no): ');
 same_area = input('Does each element have the same Area (1-yes, 0-no): ');
 
 e = [1:E];
 l1 = zeros(1,E);
 l2 = zeros(1,E);
 Modulus = zeros(1,E);
 Area = zeros(1,E);


=== Modulus, Area and Local Node Information ===
Here the program enters a ''if'' statement if "same_modulus" is equal to one, then the user is prompted to enter the Young's Modulus for all of the elements. Afterward a ''if'' statement is used again for the Area of each element. If "same_modulus" and "same_area" are equal to zero the program enters a ''for'' loop that prompts the user to enter the Young's Modulus, area and local node information of each element.
 %Receive element modulus and area if it is the same for each element
 if(same_modulus==1)
     Modulus(:)=input('\nYoungs Modulus: ');
 end
 if(same_area==1)
     Area(:)=input('\nArea: ');
 end
 
 for i=1:E
     fprintf('\nFor Element %g \n',i) 
     l1(i)=input('Global Node Number for Local Node 1: ');
     l2(i)=input('Global Node Number for Local Node 2: ');
     if(same_modulus~=1)
         Modulus(i)=input('Youngs Modulus: ');
     end
     if(same_area~=1)
         Area(i)=input('Area: ');
     end
 end

=== Element Computation ===
For the computation of the element matrix five arrays have to be defined; "L" is the length of each element, "l" is the director cosine of each element, "m" is the director sine of each element, "n" is, "k_elem" is a 3 by 3 stiffness matrix that can represent the entire 6 by 6 stiffness matrix for each element.  
 %Element Computation
 L = zeros(1,E);
 l = zeros(1,E);
 m = zeros(1,E);
 n = zeros(1,E);
 k_elem = zeros(3,3,E);
 
 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

=== Stiffness Matrix Construction ===
This section of code initializes the total stiffness matrix "k." Then several ''for'' loops are used to place the "k_elem" matrix into the total stiffness matrix. The code in these ''for'' loops places each "k_elem" matrix based on the global node and local node associated with it.  
 %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

=== 1D or 2D ===
Here the code finds rows and columns in the final stiffness matrix that contain all zeros and removes them from the total stiffness matrix.
 %remove rows of zeros if problem is 1D or 2D
 for i=1:G
     if(x(:)==0)
         k(3*i-(1+i),:)=[];
         k(:,3*i-(1+i))=[];
     end
     if(y(:)==0)
         k(3*i-i,:)=[];
         k(:,3*i-i)=[];
     end
     if(z(:)==0)
         k(3*i-(i-1),:)=[];
         k(:,3*i-(i-1))=[];       
     end
 end

== Member Forces ==

=== Eliminate Unknown Forces ===
The the rows 12, 11 and 8 correspond to the unknown forces so they are removed.
 k =
 
 1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0  -1.3404   -1.0053
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906  -1.0053   -0.7540
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0        0         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0        0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0        0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0        0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -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.6179         0         0         0         0         0    5.2358         0  -2.6179         0
          0   -3.4906         0         0         0         0         0         0         0    3.4906        0         0
    -1.3404   -1.0053         0         0         0         0         0         0   -2.6179         0   3.9583    1.0053
    -1.0053   -0.7540         0         0         0         0         0         0         0         0   1.0053    0.7540 
 
 >> k_temp=k
 
 k_temp =
 
    1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0  -1.3404   -1.0053
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906  -1.0053   -0.7540
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0        0         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0        0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0        0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0        0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -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.6179         0         0         0         0         0    5.2358         0  -2.6179         0
          0   -3.4906         0         0         0         0         0         0         0    3.4906        0         0
    -1.3404   -1.0053         0         0         0         0         0         0   -2.6179         0   3.9583    1.0053
    -1.0053   -0.7540         0         0         0         0         0         0         0         0   1.0053    0.7540 
 
 >> k_temp(12,:)=[]
 
 k_temp =
 
    1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0   -1.3404   -1.0053
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906   -1.0053   -0.7540
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0         0         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0         0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0         0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0         0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -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.6179         0         0         0         0         0    5.2358         0   -2.6179         0
          0   -3.4906         0         0         0         0         0         0         0    3.4906         0         0
    -1.3404   -1.0053         0         0         0         0         0         0   -2.6179         0    3.9583    1.0053
 
 >> k_temp(11,:)=[]
 
 k_temp =
 
    1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0   -1.3404   -1.0053
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906   -1.0053   -0.7540
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0         0         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0         0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0         0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0         0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -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.6179         0         0         0         0         0    5.2358         0   -2.6179         0
          0   -3.4906         0         0         0         0         0         0         0    3.4906         0         0
 
 
 >> k_temp(8,:)=[]
 
 k_temp =
 
    1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0   -1.3404   -1.0053
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906   -1.0053   -0.7540
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0         0         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0         0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0         0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0         0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -1.0053         0         0         0         0
          0         0   -2.6179         0         0         0         0         0    5.2358         0   -2.6179         0
          0   -3.4906         0         0         0         0         0         0         0    3.4906         0         0

=== Eliminate Known Displacements ===
The columns 12, 11 and 8 correspond to the known displacements so they are removed.
 >> k_temp(:,12)=[]
 
 k_temp =
 
   1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0   -1.3404
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906   -1.0053
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -1.0053         0         0         0
          0         0   -2.6179         0         0         0         0         0    5.2358         0   -2.6179
          0   -3.4906         0         0         0         0         0         0         0    3.4906         0
 
   
 >> k_temp(:,11)=[]
 
 k_temp =
 
    1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0         0
          0    4.9985    1.0053   -0.7540         0         0         0         0         0   -3.4906
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179         0   -2.6179         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404    1.0053         0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053   -0.7540         0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583   -1.0053         0         0
          0         0   -2.6179         0         0         0         0         0    5.2358         0
          0   -3.4906         0         0         0         0         0         0         0    3.4906
 
 >> k_temp(:,8)=[]
 
 k_temp =
 
    1.0e+03 *
 
     5.2987         0   -1.3404    1.0053   -2.6179         0         0         0         0
          0    4.9985    1.0053   -0.7540         0         0         0         0   -3.4906
    -1.3404    1.0053    6.5762   -1.0053         0         0   -2.6179   -2.6179         0
     1.0053   -0.7540   -1.0053    4.2445         0   -3.4906         0         0         0
    -2.6179         0         0         0    3.9583   -1.0053   -1.3404         0         0
          0         0         0   -3.4906   -1.0053    4.2445    1.0053         0         0
          0         0   -2.6179         0   -1.3404    1.0053    3.9583         0         0
          0         0   -2.6179         0         0         0         0    5.2358         0
          0   -3.4906         0         0         0         0         0         0    3.4906

=== Calculation of q ===
The matrix "q" is the displacements found using the reduced "k" matrix and the reduced "f" matrix. 
 >> f_temp=transpose([0 0 0 -1200 400 0 0 0 0])
 
 f_temp =
 
            0
            0
            0
        -1200
          400
            0
            0
            0
            0
 
 >> q_temp=inv(k_temp)*f_temp
 
 q_temp =
 
     0.8260
    -1.4993
     0.6112
    -2.1837
     0.5205
    -1.9258
     1.0696
     0.3056
    -1.4993
 
 >> q=transpose([transpose(q_temp(1:7)) 0 transpose(q_temp(8:9)) 0 0])
 
 q =
 
     0.8260
    -1.4993
     0.6112
    -2.1837
     0.5205
    -1.9258
     1.0696
          0
     0.3056
    -1.4993
          0
          0

=== Calculation of the Nodal Forces ===
 >> f=k*q
 
 f =
 
    1.0e+03 *
 
    -0.0000
     0.0000
     0.0000
    -1.2000
     0.4000
     0.0000
     0.0000
     0.9000
    -0.0000
    -0.0000
    -0.4000
     0.3000 
 
=== Calculation of the Member Forces === 		  
 >> q_bar=zeros(4,9);
 for i=1:E
     T= [l(i) m(i) 0    0
        -m(i) l(i) 0    0
         0    0    l(i) m(i)
         0    0   -m(i) l(i)];
     
     q_bar(:,i)=T*transpose([q(2*l1(i)-1) q(2*l1(i)) q(2*l2(i)) q(2*l2(i)-1)]);
 end
 >> q_bar
 
 q_bar =
 
    -1.7991    1.5719   -1.0696    1.4993    0.3056   -0.3056    0.8260         0    1.9258
     1.3802   -1.2284         0    0.8260   -1.4993    1.4993   -1.4993         0    0.5205
     1.6951   -0.6417    2.1837   -0.3056   -2.1837         0   -1.9258   -0.7038   -0.6112
     0.2387    0.8556   -0.6112   -1.4993    0.6112         0    0.5205    1.5604   -2.1837
 
 >> for i=1:E
     k_bar= (Modulus(i)*Area(i))/L(i)*[1  0 -1  0
                                       0  0  0  0
                                      -1  0  1  0
                                       0  0  0  0];
     
     f_bar(:,i)=k_bar*q_bar(:,i);
 end
 >> f_bar
 
 f_bar =
 
    1.0e+03 *
 
    -7.3180    4.6360   -8.5167    6.3000    6.5167   -0.8000    7.2042    1.4740    8.8556
          0         0         0         0         0         0         0         0         0
     7.3180   -4.6360    8.5167   -6.3000   -6.5167    0.8000   -7.2042   -1.4740   -8.8556
          0         0         0         0         0         0         0         0         0

=== Member Forces ===

Forces in Element 1:
'''Fx1 = -7318''',    
'''Fy1 =     0''',    
'''Fx2 =  4636''', 
'''Fy2 =     0''' 

Forces in Element 2:
'''Fx3 = -8516''', 
'''Fy3 =     0''', 
'''Fx4 =  6300''', 
'''Fy4 =     0''' 

Forces in Element 3:
'''Fx5 =  6516''', 
'''Fy5 =     0''',  
'''Fx6 =  -800''',  
'''Fy6 =     0'''  

Forces in Element 4:
'''Fx7 =  7204''',  
'''Fy7 =     0''',  
'''Fx8 =  1474''',  
'''Fy8 =     0'''  

Forces in Element 5:
'''Fx9 =  8855''',  
'''Fy9 =     0''',  
'''Fx10 =  7318''', 
'''Fy10 =     0''' 

Forces in Element 6:
'''Fx11 = -4636''', 
'''Fy11 =     0''', 
'''Fx12 =  8516''', 
'''Fy12 =     0''' 

Forces in Element 7:
'''Fx13 = -6300''', 
'''Fy13 =     0''', 
'''Fx14 = -6516''', 
'''Fy14 =     0''' 

Forces in Element 8:
'''Fx15 =   800''', 
'''Fy15 =     0''', 
'''Fx16 = -7204''', 
'''Fy16 =     0''' 

Forces in Element 9:
'''Fx17 = -1474''', 
'''Fy17 =     0''', 
'''Fx18 = -8855''', 
'''Fy18 =     0'''

== Full Matlab Code ==
<syntaxhighlight lang="matlab">
 function [k] = stiffness( G )
 
 clear
 
 %Global Node Information
 G = input('Give number of Global Nodes: ');
 
 g = [1:G];
 x = zeros(1,G);
 y = zeros(1,G);
 z = zeros(1,G);
 
 for i=1:G
     fprintf('\nFor Global Node %g \n',i) 
     x(i)=input('x coordinate: ');
     y(i)=input('y coordinate: ');
     z(i)=input('z coordinate: ');
 end
 
 %Element Information
 E = input('\nGive number of Elements: ');
 same_modulus = input('Does each element have the same Modulus (1-yes, 0-no): ');
 same_area = input('Does each element have the same Area (1-yes, 0-no): ');
 
 e = [1:E];
 l1 = zeros(1,E);
 l2 = zeros(1,E);
 Modulus = zeros(1,E);
 Area = zeros(1,E);
 
 %Receive element modulus and area if it is the same for each element
 if(same_modulus==1)
     Modulus(:)=input('\nYoung''s Modulus: ');
 end
 if(same_area==1)
     Area(:)=input('\nArea: ');
 end
 
 for i=1:E
     fprintf('\nFor Element %g \n',i) 
     l1(i)=input('Global Node Number for Local Node 1: ');
     l2(i)=input('Global Node Number for Local Node 2: ');
     if(same_modulus~=1)
         Modulus(i)=input('Young''s Modulus: ');
     end
     if(same_area~=1)
         Area(i)=input('Area: ');
     end
 end
 
 
 %Element Computation
 L = zeros(1,E);
 l = zeros(1,E);
 m = zeros(1,E);
 n = zeros(1,E);
 k_elem = zeros(3,3,E);
 
 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
 
 
 %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
 
 %remove rows of zeros if problem is 1D or 2D
 for i=1:G
     if(x(:)==0)
         k(3*i-(1+i),:)=[];
         k(:,3*i-(1+i))=[];
     end
     if(y(:)==0)
         k(3*i-i,:)=[];
         k(:,3*i-i)=[];
     end
     if(z(:)==0)
         k(3*i-(i-1),:)=[];
         k(:,3*i-(i-1))=[];       
     end
 end
</syntaxhighlight>

==References==
{{reflist}}

{{CourseCat}}