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

Class Notes

Connectivity and Location Matrix Master Arrays

The elements of a two bar truss system are shown below by Figures 1 and 2.

File:BikeElem1.jpg

Figure 1 - Element 1

File:BikeElem2.jpg

Figure 2 - Element 2

The connectivity array, denoted "conn," is shown below for the two bar truss system.

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

Each row of the connectivity array represents an element(i.e. row 1 = element 1, row 2 = element 2), and each column represents a local node number(i.e. column 1 = local node 1, column 2 = local node 2). The numbers within the connectivity array correspond to the global node numbers, such that conn(e,j) = global node number of local node j of element e.

The location matrix master array, denoted "lmm," is shown below for the two bar truss system.

lmm = $\left[ \begin{array}{cccc} {1} & {2} & {3} & {4} \\ {3} & {4} & {5} & {6} \end{array} \right]$

Somewhat similar to the connectivity array, the location matrix master array has rows that correspond to the element number and columns that represent the local degree of freedom(dof) number. The numbers within the location matrix master array correspond to the global dof number, or equation number, in the matrix K. Generally, lmm(i,j) = global dof number(or equation number) for the element stiffness coefficient corresponding to the jth local dof number.

Method 2 to derive the stiffness matrix $K_{4x4}^{e}$ is to transform a system with 2 dof's into a system with 4 dof's so that the transformed matrix is a 4x4, and hopefully invertible.

Transforming Local DOFs to a New Set of Local DOFs

The goal of the following steps is to find a transform matrix, $\tilde{\textbf{T}}^{(e)}_{4x4}$, that will change the elemental degrees of freedom matrix, $\textbf{d}^{(e)}_{4x1}$, into a new set of elemental degrees of freedom, $\tilde{\textbf{d}}^{(e)}_{4x1}$ so that the transform matrix is invertible.

$\tilde{\textbf{d}}^{(e)}_{4x1} = \tilde{\textbf{T}}^{(e)}_{4x4}\textbf{d}^{(e)}_{4x1}$

The figure below shows the current directions of the elemental degrees of freedom and the desired directions of the new degrees of freedom.

File:19.1.jpg

Local free body diagram of element e and two sets of elemental degrees of freedom.

It can be observed that the new elemental degrees of freedom are equal to the axial degrees of freedom.

$q_1^{(e)} = \tilde{d}_1^{(e)}$ and $q_2^{(e)} = \tilde{d}_2^{(e)}$

To determine which director cosines to use in determining the new elemental degree of freedom, the original degree of freedom must be computed along the $\tilde{x}$-axis.

$\tilde{d}_1^{(e)} = \textbf{d}_1^{(e)} \cdot \tilde{\textbf{i}} = (d_1^{(e)} \tilde{i} + d_2^{(e)} \tilde{j}) \cdot \tilde{\textbf{i}} = \cos\theta^{(e)}d_1^{(e)} + \sin\theta^{(e)}d_2^{(e)}$

Replacing the director cosines with their respective symbols, the equation can be entered into matrix form.

$\tilde{d}_1^{(e)} = \begin{bmatrix}l^{(e)} & m^{(e)}\end{bmatrix}\begin{bmatrix}d_1^{(e)}\\d_2^{(e)}\end{bmatrix}$

The same procedure is done with the next elemental degree of freedom but it is computed along the $\tilde{y}$-axis.

$\tilde{d}_2^{(e)} = \textbf{d}_1^{(e)} \cdot \tilde{\textbf{j}} = (d_1^{(e)}\tilde{i} + d_2^{(e)}\tilde{j}) \cdot \tilde{\textbf{j}} = -\sin\theta^{(e)}d_1^{(e)} + \cos\theta^{(e)}d_2^{(e)}$

$\tilde{d}_2^{(e)} = \begin{bmatrix}-m^{(e)} & l^{(e)}\end{bmatrix}\begin{bmatrix}d_1^{(e)}\\d_2^{(e)}\end{bmatrix}$

The two matrices are then combined and the result is shown below.

$\begin{bmatrix}\tilde{d}_1^{(e)} \\ \tilde{d}_2^{(e)}\end{bmatrix} = \begin{bmatrix}l^{(e)} & m^{(e)} \\ -m^{(e)} & l^{(e)}\end{bmatrix} \begin{bmatrix}d_1^{(e)} \\ d_2^{(e)}\end{bmatrix}$

Completing the procedure for all four elemental degrees of freedom, we now have a matrix for transforming the elemental degrees of freedom into the elemental degrees of freedom along the $\tilde{x}$ and $\tilde{y}$ axes.

$\begin{bmatrix}\tilde{d}_1^{(e)} \\ \tilde{d}_2^{(e)} \\ \tilde{d}_3^{(e)} \\ \tilde{d}_4^{(e)}\end{bmatrix} = \begin{bmatrix}\textbf{R}^{(e)}_{2x2} & \textbf{0}_{2x2} \\ \textbf{0}_{2x2} & \textbf{R}^{(e)}_{2x2}\end{bmatrix} \begin{bmatrix}d_1^{(e)} \\ d_2^{(e)} \\ d_3^{(e)} \\ d_4^{(e)}\end{bmatrix}$

$\tilde{\textbf{d}}_{4x1}^{(e)} = \tilde{\textbf{T}}_{4x4}\textbf{d}_{4x1}^{(e)}$

Transformed Forces and DOFs

Now we can examine the element as if it was horizontal and determine the forces using the degrees of freedom along the $\tilde{x}$ and $\tilde{y}$ axes.

File:19.2.jpg

Local free body diagram of spring forces and degrees of freedom.

Transverse forces do not stretch the spring, thus the matrix is limited with a majority of zeroes.

$\tilde{\textbf{f}}^{(e)} = k^{(e)} \begin{bmatrix}1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}\tilde{\textbf{d}}^{(e)}$

The resulting equation is shown below.

$\tilde{\textbf{f}}^{(e)}_{4x1} = \tilde{\textbf{k}}^{(e)}_{4x4}\tilde{\textbf{d}}^{(e)}_{4x1}$

Next, we will consider the case when

$\tilde{d}_{4}^{(e)}\neq 0$

The other degrees of freedom we are going to consider to be zero.

