User:Eml4500.f08.Lulz.Layton.Eric/Hw4

 Comment: This page is the updated version of Homework 4 for Team Lulz. The major changes from the on-time version submitted 20:59, 24 October 2008 were the additions of Sections 2 and 3 of the updated version and edits to Section 4. Changes in Section 4 include linking to Homework 2 for the 3 required MATLAB functions provided by the textbook, addition of the output from the ThreeBarTrussEx.m, and an updated MATLAB image of the Deformed and Undeformed Three-bar Truss.

The page comparing the on-time version and the updated version can be found here.

Eric Layton - Eml4500.f08.Lulz.Layton.Eric 08:13, 3 November 2008 (UTC)

Array Definitions
Consider the connectivity array "conn" for the two-bar truss system shown in Figure 1:

$$ \mathbf{conn} = \begin{bmatrix} 1 & 2\\ 2 & 3 \end{bmatrix} $$

The connectivity array describes the local nodes of a truss element in terms of the global node numbering convention. The rows correspond to individual truss elements, and the columns to local nodes 1 and 2 of those elements. In other words, conn(e,j) refers to the global node number of local node j for element e.

A related concept is the location master matrix lmm:

$$ \mathbf{lmm} = \begin{bmatrix} 1 & 2 & 3 & 4\\ 3 & 4 & 5 & 6 \end{bmatrix} $$

The location master matrix describes the local degrees of freedom for a truss element in terms of the global degrees of freedom. Like with the connectivity array, the rows correspond to individual truss elements, and the columns refer to local degrees of freedom 1-4 for those elements. In other words, lmm(i,j) refers to the global degree of freedom of local degree of freedom j of element i.

Transformation Matrix
A previous discussion centered on the inability to invert the transformation matrix $$\mathbf{T^{(e)}}$$ used in the axial displacement relationship $$\mathbf{q^{(e)}} = \mathbf{T^{(e)}}\cdot\mathbf{d^{(e)}}$$. Our goal is to convert the 2x4 transformation matrix into a 4x4 matrix $$\mathbf{\tilde{T}^{(e)}}$$ that can be easily inverted.



Figure 2 shows a general diagram of axial/transvere displacements $$\tilde{d}^{(e)}_{j}$$ on a truss element e. The relationship describing this case is $$\mathbf{\tilde{d}^{(e)}} = \mathbf{\tilde{T}^{(e)}}\cdot\mathbf{d^{(e)}}$$.

Through the same geometric process used to determine $$\mathbf{q^{(e)}}$$ in a previous homework, it can be shown that the two expressions below are true:

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

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

Combining these two expressions yields the following relationship for the axial/transverse displacements at one local node:

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

This describes a new displacement relationship given by $$\mathbf{\tilde{d}^{(e)}} = \mathbf{R^{(e)}}\cdot\mathbf{d^{(e)}}$$. Extending this to include both local nodes of a truss element yields:

$$\begin{Bmatrix} \tilde{d}^{(e)}_{1}\\\tilde{d}^{(e)}_{2}\\\tilde{d}^{(e)}_{3}\\\tilde{d}^{(e)}_{4}\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^{(e)}_{1}\\ d^{(e)}_{2}\\d^{(e)}_{3}\\ d^{(e)}_{4} \end{Bmatrix} = \begin{bmatrix} \mathbf{R^{(e)}}&\mathbf{0}\\\mathbf{0}&\mathbf{R^{(e)}}\end{bmatrix}\begin{Bmatrix} d^{(e)}_{1}\\ d^{(e)}_{2}\\d^{(e)}_{3}\\ d^{(e)}_{4} \end{Bmatrix}$$

Figure 3 shows the axial/transverse forces for a truss element. The axial/transverse force-displacement relationship is given by $$\mathbf{\tilde{f}^{(e)}} = \mathbf{\tilde{k}^{(e)}}\cdot\mathbf{\tilde{d^{(e)}}}$$. The modified element stiffness matrix $$\mathbf{\tilde{k}^{(e)}}$$is defined below:

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

Consider the case where $$\tilde{d}^{(e)}_{1} = \tilde{d}^{(e)}_{2} = \tilde{d}^{(e)}_{3}= 0$$, but $$\tilde{d}^{(e)}_{4}\neq 0$$. Because the 4th column of matrix $$\mathbf{\tilde{k}^{(e)}}$$ (corresponding to transverse displacement $$\tilde{d}^{(e)}_{4}$$) is a vector of zeros, the following is true:

$$\mathbf{\tilde{f}^{(e)}} = 0$$.

This is the interpretation of transverse degrees of freedom.

Finally, we have determined that $$\mathbf{\tilde{d}^{(e)}} = \mathbf{\tilde{T}^{(e)}}\cdot\mathbf{d^{(e)}}$$ is valid. A corresponding relationship can be defined for the axial/transverse forces:

$$\mathbf{\tilde{f}^{(e)}} = \mathbf{\tilde{T}^{(e)}}\cdot\mathbf{f^{(e)}}$$

We also know that $$\mathbf{\tilde{f}^{(e)}} = \mathbf{\tilde{k}^{(e)}}\cdot\mathbf{\tilde{d^{(e)}}}$$ is valid, leading to the following logic:

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

