# User:Eml4500.f08.ateam/hw 4

 The A-Team

A-Team Teampage

Homework Report #4

This page is intended for Homework Report #4 of the A-Team.

Due date: Friday, 24 October 08, 5 pm EDT = 17:00 EDT (or 21:00 UTC).

## Connectivity and LMM Arrays

Global/Local DOF relationship

From the diagram of the two bar truss above, we can see the relationships between the global and local nodes of elements 1 and 2. From these relationships, we can derive the connectivity array. In MATLAB, this is denoted as conn(e,j). In this array, the rows represent the elements (e) and the columns represent the local node numbers (j). The values inside the array are the global node numbers. Using the two bar truss diagram above, the connectivity array becomes:

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

For example, in this case the local node of element 1 corresponds to the global node 1. Therefore, the value of in the first row and column is 1.

Another array that we will need to derive is the location matrix master array or "LMM". This array, will represent the relationships between the local degrees of freedom (DOF's) and the global DOF's numbers. As in the connectivity array, the rows in the LMM represent the elements (e), but the columns represent the local DOF's (j). In MATLAB this command is lmm(e,j).

Element 1
Element 2

The above diagrams show the elements of the two bar truss and their local DOF's. From these diagrams, we will populate our location matrix master array. The values in the array correspond to the global DOF's. The array then becomes:

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

## Transforming Local DOFs to Another Local Set

The goal of this next step is to find $\tilde{T}_{4x4}^{(e)}$ that transforms the set of local, or elemental, degrees of freedom, ${d_{4x1}^{(e)}}$, to another set of local degrees of freedom, $\tilde{d}_{4x1}^{(e)}$, so that the matrix $\tilde{T}_{4x4}^{(e)}$ is invertible.

Transfer of local dofs

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

$\tilde{d}_1^{(e)} = [ \begin{matrix} l^{(e)} & m^{(e)} \end{matrix} ] \left[ \begin{matrix} d_1^{(e)} \\ d_2^{(e)} \end{matrix} \right]$

The right side of the equation has been seen before in previous lectures when trying to find the axial displacement, or the orthogonal projection, of the node. That means we can also write the equation as:

${q}_1^{(e)} = [ \begin{matrix} l^{(e)} & m^{(e)} \end{matrix} ] \left[ \begin{matrix} d_1^{(e)} \\ d_2^{(e)} \end{matrix} \right]$

### Derivations

Derivation of ${d_1^{(e)}}$:

$\tilde {d}_1^{(e)} = {\vec d_1^{(e)}} \cdot \vec \tilde{i}$

${d_1^{(e)}} = d_1^{(e)} \vec i + d_2^{(e)} \vec j$.

$\tilde {d_1^{(e)}} = (d_1^{(e)} \tilde{i} + d_2^{(e)} \tilde{j}) \cdot \vec \tilde{i}$

$\tilde {d_1^{(e)}} = d_1^{(e)} (\vec i \cdot \vec \tilde{i}) + d_2^{(e)} (\vec j \cdot \vec \tilde{i})$

Where:

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

$m^{(e)} = {\tilde{i}} \cdot {j} = sin(\theta^{(e)})$

Therefore:

$\tilde{d}_1^{(e)} = [ \begin{matrix} l^{(e)} & m^{(e)} \end{matrix} ] \left[ \begin{matrix} d_1^{(e)} \\ d_2^{(e)} \end{matrix} \right]$ (1)

This process is similar when finding $\tilde{d}_2^{(e)}$. The only difference is that the components of $\overrightarrow {d}_1^{(e)}$ are found along the $\tilde {j}$ or the $\tilde {y}$ axis.

Derivation of ${d_2^{(e)}}$:

${d_1^{(e)}} = d_1^{(e)} \vec i + d_2^{(e)} \vec j$.

$\tilde {d}_2^{(e)} = \vec d_1^{(e)} \cdot \vec \tilde{j}$

$\tilde {d}_2^{(e)} = (d_1^{(e)} \tilde{i} + d_2^{(e)} \tilde{j}) \cdot \vec \tilde{j}$

$\tilde {d}_2^{(e)} = d_1^{(e)} (\vec i \cdot \vec \tilde{j}) + d_2^{(e)} (\vec j \cdot \vec \tilde{j})$

$(\vec i \cdot \vec \tilde{j}) = -sin(\theta^{(e)})$ and

$(\vec j \cdot \vec \tilde{j}) = cos(\theta^{(e)})$

$\tilde {d}_2^{(e)} = -m^{(e)} d_1^{(e)} + l^{(e)} d_2^{(e)}$

$\tilde {d}_2^{(e)} = [ \begin{matrix} -m^{(e)} & l^{(e)} \end{matrix} ] \left[ \begin{matrix} d_1^{(e)} \\ d_2^{(e)} \end{matrix} \right]$ (2)

Equations (1) and (2) can now be put into a matrix form

$\begin{Bmatrix} \tilde{d}_1^{(e)} \\ \tilde{d}_2^{(e)} \end{Bmatrix} = \left[ \begin{matrix} l^{(e)} & m^{(e)} \\ -m^{(e)}& l^{(e)} \end{matrix} \right] \begin{Bmatrix} d_1^{(e)} \\ d_2^{(e)} \end{Bmatrix}$

where $\left[ \begin{matrix} l^{(e)} & m^{(e)} \\ -m^{(e)}& l^{(e)} \end{matrix} \right]$ can be written as ${R_{2x2}^{(e)}}$

Now the equation $\tilde{d}_{4x1}^{(e)} = \tilde{T}_{4x4}^{(e)} {d_{4x1}^{(e)}}$ can be formed from these equations.

$\begin{Bmatrix} \tilde{d}_1^{(e)} \\ \tilde{d}_2^{(e)}\\ \tilde{d}_3^{(e)} \\\tilde{d}_4^{(e)} \end{Bmatrix} = \left[ \begin{matrix} R_{2x2}^{(e)} & 0_{2x2} \\ 0_{2x2}& R_{2x2}^{(e)} \end{matrix} \right] \begin{Bmatrix} d_1^{(e)} \\ d_2^{(e)} \\ d_3^{(e)} \\ d_4^{(e)} \end{Bmatrix}$

### Transformed Forces and Degrees of Freedom

File:19.2.jpg
Rotated element

The element above has been rotated for an easier understanding.

$\tilde{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{d}^{(e)}$

If you try to move this node transversely, nothing will happen. That is why columns 2 and 4 are all zeros, because they do not stretch the string.

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

This gives us the final equation of:

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

## Changing global degrees of freedom to a local coordinate system aligned along the truss

The tilde vectors represent the axial and normal forces from the axis of the bar. The same equation is used to start the derivation.

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

A transformation vector will be used to get $\tilde\mathbf {d}$ as well as $\tilde {f}$.

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

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

By substituting the values for $\tilde \mathbf{T}$ it can be proven.

$\tilde \mathbf{T} = \begin{bmatrix} cos (\Theta) & sin (\Theta) & 0 & 0\\ -sin (\Theta) & cos (\Theta) & 0 & 0\\0 & 0 & cos (\Theta) & sin (\Theta)\\0 & 0 & -sin (\Theta) & cos (\Theta)\end{bmatrix}=\begin{bmatrix} \mathbf{R}^{(e)} & \mathbf{0}\\\mathbf{0} & \mathbf{R}^{(e)}\end{bmatrix}$

$\tilde \mathbf{T}^{(1)} \mathbf{f}^{(1)}= \begin{bmatrix} cos (30) & sin (30) & 0 & 0\\ -sin (30) & cos (30) & 0 & 0\\0 & 0 & cos (30) & sin (30)\\0 & 0 & -sin (30) & cos (30)\end{bmatrix}\begin{bmatrix}-4.378\\-2.5622\\4.4378\\2.5622\end{bmatrix} = \tilde\mathbf{f}^{(1)}$

$\tilde \mathbf{f}^{(1)} = \begin{bmatrix} -5.1243\\0\\5.1243\\0\end{bmatrix}$

This confirms the previous statics work for $\tilde \mathbf{f}^{(1)}$.

By substituting the previous relationships into $\tilde \mathbf{f}^{(e)}= \tilde \mathbf{k}^{(e)} \tilde \mathbf{d}^{(e)}$, this equation can be rewritten as follows.

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

## Verify transformation matrix transpose theorem

If the transformation matrix, $\tilde \mathbf{T}^{(e)}$, is invertible then the equation can be changed into the following form.

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

Now since only a square matrix can be inverted, the following process was investigated. Since $\tilde \mathbf{T}^{(e)}$ is composed of R matrices which are sin and cos functions, when the transpose of $\tilde \mathbf{R}^{(e)}$ is multiplied by $\tilde \mathbf{R}^{(e)}$, the result is the identity matrix.

$\mathbf {R}^{(e)T} \ \mathbf {R}^{(e)} = \begin{bmatrix} \sin^2 {\Theta} + \cos^2 {\Theta} & 0 \\ 0 & \sin^2 {\Theta} + \cos^2 {\Theta} \end{bmatrix}=\begin{bmatrix} 1&0\\0&1\end{bmatrix}$

Therefore if both sides are divided by $\mathbf {R}^{(e)}$, then the following relationship can be made.

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

This simplifies the math involved significantly and allows for rectangular matrices to be used instead of solely square matrices. Since $\tilde \mathbf{T}^{(e)}$ is comprised of a diagonal $\tilde \mathbf{R}^{(e)}$ matrix, the previous relationship between the transpose and inverse can be applied to $\tilde \mathbf{T}^{(e)}$ as well.

To verify that $[\tilde \mathbf{T}^{(e)-1} \ \tilde \mathbf{k}^{(e)} \ \tilde\mathbf {T}^{(e)}]$, the values from the 2-bar truss will be inserted. The first element's values will be used.

$[\tilde \mathbf{T}^{(e)-1} \ \tilde \mathbf{k}^{(e)} \ \tilde\mathbf {T}^{(e)}] = \begin{bmatrix} cos (30) & -sin (30) & 0 & 0\\ sin (30) & cos (30) & 0 & 0\\0 & 0 & cos (30) & -sin (30)\\0 & 0 & sin (30) & cos (30)\end{bmatrix}\begin{bmatrix}0\\0\\4.4352\\6.1271\end{bmatrix} \begin{bmatrix} cos (30) & sin (30) & 0 & 0\\ -sin (30) & cos (30) & 0 & 0\\0 & 0 & cos (30) & sin (30)\\0 & 0 & -sin (30) & cos (30)\end{bmatrix}$

$= \begin{bmatrix} 0.5625 & 0.325&-0.5625&-0.325\\0.325&0.1875&-0.325&-0.1875\\-0.5625&-0.325&0.5625&0.325\\-0.325&-0.1875&0.325&0.1875\end{bmatrix}= \mathbf{k}^{(1)}$

This verifies the derivation made earlier.

## Eigenvalues of the Matrix

When the eigenvalues of the 6x6 matrix $\underline{K}$ are calculated, there are 4 zero values. Three of these values correspond to the three rigid body motions and the other corresponds to the mechanism.

Eigenvalue problem: For this problem we will use the following equation: $\underline{K}\underline{v}=\lambda\underline{v}$

For this matrix, we will let {$\underline{u}_{1}, \underline{u}_{2}, \underline{u}_{3}, \underline{u}_{4}$} be the pure eigenvectors corresponding to the four zero eigenvalues. The equation above now becomes:

$\underline{K}\underline{u}_{i}=0\cdot\underline{u}_{i}=\underline{0}$

Here $\underline{K}$ is a 6x6 matrix, $\underline{u}$ is a 6x1 matrix, and i=1,...,4.

We can now compute a linear combination on the values of $\underline{u}_{i}, i=1,...,4$ with a summation:

$\sum_{i=1}^4 \alpha_{i} \underline{u}_{i} =: \underline{W}$.

In this summation, $\alpha_{i}$ is a 1x1 matrix, $\underline{u}_{i}$ is a 6x1 matrix, $\underline{W}$ is a 6x1, and $\alpha_{i}$ are real numbers. The $=$ sign here means "equal by definition" and $\underline{W}$ is a zero eigenvector corresponding to a zero eigenvalue. Our original equation for the eigenvalue problem now becomes:

$\underline{K}\underline{W}=\underline{K}(\sum_{i=1}^4\alpha_{i} \underline{u}_{i})=\sum_{i=1}^4 \alpha_{i}(\underline{K}\underline{u}_{i})=\underline{0}=0\cdot\underline{W}$

Using the equation $\underline{K}\underline{v}=\lambda\underline{v}$ we will now evaluate a three and four bar truss.

Three and Four bar trusses

Given the values in the diagram, we can compute the stiffness matrix $\underline{K}$.

HW: Plot the eigenvectors corresponding to the zero eigenvalues of the two-bar truss system.

For the two bar truss system, our matrix $\underline{K}$ is defined as:

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


With this matrix, we can now plot the eigenvalues and eigenvectors using the Matlab code [V,D]=eig(K). This will display the eigenvectors (V) and eigenvalues (D):

[V,D]=eig(K)

V =

-0.1118    0.5043   -0.0000   -0.5931    0.6174   -0.0139
-0.0814   -0.8634    0.0000   -0.3476    0.3565   -0.0080
-0.4628    0.0089   -0.0000   -0.4803   -0.5409    0.5123
0.5266   -0.0053    0.0000   -0.5429   -0.4330   -0.4904
-0.4947    0.0071   -0.7071    0.0313   -0.0765   -0.4984
0.4947   -0.0071   -0.7071   -0.0313    0.0765    0.4984

D =

-0.0000         0         0         0         0         0
0   -0.0000         0         0         0         0
0         0    0.0000         0         0         0
0         0         0    0.0000         0         0
0         0         0         0    1.4706         0
0         0         0         0         0   10.0294


HW2: Plot the eigenvectors corresponding to the zero eigenvalues of the three-bar truss system above. The matrix $\underline{K}$ has been determined to be:

A =

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    -1     0    12     0     0     0    -6
0     0    -6     0     6     0     0     0
0     0     0     0     0     0     0     0
0     0     0     0     0     0     0     0
0     0     0    -6     0     0     0     6


Using the same command from earlier ([V,D]=eig(A)) we can determine the eigenvalues and eigenvectors of this matrix:

[V,D]=eig(A)

V =

0         0         0         0         0    1.0000         0         0
-0.4538   -0.6355    0.9864         0         0         0         0         0
0         0         0    0.7071    0.7071         0         0         0
0.7670   -0.4386    0.0000         0         0         0         0         0
0         0         0   -0.7071    0.7071         0         0         0
0         0         0         0         0         0    1.0000         0
0         0         0         0         0         0         0    1.0000
-0.4538   -0.6355   -0.1644         0         0         0         0         0

D =

16.1414         0         0         0         0         0         0         0
0    1.8586         0         0         0         0         0         0
0         0    6.0000         0         0         0         0         0
0         0         0   12.0000         0         0         0         0
0         0         0         0   -0.0000         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


## Justification of Assembly of Elemental Stiffness Matrices into the Global Stiffness Matrix

k(e) = element stiffness matrix k = global stiffness matrix

Consider the example of the 2 bar truss: Recall the element force displacement (FD) relationship:

$k^{(e)}*d=f^{(e)}$

Use the Euler cut principle (Method 2) to solve for the forces. The Euler cut principle states that node 2 should be in equilibrium. We take this and solve for the forces in each of the 2 bars.

Free Body Diagrams of elements 1 and 2.

FBDs of Elements 1 and 2

The element degrees of freedom are given in the matrix: d(e). This is a 4x1 matrix.

For node 2, identify the global degrees of freedom to element degrees of freedom for both of the bars.

FBD of Node 2

Now we can use the equilibrium concept from the Euler cut method to solve for the forces.

    Sum of forces in x direction: F = 0 = -f(1)3 - f(2)1
Sum of forces in y direction: F = 0 = P -f (1)4 - f(2)2


Next we use the element force displacement relationship to obtain the displacement of node 2:

$k^{(e)}*d=f^{(e)}$

The free body diagram at Node 2 yields the following force balance equations.

$f_3^{(1)}+f_1^{(2)}=0\rightarrow (1)$

$f_4^{(1)}+f_2^{(2)}=\mathbf{P}\rightarrow (2)$

Next, the elemental Force-Displacement relation,

$k^{(e)}d^{(e)}=f^{(e)}$

can be used to rewrite Equation (1) as follows.

$(1)\rightarrow [k_{31}^{(1)}d_1^{(1)}+k_{32}^{(1)}d_2^{(1)}+k_{33}^{(1)}d_3^{(1)}+k_{34}^{(1)}d_4^{(1)}]+ [k_{11}^{(2)}d_1^{(2)}+k_{12}^{(2)}d_2^{(2)}+k_{13}^{(1)}d_3^{(2)}+k_{14}^{(1)}d_4^{(2)}]$

Similarly, Equation (2) can be written as follows.

$(2)\rightarrow [k_{41}^{(1)}d_1^{(1)}+k_{42}^{(1)}d_2^{(1)}+k_{43}^{(1)}d_3^{(1)}+k_{44}^{(1)}d_4^{(1)}]+ [k_{21}^{(2)}d_1^{(2)}+k_{22}^{(2)}d_2^{(2)}+k_{23}^{(1)}d_3^{(2)}+k_{24}^{(1)}d_4^{(2)}]$

Now, Equations (1) and (2) can be transformed into the global coordinate system as follows.

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

$(2)\rightarrow [k_{41}^{(1)}d_1+k_{42}^{(1)}d_2+k_{43}^{(1)}d_3+k_{44}^{(1)}d_4]+ [k_{21}^{(2)}d_3+k_{22}^{(2)}d_4+k_{23}^{(1)}d_5+k_{24}^{(1)}d_6]$

Now, the elemental stiffness matrices, k(e), can be assembled into the global stiffness matrix, K, through

$\mathbf{K}=\mathbf{A}\mathbf{k}$

where K is an n x n matrix, A is nel x 1 matrix and k is an ned x ned matrix. Here, n is the total number of degrees of freedom, nel is the total number of elements, ned is the total number of elemental degrees of freedom, and A is the assembly operator.

As we continue to justify the assembly of the elemental stiffness matrices into the global stiffness matrix, we must use the Principal of Virtual Work, as described previously, to obtain $\mathbf{\overline{K}}$ by eliminating the corresponding boundary conditions.

Once the Principal of Virtual Work is applied, the resulting elemental axial force equation and the elemental elemental stiffness equation can be written as follows.

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

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

Now, it is important to derive the Finite Element Method for Partial Differential Equations. Consider the Force-Displacement relation for a bar or spring with

$kd=F$

This implies that

$kd-F=0\rightarrow (3)$

or equivalently

$W(kd-F)=0\rightarrow (4)$

for all W.

Upon investigation, $(3)\Rightarrow(4)$ is trivial, but $(4)\Rightarrow(3)$ is not trivial.

Since (4) is valid for all W, select W=1. Now, (4) becomes $1(kd-F)=0\Rightarrow(3)$.

## Five-Bar Truss System: MATLAB Dissection

Analysis of the following five-bar truss system MATLAB code is important before attempting to write code for the three-bar truss system. The following code is from the FiveBarTrussEx.m file found at the website for the class textbook, Fundamental Finite Element Analysis and Applications.

Five Bar Truss System
Five-Bar Truss System Data
Element Cross-section Area (cm2) Young's Modulus (GPa) Length (m) Angle (degrees)
1

40

200

3.81

66.8

2

40

200

3.81

23.2

3

30

200

5.00

90

4

30

200

5.00

0

5

20

70

2.12

135

% 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,:)))]