$\tilde{d}_{1}^{(e)}=\tilde{d}_{2}^{(e)}=\tilde{d}_{3}^{(e)}=0$

Now, the force displacement relation for the case above is shown below.

$\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{k}}^{(e)}\mathbf{\tilde{d}}^{(e)}=\mathbf{0}$

Where $\mathbf{\tilde{f}}^{(e)}$ is a 4X1 matrix, $\mathbf{\tilde{k}}^{(e)}$ is a 4X4 matrix, $\mathbf{\tilde{d}}^{(e)}$ is a 4X1 matrix and $\mathbf{0}$ is a 4X1 matrix and the 4th column of $\mathbf{\tilde{k}}^{(e)}$.

As derived before, the following relation of the degrees of freedom is shown below.

$\mathbf{\tilde{d}}^{(e)}=\mathbf{\tilde{T}}^{(e)}\mathbf{d}^{(e)}$

Similarly, we can show the derivation for the following force relation.

$\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{T}}^{(e)}\mathbf{f}^{(e)}$

Assigned HW - Verify The Equation: $\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{T}}^{(e)}\mathbf{f}^{(e)}$

$\mathbf{\tilde{T}}^{(e)} = \begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0 \\ m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & -m^{(e)} \\ 0 & 0 & m^{(e)} & l^{(e)} \end{bmatrix}$ $\begin{Bmatrix}\ \tilde{f}_{1}^{(e)} \\ \tilde{f}_{2}^{(e)} \\ \tilde{f}_{3}^{(e)} \\ \tilde{f}_{4}^{(e)} \end{Bmatrix} = \begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0 \\ m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & -m^{(e)} \\ 0 & 0 & m^{(e)} & l^{(e)} \end{bmatrix} \begin{Bmatrix} f_{1}^{(e)} \\ f_{2}^{(e)} \\ f_{3}^{(e)} \\ f_{4}^{(e)} \end{Bmatrix}$ In order to verify the equation, the angle θ(e) = 0. This will mean that the forces in the $\tilde{x}$ will equal the forces in the $x$ direction. $\ l^{(e)}= 1$ $\ m^{(e)} =0$ $\begin{Bmatrix}\ \tilde{f}_{1}^{(e)} \\ \tilde{f}_{2}^{(e)} \\ \tilde{f}_{3}^{(e)} \\ \tilde{f}_{4}^{(e)} \end{Bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{Bmatrix} f_{1}^{(e)} \\ f_{2}^{(e)} \\ f_{3}^{(e)} \\ f_{4}^{(e)} \end{Bmatrix}$ In conclusion, $\tilde{f}_{1}^{(e)} = f_{1}^{(e)}$ $\tilde{f}_{2}^{(e)} = f_{2}^{(e)}$ $\tilde{f}_{3}^{(e)} = f_{3}^{(e)}$ $\tilde{f}_{4}^{(e)} = f_{4}^{(e)}$ The equation is verified since the forces are equal.

Looking once again at the force displacement relation

$\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{k}}^{(e)}\mathbf{\tilde{d}}^{(e)}$

The force displacement relation can be rewritten as shown, with the understanding that $\mathbf{\tilde{T}}^{(e)}$ has to be invertible.

$\left[\mathbf{\tilde{T}}^{(e)-1}\mathbf{\tilde{k}}^{(e)}\mathbf{\tilde{T}}^{(e)}\right]\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$

As shown earlier, $\mathbf{\tilde{T}}^{(e)}$ is a block diagonal matrix.

In order to see how to solve for the inverse of $\mathbf{\tilde{T}}^{(e)}$, we are going to consider the following general block diagram matrix.

$\mathbf{A}=\begin{bmatrix}\mathbf{D}_{1} & 0\\0 & \mathbf{D}_{s}\end{bmatrix}$

To show how to take the inverse of $\mathbf{A}$, a simpler example is shown here, which is a general diagonal matrix.

$\mathbf{B}=\begin{bmatrix}d_{11} & &\mathbf{0}\\ & d_{22}\\ \mathbf{0}& & d_{nn}\end{bmatrix}=Diag\left[d_{11},d_{22},...,d_{nn}\right]$

Assuming the following:

$d_{ii}\neq0$ for $i=1,2,...,n$,

The inverse of the matrix $\mathbf{B}$ is shown below

$\mathbf{B}^{-1}=Diag\left[\frac{1}{d_{11}},\frac{1}{d_{22}},...,\frac{1}{d_{nn}}\right]$.

Now, we are going back to analyzing the block diagram matrix $\mathbf{A}$.

$\mathbf{A}=Diag\left[\mathbf{D}_{1},...,\mathbf{D}_{s}\right]$.

$\mathbf{A}^{-1}=Diag\left[\mathbf{D}_{1}^{-1},...,\mathbf{D}_{s}^{-1}\right]$.

Now that the inverse of the general block diagonal matrix $\mathbf{A}$ was found, the same technique can be applied to the transfer matrix $\mathbf{\tilde{T}}^{(e)}$. The inverse is shown below.

$\mathbf{\tilde{T}}^{(e)-1}=Diag\left[\mathbf{R}^{(e)-1},\mathbf{R}^{(e)-1}\right]$.

The inverse of $\mathbf{R}^{(e)}$ in the equation above is equal to the transpose of $\mathbf{R}^{(e)}$. The proof is shown below.

The matrix $\mathbf{R}^{(e)}$ as well as its transpose are shown below.

$\mathbf{R}^{(e)}=\begin{bmatrix}l^{(e)} & m^{(e)} \\ -m^{(e)} & l^{(e)}\end{bmatrix}$

$\mathbf{R}^{(e)T}=\begin{bmatrix}l^{(e)} & -m^{(e)} \\ m^{(e)} & l^{(e)}\end{bmatrix}$

The multiplication of $\mathbf{R}^{(e)}$ and $\mathbf{R}^{(e)T}$ is shown below.

$\mathbf{R}^{(e)}\mathbf{R}^{(e)T}=\begin{bmatrix}l^{(e)} & m^{(e)} \\ -m^{(e)} & l^{(e)}\end{bmatrix}\begin{bmatrix}l^{(e)} & -m^{(e)} \\ m^{(e)} & l^{(e)}\end{bmatrix}=\begin{bmatrix}l^{(e)} + m^{(e)2}&0 \\ 0 & l^{(e)2} + m^{(e)2}\end{bmatrix}$

