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

See my comments below. Please don't remove these comments in case you amend your work; just add a comment in these comment boxes to say what you did. Eml4500.f08 12:56, 10 November 2008 (UTC)

1 Class Notes
2 Derivation of Matrix Relations and Matrices for Finite Element Analysis of 3-D Systems
3 Rerunning the MATLAB Code for 2-Bar Truss System and Interpreting results Array
- 3.1 Modified MATLAB Code
- 3.2 Interpreting results Array
  - 3.2.1 Element 1
  - 3.2.2 Element 2
4 Finite Element Analysis of Various 6-Bar Truss Systems
5 Finite Element Analysis of 3-Bar Space Truss
6 Contributing Team Members

Class Notes [edit]

Principal of Virtual Work [edit]

To justify the elimination of rows 1, 2, 5, and 6 in order to obtain K as a 2x2 matrix in a 2-bar truss, the FD relation

$\mathbf{Kd} = \mathbf{F}$

Can be rearranged to be:

$\mathbf{Kd} - \mathbf{F} = \mathbf{0}$ (1)

Which is designated to be equation 1.

Utilizing the principle of virtual work, equation 1 becomes:

$\mathbf{W}\cdot (\mathbf{Kd} - \mathbf{F}) = \mathbf{0}$ (2)

Where W is a 6x1 matrix and Kd - F is a 6x1, and 0 is a scalar number. This equation is designated as equation 2. This equation is true for all values of W and W is called the weighting matrix.

Considering this, Equation 2 being true implies Equation one is also true, and vice versa.

The proof going from Equation 1 to 2 is trivial. However, going from two to one is not trivial and to do this we must choose different values for W

1st Choice:

Choose W so that

$W_{1} = 1 , W_{2} = W_{3} = W_{4} = W_{5} = W_{6} = 0$

Therefore

$\mathbf{W^T}= \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$

and

$\mathbf{W} \cdot ( \mathbf{Kd-F} ) = 1[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]$

Thus, the result becomes

$\sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j= \mathbf{F}_1$

Which is designated to be equation (a).

2nd Choice:

Choose W such that W₂=1, W₁=W₃=W₄=W₅=W₆=0 so that:

$\mathbf{W^T}= \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 \end{bmatrix}$

Therefore:

$\mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 1[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]$

Thus, the result becomes

$\sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j= \mathbf{F}_2$

Which is designated to be equation (b)

3rd Choice: C hoose W such that W₃=1, W₁=W₂=W₄=W₅=W₆=0 so that

$\mathbf{W^T}= \begin{bmatrix} 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 1 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]$

Thus, the result becomes $\sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j= \mathbf{F}_3$

Which is designated to be equation (c)

Choice 4:

Choose W such that W₄=1, W₁=W₂=W₃=W₅=W₆=0 so that

$\mathbf{W^T}= \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 1 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]$

Thus, the result becomes $\sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j= \mathbf{F}_4$

This is designated to be equation (d)

Choice 5:

Choose W such that W₅=1, W₁=W₂=W₃=W₄=W₆=0 so that

$\mathbf{W^T}= \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 1 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]$

Thus, the result becomes $\sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j= \mathbf{F}_5$

This is designated to be equation (e)

Choice 6:

Choose W such that W₆=1, W₁=W₂=W₃=W₄=W₅=0 so that

$\mathbf{W^T}= \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 1 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]$

so that the result becomes $\sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j= \mathbf{F}_6$

This equation is designated to be equation (f)

Thus, when the equations (a) through (f) are combined, they become $\mathbf{Kd}= \mathbf{F}$ which is equation (1) from before.

The weighting coefficients must be kinetically admissible which means they cannot violate the boundary conditions so W₁ = W₂ = W₅= W₆=0 where the weighting coefficients represent the virtual displacement.

Now,

$\mathbf{W} \cdot (\mathbf{Kd}-\mathbf{F})$

Where Kd is a 6x2 matrix minus a 2x1 vector.

Thus

$= \begin{Bmatrix} W_3 \\ W_4 \end{Bmatrix} (\mathbf{\overline{K} \overline{d}- \overline{F} })$ (3)

Which is designated as Equation 3.

and is true for all values of $\begin{Bmatrix} W_3 \\ W_4 \end{Bmatrix}$

and the K matrix is reduced to $\mathbf{\overline{K}}=\begin{bmatrix} K_{33} && K_{34} \\ K_{43} && K_{44}\end{bmatrix}$

as well as:

$\mathbf{\overline{d}}= \begin{Bmatrix}d_{3} \\ d_{4} \end{Bmatrix}$

$\overline{F } = \begin{Bmatrix}F_{3} \\ F_{4} \end{Bmatrix}$

being the reduced F and d matrices.

Analysis of The Axial Force Displacement Relation Using The Principle of Virtual Work. [edit]

Recall:

$\mathbf{q}^{(i)}=\mathbf{Td}^{(i)}$

An example of the use of a weighting factor is shown in the computation of our course grade. This computation is shown below.

Course Grade = $\alpha _{o} \cdot (HW Grade)+\sum_{i=1}^{3}{\alpha _{i}}\cdot (Exam_{i})$

The weighting factors used in the above computation are $\alpha _{o}$ and $\alpha _{i}$

Deriving the following relation,

$\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$

and also recalling the following Force Displacement relation with axial degrees of freedom,

$\mathbf{\hat{k}}^{e}_{2\times 2}\mathbf{q}^{(e)}_{2\times 1}= \mathbf{p}^{(e)}_{2\times 1}$

By subtracting the force matrix $\mathbf{p}$ , the relation is now written as follows.

$\mathbf{\hat{k}}^{(e)}\mathbf{q}^{(e)}- \mathbf{p}^{(e)}=0 \qquad$ Equation (1)

Now, the Principle of Virtual Work can be applied, and it is shown below.

$\mathbf{W}\cdot (\mathbf{\hat{k}}^{e}\mathbf{q}^{(e)}-\mathbf{p}^{(e)})=0 \qquad$ Equation (2)

This equation holds true for all $\mathbf{W}$ .

Assigned HW - Prove that equation (1) $\iff$ equation (2)

The first step is to select a W such that $W_{1} = 1, W_{2}=W_{3}=W_{4}=W_{5}=W_{6}=0$

$\mathbf{W}^{T}= \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{\hat{k}q-P} ) = 1[ \sum^6_{j=1}\mathbf{\hat{k}}_{1j} \mathbf{q}_j- \mathbf{P}_1] + 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{2j} \mathbf{q}_j- \mathbf{P}_2] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{3j} \mathbf{q}_j- \mathbf{P}_3] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{4j} \mathbf{q}_j- \mathbf{P}_4] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{5j} \mathbf{q}_j- \mathbf{P}_5] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{6j} \mathbf{q}_j- \mathbf{P}_6]$

The first equation can be rewritten as shown.

$\sum^6_{j=1}\mathbf{\hat{k}}_{1j}\mathbf{q}_j= \mathbf{P}_1$

Then, this same procedure is going to be repeated five times, for each element in W being equal to one.

$W_{2} = 1, W_{1}=W_{3}=W_{4}=W_{5}=W_{6}=0$

$\mathbf{W}^{T}= \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{\hat{k}q-P} ) = 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{1j} \mathbf{q}_j- \mathbf{P}_1] + 1[ \sum^6_{j=1}\mathbf{\hat{k}}_{2j} \mathbf{q}_j- \mathbf{P}_2] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{3j} \mathbf{q}_j- \mathbf{P}_3] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{4j} \mathbf{q}_j- \mathbf{P}_4] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{5j} \mathbf{q}_j- \mathbf{P}_5] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{6j} \mathbf{q}_j- \mathbf{P}_6]$

The second equation can be rewritten as shown.

$\sum^6_{j=1}\mathbf{\hat{k}}_{2j}\mathbf{q}_j= \mathbf{P}_2$

$W_{3} = 1, W_{1}=W_{2}=W_{4}=W_{5}=W_{6}=0$

$\mathbf{W}^{T}= \begin{bmatrix} 0 & 0 & 1 & 0 & 0 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{\hat{k}q-P} ) = 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{1j} \mathbf{q}_j- \mathbf{P}_1] + 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{2j} \mathbf{q}_j- \mathbf{P}_2] + 1 [ \sum^6_{j=1}\mathbf{\hat{k}}_{3j} \mathbf{q}_j- \mathbf{P}_3] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{4j} \mathbf{q}_j- \mathbf{P}_4] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{5j} \mathbf{q}_j- \mathbf{P}_5] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{6j} \mathbf{q}_j- \mathbf{P}_6]$

The third equation can be rewritten as shown.

$\sum^6_{j=1}\mathbf{\hat{k}}_{3j}\mathbf{q}_j= \mathbf{P}_3$

$W_{4} = 1, W_{1}=W_{2}=W_{3}=W_{5}=W_{6}=0$

$\mathbf{W}^{T}= \begin{bmatrix} 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{\hat{k}q-P} ) = 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{1j} \mathbf{q}_j- \mathbf{P}_1] + 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{2j} \mathbf{q}_j- \mathbf{P}_2] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{3j} \mathbf{q}_j- \mathbf{P}_3] + 1 [ \sum^6_{j=1}\mathbf{\hat{k}}_{4j} \mathbf{q}_j- \mathbf{P}_4] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{5j} \mathbf{q}_j- \mathbf{P}_5] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{6j} \mathbf{q}_j- \mathbf{P}_6]$

The fourth equation can be rewritten as shown.

