User:Eml4500.f08.bike.mcdonald/homework report 2

Global Picture
The first step in the recipe to solve problems using the finite element method is to create a global picture of the problem itself. This global picture must contain the Global Degrees of Freedom (DOFs), also referred to as displacement DOFs and Global Forces. In order to determine what needs to be solved for, the displacement DOFs are split up between those that are known, such as constraints, and those that are unknown. The unknown parts must be solved for using the Finite Element Method.

Similarly, the Global Forces are split up into a known part, which is comprised of the applied forces to the system, and an unknown part, which are the reactions of the system.

Element Picture
Once the Global Picture has been completed, it is necessary to create an Element Picture for each of the elements of the global system. These element pictures must label all of the Element DOFs and Element Forces.

Force Displacement Relations
The first step in creating the Global Force Displacement (FD) Relation is to populate the Element Stiffness and Force Matrices. Once the Element Stiffness Matrices have been created, the Global Stiffness Matrix can be assembled. This Global Stiffness Matrix can then be inserted into the Global FD Relation.

Elimination of Known Displacement DOFs
In order to solve the Global FD Relation created in the last step, it is necessary to eliminate the known displacement DOFs in the system. As stated above, these known DOFs usually come in the form of constraints on the system. Once these constraints have been imposed on the system, it is possible to fully determine the values of the Global Displacement Matrix using the Global FD Relation.

Compute Element Forces
Once the Global Displacement Matrix has been fully determined, it is necessary to relate these global displacements with the corresponding element displacements. The fully populated Element Displacement Matrices can then be used in the Element FD Relations to solve for the forces in each element.

Compute Reactions
Because the Element Force Matrices are composed of internal element forces and reactions, computing the reactions is as simple as determining the correspondence between the Global Force Matrix and Element Force Matrices.

Global Picture
The image below shows the Global Picture of the two-member truss system.

The main reason for the creation of the Global Picture is to label each node in the two-member truss system and introduce the Global DOFs and Global Forces.

$$

\textbf{K}_{n x n}=\left[ \begin{array}{c} K_{ij} \end{array} \right]_{n x n}=\textrm{Global\ Stiffness\ Matrix}

$$

$$

\textbf{d}_{n x 1}=\begin{Bmatrix} d_{j} \end{Bmatrix}_{n x 1}=\textrm{Global\ Displacement\ Matrix}

$$

$$

\textbf{F}_{n x 1}=\begin{Bmatrix} F_{i} \end{Bmatrix}_{n x 1}=\textrm{Global\ Force\ Matrix}

$$

Numbering the Global Displacement DOFs

 * Follow the order of the global node numbers
 * For each node, follow the order of the global coordinate axes and number the displacement definitions for that node

Numbering the Global Forces

 * Follow the order of the global node numbers
 * For each node, follow the order of the global coordinate axes and number the displacement definitions for that node

Element Picture
$$

\textbf{k}_{4 x 4}^{(e)}\textbf{d}_{4 x 1}^{(e)}=\textbf{f}_{4 x 1}^{(e)}

$$

$$

\textbf{k}_{4 x 4}^{(e)}=\left[ \begin{array}{c} k_{ij}^{(e)} \end{array} \right]_{4 x 4}=\textrm{Element\ Stiffness\ Matrix}

$$

$$

\textbf{d}_{4 x 1}^{(e)}=\begin{Bmatrix} d_{j}^{(e)} \end{Bmatrix}_{4 x 1}=\textrm{Element\ Displacement\ Matrix}

$$

$$

\textbf{f}_{4 x 1}^{(e)}=\begin{Bmatrix} f_{i}^{(e)} \end{Bmatrix}_{4 x 1}=\textrm{Element\ Force\ Matrix}

$$

The Element Picture for Element 1 (e=1) is shown below.

The Element Picture for Element 2 (e=2) is shown below.

Element Force Displacement Relations
The condensed version of the Element Stiffness Matrix is shown below for an element with four different DOFs.

$$

\left[ \begin{array}{c} \textbf{k}^{(e)} \end{array} \right]_{4 x 4}\begin{Bmatrix} \textbf{d}^{(e)} \end{Bmatrix}_{4 x 1}=\begin{Bmatrix} \textbf{f}^{(e)} \end{Bmatrix}_{4 x 1}

$$


 * At the element level, the element force matrix is equal to the element stiffness matrix for element e multiplied by the element displacement matrix of element e.

From page 225 in Finite Element Analysis with Mathematica and MATLAB Computations, The second equation from the bottom is the element equation when the matrix multiplication is carried out. This equation is shown below.

$$