This means that the multiplication is equal to the identity matrix.

$\mathbf{R}^{(e)}\mathbf{R}^{(e)T}=\mathbf{I}$

Since it is equal to the identity matrix, the transpose has to be equal to the inverse of $\mathbf{R}^{(e)}.$

$\mathbf{R}^{(e)-1}=\mathbf{R}^{(e)T}$

Now that we know this, we can carry that over to the transfer matrix $\mathbf{\tilde{T}}^{(e)}$.

The inverse of $\mathbf{\tilde{T}}^{(e)}$ can be rewritten as shown.

$\mathbf{\tilde{T}}^{(e)-1} =Diag\left[\mathbf{R}^{(e)T},\mathbf{R}^{(e)T}\right]=\left(Diag\left[\mathbf{R}^{(e)},\mathbf{R}^{(e)}\right]\right)^{T}$.

Since $\mathbf{\tilde{T}}^{(e)}=Diag\left[\mathbf{R}^{(e)},\mathbf{R}^{(e)}\right]$, We can now state the following.

$\mathbf{\tilde{T}}^{(e)-1}=\mathbf{\tilde{T}}^{(e)T}$

Assigned HW: Verify the equation: $\mathbf k^{(e)} = \tilde \mathbf T^{(e)^{T}}\tilde \mathbf k^{(e)}\tilde \mathbf T^{(e)}$

$\mathbf k^{(e)} = \tilde \mathbf T^{(e)^{T}}\tilde \mathbf k^{(e)}\tilde \mathbf T^{(e)}$ By definition, $\mathbf k^{(e)} = \underset{}{k^{(e)}} \begin{bmatrix} (l^{(e)})^{2} & l^{(e)}m^{(e)} & -(l^{(e)})^{2} & -l^{(e)}m^{(e)}\\ l^{(e)}m^{(e)} & (m^{(e)})^{2} & -l^{(e)}m^{(e)} & -(m^{(e)})^{2}\\ -(l^{(e)})^{2} & -l^{(e)}m^{(e)} & (l^{(e)})^{2} & l^{(e)}m^{(e)} \\ -l^{(e)}m^{(e)} & -(m^{(e)})^{2} & l^{(e)}m^{(e)} & (m^{(e)})^{2} \\ \end{bmatrix}$ $\tilde\mathbf k^{(e)}=k^{(e)}\begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix}$ $\tilde \mathbf T^{(e)} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\ -m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \\ 0 & 0 & -m^{(e)} & l^{(e)} \end{bmatrix}$ $\tilde \mathbf T^{(e)^T} = \begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0 \\ m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & -m^{(e)} \\ 0 & 0 & m^{(e)} & l^{(e)} \end{bmatrix}$ $\mathbf k^{(e)} = k^{(e)}\begin{bmatrix} l^{(e)} & -m^{(e)} & 0 & 0 \\ m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & -m^{(e)} \\ 0 & 0 & m^{(e)} & l^{(e)} \end{bmatrix}\begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\ -m^{(e)} & l^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \\ 0 & 0 & -m^{(e)} & l^{(e)} \end{bmatrix}$ After the multiplication of the matrices, $\mathbf k^{(e)}$ is shown. Since it is the definition of $\mathbf k^{(e)}$, the equation is verified. $\mathbf k^{(e)} = \underset{}{k^{(e)}} \begin{bmatrix} (l^{(e)})^{2} & l^{(e)}m^{(e)} & -(l^{(e)})^{2} & -l^{(e)}m^{(e)} \\ l^{(e)}m^{(e)} & (m^{(e)})^{2} & -l^{(e)}m^{(e)} & -(m^{(e)})^{2} \\ -(l^{(e)})^{2} & -l^{(e)}m^{(e)} & (l^{(e)})^{2} & l^{(e)}m^{(e)} \\ -l^{(e)}m^{(e)} & -(m^{(e)})^{2} & l^{(e)}m^{(e)} & (m^{(e)})^{2} \\ \end{bmatrix}$

Eigenvalue and Eigenvector Analysis in the Finite Element Method

The mode shapes that correspond to the zero eigenvalues are to be plotted. These mode shapes may come out as a linear combination of the pure mode shapes (ie. 3 pure rigid body motions and 1 pure mechanism).

To solve an eigenvalue problem, the equation:

$\mathbf{Kv} = \lambda\mathbf{v}$

While {$\mathbf{u_{1}}, \mathbf{u_{2}}, \mathbf{u_{3}}, \mathbf{u_{4}}$} represent the pure eigenvectors corresponding to the four zero eigenvalues.

The linear combination of eigenvectors, {$\mathbf{u_{1}}, \mathbf{u_{2}}, \mathbf{u_{3}}, \mathbf{u_{4}}$}, can be compactly represented as follows:

$\sum_{i=1}^{4}{\alpha _{1}\mathbf{u}_{i}} =: \mathbf{W}$

Where =: means equal by definition, and the $\alpha _{1}$'s are real numbers.

The linear combination of eigenvectors - namely, W, is also an eigenvector itself, as proven in the following:

$\mathbf{K}\mathbf{W} = \mathbf{K}\sum_{i=1}^{4}{\alpha _{1}\mathbf{u}_{i}} = \mathbf{O}$

Due to all $\mathbf{K}\mathbf{u}_{i}$'s being equal to the zero vector.

Justification of Global Stiffness Matrix Assembly

For the justification of assembly of the element stiffness matrix into the global stiffness matrix, consider the two-bar truss system.

Let the global stiffness matrix be represented by K and the element stiffness matrix be represented by $k^{(e)}$, where e = 1,...,n and n is any number of element.

Recall the element force-displacement relationship $k^{(e)}_{4x4}d^{(e)}_{4x1}=f^{(e)}_{4x1}$ as well as the second method of the Euler cut principle shown below. This method involves an equilibrium analysis of node 2.

File:BikeEulerMethod2.jpg

Euler Method 2

Also recall the free body diagram of the two bar truss system in question. The free body diagrams of elements 1 and 2 are shown below, taking note of the naming process for the element degrees of freedom.