$\sum^6_{j=1}\mathbf{\hat{k}}_{4j}\mathbf{q}_j= \mathbf{P}_4$

$W_{5} = 1, W_{1}=W_{2}=W_{3}=W_{4}=W_{6}=0$

$\mathbf{W}^{T}= \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{\hat{k}q-P} ) = 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{1j} \mathbf{q}_j- \mathbf{P}_1] + 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{2j} \mathbf{q}_j- \mathbf{P}_2] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{3j} \mathbf{q}_j- \mathbf{P}_3] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{4j} \mathbf{q}_j- \mathbf{P}_4] + 1 [ \sum^6_{j=1}\mathbf{\hat{k}}_{5j} \mathbf{q}_j- \mathbf{P}_5] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{6j} \mathbf{q}_j- \mathbf{P}_6]$

The fifth equation can be rewritten as shown.

$\sum^6_{j=1}\mathbf{\hat{k}}_{5j}\mathbf{q}_j= \mathbf{P}_5$

$W_{6} = 1, W_{1}=W_{2}=W_{3}=W_{4}=W_{5}=0$

$\mathbf{W}^{T}= \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix}$

$\mathbf{W} \cdot ( \mathbf{\hat{k}q-P} ) = 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{1j} \mathbf{q}_j- \mathbf{P}_1] + 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{2j} \mathbf{q}_j- \mathbf{P}_2] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{3j} \mathbf{q}_j- \mathbf{P}_3] + 0 [ \sum^6_{j=1}\mathbf{\hat{k}}_{4j} \mathbf{q}_j- \mathbf{P}_4] + 0[ \sum^6_{j=1}\mathbf{\hat{k}}_{5j} \mathbf{q}_j- \mathbf{P}_5] + 1 [ \sum^6_{j=1}\mathbf{\hat{k}}_{6j} \mathbf{q}_j- \mathbf{P}_6]$

The sixth equation can be rewritten as shown.

$\sum^6_{j=1}\mathbf{\hat{k}}_{6j}\mathbf{q}_j= \mathbf{P}_6$

In conclusion, the equations that were computed above become equation (1), shown below.

$\mathbf{\hat{k}}^{(e)}\mathbf{q}^{(e)}- \mathbf{p}^{(e)}=0 \qquad$

This proves that equation (1) $\iff$ equation (2).

Principle of Virtual Work Conversion to Force Displacement Equation [edit]

$\hat{\textbf{w}}_{2x1}$ = virtual axial displacement corresponding to $\textbf{q}^{(e)}$ .
$\textbf{w}_{4x1}$ = virtual displacement in global coordinate system, corresponding to $\textbf{d}^{(e)}$ .

Substituting Equations (3) and (4) into Equation (2) yields the following Equation:
$\left(\textbf{T}^{(e)}\textbf{W} \right)\cdot \left[\mathbf{\hat{k}}^{(e)}\left(\mathbf{T}^{(e)}\mathbf{\bar{d}}^{(e)}-\mathbf{p}^{(e)} \right) \right]=0$ for all $\textbf{w}_{4x1} \; \; \; \; \; \; \;(5)$

In order to solve Equation (5), a few matrix operations must be recalled:

Assigned HW - Verification of Equation (6) Example.

$\textbf{A}_{2x3}=\begin{bmatrix} 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix}$ $\textbf{B}_{3x3}=\begin{bmatrix} 7 & 8 & 9\\ 1 & 2 & 3\\ 4 & 5 & 6 \end{bmatrix}$
$\mathbf{AB}=\begin{bmatrix} 21 & 27 & 33\\ 57 & 72 & 87 \end{bmatrix}$ $\mathbf{AB}^T=\begin{bmatrix} 21 & 57 \\ 27 & 72 \\ 33 & 87 \\ \end{bmatrix}$
$\textbf{A}^T_{2x3}=\begin{bmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \\ \end{bmatrix}$ $\textbf{B}^T_{3x3}=\begin{bmatrix} 7 & 1 & 4\\ 8 & 2 & 5\\ 9 & 3 & 6 \end{bmatrix}$ $\mathbf{B}^T\mathbf{A}^T=\begin{bmatrix} 21 & 57 \\ 27 & 72 \\ 33 & 87 \\ \end{bmatrix}$

Therefore, $\mathbf{B}^T\mathbf{A}^T = \mathbf{(AB)}^T$

To solve Equation (5), Equation (6) and (7) must be applied:

The dot product must be replaced with the transpose as shown below from Equation (7).
$\left(\textbf{T}^{(e)}\textbf{W} \right)^T\left[\mathbf{\hat{k}}^{(e)}\left(\mathbf{T}^{(e)}\mathbf{\bar{d}}^{(e)}-\mathbf{p}^{(e)} \right) \right]=0$ for all $\textbf{w}_{4x1}$

The matrix transpose is replaced with the product of individual transposes from Equation (6).
$\textbf{W}^{T}\textbf{T}^{(e)T}\left[\mathbf{\hat{k}}^{(e)}\left(\mathbf{T}^{(e)}\mathbf{\bar{d}}^{(e)}-\mathbf{p}^{(e)} \right) \right]=0$ for all $\textbf{w}_{4x1}$

After distributing the transpose matrix T the following equation is the result.
$\textbf{W}^{T}\cdot\left[\left(\textbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}\right)\mathbf{\bar{d}}^{(e)}-\left(\textbf{T}^{(e)T}\mathbf{p}^{(e)} \right) \right]=0$ for all $\textbf{w}_{4x1}$

$\mathbf{k}^{(e)}$ replaced $\textbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$ , and
$\mathbf{f}^{(e)}$ replaced $\textbf{T}^{(e)T}\mathbf{p}^{(e)}$ .

The resulting equation is:
$\mathbf{w}\cdot\left[\mathbf{k}^{(e)}\mathbf{d}^{(e)}-\mathbf{f}^{(e)}\right]=0$ for all $\mathbf{w}$

After dividing both sides by w, and moving f^(e) to the other side yielded:
$\mathbf{k}^{(e)}\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$

The previous methods have all been for discrete, non-continuous cases (matrices), and the new analytic methods will be for continuous cases (PDEs).

Motivational Model Problem:
An elastic bar with varying A(x) and E(x) is subject to a varying axial load (distributed), concentrated load, and inertia force (dynamic). Both the axial and concentrated loads are time dependent. The figure below depicts the problem.

Loaded elastic bar with varying properties.

Equation of Motion, Constitutive Relation, PDE of Motion, and Boundary Condition Derivation [edit]

Assigned HW - Literature Research of Composite Materials.

Composite materials are complex materials made from multiple substances, called constituents, that have significantly different physical or chemical properties. Although the materials making up the composite act in concert, they do not merge completely and each contains its own identity. For this reason two or more materials could be combined in order to take advantage of the good characteristics of each material.

Composites usually consist of two components, the matrix and the filler. The matrix holds the filler in an orderly pattern, and helps transfer the load among the filler. Fillers can be any material, such as ceramic, carbon fiber, or sand, and are usually added in order to increase the strength of the composite. Concrete is an example of a composite, where the cement is the matrix and the sand acts as the filler.

Composites can be classified into 4 different types depending on the filler:

Particulate - The filler particles have roughly equal dimensions in all directions, such as resin powder or gravel. The particles are very close to spherical, but do not have to be perfectly round.

Short fiber - When one dimension of the particle becomes slightly longer than the other, it can be considered a short fiber. The matrix must transfer the load very quickly when the fibers are short, and the composite can not take on a lot of the properties of the filler.

Long fiber - When the fibers become long(ratio of length to diameter approaching infinity), they are said to be continuous, and there are few breaks in the filler. This causes the composite to take on more of the filler properties, such as an increase in strength. For this reason long fibers are usually used in high performance applications such as aerospace structures.

Laminate - Composite properties are most like the filler in the direction of the fibers. Perpendicular to the fibers the composite takes on mostly matrix properties. Therefore, it is necessary to orient fibers in multiple direction to gain maximum benefit from a composite. Plies are stacked together to form a laminate in order to obtain the most efficient composite possible. It is normally most efficient to have most of the fibers oriented along the primary load direction, with just enough fibers placed in the other direction simply to hold the structure together.

References:

http://www.science.org.au/nova/059/059key.htm
http://www.acmanet.org/professionals/index.cfm

The free body diagram of section dx can be seen below in Figure 1.

File:BikefxFBD.jpg

Figure 1 - Free Body Diagram of Section dx

Summing the forces in the x direction yields

$\sum F_x = 0 = -N(x,t) + N(x+dx,t) + f(x,t)dx - m(x){\ddot{u}}dx$

$= \frac{\part N(x,t)}{\part x}dx + h.o.t + f(x,t)dx - m(x){\ddot{u}}dx$ ⇨ $Equation \left(1\right)$

The dx terms are canceled out of this equation, as well as the high order terms(h.o.t). Let this equation be call equation (1).

Why are the high order terms neglected?

Recall Taylor series expansion: $f(x+dx) = f(x) + \frac{df(x)}{dx}dx + \frac{1}{2}\frac{d^{2}f(x)}{dx^{2}}dx^{2} + ...$

where $\frac{1}{2}\frac{d^{2}f(x)}{dx^{2}}dx^{2}$ is considered a h.o.t and is therefore neglected.

From equation (1) the equation of motion(EOM) of the system can be derived. Let the EOM be called equation (2).

$\frac{\part N}{\part x} + f = m{\ddot{u}}$ ⇨ $Equation \left(2\right)$