k^{(e)} \left[ \begin{array}{cccc} (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{array} \right] \left[\begin{array}{c} u_{1} \\ v_{1} \\ u_{2} \\ v_{2}\end{array}\right] = \left[\begin{array}{c} F_{1x} \\ F_{1y} \\ F_{2x} \\ F_{2y} \end{array}\right]

$$

Where the axial stiffness of the bar element "e" is

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

When a problem with any plane assembly of bars needs to be solved, the above equation can be used. The individual elements, however, are assumed to resist axial loads.

$$

l^{(e)}\textrm{,\ }m^{(e)}=\textrm{director\ of\ cosines\ of\ the\ }\tilde {x} \textrm{\ axis\ (goes\ from\ the\ local\ node\ numbers\ 1\ to\ 2)\ with\ respect\ to\ global\ (x,y)\ coordinates.}

$$



The director of cosines are shown below as trigonometric functions.

$$ l^{(e)}=\mathbf{\vec{i}}\cdot \vec{i}=\vec{i}\cdot (\cos \theta ^{(e)}\vec{i +\sin \theta ^{(e)}}\vec{j)}=\cos \theta ^{(e)}\vec{i}\cdot \vec{i}+\sin \theta ^{(e)}\vec{j}\cdot \vec{i}=\cos \theta ^{(e)}\times 1 +\sin \theta ^{(e)}\times 0 $$

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

In general, $$k_{ij}^{(e)} = k_{ji}^{(e)}$$ or $$\textbf{k}^{(e)T} = \textbf{k}^{(e)}.$$ The lower triangle of the $$\textbf{k}^{(e)}$$ matrix is symmetric. The upper triangle part is the remaining part shown below. $$ \textbf{k}^{(e)}= \left[ \begin{array}{cccc} k_{11}^{(e)} & k_{12}^{(e)} & k_{13}^{(e)} & k_{14}^{(e)} \\ & k_{22}^{(e)} & k_{23}^{(e)} & k_{24}^{(e)} \\ & & k_{33}^{(e)} & k_{34}^{(e)} \\ & & & k_{44}^{(e)} \end{array} \right]$$

This is useful for finding the elemental stiffness matrices for any element. This allows for the reduction of the number of calculations required in order to fill in the matrices.

To assemble the elemental stiffness matrix for element 1, the director cosines and axial stiffness of the bar must be determined as shown below. The axial stiffness $$k^{(1)}$$ is found using the Elastic Modulus E, the cross-sectional area of the element A, and the length of the element L.

$$\theta^{(1)} = 30^{\circ}$$ $$l^{(1)} = \cos(\theta^{(1)}) = \cos(30^{\circ}) = \frac{\sqrt{3}}{2} $$ $$m^{(1)} = \sin(\theta^{(1)}) = \cos(30^{\circ}) = \frac{1}{2} $$ $$k^{(1)} = \frac{E^{(1)}A^{(1)}}{L^{(1)}} = \frac{(3)(1)}{4} = \frac{3}{4}$$

Once these have been determined they can be used to create the elemental stiffness matrix shown below in coefficient form.

$$ \textbf{k}^{(1)}= k^{(1)}\left[ \begin{array}{cccc} k_{11}^{(1)} & k_{12}^{(1)} & k_{13}^{(1)} & k_{14}^{(1)} \\ k_{21}^{(1)} & k_{22}^{(1)} & k_{23}^{(1)} & k_{24}^{(1)} \\ k_{31}^{(1)} & k_{32}^{(1)} & k_{33}^{(1)} & k_{34}^{(1)} \\ k_{41}^{(1)} & k_{42}^{(1)} & k_{43}^{(1)} & k_{44}^{(1)} \end{array} \right]$$ $$ = \left[ \begin{array}{c} k_{ij}^{(1)} \end{array} \right]_{4x4} $$ where $$i = 1, 2, 3, 4 $$ and $$ j = 1, 2, 3, 4$$

The coefficients can be determined using the previously found director cosines and the axial stiffness of the bar. Three sample calculations are shown below.

$$k_{11}^{(1)} = k^{(1)}(l^{(1)})^2 = (\frac{3}{4})(\frac{\sqrt{3}}{2})^2 = \frac{9}{16}$$ $$k_{12}^{(1)} = k^{(1)}(l^{(1)})(m^{(1)}) = (\frac{3}{4})(\frac{\sqrt{3}}{2})(\frac{1}{2}) = \frac{3\sqrt{3}}{16}$$ $$k_{22}^{(1)} = k^{(1)}(m^{(1)})^2 = (\frac{3}{4})(\frac{1}{2})^2 = \frac{3}{16}$$ Observations:

1) It is only necessary to compute 3 numbers: $$k_{11}^{(1)}, k_{12}^{(1)}, k_{22}^{(1)}$$. The other coefficients have the same absolute value but differ depending on whether they are positive or negative.  2) Matrix $$\textbf{k}^{(1)}$$ is symmetric; therefore, $$k_{ij}^{(1)} = k_{ji}^{(1)}$$. For example: $$k_{13}^{(1)} = k_{31}^{(1)}$$ You only need to interchange the row and column indices.