$$\mathbf{\tilde{T}^{(e)}}\cdot \mathbf{f^{(e)}} = \mathbf{\tilde{k}^{(e)}}\cdot \mathbf{\tilde{T}^{(e)}}\cdot \mathbf{d^{(e)}}$$

Assuming that $$\mathbf{\tilde{T}^{(e)}}$$ is an invertible matrix, the following is true:

$$[\mathbf{\tilde{T}^{(e)^{-1}}}\cdot\mathbf{\tilde{k}^{(e)}}\cdot\mathbf{\tilde{T}^{(e)}}]\mathbf{d^{(e)}} = \mathbf{f^{(e)}}$$

Inversion of Diagonal Matrices
Consider a general $$n\;x\;n$$ block diagonal matrix $$\mathbf{A}$$:



In order to explore the inverse, first consider an example matrix $$\mathbf{B}$$:



This matrix can alternately be described using the following notation for diagonal matrices:

$$\mathbf{B} = Diag[d_{11},d_{22},...,d_{nn}]$$

If $$d_{ii} \neq 0$$ for $$i\,=\,1,...,n$$, the inverse of this matrix can be expressed as $$\mathbf{B^{-1}} = Diag[\frac{1}{d_{11}},\frac{1}{d_{22}},...,\frac{1}{d_{nn}}]$$. Applying this procedure to $$\mathbf{A}$$ yields the following:

$$\mathbf{A} = Diag[\mathbf{D_{1}},\mathbf{D_{2}},...,\mathbf{D_{n}}]$$

$$\mathbf{A^{-1}} = Diag[\mathbf{D_{1}^{-1}},\mathbf{D_{2}^{-1}},...,\mathbf{D_{n}^{-1}}]$$

Transformation Matrix Continued
Returning to the modified transformation matrix $$\mathbf{\tilde{T}^{(e)}}$$ and applying the inversion methods discussed above, the matrix $$\mathbf{\tilde{T}^{(e)^{-1}}}$$ can be written as follows:

$$\mathbf{\tilde{T}^{(e)^{-1}}} = Diag[\mathbf{R^{(e)^{-1}}},\mathbf{R^{(e)^{-1}}}]$$, where $$\mathbf{R^{(e)}} = \begin{bmatrix} l^{(e)}&-m^{(e)}\\m^{(e)}&l^{(e)}\end{bmatrix}$$

Additionally, the transpose of matrix $$\mathbf{R^{(e)}}$$ is $$\begin{bmatrix} l^{(e)}&m^{(e)}\\-m^{(e)}&l^{(e)}\end{bmatrix}$$, leading to the following relation:

$$\mathbf{R^{(e)^{T}}}\cdot\mathbf{R^{(e)}} = \begin{bmatrix} l^{(e)}&m^{(e)}\\-m^{(e)}&l^{(e)}\end{bmatrix}\cdot \begin{bmatrix} l^{(e)}&-m^{(e)}\\m^{(e)}&l^{(e)}\end{bmatrix} = \begin{bmatrix} l^{(e)^{2}}+ m^{(e)^{2}}&0\\0&m^{(e)^{2}}+l^{(e)^{2}}\end{bmatrix}$$

Recalling that $$l^{(e)} = \cos{\mathbf{\theta}^{(e)}}$$ and $$m^{(e)} = \sin{\mathbf{\theta}^{(e)}}$$, and the trigonometric identity $$\cos{^{2}\mathbf{\theta}} + \sin{^{2}\mathbf{\theta}} = 1$$, the following is true:

$$\mathbf{R^{(e)^{T}}}\cdot\mathbf{R^{(e)}} = \begin{bmatrix} 1&0\\0&1\end{bmatrix} = \mathbf{I}$$, where $$\mathbf{I}$$ is the 2 x 2 identity matrix.

Since it is a known property that the product of a matrix and it's inverse (e.g. $$\mathbf{M}\cdot\mathbf{M^{-1}} = \mathbf{I}$$), it is apparent that $$\mathbf{R^{(e)^{T}}} = \mathbf{R^{(e)^{-1}}}$$. Thus,

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

The expression for the modified force-displacement relationship derived earlier can now be rewritten as:

$$[\mathbf{\tilde{T}^{(e)^{T}}}\cdot\mathbf{\tilde{k}^{(e)}}\cdot\mathbf{\tilde{T}^{(e)}}]\mathbf{d^{(e)}} = \mathbf{f^{(e)}}$$

Not only is the expression inside the brackets easily determined, it also resembles the form of the general force-displacement relationship, and the following can be shown:

$$\mathbf{k^{(e)}} = [\mathbf{\tilde{T}^{(e)^{T}}}\cdot\mathbf{\tilde{k}^{(e)}}\cdot\mathbf{\tilde{T}^{(e)}}]$$

Discussion of Eigenvalues and Eigenvectors
Consider the eigenvalue problem $$\mathbf{kv} = \mathbf{\lambda v}$$. Let $$[\mathbf{u_{1}},\mathbf{u_{2}},\mathbf{u_{3}},\mathbf{u_{4}}]$$ be the eigenvectors for the 4 zero eigenvalues:

$$\mathbf{ku_{i}} = 0\cdot\mathbf{u_{i}} = \mathbf{0}$$, for i = [1,2,3,4].

