User:Eml 4500.f08.delta 6.specht/HW4

 Comment: This page is a revision from the HW that was turned in before the deadline. Here is the comparison between these two versions. --Eml 4500.f08.delta 6.specht 17:47, 4 November 2008 (UTC)

Lecture Notes (October 6th - October 17th, 2008)

 * BEGINNING OF OCT 6 LECTURE

Connectivity & Location Matrix Master Arrays
Consider the two-bar truss below:



and its respective bar element:



We can notice that there is a common global node, and its corresponding dofs, when comparing both bar element diagrams. Since the FEA method can be applied to a large system with multiple common nodes an elements, then it is important to organize these element notations so as to visualize which elements share which nodes and bounded by which dofs.

The Connectivity Array, "conn"
This array is set up as follows:


 * $$\begin{matrix}

1 & 2 \end{matrix}$$   <-- Local node numbers


 * conn = $$\left.\ \begin{vmatrix}

1 & 2\\ 2 & 3 \end{vmatrix} \right\}$$ <-- Global node numbers, each row for each element

In short-hand notation, we say that conn = (e,j): global node number of local node j of element e.

The Location Matrix Master Array, "lmm"
This array is set up as follows:


 * $$\begin{matrix}

1 & 2 & 3 & 4 \end{matrix}$$   <-- Local dofs numbers


 * lmm = $$\left.\ \begin{vmatrix}

1 & 2 & 3 & 4\\ 3 & 4 & 5 & 6 \end{vmatrix} \right\}$$ <-- Global dofs numbers, each row for each element

In short-hand notation, we say that lmm = (i,j): equation number (global dof number) for the element stiffness coefficient corresponding to the jth local dof number.


 * END OF OCT 6 LECTURE

Development of Invertible Local DOF Transfer Function
Goal: Find $$ \tilde{T}_{4x4}^{(e)} $$ that transforms the set of local element d.o.f.'s $$ d^{(e)} $$ to another set of local element d.o.f.'s $$ \tilde{d}_{4x1}^{(e)} $$ such that $$ \tilde{T}^{(e)} $$ is invertible.



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

$$ (1) \hat{d}_{1}^{(e)} = \begin{bmatrix} l^{(e)}\ m^{(e)}\

\end{bmatrix}

\begin{Bmatrix}d_1^{(e)} \\ d_{2}^{(e)} \end{Bmatrix} $$

Where,

$$ \hat{d}_{1}^{(e)} = q_1^{(e)} $$

$$\hat{d}_{2}^{(e)} = \vec{d}_1^{(e)} \cdot \tilde{j}$$

Note: $$\hat{d}_{2}^{(e)}$$ is the component of $$\vec{d}_1^{(e)}$$ along $$\tilde{j}$$, i.e., $$\tilde{y}$$ axis.

Thus,

$$\tilde{d}_2^{(e)} = -\sin \theta^{(e)}d_1^{(e)} + \cos \theta^{(e)} d_2^{(e)}$$

Equation (2) is as follows:

$$(2) \hat{d}_{2}^{(e)} = \begin{bmatrix} -m^{(e)}\ l^{(e)}\

\end{bmatrix}

\begin{Bmatrix}d_1^{(e)} \\ d_{2}^{(e)} \end{Bmatrix}$$

Next, put (1) & (2) into matrix form,

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

Where,

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

Thus,

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

which represents,

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

Next, rotate the element achieving the following



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

Note: all zero values are zero because transverse force does not stretch the spring.

Next, we can conclude:

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

Note: this is the stopping point of transforming the FD relation back to original (x,y) coordinate system.

Derivation of the Global Force Vector from the Axial Displacements and Forces
$$\tilde{f}^{(e)}_{4x1} = \tilde{k}_{4x4} ^{(e)} * \tilde{d}_{4x1}^{(e)}$$ Note: Consider the case where $$ \tilde{d}_{4} \neq 0 $$ then $$ \tilde{d}_1^{(e)} = \tilde{d}_2^{(e)} = \tilde{d}_3^{(e)} = 0 $$

you will get the following relationship $$ \tilde{f}_{4x1} ^{(e)} = \tilde{k}_{4x4} ^{(e)} * \tilde{d}_{4x1} ^{(e)} = 0_{4x1} $$ where the zero matrix is the 4th column of matrix $$ k^{(e)} $$ (interpretation of the transverse dofs)