The director cosines and axial stiffness of the bar calculations have to be repeated for element 2. These calculations are shown below. $$\theta^{(2)} = -45^{\circ}$$ $$l^{(2)} = \cos(\theta^{(2)}) = \cos(-45^{\circ}) = \frac{\sqrt{2}}{2} $$ $$m^{(2)} = \sin(\theta^{(2)}) = \cos(-45^{\circ}) = \frac{-\sqrt{2}}{2} $$ $$k^{(2)} = \frac{E^{(2)}A^{(2)}}{L^{(2)}} = \frac{(5)(2)}{5} = 5$$

These are determined to find the coefficients of the stiffness matrix for element 2. That matrix is shown below in coefficient form.

$$ \textbf{k}^{(2)}= k^{(2)}\left[ \begin{array}{cccc} k_{11}^{(2)} & k_{12}^{(2)} & k_{13}^{(2)} & k_{14}^{(2)} \\ k_{21}^{(2)} & k_{22}^{(2)} & k_{23}^{(2)} & k_{24}^{(2)} \\ k_{31}^{(2)} & k_{32}^{(2)} & k_{33}^{(2)} & k_{34}^{(2)} \\ k_{41}^{(2)} & k_{42}^{(2)} & k_{43}^{(2)} & k_{44}^{(2)} \end{array} \right]$$ $$ = \left[ \begin{array}{c} k_{ij}^{(2)} \end{array} \right]_{4x4} $$ where $$i = 1, 2, 3, 4 $$ and $$ j = 1, 2, 3, 4$$

The only required calculation to determine the value of the matrix above is shown below.

$$k_{11}^{(2)} = k^{(2)}(l^{(2)})^2 = (5)(\frac{\sqrt{2}}{2})^2 = \frac{5}{2}$$ Observations: 1) The absolute value of all the coefficients in the $$\textbf{k}_{ij}^{(2)}$$ where (i,j) = 1, 2, 3, and 4 are the same. Compute one of the coefficients and all the other coefficients will have the same value but will differ depending on whether they are positive or negative.  2) It is also known that $$\textbf{k}^{(2)T} = \textbf{k}^{(2)}$$. This leads to the conclusion that $$\textbf{k}^{(2)}$$ is symmetric.

The element level displacement degrees of freedom are shown below.

$$

\textrm{Element\ 1:\ }\begin{Bmatrix} d_{1}^{(1)} & d_{2}^{(1)} & d_{3}^{(1)} & d_{4}^{(1)} \end{Bmatrix}

$$

$$

\textrm{Element\ 2:\ }\begin{Bmatrix} d_{1}^{(2)} & d_{2}^{(2)} & d_{3}^{(2)} & d_{4}^{(2)} \end{Bmatrix}

$$

Global Force Displacement Relation
The generalized Global Force Displacement Relation for a "Free-Free Structure" is shown below.

$$

\textbf{K}_{n x n}\textbf{d}_{n x 1}=\textbf{F}_{n x 1}

$$

Where $$\textbf{K}_{n x n}$$ is the Global Stiffness Matrix for the two-member truss system described above, $$\textbf{d}_{n x 1}$$ is the Global Displacement Matrix, and $$\textbf{F}_{n x 1}$$ is the Global Force Matrix.

The expanded matrix form of the Global Force-Displacement Relation for the "Free-Free Structure" above is shown below.

$$