A linear combination of {$$\mathbf{u_{i}}$$} can be defined as $$\sum_{i=1}^{4} \alpha_{i}\mathbf{u_{i}} = \mathbf{W}$$, where $$\alpha_{i}\;$$ refers to real number constants.

The vector $$\mathbf{W}$$ is also an eigenvector corresponding a zero eigenvalue, as shown below:

$$\mathbf{kW} = \mathbf{k}(\sum_{i=1}^{4} \alpha_{i}\mathbf{u_{i}}) = \sum_{i=1}^{4} \alpha_{i}\mathbf{ku_{i}} = \mathbf{0}$$

Justification of Assembly of Element Stiffness Matrices into Global Stiffness Matrix
Recalling the force displacement global relationship of $$ \mathbf{k_{4x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} = \mathbf{f_{4x1}^{(e)}} $$, the equilibrium of global node 2 can be found. The equilibrium of global node 2 will found using the Euler Cut principle, using free body diagrams of element 1 and element 2, and identifying the global dofs to the element dofs for both element 1 and 2.



Figure 4 shows the generic two-bar truss system, with global node 2 isolated for analysis. Two equations are immediately evident from the diagram:

$$(1)\;\Sigma F_{x} = -f_{3}^{(1)} - f_{1}^{(2)} = 0$$

$$(2)\;\Sigma F_{y} = \mathbf{P} - f_{4}^{(1)} - f_{2}^{(2)} = 0$$

Equation 2 can also be written as $$\mathbf{P} = f_{4}^{(1)} + f_{2}^{(2)}$$ Using the general force-displacement relationship $$\mathbf{k^{(e)}}\cdot\mathbf{d^{(e)}}= \mathbf{f^{(e)}}$$, Equation 1 can be rewritten as follows:

$$0 = [k_{31}^{(1)}d_{1}^{(1)} + k_{32}^{(1)}d_{2}^{(1)} + k_{32}^{(1)}d_{3}^{(1)} + k_{34}^{(1)}d_{4}^{(1)}] + [k_{11}^{(2)}d_{1}^{(2)} + k_{12}^{(2)}d_{2}^{(2)} + k_{13}^{(2)}d_{3}^{(2)} + k_{14}^{(2)}d_{4}^{(2)}]$$

Replacing the local displacement degrees of freedom with the global degrees of freedom produces the following:

$$0 = [k_{31}^{(1)}d_{1} + k_{32}^{(1)}d_{2} + k_{32}^{(1)}d_{3} + k_{34}^{(1)}d_{4}] + [k_{11}^{(2)}d_{3} + k_{12}^{(2)}d_{4} + k_{13}^{(2)}d_{5} + k_{14}^{(2)}d_{6}]$$

Comparison of the above with the third row of $$\mathbf{K}$$ as described in Meeting 8 shows that the two are identical.

An expression for the assembly of $$ned\;x\;ned$$ element stiffness matrices $$\mathbf{k^{(e)}}$$ for $$e\;=\; 1,...,nel$$, where $$nel$$ refers to the number of elements, into the $$n\;x\;n$$ global stiffness matrix $$\mathbf{K}$$ can be expressed as follows:

$$ \mathbf{K} = \begin{matrix} nel\\A\\e=1\end{matrix}\mathbf{k^{(e)}}$$

$$n\;$$ is the number of global degrees of freedom before the application of boundary conditions. $$ned\;$$ is the number of element degrees of freedom. $$ned\;$$ must be less than $$n\;$$. $$A\;$$ is a symbol representing the operation of assembly.

Principle of Virtual Work (PVW)
PVW was previously used as the validation for elimination of rows in the global stiffness matrix that corresponded to zero-valued boundary conditions to produce $$\mathbf{\bar{k}}$$. It was also used to derive $$\mathbf{q^{(e)}} = \mathbf{T^{(e)}} \cdot \mathbf{d^{(e)}}$$, the axial displacement relationship, and also $$\mathbf{k^{(e)}} = \mathbf{T^{(e)^{T}}} \cdot \mathbf{\hat{k}^{(e)}} \cdot \mathbf{T^{(e)}}$$, the alternate expression for the element stiffness matrix based on the transformation matrix.

Deriving FEM for PDEs
The force-displacement relationship for a bar of spring is given by $$\mathbf{kd}=\mathbf{F}$$. This relationship implies the following:

$$(3)\;\mathbf{kd} - \mathbf{F} = 0$$

An equivalent expression to Equation 3 is also implied:

$$(4)\;\mathbf{W}(\mathbf{kd} - \mathbf{F}) = 0$$ for all $$\mathbf{W}$$.

Proof:

$$(3)\leftrightarrow(4)$$ -  trivial $$(4)\leftrightarrow(3)$$ -  NOT TRIVIAL

Since $$(4)\;$$ is valid for all $$\mathbf{W}$$ selecting $$\mathbf{W}\,=\,1$$ is an acceptable operation.

$$(4)\;$$ then becomes $$\mathbf{kd} - \mathbf{F} = 0$$, which is the same as $$(3)\;$$.

Prove $$\mathbf{\tilde{f}^{(e)}}=\mathbf{\tilde{T}^{(e)}{f}^{(e)}}$$
$$\mathbf{\tilde{f}^{(e)}}=\mathbf{\tilde{T}^{(e)}{f}^{(e)}}$$

$$\mathbf{\tilde{T}^{(e)}}=\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}$$

