User:Eml4500.f08.ateam.nobrega/Homework Report Four

Consider the two bar truss from meeting 5 page 6:




The connectivity array,"conn" is featured in the figure above. Conn(e,j) is the global node number of local node j of element e. The rows in the connectivity array refer to the element number. For example, row 1 refers to element 1 and row 2 refers to element 2. However, the numbers in the matrix correspond to the global degree of freedom numbers. Each column refelects a local node number.

The location matrix master array (lmm)is featured in the figure to the right. lmm(i,j) is the equivalent number (global dof number) for the element of stiffness coefficient which corresponds to the jth local dof number. The rows in the mtarix correspond to the element numbers; while, the colums identify the local degree of freedom number. The numbers in the array are the global degree of freedom numbers NOT the element degree of freedom numbers.

Homework: Modify the MATLAB code for the 2-Bar truss system to solve the 3-Bar truss system. The solution to this homework assignment can be found in the section titled Three Bar Truss Deformation, of this homework report.

Transform a system with four dofs into a system with also four dofs (instead of only two dofs) so that the transform matric is now a 4x4 matrix and (hopefully) invertible.

Defining the Transform Matrix
The goal is to find $$ \tilde{T}^{(e)}$$, a 4x4 matrix that transforms the set of local element degrees of freedom $$ d^{(e)} $$ a 4x1 matrix to another set of local, element, degrees of freedom $$ \tilde{d}^{(e)}$$ such that $$ \tilde{T}^{(e)}$$ is invertible.

$$ \tilde{d}^{(e)}$$ is the rotated degrees of freedom matrix corresponding to the $$\tilde{x}$$ coordinate system. The relationship is described below,


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

The components that make up the $$\tilde{d}^{(e)}$$ matrix are as follows,


 * $$ \tilde{d_{1}}^{(e)} =

\begin{bmatrix} l^{(e)} & m^{(e)} \end{bmatrix} \begin{bmatrix} d_{1}^{(e)}\\ d_{2}^{(e)} \end{bmatrix} $$



Here $$\tilde{d_{1}}^{(e)}$$ directly corresponds to the axial displacement $$q_{1}^{(e)}$$. And,


 * $$ \tilde{d_{2}}^{(e)} = \overrightarrow{d_{[1]}}^{(e)} \cdot \overrightarrow{\tilde{j}}$$

where $$ \tilde{d_{2}}^{(e)}$$ is a component of $$\overrightarrow{d_{[1]}}^{(e)}$$ along the unit vector $$\overrightarrow{\tilde{j}}$$, i.e., along the $$\tilde{y}$$-axis.

Here,


 * $$ \overrightarrow{i} \cdot \overrightarrow{\tilde{j}} = \cos{(\theta + \frac{\pi}{2})} = -\sin\theta$$


 * $$ \overrightarrow{j} \cdot \overrightarrow{\tilde{j}} = \cos\theta $$

This becomes,


 * $$\tilde{d_{2}}^{(e)} = -d_{1}^{(e)}\sin\theta + d_{2}^{(e)}\cos\theta$$

Where,


 * $$-\sin\theta^{(e)} = -m^{(e)}$$
 * $$ \cos\theta^{(e)} = l^{(e)}$$

So,



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

Combining these two equations equates to



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

Here,



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

Looking at the entire element e and combining the local degrees of freedom one gets,



\begin{bmatrix} \tilde{d_{1}}^{(e)}\\ \tilde{d_{2}}^{(e)}\\ \tilde{d_{3}}^{(e)}\\ \tilde{d_{4}}^{(e)} \end{bmatrix} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ -m^{(e)} & l^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)}\\ 0 & 0 & -m^{(e)} & l^{(e)} \end{bmatrix} \begin{bmatrix} d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{1}^{(e)} \end{bmatrix} $$

The first matrix is the local element degrees of freedom corresponding to the $$\tilde{x}$$ coordinate system, the second is the transformation matrix, and the third is the matrix of the local element degrees of freedom.

The transformation matrix can be broken down as follows,





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

If the transverse displacement is taken away, the force-displacement relationship will look the same as in previous lectures.


 * $$ k^{(e)}

\begin{bmatrix} 1 & -1\\     -1 & 1 \end{bmatrix} \begin{bmatrix} q_{1}^{(e)}\\ q_{2}^{(e)} \end{bmatrix} = \begin{bmatrix} p_{1}^{(e)}\\ p_{2}^{(e)} \end{bmatrix} $$

Looking at the force-displacement relationship,



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

Here the axial displacements are present and the transverse displacements are taken to be zero because if one acts along the transverse the axial displacements are zero and the spring has no effect.

This now becomes,


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

In this equation force is represented by a 4x1 matrix, the spring stiffness is represented by a 4x4 matrix, and the displacement is a 4x1 matrix.