File:BikeElemA.jpg

Element 1

File:BikeElementB.jpg

Element 2

It is possible to relate the global degrees of freedom to element degrees of freedom for both element 1 and element 2. This is done below for node 2.

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

The equilibrium of node 2 can then be analyzed by cutting each element and relating the forces as shown below.

File:BikeNode2.jpg

Node 2

$\sum F_x = 0 = -f_3^{(1)} - f_1^{(2)} = 0 \;\;\; (1)$

$\sum F_y = 0 = P - f_4^{(1)} - f_2^{(2)} = 0 \;\;\; (2)$

Application of Force-Displacement Relation to Statics Analysis of 2-Bar Truss System

Equation 1 from above can be rewritten as follows:
$f_3^{(1)}+f_1^{(2)}=0\; \; \;(1a)$

Equation 2 from above can be rewritten as follows:
$f_4^{(1)}+f_2^{(2)}=P\; \; \;(2a)$

For Equation 1a, the forces can be written in terms of element stiffness(k) and degrees of freedom(d).$\mathbf{F=kd}$:
$\left[k_{31}^{(1)}d_1^{(1)} +k_{32}^{(1)}d_2^{(1)}+k_{33}^{(1)}d_3^{(1)}+k_{34}^{(1)}d_4^{(1)}\right]$
$+\left[k_{11}^{(2)}d_1^{(2)} +k_{12}^{(2)}d_2^{(2)}+k_{13}^{(2)}d_3^{(2)}+k_{14}^{(2)}d_4^{(2)}\right]\; \; \;(1b)$
$=0$

Assigned HW - Rewrite Equation 2a in terms of elemental stiffness(k) and degrees of freedom(d)$\mathbf{F=kd}$.
$\left[k_{41}^{(1)}d_1^{(1)} +k_{42}^{(1)}d_2^{(1)}+k_{43}^{(1)}d_3^{(1)}+k_{44}^{(1)}d_4^{(1)}\right]$ $+\left[k_{21}^{(2)}d_1^{(2)} +k_{22}^{(2)}d_2^{(2)}+k_{23}^{(2)}d_3^{(2)}+k_{24}^{(2)}d_4^{(2)}\right]\; \; \;(2b)$ $=P$

Local Degrees of Freedom to Global Degrees of Freedom Relation:
$d_1^{(1)}=d_1$
$d_2^{(1)}=d_2$
$d_3^{(1)}=d_2^{(2)}=d_3$
$d_4^{(1)}=d_1^{(2)}=d_4$
$d_3^{(2)}=d_5$
$d_4^{(2)}=d_6$

The local degrees of freedom(dof) of Equation 1b can be replaced with the global degrees of freedom using the above relations.

$\left[k_{31}^{(1)}d_1 +k_{32}^{(1)}d_2+k_{33}^{(1)}d_3+k_{34}^{(1)}d_4\right]$
$+\left[k_{11}^{(2)}d_3 +k_{12}^{(2)}d_4+k_{13}^{(2)}d_5+k_{14}^{(2)}d_6\right]$
$=0$

Assigned HW - Replace the local dofs with global dofs for Equation 2b.
$\left[k_{41}^{(1)}d_1 +k_{42}^{(1)}d_2+k_{43}^{(1)}d_3+k_{44}^{(1)}d_4\right]$ $+\left[k_{21}^{(2)}d_3 +k_{22}^{(2)}d_4+k_{23}^{(2)}d_5+k_{24}^{(2)}d_6\right]$ $=P$

Assembly of Global Stiffness Matrix from Element Stiffness Matrices

The global stiffness matrix $\mathbf{k}$ can be assembled from the element stiffness matrices $\mathbf{k^{(e)}}, e=1,...,n_{el}$ where n_el is the number of elements:

$\mathbf{k}_{nxn}=A\mathbf{k}_{n_{el}xn_{el}}^{(e)}$

New Nomenclature:
n: total number of global degrees of freedom before elimination of boundary conditions

n_ed: number of element degrees of freedom

(ned << n)

A: assembly operation

Principal of Virtual Work (PVW)

Used to eliminate the rows for corresponding boundary conditions to obtain $\mathbf{\bar{k}}_{2x2}$

Deriving Finite Element Method for Partial Differential Equations

The force displacement (FD) relation for a bar or spring is $kd=F$.

This implies that
$kd-f=0\; \; \;(3)$

Which is equivalent to
$w(kd-F)=0\; \; \;(4)$ for all values of w.

Proof:
$(3)\Rightarrow (4)$ trivial
$(4)\Rightarrow (3)$ not trivial

Since Equation 4 is valid for all w, select w=1, then Equation 4 becomes:
$1\left(kd-F \right)=0\Rightarrow (3)$

Analysis of MATLAB Code for Solving 5-Bar Truss System

In order to better understand the process of solving for the reactions, strains, stresses, and internal forces of a plane truss system using the Finite Element Method in MATLAB, the following MATLAB source code for solving a 5-bar truss system will be analyzed.

The 5-bar truss system to be analyzed is shown in the diagram below. This diagram also shows the global node numbers and element numbers of the system.

Diagram of 5-bar truss system to be solved using MATLAB

The values for Young's Modulus, area, length, and angle of inclination for each bar are shown in the table below.

Properties of 5-Bar Truss System

Element Number	Young's Modulus (GPa)	Cross-Sectional Area (cm2)	Length (m)	Angle of Inclination
1	200	40	3.81	$66.8^\circ$
2	200	40	3.81	$23.2^\circ$
3	200	30	5	$90^\circ$
4	200	30	5	$0^\circ$
5	70	20	2.12	$135^\circ$

MATLAB Code for 5-Bar Truss System

The analysis and dissection of the MATLAB code to solve the above 5-bar truss system is shown below. This code and the accompanying functions were downloaded from the online companion website for the textbook, Fundamental Finite Element Analysis and Applications: with Mathematica and MATLAB Computations , written by M. Asghar Bhatti.