Recall from 19-3 that $$ \tilde{d}^{(e)} = \tilde{T}^{(e)} * d^{(e)} $$ similarly there is the same relationship for f where $$\tilde{f}^{(e)} = \tilde{T}^{(e)} * f^{(e)} $$

Also: $$ \tilde{k}^{(e)} * \tilde{d}^{(e)}  = \tilde{f}^{(e)}$$

The lesson here is that we should use the transverse dofs instead of just using the inverse dofs therefore finding the relationship for the global matrix $$ k=d*f $$ $$\tilde{k}^{(e)} *\tilde{T}^{(e)} * d^{(e)}= \tilde{T}^{(e)} * f^{(e)}$$

If $$ T^{(e)}$$ is invertible, then the matrix that we are looking for is:

$$f^{(e)} = \begin{bmatrix} \tilde{T}^{(e)^{-1}} * \tilde{k}^{(e)} * \tilde{T}^{(e)}*d^{(e)} \end{bmatrix}$$ where $$T^{(e)}$$ is the block diagram matrix.

Considering a general block-diagram matrix A:

$$ A=

\begin{bmatrix}

D_{1} & 0 & 0 & 0\\

0 & . & 0 & 0\\

0 & 0 & . & 0\\

0 & 0 & 0 & D_{s}

\end{bmatrix} $$

so what will be the inverse of matrix A? This can be answered by observing the simple example of the diagonal matrix B:

$$ B=

\begin{bmatrix}

d_{11} & 0 & 0 & 0\\

0 & d_{22} & 0 & 0\\

0 & 0 & . & 0\\

0 & 0 & 0 & d_{mm}

\end{bmatrix} $$

$$ B= Diag \begin{bmatrix} d_{11}, d_{22}, ... , d_{mm} \end{bmatrix} $$

$$ B^{-1} = Diag\begin{bmatrix} \frac{1}{d_{11}}, \frac{1}{d_{22}}, ... , \frac{1}{d_{mm}} \end{bmatrix} $$

Assuming $$ d_{ii} \neq 0 $$  for i=1,…n

For a block diagonal matrix A

$$ A = Diag\begin{bmatrix} D_{11}, D_{22}, ... , D_{s} \end{bmatrix} $$

$$ A^{-1} = Diag\begin{bmatrix} D_{11}^{-1}, D_{22}^{-1}, ... , D_{s}^{-1} \end{bmatrix} $$

$$\tilde{T}^{(e)^{-1}} = Diag \begin{bmatrix} R^{(e)^{-1}}, ... , R^{(e)^{-1}} \end{bmatrix}$$

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

Since $$l^{2}*m^{2} = 1 $$ the matrix becomes where $$ I_{2x2} $$ is the identity matrix.

$$R_{2x2}^{(e)^{T}}*R_{2x2}^{(e)} = \begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix}_{2x2} = I_{2x2} $$

$$ R^{(e)^{T}}= R^{(e)^{-1}} $$ $$ \tilde{T}^{(e)^{-1}} = Diag\begin{bmatrix} R^{(e)^{T}}, R^{(e)^{T}} \end{bmatrix}$$ where $$ \tilde{T}^{(e)} = \begin{pmatrix} Diag \begin{bmatrix} R^{(e)}, R^{(e)} \end{bmatrix} \end{pmatrix}^{T} $$ $$ \tilde{T}^{(e)^{-1}} = \tilde{T}^{(e)^{T}} $$ $$ \begin{bmatrix} \tilde{T}^{(e)^{T}} * \tilde{k}^{(e)}* \tilde{T}^{(e)} \end{bmatrix} *d^{(e)}= f^{(e)} $$ where $$ k^{e} = \begin{bmatrix} \tilde{T}^{(e)^{T}} * \tilde{k}^{(e)}* \tilde{T}^{(e)} \end{bmatrix} $$

Mechanisms
The dynamics eigenvalue problem is given as: $$Kv=\lambda M v$$ In this equation, the mass matrix M is given as the identity matrix. This problem can then be reduced to: $$Kv=\lambda v$$

The global K matrix for the two bar truss system is:

$$K= \ \begin{bmatrix} K_{11}^{(1)} & K_{12}^{(1)} & K_{13}^{(1)} & K_{14}^{(1)} & 0 & 0  \\ K_{21}^{(1)}& K_{22}^{(1)} & K_{23}^{(1)} & K_{24}^{(1)} & 0 & 0 \\ K_{31}^{(1)} & K_{32}^{(1)} & (K_{33}^{(1)} + K_{11}^{(2)}) & (K_{34}^{(1)} + K_{12}^{(2)}) & K_{13}^{(2)} & K_{14}^{(2)} \\ K_{41}^{(1)}& K_{42}^{(1)} & (K_{43}^{(1)} + K_{21}^{(2)}) & (K_{44}^{(1)} + K_{22}^{(2)}) & K_{23}^{(2)} & K_{24}^{(2)} \\ 0 & 0 & K_{31}^{(2)} & K_{32}^{(2)} & K_{33}^{(2)} & K_{34}^{(2)} \\ 0 & 0 & K_{41}^{(2)} & K_{42}^{(2)} & K_{43}^{(2)} & K_{44}^{(2)}  \end{bmatrix}_{6X6}\,$$ $$\ \ = \begin{bmatrix} \frac{9}{16}& \frac{3\sqrt{3}}{16} & -\frac{9}{16} & -\frac{3\sqrt{3}}{16} & 0 & 0 \\  \frac{3\sqrt{3}}{16} & \frac{3}{16} & -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & 0 & 0 \\ -\frac{9}{16} & -\frac{3\sqrt{3}}{16} & (\frac{9}{16} + \frac{5}{2}) & (\frac{3\sqrt{3}}{16} + \frac{5}{2}) & -\frac{5}{2} & -\frac{5}{2} \\ -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & (\frac{3\sqrt{3}}{16} + \frac{5}{2}) & (\frac{3}{16} + \frac{5}{2}) & -\frac{5}{2} & -\frac{5}{2} \\  0 & 0 & -\frac{5}{2} & -\frac{5}{2} & \frac{5}{2} & \frac{5}{2} \\ 0 & 0 & -\frac{5}{2} & -\frac{5}{2} & \frac{5}{2} & \frac{5}{2}  \end{bmatrix}_{6X6}\,$$ $$\ \ = \begin{bmatrix} 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.5 & 2.5\\  -0.3248& -0.1875 & -2.1752 & 2.6875 & 2.5 & -2.5\\  0& 0 & -2.5 & 2.5 & 2.5 & -2.5\\ 0 & 0 & 2.5 & -2.5 & -2.5 & 2.5 \end{bmatrix}_{6X6}\,$$

When this matrix is solved for the eigenvalues, the values found are given as: $$\lambda = \begin{bmatrix} 0&0 &0  &0  & 1.4707 & 10.0294 \end{bmatrix}$$

Notice that there are four eigenvalues that are equal to zero. These correspond to the three rigid body displacements and one mechanism of the two bar truss system. Two of the three rigid body motions are the translational displacements in the two dimension cartesian coordinate system. The other rigid body motion is the rotation of the system around global node 2. The mechanism of the system can be described as the motion of each element as they pivot about node 2. The mechanism can also be described as the change in the angle between elements 1 and 2 at node 2, due to the applied force, and the displacement of node 2. If the row vector consisting of the eigenvectors that correspond to the four zero eigenvalues is given as $$\begin{Bmatrix}u_{1} & u_{2} & u_{3} & u_{4}\end{Bmatrix}$$, this equation can be expressed as (for i=1:4):

$$K_{6x6}u_{i,6x1}=0_{6x1}u_{i,6x1}$$ The linear combination of each of the eigenvectors can be equated to: $$ \sum_{i=1}^{4}{\alpha _{i}u_{i,6x1}}=:W_{6x1} $$

The value of $$\alpha_{i}$$ can be any real scalar value. The symbol used in the equation( =:) makes the right hand side equal to the left hand side by definition. The column vector $$W$$ is the eigenvector symbolizing the sum of all four of the eigenvectors corresponding to the eigenvalues ($$u_{i}$$). Thus the dynamics eigenvalue problem can be expressed as the following equation for the eigenvectors corresponding to the zero eigenvalues.

$$ K_{6x6}W_{6x1}=K_{6x6}(\sum_{i=1}^{4}{\alpha _{i}u_{i,6x1}}) $$ $$ K_{6x6}W_{6x1}=\sum_{i=1}^{4}{\alpha _{i}(K_{6x6}u_{i,6x1}}) $$ Because the values in parentheses is proven to be equal to a 6x1 column vector of zeros, the dynamics eigenvalue problem can be written as the dot product of the zero vector and the column vector $$W$$, which is the sum of the eigenvectors corresponding to the zero eigenvalues. $$ K_{6x6}W_{6x1}=0_{6x1}\cdot W_{6x1} $$