The following functions are essential for the MATLAB program above to run. They are also found at the course textbook website. These functions include PlaneTrussElement.m, NodalSoln.m, and PlaneTrussResults.m. The details of these functions are outlined below.

#### PlaneTrussElement.m

This function generates the element stiffness matrix for each element. When the function is called it calculates the X and Y coordinate positions of the element ends from the coordinates passed in. It produces the element length from this information. Next it calculates the direction cosines of the element. Using the Young's modulus and cross-sectional area of the element, the element stiffness matrix, k, is assembled. The stiffness matrix is passed back to the FiveBarTrussEx.m.

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

This function solves the force-displacement relationship for each node, resulting in the displacements and reactions at each node. This information is passed back to the FiveBarTrussEx.m.

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

This function outputs the axial strain, axial stress, and axial force of each element.

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


### MATLAB Results

The following is a printout of the exact results of the MATLAB analysis. It includes the complete Global Stiffness matrix (K), the Applied Force matrix (R), the Displacement matrix (d), the reaction forces, and the resulting stresses and strains.

K =
Columns 1 through 6
32600        76067       -32600       -76067            0            0
76067  2.9749e+005       -76067 -1.7749e+005            0    -1.2e+005
-32600       -76067  2.4309e+005  1.1914e+005       -32998        32998
-76067 -1.7749e+005  1.1914e+005  2.4309e+005        32998       -32998
0            0       -32998        32998    1.53e+005       -32998
0    -1.2e+005        32998       -32998       -32998    1.53e+005
0            0 -1.7749e+005       -76067    -1.2e+005            0
0            0       -76067       -32600            0            0