% Five bar truss example
e1=200*10^3; e2=70*10^3; a1=40*100; a2=30*100; a3=20*100;
P = 150*10^3;
nodes = 1000*[0, 0; 1.5, 3.5; 0, 5; 5, 5];
conn = [1,2; 2,4; 1,3; 3,4; 2,3];
lmm = [1,2,3,4; 3, 4, 7, 8; 1, 2, 5, 6; 5, 6, 7, 8; 3, 4, 5, 6];  
K=zeros(8);
% Generate stiffness matrix for each element and assemble it.
for i=1:2
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(e1, a1, nodes(con,:));
    K(lm, lm) = K(lm, lm) + k;
end
for i=3:4
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(e1, a2, nodes(con,:));
    K(lm, lm) = K(lm, lm) + k;
end
 
lm=lmm(5,:); con=conn(5,:);
k=PlaneTrussElement(e2, a3, nodes(con,:));
K(lm, lm) = K(lm, lm) + k
 
% Define the load vector
R = zeros(8,1); R(4)=-P
 
% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, [1,2,7,8], zeros(4,1))
results=[];
for i=1:2
    results = [results; PlaneTrussResults(e1, a1, ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
for i=3:4
    results = [results; PlaneTrussResults(e1, a2, ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
format short g
results = [results; PlaneTrussResults(e2, a3, ...
        nodes(conn(5,:),:), d(lmm(5,:)))]

Additional Functions

The following additional MATLAB functions were necessary to properly run the above MATLAB script to solve the 5-bar truss system using the Finite Element Method.

PlaneTrussElement.m

This function takes the Young's Modulus, cross-sectional area, and element end coordinates of an element as arguments in order to first determine the length of the member. Once the length has been determined, the director cosines are determined. Using the given values for Young's Modulus, cross-sectional area, the determined values for length and director cosines, the element stiffness matrix is assembled and passed back to the main MATLAB script for solving the 5-bar truss system.

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

The following MATLAB function, named NodalSoln, is responsible for determining the displacements and reactions at each node.

This function takes the global stiffness matrix, the global right hand side vector, a specified list of degrees of freedom, and the values that correspond to those specified degrees of freedom as arguments. The overall degrees of freedom of the 5-bar truss system is set to the variable dof by equating it to the length of the global right hand side vector. Because the point of this function is to reduce the size of the global force-displacement relation by eliminating the rows and columns corresponding to zero displacement degrees of freedom, the variable df is defined as the difference between the overall degrees of freedom and the specified list of degrees of freedom that have known values. This allows the function to calculate the displacements and reactions using the reduced global force-displacement relation. These values are then passed back to the main MATLAB script to allow for further calculation.

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

The following MATLAB function, named PlaneTrussResults, is responsible for computing the results for the 5-bar truss system.

This function takes the Young's Modulus, cross-sectional area, element end coordinates, and element end displacements for a given truss element as arguments and computes the axial strain, axial stress, and axial force in the specified element. The axial strain is computed by determining how much the element changed in length during the loading process, and dividing that value by the element's original length. The axial stress is computed by multiplying the axial strain by the Young's Modulus of the element. The axial force is computed by multiplying the axial stress by the cross-sectional area of the element. These values are then stored in the variable named results and passed back to the main MATLAB script used to solve the 5-bar truss system.

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];

Analysis of MATLAB Results

A printout of the results of running the MATLAB script for solving the 5-bar truss system is shown below. It is important to note that the the first column of the variable results is the axial strain of each member, the second column is the axial stress in each member, and the third column is the axial force in each member.

K =

 1.0e+005 *

   0.3260    0.7607   -0.3260   -0.7607         0         0         0         0
   0.7607    2.9749   -0.7607   -1.7749         0   -1.2000         0         0
  -0.3260   -0.7607    2.4309    1.1914   -0.3300    0.3300   -1.7749   -0.7607
  -0.7607   -1.7749    1.1914    2.4309    0.3300   -0.3300   -0.7607   -0.3260
        0         0   -0.3300    0.3300    1.5300   -0.3300   -1.2000         0
        0   -1.2000    0.3300   -0.3300   -0.3300    1.5300         0         0
        0         0   -1.7749   -0.7607   -1.2000         0    2.9749    0.7607
        0         0   -0.7607   -0.3260         0         0    0.7607    0.3260


R =

          0
          0
          0
    -150000
          0
          0
          0
          0


d =

        0
        0
   0.5390
  -0.9531
   0.2647
  -0.2647
        0
        0


reactions =

 1.0e+005 *

   0.5493
   1.5993
  -0.5493
  -0.0993


results =

  -0.0001743      -34.859 -1.3944e+005
  -3.15e-005      -6.2999       -25200
-5.2941e-005      -10.588       -31764
-5.2941e-005      -10.588       -31764
  0.00032087       22.461        44922

The important results shown in the above printout is shown in the following table.

Results for 5-Bar Truss System

Element Number	Axial Strain	Axial Stress	Axial Force
1	-0.0001743	-34.859	-1.3944 x 105
2	-3.15 x 10-5	-6.2999	-25200
3	-5.2941 x 10-5	-10.588	-31764
4	-5.2941 x 10-5	-10.588	-31764
5	0.00032087	22.461	44922

The values for the global displacement matrix are shown below.

$\begin{Bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7}\\ d_{8}\end{Bmatrix}=\begin{Bmatrix} 0 \\ 0 \\ 0.5390 \\ -0.9531 \\ 0.2647 \\ -0.2647 \\ 0 \\ 0 \end{Bmatrix}$

Plot of Undeformed and Deformed 5-Bar Truss System

In order to plot the undeformed and deformed 5-bar truss system using MATLAB, additional code was added to the 5-bar truss system MATLAB code given above. Additionally, to ensure that the deformed structure could be distinguished from the undeformed structure, a multiplication factor was applied to the deformations at nodes 2 and 3.

The plot of the undeformed and deformed 5-bar truss system is shown below with a multiplication factor on the deformations of 500.

Plot of undeformed and deformed 5-bar truss system created using MATLAB

The source code added to the original MATLAB script in order to create 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 = 4;              % total number of nodes
   n_elem = 5;              % total number of elements
   total_dof = dof * n_node;  % total dofs of system
 
% five-bar truss system data
%
   for i=1:5
        con=conn(i,:);
        L(i)=PlaneTrussElementMod(nodes(con,:));
   end
   for i=1:5
        con=conn(i,:);
        theta(i)=PlaneTrussElementMod2(nodes(con,:))
   end
 
   position(:, 1) = [ 0; 0];
   position(:, 2) = [ L(1)*cos(theta(1)); L(1)*sin(theta(1))];
   position(:, 3) = [ 0; L(3)];
   position(:, 4) = [ L(4) ; L(3) ];
 
   % 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) = 4;
 
   node_connect(1, 3) = 1;   % element 3
   node_connect(2, 3) = 3;
 
   node_connect(1, 4) = 3;   % element 4
   node_connect(2, 4) = 4;
 
   node_connect(1, 5) = 2;   % element 5
   node_connect(2, 5) = 3;
 
   % 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 6000])
      % 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,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')
   hold on
 
 
   mult_factor = 500;
   d=d*mult_factor;
 
   position_d(:, 1) = [ 0; 0];
   position_d(:, 2) = [ L(1)*cos(theta(1))+ d(3); L(1)*sin(theta(1))+ d(4) ];
   position_d(:, 3) = [ 0+d(5); L(3)+d(6)];
   position_d(:, 4) = [ L(4) ; L(3) ];
 
   % 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) = 4;
 
   node_connect_d(1, 3) = 1;   % element 3
   node_connect_d(2, 3) = 3;
 
   node_connect_d(1, 4) = 3;   % element 4
   node_connect_d(2, 4) = 4;
 
   node_connect_d(1, 5) = 2;   % element 5
   node_connect_d(2, 5) = 3;
 
   % 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 6000])
      % 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,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')
 
 
   % put the title on the figure
   title('Deformed and Undeformed States of Five-Bar Truss System')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