Justification of Assembly of Elemental Stiffness Matrix into Global Stiffness Matrix
Consider the example of the 2-bar truss system discussed earlier. If we recall the elemental Force-Discplacement relationship:

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

The equilibrium of node 2 can be found using three methods we have discussed earlier this semester. 1. The Euler Cut Principle 2. The Free Body Diagrams of each element 3. The Identity of Global D.O.F.'s to Element D.O.F.'s

At node 2, the elquilibrium can be found by simple summing the forces in each direction.



$$ \sum{F_{x}} = 0 = -f_{3}^{(1)} - f_{1}^{(2)} = 0 $$ $$ \sum{F_{y}} = 0 = P -f_{4}^{(1)} - f_{2}^{(2)} = 0 $$

When equations one and two, from the equilibrium of node two, are rewritten they become: $$f^{(1)}_{3} + f^{(2)}_{1}= 0$$ $$f^{(1)}_{4} + f^{(2)}_{2}= 0$$ The element forces can then be written in terms of the element stiffness and displacements. $$f_{3}^{(1)} = \begin{bmatrix} k_{31}^{(1)}d_{1}^{(1)} + k_{32}^{(1)}d_{2}^{(1)} + k_{33}^{(1)}d_{3}^{(1)} + k_{34}^{(1)}d_{4}^{(1)} \end{bmatrix}$$ $$f_{1}^{(2)} = \begin{bmatrix} k_{11}^{(2)}d_{1}^{(2)} + k_{12}^{(2)}d_{2}^{(2)} + k_{13}^{(2)}d_{3}^{(2)} + k_{14}^{(2)}d_{4}^{(2)} \end{bmatrix}$$ The forces written in terms of the elements stiffness and displacements are then substituted into equation one and equation one becomes. $$ \begin{bmatrix} k_{31}^{(1)}d_{1}^{(1)} + k_{32}^{(1)}d_{2}^{(1)} + k_{33}^{(1)}d_{3}^{(1)} + k_{34}^{(1)}d_{4}^{(1)} \end{bmatrix} + \begin{bmatrix} k_{11}^{(2)}d_{1}^{(2)} + k_{12}^{(2)}d_{2}^{(2)} + k_{13}^{(2)}d_{3}^{(2)} + k_{14}^{(2)}d_{4}^{(2)} \end{bmatrix}=0$$ Now recall the relationship between the global degrees of freedom and the element degrees of freedom: $$d_{1}=d_{1}^{(1)},d_{2}=d_{2}^{(1)},d_{3}=d_{3}^{(1)}=d_{1}^{(2)},d_{4}=d_{4}^{(1)}=d_{2}^{(2)},d_{5}=d_{3}^{(2)},d_{6}=d_{4}^{(2)}$$ The global degrees of freedom can now be substituted into equation one for the element degrees of freedom and equation one becomes. $$ \begin{bmatrix} k_{31}^{(1)}d_{1} + k_{32}^{(1)}d_{2} + k_{33}^{(1)}d_{3} + k_{34}^{(1)}d_{4} \end{bmatrix} + \begin{bmatrix} k_{11}^{(2)}d_{3} + k_{12}^{(2)}d_{4} + k_{13}^{(2)}d_{5} + k_{14}^{(2)}d_{6} \end{bmatrix}=0$$ Now the common terms are grouped and factored (i.e. The terms that include $$d_{3}$$ and $$d_{4}$$). The same process is repeated for equation two. Equation one results in the third row of the global stiffness matrix and equation two results in the fourth row of the global stiffness matrix. The work is shown in the homework for this section.