Since we know that $$l^{(e)} = cos(\theta)^{(e)}$$, $$m^{(e)} = sin(\theta)^{(e)}$$, and that the angle is 30 degrees, we can substitute values and compare the transformed matrix output to the values previously calculated. Confirmation will prove the relation.

$$\mathbf{\tilde{T}^{(1)}}=\begin{bmatrix} \sqrt{3}/2 & -0.5 & 0 & 0\\ 0.5 & \sqrt{3}/2 & 0 & 0\\ 0 & 0 & \sqrt{3}/2 & -0.5\\ 0 & 0 & 0.5 & \sqrt{3}/2 \end{bmatrix}$$

$$\mathbf{\tilde{T}^{(1)}f^{(1)}}=\begin{bmatrix} \sqrt{3}/2 & -0.5 & 0 & 0\\ 0.5 & \sqrt{3}/2 & 0 & 0\\ 0 & 0 & \sqrt{3}/2 & -0.5\\ 0 & 0 & 0.5 & \sqrt{3}/2 \end{bmatrix}\begin{bmatrix} -4.378\\ -2.5622\\ 4.4378\\ 2.5622 \end{bmatrix}=\mathbf{\tilde{f}^{(1)}}$$

$$\mathbf{\tilde{f}^{(1)}}=\begin{bmatrix} -5.1243\\ 0\\ 5.1243\\ 0 \end{bmatrix}$$

This confirms as expected and thus the proof is confirmed.

Prove that $$ \mathbf k^{(e)} = \tilde \mathbf T^{(e)^{T}}\tilde \mathbf k^{(e)}\tilde \mathbf T^{(e)}$$
$$ \mathbf k^{(e)} = \tilde \mathbf T^{(e)^{T}}\tilde \mathbf k^{(e)}\tilde \mathbf T^{(e)}$$

We already know the values from previous work for $$\mathbf{k^{(e)}}$$, $$\mathbf{\tilde{k}^{(e)}}$$, $$\mathbf{\tilde{T}^{(e)}}$$, and $$\mathbf{\tilde{T}^{(e)^T}}$$ is trivial to find:

$$ \mathbf{\tilde{T}^{(e)^T}} = \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}$$

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

$$\mathbf k^{(e)} = k^{(e)}\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} 1 & 0 & -1 & 0\\ 0 & 0 & 0 & 0\\ -1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0\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} (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} = \mathbf k^{(e)}$$

MATLAB Code: Five-bar and Three-bar Trusses
The first part of the required MATLAB code constructs and loads a truss system composed of five bars. The code below represents Example 1.4 (FiveBarTrussEx.m) on page 25 of the textbook. The first part defines the specifications of each bar and node including modulus of elasticity (e), cross-sectional area (a), and loading values (P). The "nodes" matrix defines the absolute locations of the four nodes. The "conn" matrix defines the global node numbers at the ends of each bar or element. Each row represents an element and nodes are numbered by the lowest first. As shown in the code, there are 5 sets of 2 numbers for the 2 nodes corresponding to each element. The "lmm" matrix represents the global degrees of freedom (dof) for each element. Each node has 2 dof's; therefore, each element has 4 corresponding dof's. Degrees of freedom are numbered from 1 to 2*(number of nodes) where the first 2 correspond to Global Node 1, the second 2 correspond to the Global Node 2, etc. From all this data, the global stiffness matrix (K) is created. The displacement of the loaded node is calculated using the relevant parts of the global stiffness matrix, inverting them, and multiplying by the load array (R). Finally, the reactions (reactions) at each node are calculated by multiplying the displacement array (d) by the global stiffness matrix. The functions "PlaneTrussElement.m", "NodalSoln.m", and "PlaneTrussResults.m" are used again from Section 6 of Homework 2.

FiveBarTrussEx.m
% Five bar truss example e1=200*10^3; e2=70*10^3; a1=40*100; a2=30*100; a3=20*100; P = 150*10^3; nodes = 1000*[0, 0; 1.5, 3.5; 0, 5; 5, 5]; conn = [1,2; 2,4; 1,3; 3,4; 2,3]; lmm = [1,2,3,4; 3, 4, 7, 8; 1, 2, 5, 6; 5, 6, 7, 8; 3, 4, 5, 6]; K=zeros(8); % Generate stiffness matrix for each element and assemble it. for i=1:2 lm=lmm(i,:); con=conn(i,:); k=PlaneTrussElement(e1, a1, nodes(con,:)); K(lm, lm) = K(lm, lm) + k; end for i=3:4 lm=lmm(i,:); con=conn(i,:); k=PlaneTrussElement(e1, a2, nodes(con,:)); K(lm, lm) = K(lm, lm) + k; end

lm=lmm(5,:); con=conn(5,:); k=PlaneTrussElement(e2, a3, nodes(con,:)); K(lm, lm) = K(lm, lm) + k

% Define the load vector R = zeros(8,1); R(4)=-P