Using MATLAB to Solve 3-Bar Truss System

The following 3-bar truss system will be solved and analyzed in order to better understand the process of solving for the reactions, strains, stresses, and internal forces of a plane truss system using the Finite Element Method in MATLAB. The diagram of the 3-bar truss system shown below also shows the global node numbers and element numbers of the system.

Diagram of 3-bar truss system to be solved using MATLAB

Properties of 3-Bar Truss System

Element Number	Young's Modulus	Cross-Sectional Area	Length	Angle of Inclination
1	2	3	5	$30^\circ$
2	4	1	5	$-30^\circ$
3	3	2	10	$-135^\circ$

MATLAB Code for Solving 3-Bar Truss System

% Three-bar truss example
clear all;
e = [2 4 3]; A = [3 1 2]; P = 30;
L=[5 5 10];
theta = [pi/6 -pi/6 -3*pi/4];
 
nodes = [
         -L(1)*cos(theta(1)), -L(1)*sin(theta(1));
         0, 0;
         L(2)*cos(theta(2)),-L(2)*sin(theta(2));
         -L(3)*cos(theta(3)),-L(3)*sin(theta(3))
        ];
 
dof=2*length(nodes);
 
conn=[1,2; 2,3; 2,4];
lmm = [1, 2, 3, 4; 3, 4, 5, 6; 3, 4, 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(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
 
%-----------------------------------------------------------------
% model with 2-D beam elements with 2 dofs per node
   dof = 2;         %  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
 
% three-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) = [ -L(1)*cos(theta(1)); -L(1)*sin(theta(1))];
   position(:, 2) = [ 0; 0];
   position(:, 3) = [ L(2)*cos(theta(2));-L(2)*sin(theta(2))];
   position(:, 4) = [ -L(3)*cos(theta(3));-L(3)*sin(theta(3)) ];
 
   % 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) = 3;
 
   node_connect(1, 3) = 2;   % element 3
   node_connect(2, 3) = 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([-8 6 -8 3])
      % 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,3)),y(node_connect(2,3)),'Node 4',...
     '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')
   hold on
 
 
   mult_factor = 1/500;
   d=d*mult_factor;
 
   position_d(:, 1) = [ -L(1)*cos(theta(1)); -L(1)*sin(theta(1))];
   position_d(:, 2) = [ 0+d(3); 0+d(4)];
   position_d(:, 3) = [ L(2)*cos(theta(2));-L(2)*sin(theta(2))];
   position_d(:, 4) = [ -L(3)*cos(theta(3));-L(3)*sin(theta(3)) ];
 
   % 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) = 3;
 
   node_connect_d(1, 3) = 2;   % element 3
   node_connect_d(2, 3) = 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([-8 6 -8 3])
      % 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,3)),y(node_connect_d(2,3)),'Node 4',...
     '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')   
 
   % put the title on the figure
   title('Deformed and Undeformed States of Three-Bar Truss System')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

Additional Functions

The following additional MATLAB functions were necessary to properly run the above MATLAB script to solve the 3-bar truss system using the Finite Element Method.

PlaneTrussElement.m

This function takes the Young's Modulus, cross-sectional area, and element end coordinates of an element as arguments in order to first determine the length of the member. Once the length has been determined, the director cosines are determined. Using the given values for Young's Modulus, cross-sectional area, the determined values for length and director cosines, the element stiffness matrix is assembled and passed back to the main MATLAB script for solving the 5-bar truss system.

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

The following MATLAB function, named NodalSoln, is responsible for determining the displacements and reactions at each node.

This function takes the global stiffness matrix, the global right hand side vector, a specified list of degrees of freedom, and the values that correspond to those specified degrees of freedom as arguments. The overall degrees of freedom of the 5-bar truss system is set to the variable dof by equating it to the length of the global right hand side vector. Because the point of this function is to reduce the size of the global force-displacement relation by eliminating the rows and columns corresponding to zero displacement degrees of freedom, the variable df is defined as the difference between the overall degrees of freedom and the specified list of degrees of freedom that have known values. This allows the function to calculate the displacements and reactions using the reduced global force-displacement relation. These values are then passed back to the main MATLAB script to allow for further calculation.

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

The following MATLAB function, named PlaneTrussResults, is responsible for computing the results for the 5-bar truss system.

This function takes the Young's Modulus, cross-sectional area, element end coordinates, and element end displacements for a given truss element as arguments and computes the axial strain, axial stress, and axial force in the specified element. The axial strain is computed by determining how much the element changed in length during the loading process, and dividing that value by the element's original length. The axial stress is computed by multiplying the axial strain by the Young's Modulus of the element. The axial force is computed by multiplying the axial stress by the cross-sectional area of the element. These values are then stored in the variable named results and passed back to the main MATLAB script used to solve the 5-bar truss system.

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];