Assembly of the global stiffness matrix
Before starting the assembly of the global stiffness matrix a quick review of the goals and the end that is being worked towards is necessary. The first goal is the use of The Principle of Virtual Work (PVW) which allows the formation of the global stiffness matrixK by eliminating rows that correspond to the boundary conditions. The ultimate goal is to use partial differential equations (PDE's) to derive the Finite Element Method (FEM). Recall now two equations that have already been derived in the process of the FEM. $$q^{(e)}_{2x1}$$ = $$T^{(e)}_{2x4}$$ $$d^{(e)}_{4x1}$$ $$K^{(e)}= T^{(e)^(T)} * \hat{K}^{(e)} * T^{(e)}$$

The assembly of the global stiffness matrix (K) is accomplished using the element stiffness matrices ($$k^{(e)}$$) and an assemble operator (A). Where e= 1,....,nel, and nel is the number of elements. The assembly operator is introduced because the sigma operator is not correct. The global stiffness matrix is assembled using a combination of element stiffness matrices. The element stiffness matrices are small when compared to the global stiffness matrix. The global stiffness matrix could be orders of magnitude larger than the element stiffness matrices. The sum of $$k_{nedxned}^{(e)}$$ doesn't equal to a $$K_{nxn}$$ matrix. Therefore the assembly operator (A) is necessary for the notation of the process to be correct. $$K=Ak^{(e)}$$ With this, new notation needs to be introduced. n is equal to the number of global degrees of freedom before eliminating the boundary conditions. ned is the number of element degrees of freedom. Where ned is very small when compared to n. Recall that the force displacement relationship for a bar or spring is kd=F. After some algebra: $$ kd-F=0$$ This will be called equation three and is equivalent to: w(kd-F)=0 This will be called equation four. This holds true for all values of w. Therefore it is trivial to imply that equation three is equal to equation four but is extremely important to understand that equation four is equal to equation three. This is true because equation four is valid for all w. So letting w=1, equation four is converted back into equation three.

Five Bar Truss Structure (MATLAB)
The code that was used to run the 5-Bar Truss System is listed as follows The first section of this code defines the material properties for each element. Since some of the values are used for more than one section, it was only necessary to declare them once. The variable nodes describes the position of each of the global nodes, which are all multiplied by a scalar value (1000) to get the actual position. The two functions that follow it are the Connectivity array (conn) and the Location Master Matrix array (lmm). The connectivity array holds the number of the global node that is being described when the following is called: conn(i,j) = global node number This function will give the global node number for the ith element, and the jth local node of that element. This is necessary to reduce the amount of storage space for the program, and to allow for the use of for loops in the body of the program. The Location Master Matrix array returns the global dof number for the ith element and the jth local dof for the element. This can be called in the same manner that the Connectivity array is called. The Connectivity array and the Location Master Matrix array can both be used to assemble the global stiffness matrix for the structure. The code in this program is divided into two sections in order to construct the global K matrix. This is done in order to use for loops to construct the matrix, since multiple matrices have the same values for the modulus of elasticity (e in the program) and for the area (a). There are three sets of values for the area, which is the maximum, so this is the limiting factor of the program. Therefore construction of this matrix must take place in at least three parts, which is done in the program for each set of corresponding area values and modulus of elasticity. The function PlaneTrussElement is discussed later. After the global stiffness matrices are constructed, the applied load that is imposed on the structure is defined. Once the load is defined, the last several lines in the program call the functions that are discussed below (PlaneTrussElement, PlaneTrussResults, NodalSoln) in order to compute the displacements and reactions at each node. There are several functions which are also called by the program in order to compute the values for the displacements, the reactions, and the results. In order to do this, the functions must be called several times, based on the fact that there are three different values for area, and two different values for the modulus of elasticity. The first function is called PlaneTrussElement, which is used to generate the element stiffness matrices for each element in the structure. In order to use this function in the body of the program values for the modulus of elasticity, area, and the coordinates at the endpoints of the element must all be included in the declaration. This then computes the length of the element, which is in turn followed by the computation of the element values ls and ms. These values are the horizontal and vertical components of the element, respectively.

The second function is the PlaneTrussResults function, and requires that the modulus, area, coordinates of the endpoints, and displacements of the endpoints be declared when the function is called. This function computes the axial displacement for the element and then uses this value in order to compute the axial force on the element. In order to do this the length is computed first, followed by the horizontal and vertical components, similar to the manner described for the first function. The transfer function (T) is then computed, which allows the displacements to be translated from the cartesian coordinate system, to the coordinate system in which the x axis is pointed from local node 1 to node 2 along the length of the element. The strain is then computed, and multiplied by the modulus in order to obtain the stress on the element. By multiplying the stress on the element by the area of the element, the axial force on the element is found, and returned to the body of the program. The last function that is called by the program is the NodalSoln function, which requires that the global stiffness matrix, the applied force vector, the element dofs, and the values for these dofs all be declared when the function is called. The number of global dofs is set to the length of the applied force vector. The value df is then set to the difference between the global dofs and the number of dofs for the element, which is iterated from 1 to the total number of global dofs. The value in the global stiffness matrix at the (df, df) position is then found. A scalar force value is established by the value in the df position of the applied force matrix minus the value of the value in the global stiffness matrix at the (df,debc) position, which is multiplied by the value of the dof before being subtracted. The displacement vector for the element is then established by setting the value in the debcth row of the specified value for the displacement, and the value in dof-debcth row to the value of the scalar force value divided by the value from the global stiffness matrix that was found earlier. The resultant force for each node is then computed by multiplying the debcth row of the global stiffness matrix by the displacement vector that was constructed, and then subtracting the value from the debcth row in the applied force vector. This value is then returned to the main body of the program.

The results from the code that is given above are listed below. First the global stiffness matrix is computed and returned. The second value given in the value of the applied force vector. After that, the displacement vector for the system is given, followed by the reactions at the fixed nodes, and the resultant forces through each element. The code that is listed below was generated in order to plot the undeformed and deformed shapes of the 5-Bar Truss Structure. This first calculates the original position of each node. Then the original position of each node is added to the vector d from the result of the code above, which is given as the displacement of each node. In order to better display the effect of the force on the global node 2, the displacement vector, d, has been multiplied by a weighting factor of 500. If this weighting factor were not present, the shape of the deformed structure would very closely coincide with the shape of the undeformed structure. The weighting factor thus allows the effect of the applied force on the system to be better understood visually, while the true structure would not deform as much as is shown. The last part of the code is necessary to include the labels at the nodes and on the elements of the structure, as well as including the key, which gives the linetype for the undeformed and deformed structures.

Three Bar Truss Example (MATLAB)
It is important to understand how to solve a bar truss system employing a finite element MATLAB code in order to compare the resulting matrices with those obtained from hand calculations.

Take this truss example with the following:

Data

MATLAB Code (Analysis of the system)

MATLAB Functions

In order to run this MATLAB code the following functions are needed:

Note: this MATLAB code was obtained from the MATLAB.zip Chapter 9 folder given by the Student Companion Site Bhatti: Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations website http://bcs.wiley.com/he-bcs/Books?action=resource&bcsId=3105&itemId=0471648078&resourceId=7629

Note: this MATLAB code was obtained from the MATLAB.zip Common folder given by the Student Companion Site Bhatti: Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations website http://bcs.wiley.com/he-bcs/Books?action=resource&bcsId=3105&itemId=0471648078&resourceId=7629

Note: this MATLAB code was obtained from the MATLAB.zip Chapter 9 folder given by the Student Companion Site Bhatti: Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations website http://bcs.wiley.com/he-bcs/Books?action=resource&bcsId=3105&itemId=0471648078&resourceId=7629

Results MATLAB Code (Plot of the system)

Derivation of $$\tilde{d}_{2}^{(e)}$$
The vector $$\tilde{d}_{2}^{(e)}$$ is the translational vector of node 1 for each element, while the vector $$\tilde{d}_{1}^{(e)}$$ is the axial displacement of the node. The value for the first vector can be found by multiplying the displacement of node 1 (given in standard cartesian coordinates) by the unit vector in the $$\tilde{y}$$ direction ($$\vec{\tilde{j}}$$). The general form of this equation is given as: $$ \tilde{d}_{2}^{(e)}=\vec{d}_{Local Node 1}^{(e)}\cdot \vec{\tilde{j}} $$

The equations for the variables on the right hand side of the equation are given as the following two equations. (1)

$$ \vec{d}_{Local Node 1}^{(e)}=d_{1}^{(e)}\vec{i}+d_{2}^{(e)}\vec{j} $$

(2)

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

Plugging the values for equations (1) and (2) into the general equation given above yields the equation for $$\tilde{d}_{2}^{(e)}$$ in terms of $$d_{1}^{(e)}$$ and $$d_{2}^{(e)}$$ for each element.

$$ \tilde{d}_{2}^{(e)}=(-sin\theta ^{(e)}\vec{i}\cdot d_{1}^{(e)}\vec{i})+(cos\theta ^{(e)}\vec{j}\cdot d_{2}^{(e)}\vec{j} ) $$

This can be reduced by noting the following:

$$ \vec{i}\cdot \vec{i}= 1 $$

$$ \vec{j}\cdot \vec{j}= 1 $$

When inputting the values for these dot products the equation reduces to: $$ \tilde{d}_{2}^{(e)}=(-sin\theta ^{(e)}\cdot d_{1}^{(e)})+(cos\theta ^{(e)}\cdot d_{2}^{(e)}) $$

Since the values for $$-sin\theta^{(e)}$$ and $$cos\theta^{(e)}$$ can be given as $$-m^{(e)}$$ and $$l^{(e)}$$ respectively, this can be reduced to the following equation:

$$ \tilde{d}_{2}^{(e)}=\begin{bmatrix} -m^{(e)} & l^{(e)} \end{bmatrix}\begin{Bmatrix} d_{1}^{(e)}\\ d_{2}^{(e)} \end{Bmatrix} $$

Verify the element stiffness matrix
Start with equation 1: $$k^{(e)}={\tilde{T}^{(e)T}}{\tilde{k}^{(e)}}{\tilde{T}^{(e)}}$$ The substitute in with: equation 2 $$R^{(e)T}=\begin{bmatrix} l^{(e)} & -m^{(e)}\\ m^{(e)} & l^{(e)} \end{bmatrix}$$ and equation 3 $$R^{(e)}=\begin{bmatrix} l^{(e)} & m^{(e)}\\ -m^{(e)} & l^{(e)} \end{bmatrix}$$ and equation 4 $${\tilde{T}^{(e)T}}=Diag\begin{bmatrix} R^{(e)T} & R^{(e)T}\end{bmatrix}$$ and equation 5 $${\tilde{T}^{(e)}}=Diag\begin{bmatrix} R^{(e)} & R^{(e)}\end{bmatrix}$$ also $${\tilde{k}^{(e)}}=k^{(e)}\begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix}$$ From previous work $$l^{(1)}=cos(30)=.866025$$ and $$m^{(1)}=sin(30)=.5$$. so $$R^{(e)T}=\begin{bmatrix} .866025 & -.5\\ .5 & .866025 \end{bmatrix}$$ and $$R^{(e)}=\begin{bmatrix} .866025 & .5\\ -.5 & .866025 \end{bmatrix}$$ and $$k^{(e)}=\begin{bmatrix} .5625 & 3248 & -.5625 & -.3248\\ .3248 & .1875 & -.3248 & -.1875\\ -.5625 & -.3248 & .5625 & .3248\\ -.3248 & -.1875 & .3248 &.1875 \end{bmatrix}$$ Now the known matrices are substituted into equations 2, 3, 4 and 5. Then equations 2, 3, 4 and 5 are substituted into equation 1 to verify that: $$k^{(e)}=\begin{bmatrix} .5625 & 3248 & -.5625 & -.3248\\ .3248 & .1875 & -.3248 & -.1875\\ -.5625 & -.3248 & .5625 & .3248\\ -.3248 & -.1875 & .3248 &.1875 \end{bmatrix}$$