% Nodal solution and reactions [d, reactions] = NodalSoln(K, R, [1,2,7,8], zeros(4,1)) results=[]; for i=1:2 results = [results; PlaneTrussResults(e1, a1, ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end for i=3:4 results = [results; PlaneTrussResults(e1, a2, ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g results = [results; PlaneTrussResults(e2, a3, ...       nodes(conn(5,:),:), d(lmm(5,:)))]

Output to the Command Window

>> FiveBarTrussEx

K =

32600       76067       -32600       -76067            0            0            0            0        76067  2.9749e+005       -76067 -1.7749e+005            0    -1.2e+005            0            0 -32600      -76067  2.4309e+005  1.1914e+005       -32998        32998 -1.7749e+005       -76067 -76067 -1.7749e+005 1.1914e+005  2.4309e+005        32998       -32998       -76067       -32600 0           0       -32998        32998    1.53e+005       -32998    -1.2e+005            0 0   -1.2e+005        32998       -32998       -32998    1.53e+005            0            0 0           0 -1.7749e+005       -76067    -1.2e+005            0  2.9749e+005        76067 0           0       -76067       -32600            0            0        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 results show the strain, stress, and axial force in each element corresponding to the three columns in the "results" matrix.

To visualize the loaded and unloaded system, the following code graphs the system in its loaded and unloaded state. In order to be able to realize the direction of node movement, the unfixed node displacements have been magnified by a factor of two (using the variable "m"). As shown, the original values from the displacement array (d3-d6) above were used below to get the right proportions of the nodal displacements.

FiveBarTrussDisplay.m
%Displaying the five bar truss clear; close;

m = 2;

d1 = 0;         %disp of node 1 in the x direction d2 = 0;         %disp of node 1 in the y direction d3 = 0.53895*m; %disp of node 2 in the x direction d4 = -0.95306*m; %disp of node 2 in the y direction d5 = 0.2647*m;  %disp of node 3 in the x direction d6 = -0.2647*m; %disp of node 3 in the y direction d7 = 0;         %disp of node 4 in the x direction d8 = 0;         %disp of node 4 in the y direction

%input nodal coordinates num_nodes = 4;              %total number of nodes num_elems = 5;              %total number of elements tot_dofs = 2 * num_nodes;   %total dofs of the system

%undeformed positions pos(1,:)=[0 1.5 0 5]; pos(2,:)=[0 3.5 5 5];

%undeformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

%element connectivity array %follows: %node_conn(local node number, element number) = global node number

node_conn(1,1)=1;    %element 1 node_conn(2,1)=2;

node_conn(1,2)=2;    %element 2 node_conn(2,2)=4;

node_conn(1,3)=1;    %element 3 node_conn(2,3)=3;

node_conn(1,4)=3;    %element 4 node_conn(2,4)=4;

node_conn(1,5)=2;    %element 5 node_conn(2,5)=3;

%plot the undeformed 2-bar truss system for i=1:num_elems node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; axis([-1 6 -1 6]) plot(xx,yy,'--rs') hold on end

n = -0.5; o = -0.25; text(0+n,0,'GN1','HorizontalAlignment','center') text(1.5+n,3.5,'GN2','HorizontalAlignment','center') text(0+n,5,'GN3','HorizontalAlignment','center') text(5+n,5,'GN4','HorizontalAlignment','center') text(0.75+o,1.75,'E1','HorizontalAlignment','center') text(3.25+o,4.25,'E2','HorizontalAlignment','center') text(0+o,2.5,'E3','HorizontalAlignment','center') text(2.5+o,5,'E4','HorizontalAlignment','center') text(0.75+o,4.25,'E5','HorizontalAlignment','center')

%deformed positions pos(1,:)=[0+d1 1.5+d3 0+d5 5+d7]; pos(2,:)=[0+d2 3.5+d4 5+d6 5+d8];

%deformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

text(1.5+d3+n,3.5+d4,'GN2','HorizontalAlignment','center') text(0+d5+n,5+d6,'GN3','HorizontalAlignment','center') text(pos(1,1)+(pos(1,2)-pos(1,1))/2+o,pos(2,1)+(pos(2,2)-pos(2,1))/2,'E1','HorizontalAlignment','center') text(pos(1,2)+(pos(1,4)-pos(1,2))/2+o,pos(2,2)+(pos(2,4)-pos(2,2))/2,'E2','HorizontalAlignment','center') text(pos(1,1)+(pos(1,3)-pos(1,1))/2+o,pos(2,1)+(pos(2,3)-pos(2,1))/2,'E3','HorizontalAlignment','center') text(pos(1,3)+(pos(1,4)-pos(1,3))/2+o,pos(2,3)+(pos(2,4)-pos(2,3))/2,'E4','HorizontalAlignment','center') text(pos(1,2)+(pos(1,3)-pos(1,2))/2+o,pos(2,2)+(pos(2,3)-pos(2,2))/2,'E5','HorizontalAlignment','center')

%plot the deformed 2-bar truss system for i=1:num_elems hold on   node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; plot(xx,yy,'-s') end

title('Deformed and Undeformed Five-bar Truss System') xlabel('X') ylabel('Y')



The dotted line represents the undeformed system shape. The solid line represents the deformed system shape. The abbreviation "GN" and "E" stand for Global Node and Element, respectively. Note only nodes 2 and 3 are displaced because they are unfixed.

For the three bar truss created in the class notes, the two bar truss system code (Two Bar Truss System.m) was modified from Section 6 of Homework 2. One node and element each were added to the system, and e, A, P, and L values were changed to coincide with the problem.

ThreeBarTrussEx.m
% Three bar truss example clear all; e = [2 4 3]; A = [3 1 2]; P = 30; L = [5 5 10]; alpha = pi/6; beta = -pi/6; gamma = -3*pi/4; nodes = [0, 0; L(1)*cos(alpha), L(1)*sin(alpha); L(1)*cos(alpha) + L(2)*cos(beta), L(1)*sin(alpha) + L(2)*sin(beta); L(1)*cos(alpha) + L(3)*cos(gamma), L(1)*sin(alpha) + L(3)*sin(gamma)]; dof = 2*length(nodes); conn = [1 2; 2 3; 2 4]; lmm = [1 2 3 4; 3 4 5 6; 3 4 7 8]; elems = size(lmm,1); K = zeros(dof); R = zeros(dof,1); debc = [1 2 5 6 7 8]; ebcVals = zeros(length(debc),1);

%load vector R = zeros(dof,1); R(4) = P;

% Assemble global stiffness matrix K = zeros(dof); for i = 1:elems lm = lmm(i,:); con = conn(i,:); k_local = e(i)*A(i)/L(i)*[1 -1; -1 1] k = PlaneTrussElement(e(i), A(i), nodes(con,:)) K(lm, lm) = K(lm, lm) + k; end K R % Nodal solution and reactions [d, reactions] = NodalSoln(K, R, debc, ebcVals) results = []; for i = 1:elems results = [results; PlaneTrussResults(e, A, nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g results

Output to the Command Window

>> ThreeBarTrussEx

k_local =

1.2000  -1.2000   -1.2000    1.2000

k =

0.9000   0.5196   -0.9000   -0.5196    0.5196    0.3000   -0.5196   -0.3000   -0.9000   -0.5196    0.9000    0.5196   -0.5196   -0.3000    0.5196    0.3000

k_local =

0.8000  -0.8000   -0.8000    0.8000

k =

0.6000  -0.3464   -0.6000    0.3464   -0.3464    0.2000    0.3464   -0.2000   -0.6000    0.3464    0.6000   -0.3464    0.3464   -0.2000   -0.3464    0.2000

k_local =

0.6000  -0.6000   -0.6000    0.6000

k =

0.3000   0.3000   -0.3000   -0.3000    0.3000    0.3000   -0.3000   -0.3000   -0.3000   -0.3000    0.3000    0.3000   -0.3000   -0.3000    0.3000    0.3000

K =

0.9000   0.5196   -0.9000   -0.5196         0         0         0         0    0.5196    0.3000   -0.5196   -0.3000         0         0         0         0   -0.9000   -0.5196    1.8000    0.4732   -0.6000    0.3464   -0.3000   -0.3000   -0.5196   -0.3000    0.4732    0.8000    0.3464   -0.2000   -0.3000   -0.3000         0         0   -0.6000    0.3464    0.6000   -0.3464         0         0         0         0    0.3464   -0.2000   -0.3464    0.2000         0         0         0         0   -0.3000   -0.3000         0         0    0.3000    0.3000         0         0   -0.3000   -0.3000         0         0    0.3000    0.3000

R =

0    0     0    30     0     0     0     0

d =

0        0  -11.6737   44.4051         0         0         0         0

reactions =

-12.5672  -7.2557   22.3866  -12.9249   -9.8194   -9.8194

results =

2.4186      4.8371       9.6742       7.2557       14.511       9.6742       14.511       6.4625       12.925        25.85       19.387       38.775        25.85       38.775       2.3145       4.6289       9.2578       6.9434       13.887       9.2578       13.887

After running this code, the output produces displacement of node 2 to be -11.6737 in the x-direction and 44.4051 in the y-direction. On the other hand, the "results" matrix has 7 columns instead of 4. This is due to the fact that vectors for the modulus of elasticity (e), cross-sectional area (A), and element length (L) were passed as variables instead of the scalar values. The first column is the strain array for elements 1-3. This is followed by two 3x3 matrices created for stresses and axial forces of the elements instead of two arrays of length 3 each. The real values for stress and axial force lie along the diagonals of each of these matrices.

The code for displaying the three bar truss was similar to the code displaying the five bar truss. However, a displacement magnification factor (m) of 0.25 was used to ease visualization of the deformation.

ThreeBarTrussDisplay.m
%Displaying the three bar truss clear; close;

m = 0.25;

d1 = 0;         %disp of node 1 in the x direction d2 = 0;         %disp of node 1 in the y direction d3 = -11.6737*m; %disp of node 2 in the x direction d4 = 44.4051*m; %disp of node 2 in the y direction d5 = 0;         %disp of node 3 in the x direction d6 = 0;         %disp of node 3 in the y direction d7 = 0;         %disp of node 4 in the x direction d8 = 0;         %disp of node 4 in the y direction

%input nodal coordinates num_nodes = 4;              %total number of nodes num_elems = 3;              %total number of elements tot_dofs = 2 * num_nodes;   %total dofs of the system

%element 1 details L1 = 5; E1 = 2; A1 = 3; theta1 = 30; %deg thetar1 = (pi * theta1)/180;

%element 2 details L2 = 5; E2 = 4; A2 = 1; theta2 = -30; %deg thetar2 = (pi * theta2)/180;

%element 3 details L3 = 10; E3 = 3; A3 = 2; theta3 = -135; %deg thetar3 = (pi * theta3)/180;

%undeformed positions pos(1,:)=[0,L1*cos(thetar1),L1*cos(thetar1)+L2*cos(thetar2),L1*cos(thetar1)+L3*cos(thetar3)]; pos(2,:)=[0,L1*sin(thetar1),L1*sin(thetar1)+L2*sin(thetar2),L1*sin(thetar1)+L3*sin(thetar3)];

%undeformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

%element connectivity array %follows: %node_conn(local node number, element number) = global node number

node_conn(1,1)=1;    %element 1 node_conn(2,1)=2;

node_conn(1,2)=2;    %element 2 node_conn(2,2)=3;

node_conn(1,3)=2;    %element 3 node_conn(2,3)=4;

%plot the undeformed 2-bar truss system for i=1:num_elems node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; axis([-5 10 -10 20]) plot(xx,yy,'--rs') hold on end

n = -1; o = -0.5; text(x(node_conn(1,1))+n,y(node_conn(1,1)),'GN1','HorizontalAlignment','center') text(x(node_conn(2,1))+n,y(node_conn(2,1)),'GN2','HorizontalAlignment','center') text(x(node_conn(2,2))+n,y(node_conn(2,2)),'GN3','HorizontalAlignment','center') text(x(node_conn(2,3))+n,y(node_conn(2,3)),'GN4','HorizontalAlignment','center') text(x(node_conn(2,1))/2+o,y(node_conn(2,1))/2,'E1','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,2))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,2))-y(node_conn(2,1)))/2,... 'E2','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,3))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,3))-y(node_conn(2,1)))/2,... 'E3','HorizontalAlignment','center')

%deformed positions pos(1,:)=[0+d1,L1*cos(thetar1)+d3,L1*cos(thetar1)+L2*cos(thetar2)+d5,L1*cos(thetar1)+L3*cos(thetar3)+d7]; pos(2,:)=[0+d2,L1*sin(thetar1)+d4,L1*sin(thetar1)+L2*sin(thetar2)+d6,L1*sin(thetar1)+L3*sin(thetar3)+d8];

%deformed nodal coordinate arrays for i=1:num_nodes x(i)=pos(1,i); y(i)=pos(2,i); end

text(x(node_conn(2,1))+n,y(node_conn(2,1)),'GN2','HorizontalAlignment','center') text(x(node_conn(2,1))/2+o,y(node_conn(2,1))/2,'E1','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,2))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,2))-y(node_conn(2,1)))/2,... 'E2','HorizontalAlignment','center') text(x(node_conn(2,1))+(x(node_conn(2,3))-x(node_conn(2,1)))/2+o,y(node_conn(2,1))+(y(node_conn(2,3))-y(node_conn(2,1)))/2,... 'E3','HorizontalAlignment','center')

%plot the deformed 2-bar truss system for i=1:num_elems hold on   node_1=node_conn(1,i); node_2=node_conn(2,i); xx=[x(node_1),x(node_2)]; yy=[y(node_1),y(node_2)]; plot(xx,yy,'-s') end

title('Deformed and Undeformed Three-bar Truss System') xlabel('X') ylabel('Y')



The dotted line represents the undeformed system shape. The solid line represents the deformed system shape. The abbreviation "GN" and "E" stand for Global Node and Element, respectively. Note only node 2 is displaced because it is unfixed.

Derivation of d2
Knowing that $$l^e = cos(\theta)$$ and $$m^e = sin(\theta)$$

$$d_2^e = d_1^e\vec i + d_2^e\vec j$$

$$\tilde d_2^e = \tilde d_1^e\vec\tilde j$$

$$\tilde d_2^e = d_1^e(\vec i*\vec\tilde j) + d_2^e(\vec j*\vec\tilde j)$$

Knowing that $$(\vec i*\vec\tilde j)= -m^e$$ and $$(\vec j*\vec\tilde j)= l^e$$ the following simplification can be preformed.

$$\tilde d_2^e = -m^ed_1^e + l^ed_2^e$$

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

Eigenvector Plotting
Now we want to plot the eigenvectors corresponding to the 0 eigenvalues of the 2-bar truss system that we have been looking at in class for the past few weeks. We want to assume no constraints on the 2-bar truss system so we will use the full 6x6 global stiffness matrix.

We know that the global stiffness matrix K is:

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

If we use the eig function in Matlab we get that the eigenvalues are as follows:

Eval=eig(K) Eval= -0.0001        -0.0000         -0.0000          0.0000          1.4706         10.0294

Next we can use the function [V,D] = EIG(X)which gives us a diagonal matrix D of eigenvalues and a full matrix V whose columns contain the corresponding eigenvectors that corresponds to the equation:

$$\textbf{X}\times\textbf{V}=\textbf{D}\times\textbf{V}$$

For our particular case, we set X equivalent to K and the diagonal matrix D equivalent to $$\lambda$$ to get the following:

$$\textbf{K}\times\textbf{V}=\lambda\times\textbf{V}$$

Therefore, using Matlab we get:

[V,Lam]=eig(K) Lam = -0.0001        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  V = -0.4494   0.3877    0.5158    0.0160    0.6174   -0.0139    0.6754   -0.4281    0.4825    0.0244    0.3565   -0.0080    0.1682    0.3877    0.5158    0.0160   -0.5409    0.5123   -0.3942   -0.4281    0.4825    0.0244   -0.4330   -0.4904    0.2812    0.4163   -0.0116    0.7023   -0.0765   -0.4984   -0.2812   -0.3995   -0.0449    0.7107    0.0765    0.4984

Where Lam contains the previously found eigenvalues and V contains the corresponding eigenvectors. It is important to note now that while Matlab calculated the eigenvalues to precise decimals, we can assume that the first four values in the diagonal Lam matrix are 0.

Now we would like to plot the eigenvectors that correspond to these four 0 eigenvalues so we may get a look at their mode shapes. Using the following Matlab function allows us to plot the four mode shapes at the same time.

for i=1:4 figure plot(V(:,i)); title(['Mode Shape ',num2str(i)]) end



Now let us consider the following three-bar truss system.



Assume that all three bars have the same material properties:

L=1    E=2      A=3

$$k=\frac{EA}L=\frac{(2)(3)}1=6$$

For element 1 and element 3:

$$l^{(1,3)}=0$$

$$m^{(1,3)}=1$$

For elemnt 2:

$$l^{(2)}=1$$

$$m^{(2)}=0$$

$$\textbf{k}^{(1)}=\textbf{k}^{(3)}= \left [\begin{array}{c c c c} 0 & 0 & 0 & 0\\ 0 & 36 & 0 & -36\\ 0 & 0 & 0 & 0\\ 0 & -36 & 0 &-36 \end{array}\right]$$

$$\textbf{k}^{(2)}= \left [\begin{array}{c c c c} 36 & 0 & -36 & 0\\ 0 & 0 & 0 & 0\\ -36 & 0 & 36 & 0\\ 0 & 0 & 0 & 0 \end{array} \right]$$

$$\textbf{K}= \left [\begin{array}{c c c c c c c c} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 36 & 0 & -36 & 0 & 0 & 0 & 0\\ 0 & 0 & 36 & 0 & -36 & 0 & 0 & 0\\ 0 & -36 & 0 & 36 & 0 & 0 & 0 & 0\\ 0 & 0 & -36 & 0 & 36 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 36 & 0 & -36\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & -36 & 0 & 36 \end{array} \right] $$

In this example we are assuming that global nodes 1 and 4 are constrained so that they don't move. This means that $$d_1=d_2=d_7=d_8=0$$. Therefore we only have to consider the following modified global stiffness matrix.

$$\overline{\textbf{K}}= \left[ \begin{array} {c c c c} 36 & 0 & -36 & 0\\ 0 & 36 & 0 & 0\\ -36 & 0 & 36 & 0\\ 0 & 0 & 0 & 36 \end{array} \right]$$

If we follow the same steps as above we get:

Kbar = 36    0   -36     0     0    36     0     0   -36     0    36     0     0     0     0    36 Eval=eig(Kbar) Eval = 0   36    36    72 [V,Lam]=eig(Kbar) V = -0.7071        0         0   -0.7071         0    1.0000         0         0   -0.7071         0         0    0.7071         0         0    1.0000         0 Lam = 0    0     0     0     0    36     0     0     0     0    36     0     0     0     0    72

In this example we have only one zero eigenvalue and it's eigenvectors can be plotted as follows:

plot(V(:,1)); title(['Mode Shape',num2str(1)])



Derivation and Assembly of Row 4
Equation 2 is:

$$f_4^1+f_2^2=0$$

Which can also be written as:

$$\begin{bmatrix} k_{41}^1d_1^1+k_{42}^1d_2^1+k_{43}^1d_3^1+k_{44}^1d_4^1 \end{bmatrix}+\begin{bmatrix} k_{21}^2d_1^2+k_{22}^2d_2^2+k_{23}^2d_3^2+k_{24}^2d_4^2 \end{bmatrix}=0$$

The previous equation rewritten using global dof rather then local dof is:

$$\begin{bmatrix} k_{41}^1d_1^1+k_{42}^1d_2^1+k_{43}^1d_3^1+k_{44}^1d_4^1 \end{bmatrix}+\begin{bmatrix} k_{21}^2d_3^2+k_{22}^2d_4^2+k_{23}^2d_5^2+k_{24}^2d_6^2 \end{bmatrix}=0$$

Contributing Team Members
Aaron Fisher - Eml4500.f08.Lulz.fisher 00:37, 24 October 2008 (UTC)

Eric Layton - Eml4500.f08.Lulz.Layton.Eric 02:19, 22 October 2008 (UTC)

John Saxon - Eml4500.f08.Lulz.js 05:23, 24 October 2008 (UTC)

Andrew Strack - Eml4500.f08.lulz.strack 10:25, 24 October 2008 (UTC)

Benjamin Mitchell - Eml4500.f08.Lulz.mitchell.bm 20:09, 24 October 2008 (UTC)

Sam Miorelli - Eml4500.f08.lulz.abcd 21:01, 24 October 2008 (UTC)