Equation (3), the constitutive relation, can then be determined.

$N\left(x,t\right) = A(x)\sigma(x,t)$ ⇨ $Equation \left(3\right)$

where $\sigma\left(x,t\right) = E(x)\epsilon(x,t)$

and $\epsilon(x,t) = \frac{\part u(x,t)}{\part x}$

Substituting equation (3) into equation (2) yields the partial differential equation of motion, equation (4) shown below.

$\frac{\part}{\part x}[A(x)E(x)\frac{\part u}{\part x}] + f(x,t) = m(x){\ddot{u}}$ ⇨ $Equation \left(4\right)$

There are 2 initial conditions(displacement and velocity), and therefore 2 boundary conditions(2nd order derivative with respect to x).

File:Bikebc1.jpg

Figure 2

For Figure 2 the boundary conditions are $u\left(0,t\right) = 0 = u(L,t)$

File:Bikebc2.jpg

Figure 3

For Figure 3 the boundary condition at 0 is also $u\left(0,t\right) = 0$ , but the boundary condition at L can be derived as follows

$N\left(L,t\right) = F(t)$

where $N = A\left(L\right)\sigma(L,t)$

and $\sigma\left(L,t\right) = E(L)\epsilon(L,t)$

and $\epsilon(L,t) = \frac{\part u(L,t)}{\part x}$

The boundary condition at L can therefore be simplified to

$\frac{\part u(L,t)}{\part x} = \frac{F(t)}{A(L)E(L)}$

Derivation of Matrix Relations and Matrices for Finite Element Analysis of 3-D Systems [edit]

Axial FD Relation [edit]

When an additional displacement DOF is added to the system, the Axial Force-Displacement Relation does not change. This is because a member can only have two-reaction forces, and those forces are constrained to be along the direction of the member. Because deformation in a member is caused by these forces, the only axial displacements in a space truss member are found at the bar ends in the direction of the member.

Thus, the Axial Force-Displacement Relation is as shown below, where P represents the axial forces at each node in the member and q represents the axial displacement of each node.

$\mathbf{P}^{(e)} = \mathbf{\hat{k}}^{(e)}\mathbf{q}^{(e)}$

$\begin{Bmatrix}P_{1}^{(e)} \\ P_{2}^{(e)}\end{Bmatrix} = k^{(e)}\begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}\begin{Bmatrix}q_{1}^{(e)}\\ q_{2}^{(e)}\end{Bmatrix}$

Element DOFs and Element Forces in Global xyz Coordinates [edit]

The following image shows the proper setup of the element displacement DOFs and element forces for a space truss element.

Proper setup of the element displacement DOFs and element forces for a space truss element

Since each element in a space truss system has 6 displacement DOFs and forces, the element displacement DOF and force matrices are of size 6 x 1. These matrices are shown below.

$\mathbf{d}^{(e)} = \begin{Bmatrix}d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{4}^{(e)}\\ d_{5}^{(e)}\\ d_{6}^{(e)}\end{Bmatrix} \qquad \mathbf{f}^{(e)} = \begin{Bmatrix}f_{1}^{(e)}\\ f_{2}^{(e)}\\ f_{3}^{(e)}\\ f_{4}^{(e)}\\ f_{5}^{(e)}\\ f_{6}^{(e)}\end{Bmatrix}$

Transformation Matrix [edit]

For plane truss members, the transformation matrix is as shown below, where $l^{(e)}$ and $m^{(e)}$ are the director cosines for a 2D system.

$\mathbf{T}^{(e)}_{2 x 4} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix}$

For space truss members, the transformation matrix is slightly different, due to the fact that an additional coordinate, z, is added to the system. The director cosine corresponding to the z-coordinate is represented by $n^{(e)}$ . The transformation matrix for a space truss system is shown below.

$\mathbf{T}^{(e)}_{2 x 6} = \begin{bmatrix} l^{(e)} & m^{(e)} & n^{(e)} & 0 & 0 & 0 \\ 0 & 0 & 0 & l^{(e)} & m^{(e)} & n^{(e)} \end{bmatrix}$

Relation between Axial DOFs and Element DOFs in Global Coordinates [edit]

The relation between the axial displacement DOFs and element DOFs using the transformation matrix is shown below.

$\mathbf{q}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{d}_{6\times 1}^{(e)}$

$\begin{Bmatrix} q_{1}^{(e)}\\ q_{2}^{(e)} \end{Bmatrix} = \begin{bmatrix} l^{(e)} &m^{(e)} &n^{(e)} & 0 & 0& 0\\ 0& 0 & 0 &l^{(e)} &m^{(e)} &n^{(e)} \end{bmatrix}\begin{Bmatrix} d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{4}^{(e)}\\ d_{5}^{(e)}\\ d_{6}^{(e)} \end{Bmatrix}$

Relation between Axial Forces and Element Forces in Global Coordinates [edit]

The relation between the axial forces and element forces using the transformation matrix is shown below.

$\mathbf{P}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)}$

$\begin{Bmatrix} P_{1}^{(e)}\\ P_{2}^{(e)} \end{Bmatrix} = \begin{bmatrix} l^{(e)} &m^{(e)} &n^{(e)} & 0 & 0& 0\\ 0& 0 & 0 &l^{(e)} &m^{(e)} &n^{(e)} \end{bmatrix}\begin{Bmatrix} f_{1}^{(e)}\\ f_{2}^{(e)}\\ f_{3}^{(e)}\\ f_{4}^{(e)}\\ f_{5}^{(e)}\\ f_{6}^{(e)} \end{Bmatrix}$

Derivations Using the Principal of Virtual Work [edit]

In this section, various relations for space truss elements are derived using the Principal of Virtual Work.

Global FD Relation [edit]

The element axial force-displacement relation is shown below.

$\mathbf{P}^{(e)} = \mathbf{\hat{k}}^{(e)}\mathbf{q}^{(e)}$

Relations for the element axial displacement DOFs and element axial forces in terms of the transformation matrix and the element displacement DOFs and forces in global coordinates are shown below.

$\mathbf{q}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{d}_{6\times 1}^{(e)}$

$\mathbf{P}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)}$

Substituting these relations into the element axial force-displacement relation, the following relation is developed.

$\mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)} = \mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6}\mathbf{d}_{6\times 1}^{(e)}$

Applying the Principal of Virtual Work, the above relation can be written as shown below, greatly simplifying the problem.

$\mathbf{f}_{6\times 1}^{(e)} = [\mathbf{T}^{(e)^{T}}_{2\times 6}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6}]\mathbf{d}_{6\times 1}^{(e)}$

Error: No, you cannot "solve" the equation further above to obtain the above equation; you need to use the PVW. Eml4500.f08 11:06, 10 November 2008 (UTC)

The above derivation has been modified to explain how the PVW is used. Eml4500.f08.bike.mcdonald 21:55, 11 November 2008 (UTC)

Multiplying out the term inside the brackets, it can be shown that it is equal to the element stiffness matrix, as shown below.

$\mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} = \begin{bmatrix} l^{(e)}&0 \\ m^{(e)}&0 \\ n^{(e)}&0 \\ 0&l^{(e)} \\ 0&m^{(e)} \\ 0&n^{(e)} \end{bmatrix}k^{(e)}\begin{bmatrix} 1 &-1 \\ -1&1 \end{bmatrix}\begin{bmatrix} l^{(e)}& m^{(e)} &n^{(e)} & 0 & 0 &0 \\ 0& 0 & 0 & l^{(e)} &m^{(e)} &n^{(e)} \end{bmatrix}$

$\mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} = k^{(e)}\begin{bmatrix} l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} & -l^{(e)^{2}} & -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} \\ m^{(e)}l^{(e)}& m^{(e)^{2}} & m^{(e)}n^{(e)} & -m^{(e)}l^{(e)} & -m^{(e)^{2}} & -m^{(e)}n^{(e)} \\ n^{(e)}l^{(e)}& m^{(e)}n^{(e)} & n^{(e)^{2}} & -n^{(e)}l^{(e)} & -m^{(e)}n^{(e)} & -n^{(e)^{2}} \\ -l^{(e)^{2}}& -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} & l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} \\ -m^{(e)}l^{(e)}& -m^{(e)^{2}} & -m^{(e)}n^{(e)} & m^{(e)}l^{(e)} & m^{(e)^{2}} & m^{(e)}n^{(e)} \\ -n^{(e)}l^{(e)}& -m^{(e)}n^{(e)} & -n^{(e)^{2}} & n^{(e)}l^{(e)} & m^{(e)}n^{(e)} & n^{(e)^{2}} \end{bmatrix}$

$\mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} = \mathbf{k}^{(e)}_{6\times 6}$

Substituting this relation into the above relation for the element force matrix in global coordinates gives the desired result, the element force-displacement relation, proved using the Principle of Virtual Work.

$\mathbf{f}_{6\times 1}^{(e)} = \mathbf{k}^{(e)}_{6\times 6}\mathbf{d}_{6\times 1}^{(e)}$

Element Stiffness Matrix in Global Coordinates in Terms of Transformation Matrix and Element Stiffness Matrix from Axial FD Relation [edit]

The following relation shows the element stiffness matrix of a space truss element in terms of the transformation matrix, $T^{(e)}$ , and the element stiffness matrix from the axial force-displacement relation, $\hat{k}^{(e)}$ . This relation was derived in the above derivation of the global force-displacement relation.

$\mathbf{k}^{(e)}_{6\times 6} = \mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6}$

Element Forces in Global Coordinates in Terms of Transformation Matrix and Element Axial Forces [edit]