Results of MATLAB Analysis on 3-Bar Truss System

Running the MATLAB code for solving a 3-bar truss system produced the following results.

k =

         0.9      0.51962         -0.9     -0.51962
     0.51962          0.3     -0.51962         -0.3
        -0.9     -0.51962          0.9      0.51962
    -0.51962         -0.3      0.51962          0.3


k =

         0.6      0.34641         -0.6     -0.34641
     0.34641          0.2     -0.34641         -0.2
        -0.6     -0.34641          0.6      0.34641
    -0.34641         -0.2      0.34641          0.2


k =

         0.3          0.3         -0.3         -0.3
         0.3          0.3         -0.3         -0.3
        -0.3         -0.3          0.3          0.3
        -0.3         -0.3          0.3          0.3


K =

 Columns 1 through 7

         0.9      0.51962         -0.9     -0.51962            0            0            0
     0.51962          0.3     -0.51962         -0.3            0            0            0
        -0.9     -0.51962          1.8        1.166         -0.6     -0.34641         -0.3
    -0.51962         -0.3        1.166          0.8     -0.34641         -0.2         -0.3
           0            0         -0.6     -0.34641          0.6      0.34641            0
           0            0     -0.34641         -0.2      0.34641          0.2            0
           0            0         -0.3         -0.3            0            0          0.3
           0            0         -0.3         -0.3            0            0          0.3

 Column 8

           0
           0
        -0.3
        -0.3
           0
           0
         0.3
         0.3


R =

    0
    0
    0
   30
    0
    0
    0
    0


d =

           0
           0
     -435.17
      671.77
           0
           0
           0
           0


reactions =

      42.588
      24.588
      28.392
      16.392
     -70.981
     -70.981


results =

     -8.1962      -16.392      -32.785      -24.588      -49.177      -32.785      -49.177
      8.1962       16.392       32.785       24.588       49.177       32.785       49.177
      -16.73      -33.461      -66.921      -50.191      -100.38      -66.921      -100.38

The values for the global displacement matrix are shown below.

$\begin{Bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7}\\ d_{8}\end{Bmatrix}=\begin{Bmatrix} 0 \\ 0 \\ -435.17 \\ 671.77 \\ 0 \\ 0 \\ 0 \\ 0 \end{Bmatrix}$

Plot of Undeformed and Deformed 3-Bar Truss System

The plot of the undeformed and deformed 3-bar truss system is shown below with a multiplication factor on the deformations of 1/500.

Plot of undeformed and deformed 3-bar truss system

Analysis of Unstable 3-Bar Truss System

Finite Element Method Analysis

The setup for the unstable 3-bar truss system presented by Professor Loc Vu Quoc is shown below.

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

Properties of Unstable 3-Bar Truss System

Element Number	Young's Modulus	Cross-Sectional Area	Length	Angle of Inclination
1	2	3	1	$90^\circ$
2	2	3	1	$0^\circ$
3	2	3	1	$-90^\circ$

In order to analyze this 3-bar truss system, the following MATLAB code was developed by modifying the 2-bar truss system MATLAB code, shown in Homework Report 3.

% Three-bar truss example
clear all;
e = [2 2 2]; A = [3 3 3]; P = 1;
L=[1 1 1];
theta = [pi/2 0 -pi/2];
 
nodes = [
         0 , 0;
         0 , L(1);
         L(2), L(1);
         L(2) , 0
        ];
 
dof=2*length(nodes);
 
conn=[1,2; 2,3; 3,4];
lmm = [1, 2, 3, 4; 3, 4, 5, 6; 5, 6, 7, 8];
elems=size(lmm,1);
K=zeros(dof); R = zeros(dof,1);
debc = [1, 2, 7, 8];
ebcVals = zeros(length(debc),1);
 
%load vector
R = zeros(dof,1); R(4) = P; R(5) = 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

The results obtained from running the code given above is shown below.

k =

    0     0     0     0
    0     6     0    -6
    0     0     0     0
    0    -6     0     6


k =

    6     0    -6     0
    0     0     0     0
   -6     0     6     0
    0     0     0     0


k =

    0     0     0     0
    0     6     0    -6
    0     0     0     0
    0    -6     0     6


K =

    0     0     0     0     0     0     0     0
    0     6     0    -6     0     0     0     0
    0     0     6     0    -6     0     0     0
    0    -6     0     6     0     0     0     0
    0     0    -6     0     6     0     0     0
    0     0     0     0     0     6     0    -6
    0     0     0     0     0     0     0     0
    0     0     0     0     0    -6     0     6


R =
 
    0
    0
    0
    1
    1
    0
    0
    0

Warning: Matrix is singular to working precision.
> In NodalSoln at 12
  In ThreeBarTrussMod at 38

d =

    0
    0
  Inf
  NaN
  NaN
  NaN
    0
    0


reactions =

  NaN
  NaN
  NaN
  NaN


results =

  NaN   NaN   NaN   NaN   NaN   NaN   NaN
  NaN   NaN   NaN   NaN   NaN   NaN   NaN
  NaN   NaN   NaN   NaN   NaN   NaN   NaN

The results above show that the Finite Element Method was not able to determine the solution to the 3-bar truss system introduced in this section. This is due to the fact that the global stiffness matrix for this 3-bar truss system is singular, making the force-displacement relation impossible to solve. The significance of the global stiffness being singular is that it shows that the system is unstable. This means that any force applied to this structure will result in an infinite deformation, which is why it cannot be solved for using the Finite Element Method.

Analysis of Eigenvalues and Eigenvectors

The eigenvalues (D) and eigenvectors (V) of the K matrix were found and are shown below.

As seen above, the first five columns are the zero eigenvectors. Thus, the first 5 columns of the eigenvector matrix correspond to the zero eigenvectors.

The following code was used to plot these mode shapes.

EDU>> zeroe3bar = V(:,1) + V(:,2) + V(:,3) + V(:,4);
EDU>> plot(zeroe3bar,':')
EDU>> hold
Current plot held
EDU>> stem(zeroe3bar)
EDU>> xlabel('Elements');ylabel('Mode Amplitude')