\left[ \begin{array}{cccccc} K_{11} & K_{12} & K_{13} & K_{14} & K_{15} & K_{16} \\ K_{21} & K_{22} & K_{23} & K_{24} & K_{25} & K_{26} \\ K_{31} & K_{32} & K_{33} & K_{34} & K_{35} & K_{36} \\ K_{41} & K_{42} & K_{43} & K_{44} & K_{45} & K_{46} \\ K_{51} & K_{52} & K_{53} & K_{54} & K_{55} & K_{56} \\ K_{61} & K_{62} & K_{63} & K_{64} & K_{65} & K_{66} \end{array} \right]\begin{Bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \end{Bmatrix}=\begin{Bmatrix} F_{1} \\ F_{2} \\ F_{3} \\ F_{4} \\ F_{5} \\ F_{6} \end{Bmatrix}

$$

The global level displacement degrees of freedom are shown below.

$$

\begin{Bmatrix} d_{1} & d_{2} & d_{3} & d_{4} & d_{5} & d_{6} \end{Bmatrix}

$$

In order to go from the Element Matrices (Stiffness, Displacement, and Force) to the Global Matrices, it is necessary to go through an assembly process. This assembly process starts with identifying the correspondence between element displacement degrees of freedom and global displacement degrees of freedom as shown below.

$$

\textrm{Node\ 1:\ }d_{1}=d_{1}^{(1)}\textrm{,\ }d_{2}=d_{2}^{(1)}

$$

$$

\textrm{Node\ 2:\ }d_{3}=d_{3}^{(1)}=d_{1}^{(2)}\textrm{,\ }d_{4}=d_{4}^{(1)}=d_{2}^{(2)}

$$

$$

\textrm{Node\ 3:\ }d_{5}=d_{3}^{(2)}\textrm{,\ }d_{6}=d_{4}^{(2)}

$$

The conceptual step of the assembly process is shown below.



The section in the above matrix where the red shaded section and blue shaded section overlap is caused by the fact that $$d_{3}$$ and $$d_{4}$$ are common to both elements of the two-member truss system.

The above figure from the previous lecture shows the combination of the two local 4x4 k matrices into a 6x6 global k matrix. The top right corner is the k matrix corresponding to the first element and the bottom left corner corresponds to the k matrix of element 2. The 4x4 section in the middle of this global matrix is a combination of the 2 local matrices. The global matrix in terms of the local matrix elements is as follows:

$$

\left[ \begin{array}{cccccc} K^{(1)}_{11} & K^{(1)}_{12} & K^{(1)}_{13} & K^{(1)}_{14} & 0 & 0 \\ K^{(1)}_{21} & K^{(1)}_{22} & K^{(1)}_{23} & K^{(1)}_{24} & 0 & 0 \\ K^{(1)}_{31} & K^{(1)}_{32} & (K^{(1)}_{33} + K^{(2)}_{11} & (K^{(1)}_{34} + K^{(2)}_{12})& K^{(2)}_{13} & K^{(2)}_{14} \\ K^{(1)}_{41} & K^{(1)}_{42} & (K^{(1)}_{43} + K^{(2)}_{21}) & (K^{(1)}_{44} + K^{(2)}_{22}) & K^{(2)}_{23} & K^{(2)}_{24} \\ 0 & 0 & K^{(2)}_{31} & K^{(2)}_{32} & K^{(2)}_{33} & K^{(2)}_{34} \\ 0 & 0 & K^{(2)}_{41} & K^{(2)}_{42} & K^{(2)}_{43} & K^{(2)}_{44} \end{array} \right]\begin{Bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \end{Bmatrix}=\begin{Bmatrix} F_{1} \\ F_{2} \\ F_{3} \\ F_{4} \\ F_{5} \\ F_{6} \end{Bmatrix}

$$

Using the K values given earlier in class, some sample element calculations are as follows, with the rest able to be found using the symmetry of the matrix:

$$ K_{11} = \frac{9}{16}$$

$$ K_{12} = \frac{3\sqrt{3}}{16}$$

$$ K_{33} = \frac{49}{16}$$

$$ K_{34} = \textbf{-2.175}$$

$$ K_{43} = \textbf{-2.175}$$

$$ K_{44} = \textbf{2.6875}$$

The notation used above is for the global elements, rather than the local elements.

Elimination of Known Displacement DOFs
This matrix in compact notation is shown below.

$$

\left[ \begin{array}{c} K_{ij} \end{array} \right]_{6 x 6}\begin{Bmatrix} d_{j} \end{Bmatrix}_{6 x 1}=\begin{Bmatrix} F_{i} \end{Bmatrix}_{6 x 1}

$$

The summation notation form of the above compact matrix notation is shown below.

$$

\displaystyle\sum_{j=1}^6 K_{ij}d_{j}=F_{i}, i=1,...,6

$$

Also, due to the symmetry of the matrix, $$ \textbf{K}_{34} = \textbf{K}_{43}$$, therefore only one needs to be calculated when computing the values of the elements.

The green shaded top right and bottom left corners of the figure correspond to the $$ \textbf{0}$$ elements of the matrix.

After constructing this global matrix, the next step is to reduce the equation by utilizing the boundary conditions. In this particular problem, global nodes 1 and 2 are fixed, thus there is no displacement involved, and:

$$ \textbf{d} = \begin{Bmatrix} 0 \\ 0 \\ d_{3} \\ d_{4} \\ 0 \\ 0 \end{Bmatrix} $$

Utilizing the Principle of Virtual Work allows us to ignore the first, second, fifth and sixth column of the $$ \textbf{k}$$ matrix, as all $$ \textbf{F}$$ elements that it would produce would not be used in the calculation. Thus:

$$

\left[ \begin{array}{cc} K_{13} & K_{14} \\ K_{23} & K_{24} \\ K_{33} & K_{34} \\ K_{43} & K_{44} \\ K_{53} & K_{54}\\K_{63} & K_{64}\end{array} \right]\begin{Bmatrix} d_{3} \\ d_{4}\end{Bmatrix}=\begin{Bmatrix} F_{3} \\ F_{4}\end{Bmatrix}

$$

Step 4 in the finite element analysis process begins with the reduction of the global force-displacement relations. This reduction is done by first eliminating the known degrees of freedom(DOF’s). $$ d_{1}=d_{2}=d_{5}=d_{6}=0 $$

Deleting the corresponding columns in the global stiffness matrix and the force matrix results in the following relationship.

$$ \left[ \begin{array}{cc} k_{33} & k_{34}\\ k_{43}& k_{44} \end{array} \right]_{2x2} $$ $$\begin{Bmatrix} d_{3}\\d_{4} \end{Bmatrix}_{2x1}=\begin{Bmatrix}F_{3}\\F_{4} \end{Bmatrix}_{2x1}$$

In order to find the inverse of the stiffness matrix one must first calculate the determinant as follows. det $$k = k_{33}k_{44}-k_{34}k_{43} $$

The inverse of the stiffness matrix can then be calculated as shown below.

$$k^{-1} = \frac {1}{\text{det(k)}}$$$$ \left[ \begin{array}{cc} k_{44} & -k_{34}\\ -k_{43}& k_{33} \end{array} \right] $$

If done correctly the stiffness matrix multiplied by its inverse will yield the identity matrix, I. $$I= \left[ \begin{array}{cc} 1 & 0\\ 0& 1 \end{array} \right] $$

For the two beam truss system shown below, the displacement matrix is calculated by multiplying the inverse matrix by $$\begin{Bmatrix}{0} \\ {P} \end{Bmatrix}$$, where P=7.



$$\begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}=k^{-1}\begin{Bmatrix}{0}\\{7}\end{Bmatrix}$$

det $$k=(3.0625)(2.6875)-(-2.175)(-2.175)=3.50$$

Therefore, $$ k^{-1}=\frac {1}{3.50}$$$$ \left[ \begin{array}{cc}2.6875&2.175\\ 2.175&3.0625\end{array} \right]=\left[ \begin{array}{cc}.7678&.6214\\ .6214&.8750\end{array} \right]$$

so $$\begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}=\left[ \begin{array}{cc}.7678&.6214\\ .6214&.8750\end{array} \right]\begin{Bmatrix}{0}\\{7}\end{Bmatrix}=\begin{Bmatrix}{4.352}\\{6.1271}\end{Bmatrix}$$

Step 5 consists of computing the reactions, which can be done by two different methods. Method 1 is to use the element force-displacement relations and method 2 is the Euler cut principle. The force-displacement relation corresponding to the two beam truss system is shown below.

$$k_{4x4}^{(e)}d_{4x1}^{(e)}=f_{4x1}^{(e)}$$ where $$e=1,2$$

$$d^{(1)}=\begin{Bmatrix}{0}\\{0}\\{4.352}\\{6.1271}\end{Bmatrix}$$ and $$f^{(1)}=\begin{Bmatrix}{f_{1}^{(1)}}\\{f_{2}^{(1)}}\\{f_{3}^{(1)}}\\{f_{4}^{(1)}}\end{Bmatrix}$$

$$d^{(2)}=\begin{Bmatrix}{4.352}\\{6.1271}\\{0}\\{0}\end{Bmatrix}$$ and $$f^{(2)}=\begin{Bmatrix}{f_{1}^{(2)}}\\{f_{2}^{(2)}}\\{f_{3}^{(2)}}\\{f_{4}^{(2)}}\end{Bmatrix}$$

Compute Element Forces
The resulting element displacement and stiffness matrices can be used to determine the element forces. The forces on Element 1 are calculated below.

Element 1 :

The element stiffness matrix multiplied by the element displacement matrix is k (1) d (1)

$$

=\left[ \begin{array}{cccc} K_{11} & K_{12} & K_{13} & K_{14} \\ K_{21} & K_{22} & K_{23} & K_{24} \\ K_{31} & K_{32} & K_{33} & K_{34} \\ K_{41} & K_{42} & K_{43} & K_{44} \end{array} \right]\begin{Bmatrix} 0 \\ 0 \\ 4.352 \\ 6.127  \end{Bmatrix}

$$

The matrices can be reduced due to the zeros in the displacement matrix and set equal to the element force matrix k (1) d (1) = f (1)

$$

=\left[ \begin{array}{cc} K_{13} & K_{14} \\ K_{23} & K_{24} \\ K_{33} & K_{34} \\ K_{43} & K_{44} \end{array} \right]\begin{Bmatrix} 4.352 \\ 6.127  \end{Bmatrix} = \left[ \begin{array}{c} -4.4378 \\ -2.5622 \\ 4.4378 \\ 2.5622 \end{array} \right]  =  \left[ \begin{array}{c} f_{1}^{(1)} \\ f_{2}^{(1)} \\ f_{3}^{(1)} \\ f_{4}^{(1)} \end{array} \right]

$$

Observation:

Element 1 is in equilibrium therefore:

(1) $$\sum F_x = f_1^{(1)} + f_3^{(1)} = 0 $$

(2) $$\sum F_y = f_2^{(1)} + f_4^{(1)} = 0 $$

(3) $$\sum M_{any pt} = 0 $$



Element 2 :

The figure below shows element 2 labeled with axial forces P.



Method 2 to solve for a 2-bar truss system using statics method:



P can be brought back into the big picture by performing an equilibrium balance of node 2.



Compute Reactions
The final step in solving for the reactions is to relate the element forces found in the previous step with the global forces that correspond to the reactions of the system. The reactions of the system are shown below.

$$

\begin{Bmatrix} f_{1} \\ f_{2} \\ f_{5} \\ f_{6} \end{Bmatrix}=\begin{Bmatrix} -4.438 \\ -2.562 \\ 4.438 \\ -4.438 \end{Bmatrix}

$$

Analysis and Verification of MATLAB Code for Two-Member Truss System
In this section, the two-member truss problem covered during class will be solved manually using the Finite Element Method and compared with the execution of the Finite Element Method by means of a MATLAB code written by Professor Loc Vu Quoc from the University of Florida Department of Mechanical and Aerospace Engineering.

Finite Element Method Manual Calculations


The next step to solving the two-member truss system shown above is to populate the Element Stiffness Matrices, $$\textbf{k}^{(1)}$$ and $$\textbf{k}^{(2)}$$. This step is shown below using the equation for the Element Stiffness Matrices given in the course notes above.

$$

\textbf{k}^{(1)}=k^{(1)} \left[ \begin{array}{cccc} (l^{(1)})^{2} & l^{(1)}m^{(1)} & -(l^{(1)})^{2} & -l^{(1)}m^{(1)} \\ l^{(1)}m^{(1)} & (m^{(1)})^{2} & -l^{(1)}m^{(1)} & -(m^{(1)})^{2} \\ -(l^{(1)})^{2} & -l^{(1)}m^{(1)} & (l^{(1)})^{2} & l^{(1)}m^{(1)} \\ -l^{(1)}m^{(1)} & -(m^{(1)})^{2} & l^{(1)}m^{(1)} & (m^{(1)})^{2} \end{array}\right]=\left[\begin{array}{cccc} 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 \end{array} \right]

$$

$$

\textbf{k}^{(2)}=k^{(2)} \left[ \begin{array}{cccc} (l^{(2)})^{2} & l^{(2)}m^{(2)} & -(l^{(2)})^{2} & -l^{(2)}m^{(2)} \\ l^{(2)}m^{(2)} & (m^{(2)})^{2} & -l^{(2)}m^{(2)} & -(m^{(2)})^{2} \\ -(l^{(2)})^{2} & -l^{(2)}m^{(2)} & (l^{(2)})^{2} & l^{(2)}m^{(2)} \\ -l^{(2)}m^{(2)} & -(m^{(2)})^{2} & l^{(2)}m^{(2)} & (m^{(2)})^{2} \end{array}\right]=\left[\begin{array}{cccc} 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{array} \right]

$$

Once both Element Stiffness Matrices have been populated, it is necessary to assemble the Global Stiffness Matrix using an assembly process involving both $$\textbf{k}^{(1)}$$ and $$\textbf{k}^{(2)}$$. This process is shown below.

$$

\textbf{K}=\left[ \begin{array}{cccccc} K^{(1)}_{11} & K^{(1)}_{12} & K^{(1)}_{13} & K^{(1)}_{14} & 0 & 0 \\ K^{(1)}_{21} & K^{(1)}_{22} & K^{(1)}_{23} & K^{(1)}_{24} & 0 & 0 \\ K^{(1)}_{31} & K^{(1)}_{32} & (K^{(1)}_{33} + K^{(2)}_{11} & (K^{(1)}_{34} + K^{(2)}_{12})& K^{(2)}_{13} & K^{(2)}_{14} \\ K^{(1)}_{41} & K^{(1)}_{42} & (K^{(1)}_{43} + K^{(2)}_{21}) & (K^{(1)}_{44} + K^{(2)}_{22}) & K^{(2)}_{23} & K^{(2)}_{24} \\ 0 & 0 & K^{(2)}_{31} & K^{(2)}_{32} & K^{(2)}_{33} & K^{(2)}_{34} \\ 0 & 0 & K^{(2)}_{41} & K^{(2)}_{42} & K^{(2)}_{43} & K^{(2)}_{44} \end{array} \right]=\left[ \begin{array}{cccccc} 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.0625 & -2.1752 & -2.5 & 2.5 \\ -0.3248 & -0.1875 & -2.1752 & 2.6875 & 2.5 & -2.5 \\ 0 & 0 & -2.5 & 2.5 & 2.5 & -2.5 \\ 0 & 0 & 2.5 & -2.5 & -2.5 & 2.5 \end{array} \right]

$$

Next, the constraints imposed on the Global Displacement Matrix due to the pinned nodes are taken into account. These constraints are: $$d_{1}=0$$, $$d_{2}=0$$, $$d_{5}=0$$, $$d_{6}=0$$.

$$

\textbf{d} = \begin{Bmatrix} 0 \\ 0 \\ d_{3} \\ d_{4} \\ 0 \\ 0 \end{Bmatrix}

$$

Because the first and last two rows of the Global Displacement Matrix are all 0, the corresponding rows and columns can be deleted from the Global Stiffness Matrix and Global Force Matrix in order to reduce the Global Force Displacement Relation to a 2 x 2 linear system. This system is shown below.

$$

\left[ \begin{array}{cc} K_{33} & K_{34} \\ K_{43} & K_{44}\end{array} \right]\begin{Bmatrix} d_{3} \\ d_{4}\end{Bmatrix}=\begin{Bmatrix} F_{3} \\ F_{4}\end{Bmatrix}

$$

The process of solving this linear system for the global displacements $$d_{3}$$ and $$d_{4}$$ is shown below.

$$

\left[ \begin{array}{cc} 3.0625 & -2.1752 \\ -2.7152 & 2.6875\end{array} \right]\begin{Bmatrix} d_{3} \\ d_{4}\end{Bmatrix}=\begin{Bmatrix} 0 \\ 7\end{Bmatrix}

$$

$$

\begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}=\textbf{K}^{-1}\begin{Bmatrix}{0}\\{7}\end{Bmatrix}

$$

$$

\textbf{K}^{-1}=\dfrac{1}{\text{det} \textbf{K}}\left[ \begin{array}{cc} K_{44} & -K_{34} \\ -K_{43} & K_{33}\end{array} \right]=\dfrac{1}{(3.0625)(2.6875)-(-2.1752)(-2.1752)}\left[ \begin{array}{cc} 2.6875 & 2.1752 \\ 2.1752 & 3.0625\end{array} \right]=\left[ \begin{array}{cc} 0.7681 & 0.6217 \\ 0.6217 & 0.8753\end{array} \right]

$$

$$

\begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}=\left[ \begin{array}{cc} 0.7681 & 0.6217 \\ 0.6217 & 0.8753\end{array} \right]\begin{Bmatrix}{0}\\{7}\end{Bmatrix}