The following relation was derived and shown in the above derivation of the space truss element force-displacement relation.

$\mathbf{P}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)}$

Multiplying both sides of this relation by $\mathbf{T}^{(e)^{T}}_{6\times 2}$ gives the relation shown below for the element forces in global coordinates in terms of the transformation matrix and the element axial forces.

$\mathbf{f}_{6\times 1}^{(e)} = \mathbf{P}_{2\times 1}^{(e)}\mathbf{T}^{(e)^{T}}_{6\times 2}$

Element Stiffness Matrix in Global Coordinates in Terms of Director Cosines [edit]

The following relation shows the element stiffness matrix of a space truss element in terms of the director cosines, $l^{(e)}$ , $m^{(e)}$ , and $n^{(e)}$ . This relation was derived in the above derivation of the global force-displacement relation.

$\mathbf{k}_{6\times 6}^{(e)} = k^{(e)}\begin{bmatrix} l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} & -l^{(e)^{2}} & -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} \\ m^{(e)}l^{(e)}& m^{(e)^{2}} & m^{(e)}n^{(e)} & -m^{(e)}l^{(e)} & -m^{(e)^{2}} & -m^{(e)}n^{(e)} \\ n^{(e)}l^{(e)}& m^{(e)}n^{(e)} & n^{(e)^{2}} & -n^{(e)}l^{(e)} & -m^{(e)}n^{(e)} & -n^{(e)^{2}} \\ -l^{(e)^{2}}& -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} & l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} \\ -m^{(e)}l^{(e)}& -m^{(e)^{2}} & -m^{(e)}n^{(e)} & m^{(e)}l^{(e)} & m^{(e)^{2}} & m^{(e)}n^{(e)} \\ -n^{(e)}l^{(e)}& -m^{(e)}n^{(e)} & -n^{(e)^{2}} & n^{(e)}l^{(e)} & m^{(e)}n^{(e)} & n^{(e)^{2}} \end{bmatrix}$

Rerunning the MATLAB Code for 2-Bar Truss System and Interpreting `results` Array [edit]

The setup of the 2-bar truss system to be analyzed is shown below. This system was previously analyzed and verified using statics in Homework Report 2, however there was a bug in the MATLAB code that needed to be corrected. The steps taken to correct the bug will be shown in this section, and the new results that come from the modified MATLAB code will be presented.

Given truss arrangement and loading of two-member truss system

Free body diagram of two-member truss system

The properties of the 2-bar truss system to be analyzed is shown in the table below.

**Properties of 2-Bar Truss System**
Element Number	Young's Modulus	Cross-Sectional Area	Length	Angle of Inclination
1	3	1	4	$30^\circ$
2	5	2	2	$45^\circ$

Modified MATLAB Code [edit]

When the MATLAB code for solving the 2-bar truss system using the Finite Element Method was first analyzed in Homework Report 2, there was a bug that showed up when calculating the results array. This bug was due to the fact that the MATLAB code passed the entire e and A arrays into the PlaneTrussResults(e,A,coord,disps) function, instead of passing only the proper array index in.

The section of code with the error is shown below prior to fixing the bug.

for i=1:elems
    results = [results; PlaneTrussResults(e, A, ...
    nodes(conn(i,:),:), d(lmm(i,:)))];
end

In order to fix the bug, the variables e and A were replaced with e(i) and A(i), which allowed for the values at each index in these arrays to be properly used when calculating the results. This is because the function PlaneTrussResults is located inside a for loop from the value 1 to the variable elems, which in this case has a value of 2.

The modified section of code is shown below.

for i=1:elems
    results = [results; PlaneTrussResults(e(i), A(i), ...
    nodes(conn(i,:),:), d(lmm(i,:)))];
end

The overall modified MATLAB code used to solve the 2-bar truss system is shown below.

% MATLAB Code to perform Finite Element Analysis on
% Two-Bar Truss System
% NOTE: Modification made to fix bug in calculating results
 
clear all;
e = [3 5]; 
A = [1 2]; 
P = 7;
L=[4 2];
alpha = pi/3;
beta = pi/4;
 
nodes = [
         0, 0; % x,y position of node 1
         L(1)*cos(pi/2-alpha), L(1)*sin(pi/2-alpha); % x,y position of node 2
         L(1)*cos(pi/2-alpha)+L(2)*sin(beta),L(1)*sin(pi/2-alpha)-L(2)*cos(beta) % x,y position of node 3
         ];
 
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(i), A(i), ...
    nodes(conn(i,:),:), d(lmm(i,:)))];
end
format short g
results

Interpreting `results` Array [edit]

The MATLAB command window output after running the MATLAB code shown in the above section is shown below.

k_local =

   0.7500   -0.7500
  -0.7500    0.7500


k =

   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


k_local =

    5    -5
   -5     5


k =

   2.5000   -2.5000   -2.5000    2.5000
  -2.5000    2.5000    2.5000   -2.5000
  -2.5000    2.5000    2.5000   -2.5000
   2.5000   -2.5000   -2.5000    2.5000


K =

   0.5625    0.3248   -0.5625   -0.3248         0         0
   0.3248    0.1875   -0.3248   -0.1875         0         0
  -0.5625   -0.3248    3.0625   -2.1752   -2.5000    2.5000
  -0.3248   -0.1875   -2.1752    2.6875    2.5000   -2.5000
        0         0   -2.5000    2.5000    2.5000   -2.5000
        0         0    2.5000   -2.5000   -2.5000    2.5000


R =

    0
    0
    0
    7
    0
    0


d =

        0
        0
   4.3520
   6.1271
        0
        0


reactions =

  -4.4378
  -2.5622
   4.4378
  -4.4378


results =

      1.7081       5.1244       5.1244
      0.6276        3.138        6.276

The last value shown in the MATLAB command window output is the results array. The first row of the array is associated with Element 1, and the second row is associated with Element 2. The first column gives the value of axial strain in the corresponding elements, the second column gives the value of axial stress in the corresponding elements, and the third column gives the value of axial force in the corresponding elements.

A table showing the results of the Finite Element Analysis on the 2-bar truss system is given below.

**Results for 2-Bar Truss System**
Element Number	Axial Strain	Axial Stress	Axial Force
1	1.7081	5.1244	5.1244
2	0.6276	3.138	6.276

In order to verify that the results obtained from the MATLAB code shown above, the calculation of axial force, axial stress, and axial strain using statics is shown below for both element 1 and 2 of the 2-bar truss system.

The axial forces of both elements were obtained by summing forces in both the x and y directions, as shown below.

$\sum{F_{x}}= F_{1}\cos{30^\circ} - F_{2}\cos{45^\circ}$

$\sum{F_{y}}= F_{1}\sin{30^\circ} + F_{2}\sin{45^\circ}$

$F_{1} = 5.124$

$F_{2} = 6.277$

Element 1 [edit]

Knowing the axial force in element 1, the axial stress and axial strain can be calculated as shown below.

$F_{1} = 5.124$

$\sigma_{1} = \frac{F_{1}}{A_{1}} = \frac{5.124}{1} = 5.124$

$\epsilon_{1} = \frac{\sigma_{1}}{E_{1}} = \frac{5.124}{3} = 1.708$

Comparing these values with the values obtained from the Finite Element Method solution above, it can be seen that both methods are in complete agreement.

Element 2 [edit]

Knowing the axial force in element 2, the axial stress and axial strain can be calculated as shown below.

$F_{2} = 6.277$

$\sigma_{2} = \frac{F_{2}}{A_{2}} = \frac{6.277}{2} = 3.139$

$\epsilon_{2} = \frac{\sigma_{2}}{E_{2}} = \frac{3.139}{5} = 0.628$

Comparing these values with the values obtained from the Finite Element Method solution above, it can be seen that both methods are in complete agreement.

Finite Element Analysis of Various 6-Bar Truss Systems [edit]

The setup of the 6-bar truss system as given on page 226 of FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti is shown in the figure below.

Setup of 6-bar truss system to be solved using Finite Element Method

6-Bar Truss System With Same Young's Modulus and Cross-Sectional Area Values [edit]

The properties of the 6-bar truss system with the same values for Young's Modulus and cross-sectional area are shown in the table below.

**Properties of 6-Bar Truss System with Same Values for Young's Modulus and Cross-Sectional Area**
Element Number	Young's Modulus (GPa)	Cross-Sectional Area (m²)	Length (m)	Angle of Inclination
1	200	0.001	4	$0^\circ$
2	200	0.001	2.83	$135^\circ$
3	200	0.001	2.24	$153.4^\circ$
4	200	0.001	3	$90^\circ$
5	200	0.001	2.83	$45^\circ$
6	200	0.001	2.24	$26.6^\circ$

MATLAB Code [edit]

The following MATLAB code was provided as an online resource from FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti.

% Six bar truss example
e = 200*10^3; A = 0.001*1000^2; P = 20000.; 
alpha = pi/6; 
nodes = 1000*[0, 0; 4, 0; 0, 3; 4, 3; 2, 2]; 
dof=2*length(nodes);
conn=[1,2; 2,5; 5,3; 2,4; 1,5; 5,4];
lmm = [1, 2, 3, 4; 3, 4, 9, 10; 9, 10, 5, 6;
    3, 4, 7, 8; 1, 2, 9, 10; 9, 10, 7, 8];
elems=size(lmm,1);
K=zeros(dof); R = zeros(dof,1);
debc = [1, 2, 5, 6, 7, 8];
ebcVals = zeros(length(debc),1);
 