Columns 7 through 8
0            0
0            0
-1.7749e+005       -76067
-76067       -32600
-1.2e+005            0
0            0
2.9749e+005        76067
76067        32600

R =
0
0
0
-150000
0
0
0
0

d =
0
0
0.53895
-0.95306
0.2647
-0.2647
0
0

reactions =
54927
1.5993e+005
-54927
-9926.7

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 d matrix represents the displacements at each node.

$\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.53895 \\ -0.95306 \\ 0.2647 \\ -0.2647 \\ 0 \\ 0 \\ \end{Bmatrix}$

Reaction forces at Global Node 1 and Global Node 4:

Five-Bar Truss System Reaction Forces
Reaction Force
F(1)1

54927

F(1)2

1.5993x105

F(4)1

-54927

F(4)1

-9926.7

Resulting Stresses and Strains:

Five-Bar Truss System Results
Element Axial Strain Axial Stress Axial Force
1

-1.74x10-4

-34.859

-1.39x105

2

-1.39x105

-6.2999

-25200

3

-5.29x105

-10.588

5.00

4

-5.294x10-5

-10.588

-31764

5

3.21x10-4

22.461

44922

### Plot of Undeformed and Deformed Five-Bar Truss System

% Matlab Script
%*********************************************************************
% filename: fivebar.m
%
% PURPOSE:
%    Example of Five-Bar Truss deformation due to applied forces
%
% AUTHOR: Daniel DiDomenico
%
% Modified on: Sun, 19 Oct 2008, 16:08:58 EDT
% Created on : Sun, 19 Oct 2008, 09:40:47 EST
%
% DEPENDENCIES:
%   call:
%
% REMARKS:  This code is added to the FiveBarTrussEx.m file to plot
%           the system deformation.
%
%*********************************************************************