$$

$$

\begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}=\begin{Bmatrix}4.3519 \\ 6.1271 \end{Bmatrix}

$$

In order determine the internal forces and reactions in both members, it is necessary to use the newly determined values from the Global Displacement Matrix in the Element Force Displacement Relation. The process of solving for the internal forces and reactions for Element 1 is shown below.

$$

\textbf{k}^{(1)}\textbf{d}^{(1)}=\textbf{f}^{(1)}

$$

$$

\textbf{d}^{(1)}=\begin{Bmatrix} 0 \\ 0 \\ 4.3519 \\ 6.1271 \end{Bmatrix}

$$

$$

\left[\begin{array}{cccc} 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 \end{array} \right]\begin{Bmatrix} 0 \\ 0 \\ 4.3519 \\ 6.1271 \end{Bmatrix}=\begin{Bmatrix} f_{1}^{(1)} \\ f_{2}^{(1)} \\ f_{3}^{(1)} \\ f_{4}^{(1)} \end{Bmatrix}

$$

$$

\begin{Bmatrix} f_{1}^{(1)} \\ f_{2}^{(1)} \\ f_{3}^{(1)} \\ f_{4}^{(1)} \end{Bmatrix}=\begin{Bmatrix} -4.438 \\ -2.562 \\ 4.438 \\ 2.562 \end{Bmatrix}