Which produced a graph as follows:

File:Feamodegraph.jpg

Adding an Additional Member to the Unstable 3-Bar Truss System

Finite Element Method Analysis

The setup for a truss system that will fix the problems encountered with the unstable 3-bar truss system is shown below. This setup has an additional member added as a cross-member to increase structural stability under loading.

Setup of unstable 3-bar truss system with additional member added for stability

Properties of Modified Unstable 3-Bar Truss System

Element Number	Young's Modulus	Cross-Sectional Area	Length	Angle of Inclination
1	2	3	1	$90^\circ$
2	2	3	1	$0^\circ$
3	2	3	1	$-90^\circ$
4	2	3	$\sqrt{2}$	$45^\circ$

The MATLAB code used to solve and plot the modified unstable three-bar truss system is shown below.

% Three bar truss example
E = [2 2 2 2];
A = [3 3 3 3];
L = [1 1 1 sqrt(2)]
P = 0.5;
nodes = [0, 0; 0, L(1); L(2), L(1); L(2), 0];
conn = [1,2; 2,3; 3,4; 1,3];
lmm = [1, 2, 3, 4; 3, 4, 5, 6; 5, 6, 7, 8; 1, 2, 5, 6];  
K=zeros(8);
% Generate stiffness matrix for each element and assemble it.
for i=1:4
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(E(i), A(i), nodes(con,:))
    K(lm, lm) = K(lm, lm) + k;
end
K
% Define the load vector
R = zeros(8,1); R(5)= P
 
% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, [1,2,7,8], zeros(4,1))
results=[];
for i=1:4
    results = [results; PlaneTrussResults(E(i), A(i), ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
 
format short g
 
 
%-----------------------------------------------------------------
% model with 2-D beam elements with 2 dofs per node
   dof = 2;         %  dof per node
 
%
% input the nodal coordinates
%
   n_node = 4;              % total number of nodes
   n_elem = 4;              % total number of elements
   total_dof = dof * n_node;  % total dofs of system
 
% three-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) = [ 0; L(1)];
   position(:, 3) = [ L(2); L(1)];
   position(:, 4) = [ L(2) ; 0 ];
 
   % 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) = 3;
 
   node_connect(1, 3) = 3;   % element 3
   node_connect(2, 3) = 4;
 
   node_connect(1, 4) = 1;   % element 4
   node_connect(2, 4) = 3;
 
   % 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([-.5 1.5 -.5 1.5])
      % 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,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')
 
   mult_factor = 1;
   d=d*mult_factor;
 
   position_d(:, 1) = [ 0; 0];
   position_d(:, 2) = [ 0+d(3); L(1)+d(4)];
   position_d(:, 3) = [ L(2)+d(5); L(1)+d(6)];
   position_d(:, 4) = [ L(2) ; 0 ];
 
   % 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) = 3;
 
   node_connect_d(1, 3) = 3;   % element 3
   node_connect_d(2, 3) = 4;
 
   node_connect_d(1, 4) = 1;   % element 34
   node_connect_d(2, 4) = 3;
 
   % 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([-.5 1.5 -.5 1.5])
      % 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,3)),y(node_connect_d(2,3)),'Node 4',...
     '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')
 
   % put the title on the figure
   title('Deformed and Undeformed States of Modified Unstable Three-Bar Truss System')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

The results obtained from running the above MATLAB code developed to solve the modified unstable three-bar truss system using the Finite Element Method are shown below.

k =

    0     0     0     0
    0     6     0    -6
    0     0     0     0
    0    -6     0     6


k =

    6     0    -6     0
    0     0     0     0
   -6     0     6     0
    0     0     0     0


k =

    0     0     0     0
    0     6     0    -6
    0     0     0     0
    0    -6     0     6


k =

      2.1213       2.1213      -2.1213      -2.1213
      2.1213       2.1213      -2.1213      -2.1213
     -2.1213      -2.1213       2.1213       2.1213
     -2.1213      -2.1213       2.1213       2.1213


K =

 Columns 1 through 7

      2.1213       2.1213            0            0      -2.1213      -2.1213            0
      2.1213       8.1213            0           -6      -2.1213      -2.1213            0
           0            0            6            0           -6            0            0
           0           -6            0            6            0            0            0
     -2.1213      -2.1213           -6            0       8.1213       2.1213            0
     -2.1213      -2.1213            0            0       2.1213       8.1213            0
           0            0            0            0            0            0            0
           0            0            0            0            0           -6            0

 Column 8

           0
           0
           0
           0
           0
          -6
           0
           6


R =

           0
           0
           0
           0
         0.5
           0
           0
           0


d =

           0
           0
     0.31904
           0
     0.31904
   -0.083333
           0
           0


reactions =

        -0.5
        -0.5
           0
         0.5

The plot developed using the above code showing the deformed and undeformed structure of the modified unstable 3-bar truss system is shown below.

Plot of deformed and undeformed shapes of unstable 3-bar truss system with additional member added for stability

Adding the additional member to the unstable 3-bar truss system caused the system to become stable. This is shown by the fact that the global stiffness matrix is not singular for this problem, allowing the global force-displacement relation to be solved.

Eigenvalue and Eigenvector Analysis

The same analysis as was conducted for the three bar system is now performed on the four bar system and the eigenvalues and eigenvectors are found to be as follows:

Columns 2 through 5 correspond to the zero eigenvalues and are the eigenvectors which correspond to the mode shapes. They are plotted using the same code as before and the graph can be seen below:

File:Fea4modegraph.jpg

Contributing Team Members

Andrew McDonald - Eml4500.f08.bike.mcdonald 14:53, 17 October 2008 (UTC)
Garrett Pataky - Eml4500.f08.bike.pataky 17:42, 24 October 2008 (UTC)
Sam Bernal - Eml4500.f08.bike.bernal 19:51, 24 October 2008 (UTC)
Bobby Sweeting - Eml4500.f08.bike.sweeting 19:35, 13 October 2008 (UTC)
Shawn Gravois - Eml4500.f08.bike.gravois 22:25, 17 October 2008 (UTC)
Eric Viale - Eml4500.f08.bike.viale 19:47, 24 October 2008 (UTC)