% 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

% two-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

%Nodal position
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

d = d*500

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

% 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')

Five-Bar Truss System Deformation

## Three-Bar Truss System: MATLAB Dissection

The following code solves the three-bar truss system presented in class and shown below. The two-bar truss system code provided by Dr. Vu-Quoc on the EML4500 website was modified to account for the additional element. The original code can be found here. After analysis, the resulting deformation was scaled by a factor of 1/500 for presentation.

Three-Bar Truss System
% Matlab Script
%*********************************************************************
% filename: threebar.m
%
% PURPOSE:
%    Example of a Three-Bar Truss system deformation due to applied loads
%
% AUTHOR: Daniel DiDomenico
%
% Modified on: Sun, 19 Oct 2008, 16:08:58 EDT
% Created on : Sun, 19 Oct 2008, 18:08:58 EST
%
% DEPENDENCIES:
%   call:
%
% REMARKS:
%
%*********************************************************************

% Three bar truss example
clear all;
e = [ 2, 4, 3 ];
A = [ 3, 1, 2 ];
L = [ 5, 5, 10 ];
P = 30;

%dof = 2;
n_nodes = 4
n_elem = 3
%total_dof = dof*n_nodes;

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

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

% two-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

%Nodal position
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_nodes
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([-9 6 -9 4])
% 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