Note: Consider the case where d4(e) &ne; 0, where d1(e)=d2(e)=d3(e)=0

Multiplying out the equation out results in a 0 matrix as listed below:

f4(e) = k4(e) d4(e) = 0 where f is a 4x1 matrix, k is a 4x4 matrix, d is a 4x1 matrix which creates a 4x1 0 matrix.

The 0 matrix represents the 4th column of the k(e) which is the interpretation of the transverse dof's.

From meeting 14 page 3: d4(e) = T4(e) d4(e)

Similarly using the same argument f4(e) = T4(e) d4(e)

Also k4(e) * d4(e) = f4(e)
 * which leads to k4(e) T4(e) d4(e) = T4(e) f4(e)

In the equation above the T4(e) is a block diagonal matrix.



\begin{bmatrix} R_{2x2}^{(e)} & 0_{2x2}\\ 0_{2x2} & R_{2x2}^{(e)}\\ \end{bmatrix} $$

If T(e) is invertible, then [T(e)-1 k(e) T(e)]d(e) = f(e).

Consider a general block diagonal matrix with:



A= \begin{bmatrix} D_1& & &  & 0\\ &  .& &  & &\\      &   &  .& & &\\      &   &   & .&&\\      0&   &  &  & D_5\\ \end{bmatrix} $$

Question: What is the inverse of A (A-1? Simple example: Diagonal Matrix



B= \begin{bmatrix} d_{11}& & &  & 0\\ &  d_{12}& &  & &\\ &  &  .& & &\\      &   &   & .&&\\      0&   &  &  & d_{mn}\\ \end{bmatrix} = diag[d_{11}, d_{22}, .... , d_{mn}] $$

Thus,



B^{-1} = Diag \begin{bmatrix} \tfrac{1}{d_{11}},& \tfrac{1}{d_{22}},& .&  .& .& .,&\tfrac{1}{d_{mn}}\\ \end{bmatrix} $$

Assuming that all values of the diagonal coefficients are not equal to 0, i.e. dii &ne; 0 for all values of i=1......n.

For a block diagonal matrix A,



A = Diag \begin{bmatrix} D_1,& .& .& .&, D_s \end{bmatrix} $$ and



A^{-1} = Diag \begin{bmatrix} D_1^{-1}, & .& .& .& .&, D_s^{-1}\\ \end{bmatrix} $$ where D1-1, is the inverse of the matrix D1.

The transform matrix is defined to be
 * T(e)-1=Diag[R(e)-1, R(e)-1]

On meeting 19 page 2 R(e) is defined as



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

R^{(e)T} = \begin{bmatrix} l^{(e)} & -m^{(e)}\\ m^{(e)} & l^{(e)} \end{bmatrix} $$ Thus multiplying R(e) by its transpose will result in



R^{(e)T}R^{(e)} = \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} 1 & 0\\     0 & l \end{bmatrix} $$ The resulting matrix represents a 2x2 I or identity matrix.

Thus,
 * R(e)-1=R(e)T
 * T(e)-1=Diag[R(e)T, R(e)T]=(Diag[R(e), R(e)])T

Also,
 * T(e)-1=T(e)T

Recall Equations:

f_{3}^{(1)} + f_{1}^{(2)} = 0 $$
 * Equation 1: $$

$$ [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}^{(1)}*d_{1}^{(2)} +            k_{12}^{(1)}*d_{2}^{(2)} + k_{13}^{(1)}*d_{3}^{(2)} + k_{14}^{(1)}*d_{4}^{(2)}] = 0 $$
 * Equation 1: Expanded

$$ [k_{31}^{(1)}*d_{1} + k_{32}^{(1)}*d_{2} + k_{33}^{(1)}*d_{3} + k_{34}^{(1)}*d_{4}] + [k_{11}^{(1)}*d_{3} + k_{12}^{(1)}*d_{4} + k_{13}^{(1)}*d_{5} + k_{14}^{(1)}*d_{6}] = 0 $$
 * Equation 1: Expanded using global degrees of freedom

f_{4}^{(1)} + f_{2}^{(2)} = P $$
 * Equation 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}^{(1)}*d_{1}^{(2)} +           k_{22}^{(1)}*d_{2}^{(2)} + k_{23}^{(1)}*d_{3}^{(2)} + k_{24}^{(1)}*d_{4}^{(2)}] = P $$
 * Equation 2: Expanded

$$ [k_{41}^{(1)}*d_{1} + k_{42}^{(1)}*d_{2} + k_{43}^{(1)}*d_{3} + k_{44}^{(1)}*d_{4}] + [k_{21}^{(1)}*d_{3} + k_{22}^{(1)}*d_{4} + k_{23}^{(1)}*d_{5} + k_{24}^{(1)}*d_{6}] = P $$
 * Equation 2: Expanded using global degrees of freedom