Proof of Global Transfer Function Use with Element Force Matrix
Known formulas: equation 2 $${\tilde{f}^{(e)}}={\tilde{k}^{(e)}}{\tilde{d}^{(e)}}$$ equation 3 $${\tilde{d}^{(e)}}={\tilde{T}^{(e)}}d^{(e)}$$ equation 4, where I is the identity matrix $${\tilde{k}^{(e)}}=k^{(e)}I$$ equation 5 $$f^{(e)}=k^{(e)}d^{(e)}$$ Substituting equation 2 into equation 1: $${\tilde{f}^{(e)}}={\tilde{k}^{(e)}}{\tilde{T}^{(e)}}d^{(e)}$$ The substituting equation 3 and equation 4 into equation 1: $${\tilde{f}^{(e)}}=k^{(e)}I{\tilde{T}^{(e)}}d^{(e)}$$ Finally, equation 5 is substituted into equation 1 and the Identity matrix is dropped. This gives the wanted results. $${\tilde{f}^{(e)}}={\tilde{T}^{(e)}}f^{(e)}$$

Plot of the Eigenvectors Corresponding to the Zero Eigenvalues
Using MATLAB, the output of the Global Stiffness Matrix is:

The V and D matrices below correspond to the eigenvectors and eignevalues, respectively.