d = d/500;

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_nodes
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([-9 6 -9 4])
% solid line
plot(xx,yy,'-r')
% hold the current figure for the next element
hold on
end

% 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')


Three-Bar Truss System Deformation

### MATLAB Results

n_nodes =

4

n_elem =

3

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 6

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

Columns 7 through 8

0            0
0            0
-0.3         -0.3
-0.3         -0.3
0            0
0            0
0.3          0.3
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


## Unstable Three-Bar Truss: MATLAB Dissection

The following code is the three-bar truss system code modified for the unstable three-bar truss system shown at the left of the figure below. To stabilize this truss, an additional member is added at shown at the right of the figure. A force, P, of magnitude 0.5 is applied as shown to evaluate the truss system stability.

Unstable Truss system stabilized by additional member

### Unstable Truss System

% Matlab Script
%*********************************************************************
% filename: three_unstable.m
%
% PURPOSE:
%    Example of an Unstable Three-Bar Truss system deformation due to applied loads
%
% AUTHOR: Daniel DiDomenico
%
% Modified on:
% Created on : Wed, 22 Oct 2008, 16:08:58 EST
%
% DEPENDENCIES:
%   call:
%
% REMARKS:
%
%*********************************************************************

% Three bar truss example
clear all;
e = [ 2, 2, 2 ];
A = [ 3, 3, 3 ];
L = [ 1, 1, 1 ];
P = 1;