Assembly of   $$  \displaystyle \mathbf     k^{(e)} $$, e = 1,. . ., nel (number of elements) into the global stiffness matrix

$$\displaystyle \mathbf K = \overset{e=nel}{\underset{e=1}{\mathbb A}} \mathbf k^{(e)}$$

$$\displaystyle \overset{e=nel}{\underset{e=1}{\mathbb A}} $$ = Assembly Operator

n : total number of global degrees of freedom before eliminating the boundry conditions

ned : number of element degrees of freedom

ned << n : ned is very small when compared to n

Further Evaluation and Understanding of the Zero Eigen Values
Meeting 21 on Monday October 13th 2008 was a review and further understanding of a previously assigned homework problem dealing with the Eigen values of the Global K Matrix for a two bar truss system. The solution to the Eigen value problem actually yields four zero Eigen values which corresponds to rigid body motion. The four zero values correspond to the global dof constraints: two translations, one rotation, and one mechanism which relates the idea of zero stored energy in the truss system if you remove both supports from the two bar truss system.

The next section is a plot of the Eigen vectors which correspond to the zero Eigen values of the two bar truss system.

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

The corresponding Eigen values in a column vector are as follows:

eig(K)

ans =

-0.0000  -0.0000    0.0000    0.0000    1.4706   10.0294

The command [U,D]=eig(K) produces a matrix U whose columns are the eigenvectors of K and a diagonal matrix D with the eigenvalues of K on its diagonal.

[U,D]=eig(K)

U =

-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

From these matrices (U and D) the four Eigenvectors corresponding to the zero Eigen values can be pulled from matrix U and are as follows:

E1 =

-0.1118  -0.0814   -0.4628    0.5266   -0.4947    0.4947

E2 =

0.5043  -0.8634    0.0089   -0.0053    0.0071   -0.0071

E3 =

0        0         0         0   -0.7071   -0.7071

E4 =

-0.5931  -0.3476   -0.4803   -0.5429    0.0313   -0.0313





The results of the plots are the Mode Shapes which are a linear combination of the pure mode shapes (pure rigid body motion two translations, one rotation, and one pure mechanism). The results explained above can also be represented mathematically with the following arguments.

$$Kv = \lambda v$$

Let $$\begin{Bmatrix} U_1 & U_2 & U_3 & U_4 \end{Bmatrix}$$ be the pure eigenvetors corresponding to the four zero eigenvalues. $$K*U_i = 0*U_i$$ for i = 1,...,4.

To get W for each zero eigenvalue, we take the sum of the four contributors, each multiplied by a constant $$\alpha _{1}$$, which gives us the following relation:

$$ W = \alpha_1 u_1 + \alpha_2 u_2 + \alpha_3 u_3 + \alpha_4 u_4$$

W consists of the combination of eigenvectors for zero eigenvalues, therefore W also represents an eigenvector corresponding to a zero eigenvalue.

Eigenvalue Analysis - Continuation
The links below illustrate the application of Eigenvalue Modal Analysis to real world examples:


 * Frequency analysis of shell structures with contact
 * Frequency Solutions with Contact

The first link demonstrates various modes of the oscillations in a shell structure. The second link is an example of rigid body modes with zero eigenvalues.

In order to generate the plots of modal analysis, the Eigenvalue is simply multiplied by a sinusoidal function. For example:

$$d(t)=\lambda cos(\omega t)$$

where

$$ \omega= \frac{2\pi}{f}$$

f is the frequency of vibration.undefined

Assembly of Element Stiffness Matrix into Global Stiffness Matrix
Consider the previous example of the two bar truss and recall the element force-displacement relationship:

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

We have previously proved the equilibrium of Node 2 from the free body diagrams of Element 1 and Element 2.

$$d_{3}=d_3^{(1)}=d_1^{(2)}$$

$$d_{4}=d_4^{(1)}=d_2^{(2)}$$

For Node 2, identify the global degree of freedom to element degree of freedom for both Element 1 and 2.

Equilibrium of Node 2


$$ \sum F_x=0=-f_{3}^{(1)}-f_{1}^{(2)} $$

$$ \sum F_y=0=P-f_{4}^{(1)}-f_{2}^{(2)} $$

Next replace each f above with the individual element kd relationship.

$$ \left(k^{(1)}_{31}d^{(1)}_{1}+k^{(1)}_{32}d^{(1)}_{2}+k^{(1)}_{33}d^{(1)}_{3}+k^{(1)}_{34}d^{(1)}_{4} \right)+\left(k^{(2)}_{11}d^{(2)}_{1}+k^{(2)}_{12}d^{(2)}_{2}+k^{(2)}_{13}d^{(2)}_{3}+k^{(2)}_{14}d^{(2)}_{4} \right)=0$$

