User:Eml4500.f08.ateam/hw 4

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


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



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.



$$ \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
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.



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:

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

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:

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

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.



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.



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  file found at the website for the class textbook, Fundamental Finite Element Analysis and Applications.



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

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.

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.

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

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:

Resulting Stresses and Strains:

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.





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.



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.

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



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)