n_nodes = 4
n_elem = 3

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

R = zeros(dof,1); R(5) = P;

% Assemble global stiffness matrix
K = zeros(dof);

for i = 1:elems
lm = lmm(i,:);
con = conn(i,:);
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

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

%Nodal position

position(:, 1) = [  0; 0 ];
position(:, 2) = [ L(1)*cos(theta(1)); L(1)*sin(theta(1)) ];
position(:, 3) = [ L(2)*cos(theta(2))+ L(1)*cos(theta(1)); L(2)*sin(theta(1))+L(1)*cos(theta(1)) ];
position(:, 4) = [ L(1); 0 ];

% set up the nodal coordinate arrays x, y
for i = 1 : n_nodes
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;

% 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,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

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(1); 0 ];

% set up the nodal coordinate arrays x, y
for i = 1 : n_nodes
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;

% 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

% put the title on the figure
title('Deformed and Undeformed States of An Unstable Three-Bar Truss System')
% label x axis
xlabel('x')
% label y axis
ylabel('y')


#### Breakdown of FEM

Running this MATLAB code results in a breakdown of the FEM. The system above is unstable and results in a singular stiffness matrix. FEM cannot be used to solve for the displacements and forces at each node. This breakdown is shown below.

EDU>> three_unstable

n_nodes =