As can be seen by the D matrix, the first four columns of the V matrix correspond to the eigenvectors of the zero eigenvalues.

Columns 1 & 2

Columns 3 & 4

Completion of $$f_1^{(2)} $$
The force component in the x direction of the local node 1 of element 2, $$f_1^{(2)}$$ can be expressed as:

$$f_1^{(2)} = [k_{11}^{(2)}d_1^{(2)} + k_{12}^{(2)}d_2^{(2)} + k_{13}^{(2)}d_3^{(2)} + k_{14}^{(2)}d_4^{(2)}]$$

Derivation of Fourth Row of Global Force Vector
$$f_4^{(1)} + f_2^{(2)} = P$$ =F_4

Where,

1. $$f_4^{(1)} = [k_{41}^{(1)}d_1^{(1)} + k_{42}^{(1)}d_2^{(1)} + k_{43}^{(1)}d_3^{(1)} + k_{44}^{(1)}d_4^{(1)}]$$

2. $$f_2^{(2)} = [k_{21}^{(2)}d_1^{(2)} + k_{22}^{(2)}d_2^{(2)} + k_{23}^{(2)}d_3^{(2)} + k_{24}^{(2)}d_4^{(2)}]$$

Thus,

$$f_{3}^{(1)}+f_{1}^{(2)}=[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}^{(2)}d_3^{(2)} + k_{24}^{(2)}d_4^{(2)}]$$ $$f_{3}^{(1)}+f_{1}^{(2)}=[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}^{(2)}d_5 + k_{24}^{(2)}d_6]$$