$$

In the above Element Force Matrix for Element 1, $$f_{1}^{(1)}$$ and $$f_{2}^{(1)}$$ represent the reactions at Global Node 1, whereas $$f_{3}^{(1)}$$ and $$f_{4}^{(1)}$$ represent the internal forces in Element 1.

The process shown above is repeated below in order to solve for the internal forces and reactions of Element 2.

$$

\textbf{k}^{(2)}\textbf{d}^{(2)}=\textbf{f}^{(2)}

$$

$$

\textbf{d}^{(2)}=\begin{Bmatrix} 4.3519 \\ 6.1271 \\ 0 \\ 0 \end{Bmatrix}

$$

$$

\left[\begin{array}{cccc} 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{array} \right]\begin{Bmatrix} 4.3519 \\ 6.1271 \\ 0 \\ 0 \end{Bmatrix}=\begin{Bmatrix} f_{1}^{(2)} \\ f_{2}^{(2)} \\ f_{3}^{(2)} \\ f_{4}^{(2)} \end{Bmatrix}

$$

$$

\begin{Bmatrix} f_{1}^{(2)} \\ f_{2}^{(2)} \\ f_{3}^{(2)} \\ f_{4}^{(2)} \end{Bmatrix}=\begin{Bmatrix} -4.438 \\ 4.438 \\ 4.438 \\ -4.438 \end{Bmatrix}