%load vector
R = zeros(dof,1); R(3) = P*sin(alpha); R(4) = P*cos(alpha); 
 
% Assemble global stiffness matrix
K=zeros(dof);
for i=1:elems
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(e, A, 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

In order for this code to run properly, it needs three additional MATLAB functions. These functions are listed and shown in detail below.

PlaneTrussElement.m [edit]

function k = PlaneTrussElement(e, A, coord)
% PlaneTrussElement(e, A, coord)
% Generates stiffness matrix for a plane truss element
% e = modulus of elasticity
% A = area of cross-section
% coord = coordinates at the element ends
 
x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms;
    ls*ms, ms^2,-ls*ms,-ms^2;
    -ls^2,-ls*ms,ls^2,ls*ms;
    -ls*ms,-ms^2,ls*ms,ms^2];

NodalSoln.m [edit]

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

PlaneTrussResults.m [edit]

function results = PlaneTrussResults(e, A, coord, disps)
% results = PlaneTrussResults(e, A, coord, disps)
% Compute plane truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
% The output quantities are eps = axial strain
% sigma = axial stress and force = axial force.
 
x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
T=[ls,ms,0,0; 0,0,ls,ms];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;
results=[eps, sigma, force];

Plot of Undeformed and Deformed Structure [edit]

The plot of the undeformed and deformed shapes of the 6-bar truss system with the same values for Young's Modulus and cross-sectional area is shown below with a multiplication factor on the deformation of 3000.

Plot of undeformed and deformed shapes of the 6-bar truss system analyzed using the Finite Element Method

The source code added to the original MATLAB code to produce the above plot is given below.

%-----------------------------------------------------------------
% model with 2-D beam elements with 2 dofs per node
   dof = 2;         %  dof per node
 
%
% input the nodal coordinates
%
   n_node = 5;              % total number of nodes
   n_elem = 6;              % total number of elements
   total_dof = dof * n_node;  % total dofs of system
 
% five-bar truss system data
%
   for i=1:n_elem
        con=conn(i,:);
        L(i)=PlaneTrussElementMod(nodes(con,:));
   end
   for i=1:n_elem
        con=conn(i,:);
        theta(i)=PlaneTrussElementMod2(nodes(con,:));
   end
 
   position(:, 1) = [ 0; 0];
   position(:, 2) = [ L(1); 0];
   position(:, 3) = [ 0; L(4)];
   position(:, 4) = [ L(1) ; L(4) ];
   position(:, 5) = [ L(5)*cos(theta(5)) ; L(5)*sin(theta(5)) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position(1,i);
      y(i) = position(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect(1, 1) = 1;   % element 1
   node_connect(2, 1) = 2;
 
   node_connect(1, 2) = 2;   % element 2
   node_connect(2, 2) = 5;
 
   node_connect(1, 3) = 5;   % element 3
   node_connect(2, 3) = 3;
 
   node_connect(1, 4) = 2;   % element 4
   node_connect(2, 4) = 4;
 
   node_connect(1, 5) = 1;   % element 5
   node_connect(2, 5) = 5;
 
   node_connect(1, 6) = 5;   % element 6
   node_connect(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem  
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 5000 -1000 4000])
      % solid line
      plot(xx,yy,'--')
      % hold the current figure for the next element
      hold on
   end
 
 
   text(x(node_connect(1,1)),y(node_connect(1,1)),'Node 1',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(1,2)),y(node_connect(1,2)),'Node 2',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(2,3)),y(node_connect(2,3)),'Node 3',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(2,4)),y(node_connect(2,4)),'Node 4',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(2,4)),y(node_connect(2,4)),'Node 5',...
     'HorizontalAlignment', 'center')
 
   text(x(node_connect(2,1))/2,y(node_connect(2,1))/2,'Element 1',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 2',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 3',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 4',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 5',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 6',...
    'HorizontalAlignment', 'center')
   hold on
 
   mult_factor = 3000;
   d=d*mult_factor;
 
   position_d(:, 1) = [ 0; 0];
   position_d(:, 2) = [ L(1)+d(3); 0+d(4)];
   position_d(:, 3) = [ 0; L(4)];
   position_d(:, 4) = [ L(1) ; L(4) ];
   position_d(:, 5) = [ L(5)*cos(theta(5))+d(9) ; L(5)*sin(theta(5))+d(10) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position_d(1,i);
      y(i) = position_d(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect_d(1, 1) = 1;   % element 1
   node_connect_d(2, 1) = 2;
 
   node_connect_d(1, 2) = 2;   % element 2
   node_connect_d(2, 2) = 5;
 
   node_connect_d(1, 3) = 5;   % element 3
   node_connect_d(2, 3) = 3;
 
   node_connect_d(1, 4) = 2;   % element 4
   node_connect_d(2, 4) = 4;
 
   node_connect_d(1, 5) = 1;   % element 5
   node_connect_d(2, 5) = 5;
 
   node_connect_d(1, 6) = 5;   % element 6
   node_connect_d(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem 
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect_d(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect_d(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 5000 -1000 4000])
      % solid line
      plot(xx,yy,'-r')
      % hold the current figure for the next element
      hold on
   end
 
 
   text(x(node_connect_d(1,1)),y(node_connect_d(1,1)),'Node 1',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(1,2)),y(node_connect_d(1,2)),'Node 2',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(2,3)),y(node_connect_d(2,3)),'Node 3',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(2,4)),y(node_connect_d(2,4)),'Node 4',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(2,4)),y(node_connect_d(2,4)),'Node 5',...
     'HorizontalAlignment', 'center')
 
   text(x(node_connect_d(2,1))/2,y(node_connect_d(2,1))/2,'Element 1',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 2',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 3',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 4',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 5',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 6',...
    'HorizontalAlignment', 'center')
 
 
   % put the title on the figure
   title('Deformed and Undeformed States of 6-Bar Truss System')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

6-Bar Truss System Where Young's Modulus Values Are Not Equal [edit]

The properties of the same 6-bar truss system analyzed above but different and varying values for Young's Modulus are shown in the table below.

**Properties of 6-Bar Truss System with Different Values for Young's Modulus**
Element Number	Young's Modulus (GPa)	Cross-Sectional Area (m²)	Length (m)	Angle of Inclination
1	150	0.001	4	$0^\circ$
2	180	0.001	2.83	$135^\circ$
3	200	0.001	2.24	$153.4^\circ$
4	200	0.001	3	$90^\circ$
5	220	0.001	2.83	$45^\circ$
6	250	0.001	2.24	$26.6^\circ$

MATLAB Code [edit]

The following MATLAB code was provided as an online resource from FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti.

% Six bar truss example
e = [150*10^3;180*10^3;200*10^3;200*10^3;220*10^3;250*10^3]; A = 0.001*1000^2; P = 20000.; 
alpha = pi/6; 
nodes = 1000*[0, 0; 4, 0; 0, 3; 4, 3; 2, 2]; 
dof=2*length(nodes);
conn=[1,2; 2,5; 5,3; 2,4; 1,5; 5,4];
lmm = [1, 2, 3, 4; 3, 4, 9, 10; 9, 10, 5, 6;
    3, 4, 7, 8; 1, 2, 9, 10; 9, 10, 7, 8];
elems=size(lmm,1);
K=zeros(dof); R = zeros(dof,1);
debc = [1, 2, 5, 6, 7, 8];
ebcVals = zeros(length(debc),1);
 
%load vector
R = zeros(dof,1); R(3) = P*sin(alpha); R(4) = P*cos(alpha); 
 
% Assemble global stiffness matrix
K=zeros(dof);
for i=1:elems
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(e(i), A, 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(i), A, ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
format short g
results

In order for this code to run properly, it needs three additional MATLAB functions. These functions are listed and shown in detail below.

PlaneTrussElement.m [edit]

function k = PlaneTrussElement(e, A, coord)
% PlaneTrussElement(e, A, coord)
% Generates stiffness matrix for a plane truss element
% e = modulus of elasticity
% A = area of cross-section
% coord = coordinates at the element ends
 
x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms;
    ls*ms, ms^2,-ls*ms,-ms^2;
    -ls^2,-ls*ms,ls^2,ls*ms;
    -ls*ms,-ms^2,ls*ms,ms^2];

NodalSoln.m [edit]

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

PlaneTrussResults.m [edit]

function results = PlaneTrussResults(e, A, coord, disps)
% results = PlaneTrussResults(e, A, coord, disps)
% Compute plane truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
% The output quantities are eps = axial strain
% sigma = axial stress and force = axial force.
 
x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
T=[ls,ms,0,0; 0,0,ls,ms];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;
results=[eps, sigma, force];

Plot of Undeformed and Deformed Structure [edit]

The plot of the undeformed and deformed shapes of the 6-bar truss system with the same values for Young's Modulus and cross-sectional area is shown below with a multiplication factor on the deformation of 3000.

Plot of undeformed and deformed shapes of the 6-bar truss system with different Young's Modulus values analyzed using the Finite Element Method

The source code added to the original MATLAB code to produce the above plot is given below.

%-----------------------------------------------------------------
% model with 2-D beam elements with 2 dofs per node
   dof = 2;         %  dof per node
 
%
% input the nodal coordinates
%
   n_node = 5;              % total number of nodes
   n_elem = 6;              % total number of elements
   total_dof = dof * n_node;  % total dofs of system
 