$$ \left(k^{(1)}_{41}d^{(1)}_{1}+k^{(1)}_{42}d^{(1)}_{2}+k^{(1)}_{43}d^{(1)}_{3}+k^{(1)}_{44}d^{(1)}_{4} \right)+\left(k^{(2)}_{21}d^{(2)}_{1}+k^{(2)}_{22}d^{(2)}_{2}+k^{(2)}_{23}d^{(2)}_{3}+k^{(2)}_{24}d^{(2)}_{4} \right)=P $$

Principle of Virtual Work
Motivation:


 * The principle of virtual work allows us to eliminate the corresponding rows to the boundry conditions to obtain the global stiffness matrix K. Reducing the global stiffness matrix from a (6x2) matrix to a (2x2) matrix.

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



k^{(e)} = T^{(e)T}\cdot k^{(e)}\cdot T^{(e)} $$


 * The ultimate goal is that the principal of virtual work allows us to derive the Finite Element Method for Partial Differential Equations which are heavily used in the world of Heat Transfer, Vibrations, Finite Element Analysis Commercial Codes and many more.

Force Displacement Relation for a bar or spring is expressed in the form $$ kd=F $$, which implies

$$ kd-F=0 $$ Equation (3)

$$ w(kd-F)=0 $$ for all w Equation (4)

Proof:

(3) => (4) Trivial

(4) => (3) Not Trivial

The reason why going from equation 3 to 4 is trivial is because since equation 4 is valid for all w, select w = 1 and then equation 4 becomes....

$$ 1(kd-F) = 0 $$ which is equal to equation 3

Evaluating the Five Bar Truss
The extra functions that are needed to finish the code include the functions PlaneTrussElement.m, NodalSoln.m, and PlaneTrussResults.m. They are as follows,

PlaneTrussElement.m NodalSoln.m PlaneTrussResults.m Once the extra functions have been created this code gives,

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 Connectivity Array


The connectivity array, con, uses the local node numbers of each element to locate the global node numbers. The matrix shows which global node numbers are shared by what element. For the five-bar truss the connectivity matrix has the dimensions 5x2 based on the five elements of the truss and the two nodes per each element.

The Master Location Matrix


The master location matrix, lmm, uses the local degrees of freedom numbers of each element to place the global degrees of freedom numbers of each element. The matrix illustrates which global dof numbers are shared by what element. For the five-bar truss the master location matrix has the dimensions 5x4 based on the five elements of the truss system with four degrees of freedom per each element.

The Element Stiffness Coefficient
The element stiffness coefficient for the five-bar truss is computed using the function PlaneTrussElement.m. This function brings in the element modulus of elasticity, the element area, and the element node coordinate system. From this coordinate system the function calculates the cosine directors l and m which are defined as follows,



l^{(e)} = \cos{\theta^{(e)}} = \frac{x_{2} - x_{1}}{length} $$



m^{(e)} = \sin{\theta^{(e)}} = \frac{y_{2} - y_{1}}{length} $$

The stiffness coefficient matrix is then calculate based on,



K^{(e)}= k^{(e)} \cdot \begin{bmatrix} (l^{(e)})^2 & l^{(e)}m^{(e)} & -(l^{(e)})^2 & -l^{(e)}m^{(e)}\\ l^{(e)}m^{(e)} & (m^{(e)})^2 & -l^{(e)}m^{(e)} & -(m^{(e)})^2\\ -(l^{(e)})^2 & -l^{(e)}m^{(e)} & (l^{(e)})^2 & l^{(e)}m^{(e)}\\ -l^{(e)}m^{(e)} & -(m^{(e)})^2 & l^{(e)}m^{(e)} & (m^{(e)})^2 \end{bmatrix} $$

Where,



k^{(e)} = \frac{E^{(e)}A^{(e)}}{L^{(e)}} $$

Plot of Five Bar Truss Deformation
Based on the above computation and referencing the two-bar truss deformation plot the five-bar truss can now be plotted as both un-deformed and deformed. To make the deformation noticeable a magnifying factor of 1000 was used. The code is as follows,



Contributing Writers
Eml4500.f08.ateam.boggs.t 17:02, 24 October 2008 (UTC)

Eml4500.f08.ateam.carr 17:16, 24 October 2008 (UTC)

Eml4500.f08.ateam.shah 17:36, 24 October 2008 (UTC)

Eml4500.f08.ateam.nobrega 18:10, 24 October 2008 (UTC)

Eml4500.f08.ateam.rivero 18:10, 24 October 2008 (UTC)

Eml4500.f08.ateam.mcnally 18:15, 24 October 2008 (UTC