$$

In the above Element Force Matrix for Element 2, $$f_{1}^{(2)}$$ and $$f_{2}^{(2)}$$ represent the internal forces in Element 2, whereas $$f_{3}^{(2)}$$ and $$f_{4}^{(2)}$$ represent the reactions at Global Node 3.

Relating the element forces solved for with the global forces, an array possessing only the reactions at Global Node 1 and Global Node 3 is shown below.

$$

\begin{Bmatrix} f_{1} \\ f_{2} \\ f_{5} \\ f_{6} \end{Bmatrix}=\begin{Bmatrix} -4.438 \\ -2.562 \\ 4.438 \\ -4.438 \end{Bmatrix}

$$

MATLAB Code
% Two bar truss example clear all; e = [3 5]; A = [1 2]; P = 7; L=[4 2]; alpha = pi/3; beta = pi/4; nodes = [0, 0; L(1)*cos(pi/2-alpha), L(1)*sin(pi/2-alpha); L(1)*cos(pi/2-alpha)+L(2)*sin(beta),L(1)*sin(pi/2-alpha)-L(2)*cos(beta)]; dof=2*length(nodes); conn=[1,2; 2,3]; lmm = [1, 2, 3, 4; 3, 4, 5, 6]; elems=size(lmm,1); K=zeros(dof); R = zeros(dof,1); debc = [1, 2, 5, 6]; ebcVals = zeros(length(debc),1); %load vector R = zeros(dof,1); R(4) = P; % Assemble global stiffness matrix K=zeros(dof); for i=1:elems lm=lmm(i,:); con=conn(i,:); k_local=e(i)*A(i)/L(i)*[1 -1; -1 1] k=PlaneTrussElement(e(i), A(i), nodes(con,:)) K(lm, lm) = K(lm, lm) + k; end K R % Nodal solution and reactions [d, reactions] = NodalSoln(K, R, debc, ebcVals) results=[]; for i=1:elems results = [results; PlaneTrussResults(e, A, ... nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g results k_local = 0.75       -0.75      -0.75         0.75 k = 0.5625     0.32476      -0.5625     -0.32476      0.32476       0.1875     -0.32476      -0.1875      -0.5625     -0.32476       0.5625      0.32476     -0.32476      -0.1875      0.32476       0.1875 k_local = 5   -5      -5     5 k = 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 K = 0.5625     0.32476      -0.5625     -0.32476            0            0      0.32476       0.1875     -0.32476      -0.1875            0            0      -0.5625     -0.32476       3.0625      -2.1752         -2.5          2.5     -0.32476      -0.1875      -2.1752       2.6875          2.5         -2.5            0            0         -2.5          2.5          2.5         -2.5            0            0          2.5         -2.5         -2.5          2.5 R = 0         0          0          7          0          0 d = 0         0      4.352     6.1271          0          0 reactions = -4.4378      -2.5622        4.4378       -4.4378 results = 1.7081      5.1244       8.5406       5.1244       17.081        0.6276       1.8828        3.138       1.8828        6.276

Conclusions
As shown in the previous sections, the results of the manual calculations match the results determined from the MATLAB Code written by Professor Loc Vu Quoc. This shows that the proper methodology was taken and accurate calculations were made throughout the process.

Contributing Team Members
Andrew McDonald - Eml4500.f08.bike.mcdonald 22:01, 25 September 2008 (UTC) Garrett Pataky - Eml4500.f08.bike.pataky 22:02, 25 September 2008 (UTC) Samuel Bernal - Eml4500.f08.bike.bernal 13:07, 26 September 2008 (UTC) Bobby Sweeting - Eml4500.f08.bike.sweeting 18:33, 25 September 2008 (UTC) Shawn Gravois - Eml4500.f08.bike.gravois 18:22, 25 September 2008 (UTC) Eric Viale - Eml4500.f08.bike.viale 21:41, 25 September 2008 (UTC)