% five-bar truss system data
%
   for i=1:n_elem
        con=conn(i,:);
        L(i)=PlaneTrussElementMod(nodes(con,:));
   end
   for i=1:n_elem
        con=conn(i,:);
        theta(i)=PlaneTrussElementMod2(nodes(con,:));
   end
 
   position(:, 1) = [ 0; 0];
   position(:, 2) = [ L(1); 0];
   position(:, 3) = [ 0; L(4)];
   position(:, 4) = [ L(1) ; L(4) ];
   position(:, 5) = [ L(5)*cos(theta(5)) ; L(5)*sin(theta(5)) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position(1,i);
      y(i) = position(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect(1, 1) = 1;   % element 1
   node_connect(2, 1) = 2;
 
   node_connect(1, 2) = 2;   % element 2
   node_connect(2, 2) = 5;
 
   node_connect(1, 3) = 5;   % element 3
   node_connect(2, 3) = 3;
 
   node_connect(1, 4) = 2;   % element 4
   node_connect(2, 4) = 4;
 
   node_connect(1, 5) = 1;   % element 5
   node_connect(2, 5) = 5;
 
   node_connect(1, 6) = 5;   % element 6
   node_connect(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem  
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 10000 -1000 4000])
      % solid line
      plot(xx,yy,'--')
      % hold the current figure for the next element
      hold on
   end
 
 
   text(x(node_connect(1,1)),y(node_connect(1,1)),'Node 1',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(1,2)),y(node_connect(1,2)),'Node 2',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(2,3)),y(node_connect(2,3)),'Node 3',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(2,4)),y(node_connect(2,4)),'Node 4',...
     'HorizontalAlignment', 'center')
   text(x(node_connect(2,4)),y(node_connect(2,4)),'Node 5',...
     'HorizontalAlignment', 'center')
 
   text(x(node_connect(2,1))/2,y(node_connect(2,1))/2,'Element 1',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 2',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 3',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 4',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 5',...
    'HorizontalAlignment', 'center')
   text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2 , y(node_connect(2,1)) + (y(node_connect(2,2))-y(node_connect(2,1)))/2 ,'Element 6',...
    'HorizontalAlignment', 'center')
   hold on
 
   mult_factor = 3000;
   d=d*mult_factor;
 
   position_d(:, 1) = [ 0; 0];
   position_d(:, 2) = [ L(1)+d(3); 0+d(4)];
   position_d(:, 3) = [ 0; L(4)];
   position_d(:, 4) = [ L(1) ; L(4) ];
   position_d(:, 5) = [ L(5)*cos(theta(5))+d(9) ; L(5)*sin(theta(5))+d(10) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position_d(1,i);
      y(i) = position_d(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect_d(1, 1) = 1;   % element 1
   node_connect_d(2, 1) = 2;
 
   node_connect_d(1, 2) = 2;   % element 2
   node_connect_d(2, 2) = 5;
 
   node_connect_d(1, 3) = 5;   % element 3
   node_connect_d(2, 3) = 3;
 
   node_connect_d(1, 4) = 2;   % element 4
   node_connect_d(2, 4) = 4;
 
   node_connect_d(1, 5) = 1;   % element 5
   node_connect_d(2, 5) = 5;
 
   node_connect_d(1, 6) = 5;   % element 6
   node_connect_d(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem 
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect_d(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect_d(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 10000 -1000 4000])
      % solid line
      plot(xx,yy,'-r')
      % hold the current figure for the next element
      hold on
   end
 
 
   text(x(node_connect_d(1,1)),y(node_connect_d(1,1)),'Node 1',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(1,2)),y(node_connect_d(1,2)),'Node 2',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(2,3)),y(node_connect_d(2,3)),'Node 3',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(2,4)),y(node_connect_d(2,4)),'Node 4',...
     'HorizontalAlignment', 'center')
   text(x(node_connect_d(2,4)),y(node_connect_d(2,4)),'Node 5',...
     'HorizontalAlignment', 'center')
 
   text(x(node_connect_d(2,1))/2,y(node_connect_d(2,1))/2,'Element 1',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 2',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 3',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 4',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 5',...
    'HorizontalAlignment', 'center')
   text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1)) + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 6',...
    'HorizontalAlignment', 'center')
 
 
   % put the title on the figure
   title('Deformed and Undeformed States of 6-Bar Truss System')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

Plot Comparing Deformed Structures of Both 6-Bar Truss Systems [edit]

The following plot shows a comparison of the deformed structures of the two 6-bar truss systems analyzed previously using the Finite Element Method. The blue dotted line represents the undeformed structure of the 6-bar truss system, the solid red line represents the deformed structure of the 6-bar truss system with the same values of Young's Modulus, and the solid green line represents the deformed structure of the 6-bar truss system with different values of Young's Modulus.

Plot comparing the deformed structures of the two 6-bar truss systems analyzed using the Finite Element Method

The MATLAB code used to produce the above plot is shown below.

% model with 2-D beam elements with 2 dofs per node
   dof = 2;         %  dof per node
 
% input the nodal coordinates
 
   n_node = 5;              % total number of nodes
   n_elem = 6;              % total number of elements
   total_dof = dof * n_node;  % total dofs of system
 
% five-bar truss system data
 
   for i=1:n_elem
        con=conn(i,:);
        L(i)=PlaneTrussElementMod(nodes(con,:));
   end
   for i=1:n_elem
        con=conn(i,:);
        theta(i)=PlaneTrussElementMod2(nodes(con,:));
   end
 
mult_factor = 3000;
 
% displacement dofs of system with same values for Young's Modulus
d = mult_factor*[ 0; 0; 0.21311; 0.24998; 0; 0; 0; 0; -0.0060971; 0.012242]
 
% displacement dofs of system with different values for Young's Modulus 
d_diff_e = mult_factor*[ 0; 0; 0.26485; 0.26083; 0; 0; 0; 0; 0.00063864; -0.001246]
 