Observation: Last two items in 1. & first two items in 2. correspond to the 4th row in the global stiffness matrix K6x6.

Hence,

$$[(k_{43}^{(1)} + k_{21}^{(2)})(k_{44}^{(1)} + k_{22}^{(2)})]$$

Three-Bar Truss Eigenvalues
We have two truss systems, one 3-bar and 5-bar. They are illustrated below:



Where $$E^{(1)}=E^{(2)}=E^{(3)}=1$$; $$A^{(1)}=A^{(2)}=A^{(3)}=1$$; $$L^{(1)}=L^{(2)}=L^{(3)}=1$$. It is possible to analyze the eigenvalues and eigen vectors according to $$\mathbf{K}v=\mathbf{\lambda} v$$ for the 3-bar truss case.

We construct the stffness matrix from each individual element stiffness matrix. It was found to be:
 * K=$$\begin{bmatrix}

0&0&0&0&0&0&0&0\\ 0&1&0&-1&0&0&0&0\\ 0&0&1&0&-1&0&0&0\\ 0&-1&0&2&0&0&0&-1\\ 0&0&-1&0&1&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&-1&0&0&0&1 \end{bmatrix}$$ This matrix was then inputed into MATLAB and its eigenvalues and eigenvectors were determined. The MATLAB results are as follow:

V = -0.0000   0.0000    0.0000   -0.4740    0.7103   -0.2504    0.4082    0.0000   -0.7071   -0.4674   -0.1769    0.2838   -0.0000   -0.7071         0   -0.2136   -0.3585   -0.5702   -0.8165   -0.0000    0.0000   -0.4674   -0.1769    0.2838    0.0000    0.7071    0.0000   -0.2136   -0.3585   -0.5702   -0.0000    0.0000   -0.0000   -0.0556    0.1402   -0.1014   -0.0000    0.0000    0.0000   -0.1604    0.3534   -0.1873    0.4082    0.0000    0.7071   -0.4674   -0.1769    0.2838 D = 3.0000        0         0         0         0         0         0    2.0000         0         0         0         0         0         0    1.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

Contributing Team Members
Iain Specht--Eml 4500.f08.delta 6.specht 23:57, 16 October 2008 (UTC) Micheal Krueger--Eml4500.f08.delta 6.krueger 20:29, 19 October 2008 (UTC) Nicolas Castrillon--Eml4500.f08.delta 6.castrillon 15:56, 21 October 2008 (UTC) Diana Guzman--Eml4500.f08.delta 6.guzman 06:26, 24 October 2008 (UTC) Sean Lopez--Eml4500.f08.delta 6.lopez 07:51, 24 October 2008 (UTC) Eric Ramirez--Eml4500.f08.delta 6.ramirez 21:03, 24 October 2008 (UTC)