4

n_elem =

3

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
0
1
0
0
0

Warning: Matrix is singular to working precision.
> In NodalSoln at 12
In three_unstable at 67

d =

0
0
Inf
NaN
NaN
NaN
0
0

reactions =

NaN
NaN
NaN
NaN


### Stabilization Due to Additional Member

To add the cross-bar member to the truss the MATLAB code is modified to the following.

% Matlab Script
%*********************************************************************
% filename: three_stable.m
%
% PURPOSE:
%    Example of an Unstable Three-Bar Truss system stabilized by an
%
% AUTHOR: Daniel DiDomenico
%
% Modified on:
% Created on : Wed, 22 Oct 2008, 17:08:58 EST
%
% DEPENDENCIES:
%   call:
%
% REMARKS:
%
%*********************************************************************

% Three bar truss example
clear all;
e = [ 2, 2, 2, 2 ];
A = [ 3, 3, 3, 3];
L = [ 1, 1, 1, sqrt(2) ];
P = .5;

n_nodes = 4
n_elem = 4

theta = [pi/2 0 -pi/2 pi/4];

nodes = [0,0;
L(1)*cos(theta(1)), L(1)*sin(theta(1));
L(2)*cos(theta(2)), L(1)*sin(theta(1));
L(2)*cos(theta(2)), 0;
];

