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

Connectivity and Location Matrix Master Arrays
The elements of a two bar truss system are shown below by Figures 1 and 2.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

The resulting equation is shown below.

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

Next, we will consider the case when $$ \tilde{d}_{4}^{(e)}\neq 0$$ The other degrees of freedom we are going to consider to be zero. $$\tilde{d}_{1}^{(e)}=\tilde{d}_{2}^{(e)}=\tilde{d}_{3}^{(e)}=0$$ Now, the force displacement relation for the case above is shown below. $$\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{k}}^{(e)}\mathbf{\tilde{d}}^{(e)}=\mathbf{0}$$ Where $$\mathbf{\tilde{f}}^{(e)}$$ is a 4X1 matrix, $$\mathbf{\tilde{k}}^{(e)}$$ is a 4X4 matrix, $$\mathbf{\tilde{d}}^{(e)}$$ is a 4X1 matrix and $$\mathbf{0}$$ is a 4X1 matrix and the 4th column of $$\mathbf{\tilde{k}}^{(e)}$$.

As derived before, the following relation of the degrees of freedom is shown below. $$\mathbf{\tilde{d}}^{(e)}=\mathbf{\tilde{T}}^{(e)}\mathbf{d}^{(e)}$$ Similarly, we can show the derivation for the following force relation. $$\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{T}}^{(e)}\mathbf{f}^{(e)}$$ Looking once again at the force displacement relation $$\mathbf{\tilde{f}}^{(e)}=\mathbf{\tilde{k}}^{(e)}\mathbf{\tilde{d}}^{(e)}$$

The force displacement relation can be rewritten as shown, with the understanding that $$\mathbf{\tilde{T}}^{(e)}$$ has to be invertible. $$\left[\mathbf{\tilde{T}}^{(e)-1}\mathbf{\tilde{k}}^{(e)}\mathbf{\tilde{T}}^{(e)}\right]\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$$

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

In order to see how to solve for the inverse of $$\mathbf{\tilde{T}}^{(e)}$$, we are going to consider the following general block diagram matrix. $$ \mathbf{A}=\begin{bmatrix}\mathbf{D}_{1} & 0\\0 & \mathbf{D}_{s}\end{bmatrix}$$

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

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

Assuming the following: $$d_{ii}\neq0$$ for $$ i=1,2,...,n$$,

The inverse of the matrix $$\mathbf{B}$$ is shown below $$\mathbf{B}^{-1}=Diag\left[\frac{1}{d_{11}},\frac{1}{d_{22}},...,\frac{1}{d_{nn}}\right] $$.

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

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

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

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

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

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

The matrix $$\mathbf{R}^{(e)}$$ as well as its transpose are shown below. $$\mathbf{R}^{(e)}=\begin{bmatrix}l^{(e)} & m^{(e)} \\ -m^{(e)} & l^{(e)}\end{bmatrix} $$

$$\mathbf{R}^{(e)T}=\begin{bmatrix}l^{(e)} & -m^{(e)} \\ m^{(e)} & l^{(e)}\end{bmatrix} $$ The multiplication of $$\mathbf{R}^{(e)}$$ and $$\mathbf{R}^{(e)T}$$ is shown below. $$\mathbf{R}^{(e)}\mathbf{R}^{(e)T}=\begin{bmatrix}l^{(e)} & m^{(e)} \\ -m^{(e)} & l^{(e)}\end{bmatrix}\begin{bmatrix}l^{(e)} & -m^{(e)} \\ m^{(e)} & l^{(e)}\end{bmatrix}=\begin{bmatrix}l^{(e)} + m^{(e)2}&0 \\ 0 & l^{(e)2} + m^{(e)2}\end{bmatrix}$$ This means that the multiplication is equal to the identity matrix. $$\mathbf{R}^{(e)}\mathbf{R}^{(e)T}=\mathbf{I}$$ Since it is equal to the identity matrix, the transpose has to be equal to the inverse of $$\mathbf{R}^{(e)}.$$

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

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

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

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

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

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

To solve an eigenvalue problem, the equation:

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

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

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

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

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

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

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

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

Justification of Global Stiffness Matrix Assembly
For the justification of assembly of the element stiffness matrix into the global stiffness matrix, consider the two-bar truss system.

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

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

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

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

$$

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

$$

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

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

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

Application of Force-Displacement Relation to Statics Analysis of 2-Bar Truss System
Equation 1 from above can be rewritten as follows: $$f_3^{(1)}+f_1^{(2)}=0\; \; \;(1a)$$

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

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

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

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

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

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

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

New Nomenclature: n: total number of global degrees of freedom before elimination of boundary conditions ned: number of element degrees of freedom
 * (ned << n)

A: assembly operation
 * Principal of Virtual Work (PVW)
 * Used to eliminate the rows for corresponding boundary conditions to obtain $$\mathbf{\bar{k}}_{2x2}$$

Deriving Finite Element Method for Partial Differential Equations
The force displacement (FD) relation for a bar or spring is $$kd=F$$. This implies that $$kd-f=0\; \; \;(3)$$ Which is equivalent to $$w(kd-F)=0\; \; \;(4)$$ for all values of w.

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

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

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

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



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

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

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

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

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

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

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

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

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

K =

1.0e+005 *

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

R =

0          0           0     -150000           0           0           0           0

d =

0        0    0.5390   -0.9531    0.2647   -0.2647         0         0

reactions =

1.0e+005 *

0.5493   1.5993   -0.5493   -0.0993

results =

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

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

The values for the global displacement matrix are shown below.

$$

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

$$

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

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



The source code added to the original MATLAB script in order to create the above plot is shown below.

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



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

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

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

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

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

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

Results of MATLAB Analysis on 3-Bar Truss System
Running the MATLAB code for solving a 3-bar truss system produced the following results.

k =

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

k =

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

k =

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

K =

Columns 1 through 7

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

Column 8

0           0         -0.3         -0.3            0            0          0.3          0.3

R =

0    0     0    30     0     0     0     0

d =

0           0      -435.17       671.77            0            0            0            0

reactions =

42.588      24.588       28.392       16.392      -70.981      -70.981

results =

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

The values for the global displacement matrix are shown below.

$$

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

$$

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



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



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

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

k =

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

k =

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

k =

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

K =

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

R =

0    0     0     1     1     0     0     0

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

d =

0    0   Inf NaN NaN NaN 0    0

reactions =

NaN NaN NaN NaN

results =

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

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

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

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

The following code was used to plot these mode shapes.

Which produced a graph as follows:

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



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

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

k =

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

k =

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

k =

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

k =

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

K =

Columns 1 through 7

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

Column 8

0           0            0            0            0           -6            0            6

R =

0           0            0            0          0.5            0            0            0

d =

0           0      0.31904            0      0.31904    -0.083333            0            0

reactions =

-0.5        -0.5            0          0.5

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



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

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

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

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