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

1 Connectivity and LMM Arrays
2 Transforming Local DOFs to Another Local Set
- 2.1 Derivations
- 2.2 Transformed Forces and Degrees of Freedom
3 Changing global degrees of freedom to a local coordinate system aligned along the truss
4 Verify transformation matrix transpose theorem
5 Eigenvalues of the Matrix
6 Justification of Assembly of Elemental Stiffness Matrices into the Global Stiffness Matrix
7 Five-Bar Truss System: MATLAB Dissection
8 Three-Bar Truss System: MATLAB Dissection
- 8.1 MATLAB Results
9 Unstable Three-Bar Truss: MATLAB Dissection
- 9.1 Unstable Truss System
  - 9.1.1 Breakdown of FEM
- 9.2 Stabilization Due to Additional Member
10 Contributing Team Members

Connectivity and LMM Arrays [edit]

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 [edit]

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 [edit]

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 [edit]

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 [edit]

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 [edit]

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 [edit]

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 [edit]

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⁽¹⁾₃ - f⁽²⁾₁
    Sum of forces in y direction: F = 0 = P -f ⁽¹⁾₄ - f⁽²⁾₂

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 [edit]

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 (cm²)	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,:)))]

Additional Functions [edit]

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` [edit]

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` [edit]

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` [edit]

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 [edit]

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⁽¹⁾₁	54927
F⁽¹⁾₂	1.5993x10⁵
F⁽⁴⁾₁	-54927
F⁽⁴⁾₁	-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.39x10⁵
2	-1.39x10⁵	-6.2999	-25200
3	-5.29x10⁵	-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 [edit]

% 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 [edit]

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);
 
%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
 
 
 % 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 [edit]

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 [edit]

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 [edit]

% 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);
 
%load vector
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 [edit]

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 [edit]

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 
%        additional member
%
% 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);
 
%load vector
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 [edit]

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)

User:Eml4500.f08.ateam/hw 4

Contents

Connectivity and LMM Arrays [edit]

Transforming Local DOFs to Another Local Set [edit]

Derivations [edit]

Transformed Forces and Degrees of Freedom [edit]

Changing global degrees of freedom to a local coordinate system aligned along the truss [edit]

Verify transformation matrix transpose theorem [edit]

Eigenvalues of the Matrix [edit]

Justification of Assembly of Elemental Stiffness Matrices into the Global Stiffness Matrix [edit]

Five-Bar Truss System: MATLAB Dissection [edit]

Additional Functions [edit]

`PlaneTrussElement.m` [edit]

`NodalSoln.m` [edit]

`PlaneTrussResults.m` [edit]

MATLAB Results [edit]

Plot of Undeformed and Deformed Five-Bar Truss System [edit]

Three-Bar Truss System: MATLAB Dissection [edit]

MATLAB Results [edit]

Unstable Three-Bar Truss: MATLAB Dissection [edit]

Unstable Truss System [edit]

Breakdown of FEM [edit]

Stabilization Due to Additional Member [edit]

Contributing Team Members [edit]

Navigation menu

Personal tools

Namespaces

Variants

Views

Actions

Search

Navigation

Community

Toolbox

Wikimedia projects

Print/export

User:Eml4500.f08.ateam/hw 4

Contents

Connectivity and LMM Arrays [edit]

Transforming Local DOFs to Another Local Set [edit]

Derivations [edit]

Transformed Forces and Degrees of Freedom [edit]

Changing global degrees of freedom to a local coordinate system aligned along the truss [edit]

Verify transformation matrix transpose theorem [edit]

Eigenvalues of the Matrix [edit]

Justification of Assembly of Elemental Stiffness Matrices into the Global Stiffness Matrix [edit]

Five-Bar Truss System: MATLAB Dissection [edit]

Additional Functions [edit]

PlaneTrussElement.m [edit]

NodalSoln.m [edit]

PlaneTrussResults.m [edit]

MATLAB Results [edit]

Plot of Undeformed and Deformed Five-Bar Truss System [edit]

Three-Bar Truss System: MATLAB Dissection [edit]

MATLAB Results [edit]

Unstable Three-Bar Truss: MATLAB Dissection [edit]

Unstable Truss System [edit]

Breakdown of FEM [edit]

Stabilization Due to Additional Member [edit]

Contributing Team Members [edit]

Navigation menu

Search

`PlaneTrussElement.m` [edit]

`NodalSoln.m` [edit]

`PlaneTrussResults.m` [edit]