dof = 2*length(nodes);

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 ];
elems = size(lmm,1);
K = zeros(dof);
R = zeros(dof,1);
debc = [ 1, 2, 7, 8 ];
ebcVals = zeros(length(debc),1);

R = zeros(dof,1); R(3) = P;

% Assemble global stiffness matrix
K = zeros(dof);

for i = 1:elems
lm = lmm(i,:);
con = conn(i,:);
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

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

%Nodal position

position(:, 1) = [  0; 0 ];
position(:, 2) = [ L(1)*cos(theta(1)); L(1)*sin(theta(1)) ];
position(:, 3) = [ L(2)*cos(theta(2))+ L(1)*cos(theta(1)); L(2)*sin(theta(1))+L(1)*cos(theta(1)) ];
position(:, 4) = [ L(1); 0 ];

% set up the nodal coordinate arrays x, y
for i = 1 : n_nodes
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,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

position_d(:, 1) = [  0; 0 ];
position_d(:, 2) = [ 0 + d(3); L(1)+ d(4) ];
position_d(:, 3) = [ L(4)*cos(theta(4))+ d(5); L(4)*sin(theta(4))+ d(6) ];
position_d(:, 4) = [ L(1); 0 ];

% set up the nodal coordinate arrays x, y
for i = 1 : n_nodes
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 4
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

% put the title on the figure
title('Deformed and Undeformed States of An Unstable Three-Bar Truss System Stabilized by Additional Member')
% label x axis
xlabel('x')
% label y axis
ylabel('y')


This code produces the following results.

EDU>> three_stable

n_nodes =

4

n_elem =

4

k =

2.2496e-032  3.6739e-016 -2.2496e-032 -3.6739e-016
3.6739e-016            6 -3.6739e-016           -6
-2.2496e-032 -3.6739e-016  2.2496e-032  3.6739e-016
-3.6739e-016           -6  3.6739e-016            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 -2.2496e-032 -3.6739e-016      -2.1213      -2.1213            0
2.1213       8.1213 -3.6739e-016           -6      -2.1213      -2.1213            0
-2.2496e-032 -3.6739e-016            6  3.6739e-016           -6            0            0
-3.6739e-016           -6  3.6739e-016            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.5
0
0
0
0
0

d =

0
0
0.40237
-2.4638e-017
0.31904
-0.083333
0
0

reactions =

-0.5
-0.5
0
0.5


Stabilization due to additional member

## Contributing Team Members

The following students contributed to this report:

Mark Barry Eas4200c.f08.blue.a 16:16, 24 October 2008 (UTC)

Daniel DiDomenico Eml4500.f08.ateam.didomenico 23:02, 19 October 2008 (UTC)

Sean Miller EML4500.f08.Ateam.Miller 03:49, 24 October 2008 (UTC)

Jil Paquette Eml4500.f08.Ateam.paquette 17:09, 24 October 2008 (UTC)

Matthew Philbin Eml4500.f08.Ateam.philbin 22:47, 23 October 2008 (UTC)

Ian Slevinski Eml4500.f08.Ateam.Slevinski 21:21, 23 October 2008(UTC)