% plot the undeformed shape
   position(:, 1) = [ 0; 0];
   position(:, 2) = [ L(1); 0];
   position(:, 3) = [ 0; L(4)];
   position(:, 4) = [ L(1) ; L(4) ];
   position(:, 5) = [ L(5)*cos(theta(5)) ; L(5)*sin(theta(5)) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position(1,i);
      y(i) = position(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect(1, 1) = 1;   % element 1
   node_connect(2, 1) = 2;
 
   node_connect(1, 2) = 2;   % element 2
   node_connect(2, 2) = 5;
 
   node_connect(1, 3) = 5;   % element 3
   node_connect(2, 3) = 3;
 
   node_connect(1, 4) = 2;   % element 4
   node_connect(2, 4) = 4;
 
   node_connect(1, 5) = 1;   % element 5
   node_connect(2, 5) = 5;
 
   node_connect(1, 6) = 5;   % element 6
   node_connect(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem  
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 6000 -1000 4000])
      % solid line
      plot(xx,yy,'--')
      % hold the current figure for the next element
      hold on
   end
 
% plot the deformed shape of the system with same Young's Modulus values
   position_d(:, 1) = [ 0; 0];
   position_d(:, 2) = [ L(1)+d(3); 0+d(4)];
   position_d(:, 3) = [ 0; L(4)];
   position_d(:, 4) = [ L(1) ; L(4) ];
   position_d(:, 5) = [ L(5)*cos(theta(5))+d(9) ; L(5)*sin(theta(5))+d(10) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position_d(1,i);
      y(i) = position_d(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect_d(1, 1) = 1;   % element 1
   node_connect_d(2, 1) = 2;
 
   node_connect_d(1, 2) = 2;   % element 2
   node_connect_d(2, 2) = 5;
 
   node_connect_d(1, 3) = 5;   % element 3
   node_connect_d(2, 3) = 3;
 
   node_connect_d(1, 4) = 2;   % element 4
   node_connect_d(2, 4) = 4;
 
   node_connect_d(1, 5) = 1;   % element 5
   node_connect_d(2, 5) = 5;
 
   node_connect_d(1, 6) = 5;   % element 6
   node_connect_d(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem 
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect_d(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect_d(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 6000 -1000 4000])
      % solid line
      plot(xx,yy,'-r')
      % hold the current figure for the next element
      hold on
   end
 
% plot deformed shape for the system with different Young's Modulus values
   position_d(:, 1) = [ 0; 0];
   position_d(:, 2) = [ L(1)+d_diff_e(3); 0+d_diff_e(4)];
   position_d(:, 3) = [ 0; L(4)];
   position_d(:, 4) = [ L(1) ; L(4) ];
   position_d(:, 5) = [ L(5)*cos(theta(5))+d_diff_e(9) ; L(5)*sin(theta(5))+d_diff_e(10) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position_d(1,i);
      y(i) = position_d(2,i);
   end
 
% set up the element connectivity array node_connect, defined as
% follows:
% node_connect(local node number, element number) = global node number
 
   node_connect_d_diff_e(1, 1) = 1;   % element 1
   node_connect_d_diff_e(2, 1) = 2;
 
   node_connect_d_diff_e(1, 2) = 2;   % element 2
   node_connect_d_diff_e(2, 2) = 5;
 
   node_connect_d_diff_e(1, 3) = 5;   % element 3
   node_connect_d_diff_e(2, 3) = 3;
 
   node_connect_d_diff_e(1, 4) = 2;   % element 4
   node_connect_d_diff_e(2, 4) = 4;
 
   node_connect_d_diff_e(1, 5) = 1;   % element 5
   node_connect_d_diff_e(2, 5) = 5;
 
   node_connect_d_diff_e(1, 6) = 5;   % element 6
   node_connect_d_diff_e(2, 6) = 4;
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem 
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect_d_diff_e(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect_d_diff_e(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-1000 6000 -1000 4000])
      % solid line
      plot(xx,yy,'-g')
      % hold the current figure for the next element
      hold on
   end
 
   % put the title on the figure
   title('Comparison of Deformed States of 6-Bar Truss Systems')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

Finite Element Analysis of 3-Bar Space Truss [edit]

The setup of the 3-bar space truss system as given on page 230 of FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti is shown in the figure below.

Setup of 3-bar space truss system to be solved using Finite Element Method

The element and node properties of the 3-bar space truss system are shown in the following two tables.

**Element Properties of 3-Bar Space Truss System**
Element Number	Young's Modulus (GPa)	Cross-Sectional Area (mm²)
1	200	200
2	200	200
3	200	600

**Node Properties of 3-Bar Space Truss System**
Node Number	x (m)	y (m)	z (m)
1	0.96	1.92	0
2	-1.44	1.44	0
3	0	0	0
4	0	0	2

In order to properly create MATLAB code to solve the 3-bar space truss system using the Finite Element Method, the conn and lmm arrays must be populated. The proper assembly of the conn array is shown below.

$\texttt{conn}=\left[ \begin{array}{cc} 1 & 4 \\ 2 & 4 \\ 3 & 4 \end{array} \right]$

The proper assembly of the lmm array is shown below.

$\texttt{lmm}=\left[ \begin{array}{cccccc} 1 & 2 & 3 & 10 & 11 & 12 \\ 4 & 5 & 6 & 10 & 11 & 12 \\ 7 & 8 & 9 & 10 & 11 & 12 \end{array} \right]$

Analysis 3-Bar Space Truss Using MATLAB and the Finite Element Method [edit]

MATLAB Code [edit]

The MATLAB code, shown below, and used to solve the 3-bar space truss system was provided as an online resource by FUNDAMENTAL Finite Element Analysis and Applications with Mathematica and MATLAB Computations, written by M. Asghar Bhatti.

% Three bar space truss example
a1 = 200; a2 = 600; e = 200000; P = 20000;
nodes = 1000*[.96, 1.92, 0; -1.44, 1.44, 0; 0, 0, 0; 0, 0, 2]; 
dof=3*length(nodes);
conn=[1,4; 2,4; 3,4];
lmm = [1, 2, 3, 10, 11, 12; 
    4, 5, 6, 10, 11, 12; 
    7, 8, 9, 10, 11, 12];
debc = [1:9];
ebcVals = zeros(length(debc),1);
 
%load vector
R = zeros(dof,1); R(11) = -P; 
 
% Assemble global stiffness matrix
K=zeros(dof);
for i=1:2
    lm=lmm(i,:);
    con=conn(i,:);
    k=SpaceTrussElement(e, a1, nodes(con,:));
    K(lm, lm) = K(lm, lm) + k;
end
lm=lmm(3,:);
con=conn(3,:);
k=SpaceTrussElement(e, a2, nodes(con,:));
K(lm, lm) = K(lm, lm) + k
 
% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, debc, ebcVals)
results=[];
for i=1:2
    results = [results; SpaceTrussResults(e, a1, ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
results = [results; SpaceTrussResults(e, a2, ...
        nodes(conn(3,:),:), d(lmm(3,:)))];
format short g
results

The following additional functions were called in the above code and were necessary to properly solve the 3-bar truss system using MATLAB.

SpaceTrussElement.m [edit]

function k = SpaceTrussElement(e, A, coord)
% k = SpaceTrussElement(e, A, coord)
% Generates stiffness matrix of a space truss element
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
 
x1=coord(1,1); y1=coord(1,2); z1=coord(1,3);
x2=coord(2,1); y2=coord(2,2); z2=coord(2,3);
L=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L; ns=(z2-z1)/L;
k = e*A/L*[ls^2, ls*ms, ls*ns, -ls^2, -(ls*ms), -(ls*ns);
    ls*ms, ms^2, ms*ns, -(ls*ms), -ms^2, -(ms*ns);
    ls*ns, ms*ns, ns^2, -(ls*ns), -(ms*ns), -ns^2;
    -ls^2, -(ls*ms), -(ls*ns), ls^2, ls*ms, ls*ns;
    -(ls*ms), -ms^2, -(ms*ns), ls*ms, ms^2, ms*ns;
    -(ls*ns), -(ms*ns), -ns^2, ls*ns, ms*ns, ns^2];

NodalSoln.m [edit]

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

SpaceTrussResults.m [edit]

function results = SpaceTrussResults(e, A, coord, disps)
% results = SpaceTrussResults(e, A, coord, disps)
% Compute space truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
 
x1=coord(1,1); y1=coord(1,2); z1=coord(1,3);
x2=coord(2,1); y2=coord(2,2); z2=coord(2,3);
L=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L; ns=(z2-z1)/L;
T=[ls,ms,ns,0,0,0; 0,0,0,ls,ms,ns];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;
results=[eps, sigma, force];

MATLAB Results [edit]

The output in the MATLAB command window after running the above code is shown below.

K =

 1.0e+004 *

 Columns 1 through 11

   0.1460    0.2919   -0.3041         0         0         0         0         0         0   -0.1460   -0.2919
   0.2919    0.5839   -0.6082         0         0         0         0         0         0   -0.2919   -0.5839
  -0.3041   -0.6082    0.6335         0         0         0         0         0         0    0.3041    0.6082
        0         0         0    0.3567   -0.3567    0.4954         0         0         0   -0.3567    0.3567
        0         0         0   -0.3567    0.3567   -0.4954         0         0         0    0.3567   -0.3567
        0         0         0    0.4954   -0.4954    0.6880         0         0         0   -0.4954    0.4954
        0         0         0         0         0         0         0         0         0         0         0
        0         0         0         0         0         0         0         0         0         0         0
        0         0         0         0         0         0         0         0    6.0000         0         0
  -0.1460   -0.2919    0.3041   -0.3567    0.3567   -0.4954         0         0         0    0.5026   -0.0647
  -0.2919   -0.5839    0.6082    0.3567   -0.3567    0.4954         0         0         0   -0.0647    0.9405
   0.3041    0.6082   -0.6335   -0.4954    0.4954   -0.6880         0         0   -6.0000    0.1913   -1.1036

 Column 12

   0.3041
   0.6082
  -0.6335
  -0.4954
   0.4954
  -0.6880
        0
        0
  -6.0000
   0.1913
  -1.1036
   7.3216


d =

        0
        0
        0
        0
        0
        0
        0
        0
        0
  -0.1871
  -2.5920
  -0.3858

 
reactions =
 
 1.0e+004 *

   0.6667
   1.3333
  -1.3889
  -0.6667
   0.6667
  -0.9259
        0
        0
   2.3148


results =

  0.00050936       101.87        20375
  0.00033036       66.072        13214
  -0.0001929       -38.58       -23148

The computed reactions of the 3-bar space truss system are shown below.

$\mathbf{R}= \left[ \begin{array}{c} R_{1} \\ R_{2} \\ R_{3} \\ R_{4} \\ R_{5} \\ R_{6} \\ R_{7} \\ R_{8} \\ R_{9} \\ \end{array} \right] = \left[ \begin{array}{c} 6667 \\ 13333 \\ -13889 \\ -6667 \\ 6667 \\ -9259 \\ 0 \\ 0 \\ 23148 \\ \end{array} \right]$

**Results for 2-Bar Truss System**
Element Number	Axial Strain	Axial Stress	Axial Force
1	0.00050936	101.87	20375
2	0.00033036	66.072	13214
3	-0.0001929	-38.58	-23148

Determination of Statically Determinate or Indeterminate System [edit]

The 3-bar space truss system analyzed above is a statically determinate system. This is because each of the three members of the system have an axial force directed along the member, which gives the system 3 unknown forces. Because the space truss system has 3 displacement DOFs at each node, the summation of forces can be performed in the x-direction, y-direction, and z-direction. This gives 3 unique equations. Because there are 3 equations and 3 unknowns, this system can be solved using only the methods of statics.

In order to verify that the results obtained from the MATLAB code shown above, the calculation of axial force, axial stress, and axial strain using statics is shown below for the 3-bar space truss system.

Three unique equations were obtained by summing forces in the x, y, and z directions, as shown below.

$\sum{F_{x}}= R_{1}l^{(1)} + R_{2}l^{(2)} + R_{3}l^{(3)} = 0$

$\sum{F_{y}}= R_{1}m^{(1)} + R_{2}m^{(2)} + R_{3}m^{(3)} = -P$

$\sum{F_{z}}= R_{1}n^{(1)} + R_{2}n^{(2)} + R_{3}n^{(3)} = 0$

$R_{1} = 20375$

$R_{2} = 13214$

$R_{3} = -23149$

Element 1 [edit]

Knowing the axial force in element 1, the axial stress and axial strain can be calculated as shown below.

$R_{1} = 20375$

$\sigma_{1} = \frac{R_{1}}{A_{1}} = \frac{20375}{200} = 101.875$

$\epsilon_{1} = \frac{\sigma_{1}}{E_{1}} = \frac{101.875}{200\times10^{3}} = 5.094\times10^{-4}$

Element 2 [edit]

Knowing the axial force in element 2, the axial stress and axial strain can be calculated as shown below.

$R_{2} = 13214$

$\sigma_{2} = \frac{R_{2}}{A_{2}} = \frac{13214}{200} = 66.07$

$\epsilon_{2} = \frac{\sigma_{2}}{E_{2}} = \frac{66.07}{200\times10^{3}} = 3.3035\times10^{-4}$

Element 3 [edit]

Knowing the axial force in element 3, the axial stress and axial strain can be calculated as shown below.

$R_{3} = -23149$

$\sigma_{3} = \frac{R_{3}}{A_{3}} = \frac{-23149}{600} = -38.58$

$\epsilon_{3} = \frac{\sigma_{3}}{E_{3}} = \frac{-38.58}{200\times10^{3}} = -1.929\times10^{-4}$

Comparing these values with the values obtained from the Finite Element Method solution above, it can be seen that both methods are in complete agreement.

Plots of Undeformed and Deformed Structure using MATLAB [edit]

The following sections show plots of the undeformed and deformed shapes of the 3-bar space truss system described above. These plots were created using the following MATLAB code, which was added to the end of the MATLAB code used to solve the 3-bar space truss system using the Finite Element Method.

%-----------------------------------------------------------------
% model with 2-D beam elements with 2 dofs per node
   dof = 3;         %  dof per node
 
%
% input the nodal coordinates
%
   n_node = 4;              % total number of nodes
   n_elem = 3;              % total number of elements
   total_dof = dof * n_node;  % total dofs of system
 
% five-bar truss system data
%
   for i=1:n_elem
        con=conn(i,:);
        [L(i) ls(i) ms(i) ns(i)] = SpaceTrussElementMod(nodes(con,:));
   end
 
   position(:, 1) = 1000*[ 0.96 ; 1.92 ; 0 ];
   position(:, 2) = 1000*[ -1.44 ; 1.44 ; 0 ];
   position(:, 3) = 1000*[ 0 ; 0 ; 0 ];
   position(:, 4) = 1000*[ 0 ; 0 ; 2 ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position(1,i);
      y(i) = position(2,i);
      z(i) = position(3,i);
   end
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem  
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = conn(i,1);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = conn(i,2);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % zz : 1x2 array containing z coordinates of node_1 and node_2
      zz = [z(node_1),z(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-2000 1500 -3000 2500 -500 2500])
      % solid line
      view([1000,1000,1000])
      plot3(xx,yy,zz,'--')
      % hold the current figure for the next element
      hold on
   end
 
 
   text(x(conn(1,1)),y(conn(1,1)),z(conn(1,1)),'Node 1',...
     'HorizontalAlignment', 'center')
   text(x(conn(2,1)),y(conn(2,1)),z(conn(2,1)),'Node 2',...
     'HorizontalAlignment', 'center')
   text(x(conn(3,1)),y(conn(3,1)),z(conn(3,1)),'Node 3',...
     'HorizontalAlignment', 'center')
   text(x(conn(3,2)),y(conn(3,2)),z(conn(3,2)),'Node 4',...
     'HorizontalAlignment', 'center')
 
   text(x(conn(1,2))/2,y(conn(1,2))/2,z(conn(1,2))/2,'Element 1',...
    'HorizontalAlignment', 'center')
   text( x(conn(1,2)) + (x(conn(2,2))-x(conn(1,2)))/2 , y(conn(1,2)) + (y(conn(2,2))-y(conn(1,2)))/2, z(conn(1,2)) + (z(conn(2,2))-z(conn(1,2)))/2 ,'Element 2',...
    'HorizontalAlignment', 'center')
   text( x(conn(1,2)) + (x(conn(2,2))-x(conn(1,2)))/2 , y(conn(1,2)) + (y(conn(2,2))-y(conn(1,2)))/2, z(conn(1,2)) + (z(conn(2,2))-z(conn(1,2)))/2 ,'Element 3',...
    'HorizontalAlignment', 'center')
 
   hold on
 
   mult_factor = 1;
   d=d*mult_factor;
 
   position_d(:, 1) = 1000*[ 0.96 ; 1.92 ; 0 ];
   position_d(:, 2) = 1000*[ -1.44 ; 1.44 ; 0 ];
   position_d(:, 3) = 1000*[ 0 ; 0 ; 0 ];
   position_d(:, 4) = 1000*[ 0+d(10) ; 0+d(11) ; 2+d(12) ];
 
   % set up the nodal coordinate arrays x, y
   for i = 1 : n_node
      x(i) = position_d(1,i);
      y(i) = position_d(2,i);
      z(i) = position_d(3,i);
   end
 
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem  
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = conn(i,1);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = conn(i,2);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % zz : 1x2 array containing z coordinates of node_1 and node_2
      zz = [z(node_1),z(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([-2000 1500 -3000 2500 -500 2500])
      % solid line
      view([1000,1000,1000])
      plot3(xx,yy,zz,'-r')
      % hold the current figure for the next element
      hold on
   end
 
   % put the title on the figure
   title('Deformed and Undeformed States of 3-Bar Space Truss System')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')
   % label z axis
   zlabel('z')

Plane View Down x-axis [edit]

Plot of the undeformed and deformed shapes of the 3-bar space truss system from a plane view down the x-axis

Plane View Down y-axis [edit]

Plot of the undeformed and deformed shapes of the 3-bar space truss system from a plane view down the y-axis

Plane View Down z-axis [edit]

Plot of the undeformed and deformed shapes of the 3-bar space truss system from a plane view down the z-axis

Perspective Views [edit]

A perspective view of the undeformed and deformed shapes of the 3-bar space truss analyzed using the Finite Element Method is shown below. This perspective view is from the coordinates: x = -2 m, y = -2 m, z = 3 m

Plot of the undeformed and deformed shapes of the 3-bar space truss system from a perspective view with coordinates of x = -2 m, y = -2 m, z = 3 m

In order to obtain a more clear image of the deformation of the 3-bar space truss system, a perspective view from the coordinates of x = 1 m, y = 1 m, and z = 1 m is shown below.

Plot of the undeformed and deformed shapes of the 3-bar space truss system from a perspective view with coordinates of x = 1 m, y = 1 m, z = 1 m

Contributing Team Members [edit]

Andrew McDonald - Eml4500.f08.bike.mcdonald 21:25, 2 November 2008 (UTC)
Garrett Pataky - Eml4500.f08.bike.pataky 20:25, 3 November 2008 (UTC)
Sam Bernal - Eml4500.f08.bike.bernal 03:32, 7 November 2008 (UTC)
Bobby Sweeting - Eml4500.f08.bike.sweeting 22:10, 2 November 2008 (UTC)
Shawn Gravois - Eml4500.f08.bike.gravois 06:09, 5 November 2008 (UTC)
Eric Viale - Eml4500.f08.bike.viale 21:44, 6 November 2008 (UTC)

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

Contents

Class Notes [edit]

Principal of Virtual Work [edit]

Analysis of The Axial Force Displacement Relation Using The Principle of Virtual Work. [edit]

Principle of Virtual Work Conversion to Force Displacement Equation [edit]

Equation of Motion, Constitutive Relation, PDE of Motion, and Boundary Condition Derivation [edit]

Derivation of Matrix Relations and Matrices for Finite Element Analysis of 3-D Systems [edit]

Axial FD Relation [edit]

Element DOFs and Element Forces in Global xyz Coordinates [edit]

Transformation Matrix [edit]

Relation between Axial DOFs and Element DOFs in Global Coordinates [edit]

Relation between Axial Forces and Element Forces in Global Coordinates [edit]

Derivations Using the Principal of Virtual Work [edit]

Global FD Relation [edit]

Element Stiffness Matrix in Global Coordinates in Terms of Transformation Matrix and Element Stiffness Matrix from Axial FD Relation [edit]

Element Forces in Global Coordinates in Terms of Transformation Matrix and Element Axial Forces [edit]

Element Stiffness Matrix in Global Coordinates in Terms of Director Cosines [edit]

Rerunning the MATLAB Code for 2-Bar Truss System and Interpreting results Array [edit]

Modified MATLAB Code [edit]

Interpreting results Array [edit]

Element 1 [edit]

Element 2 [edit]

Finite Element Analysis of Various 6-Bar Truss Systems [edit]

6-Bar Truss System With Same Young's Modulus and Cross-Sectional Area Values [edit]

MATLAB Code [edit]

PlaneTrussElement.m [edit]

NodalSoln.m [edit]

PlaneTrussResults.m [edit]

Plot of Undeformed and Deformed Structure [edit]

6-Bar Truss System Where Young's Modulus Values Are Not Equal [edit]

MATLAB Code [edit]

PlaneTrussElement.m [edit]

NodalSoln.m [edit]

PlaneTrussResults.m [edit]

Plot of Undeformed and Deformed Structure [edit]

Plot Comparing Deformed Structures of Both 6-Bar Truss Systems [edit]

Finite Element Analysis of 3-Bar Space Truss [edit]

Analysis 3-Bar Space Truss Using MATLAB and the Finite Element Method [edit]

MATLAB Code [edit]

SpaceTrussElement.m [edit]

NodalSoln.m [edit]

SpaceTrussResults.m [edit]

MATLAB Results [edit]

Determination of Statically Determinate or Indeterminate System [edit]

Element 1 [edit]

Element 2 [edit]

Element 3 [edit]

Plots of Undeformed and Deformed Structure using MATLAB [edit]

Plane View Down x-axis [edit]

Plane View Down y-axis [edit]

Plane View Down z-axis [edit]

Perspective Views [edit]

Contributing Team Members [edit]

Navigation menu

Search

Rerunning the MATLAB Code for 2-Bar Truss System and Interpreting `results` Array [edit]

Interpreting `results` Array [edit]