User:Eml4500.f08.FEABBQ/HW3/Notes

Class Notes: September 22, 2008
Derivation of element FD cont. global coordinate system

Equation [1]

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



Equation [2]

$$ \mathbf{k^{(e)}} \mathbf{q^{(e)}} = \mathbf{P^{(e)}} $$

$$ k^{(e)} \left[ \begin{matrix} 1 & -1 \\ -1 & 1 \end{matrix} \right] \begin{Bmatrix} q_{1}^{(e)} \\ q_{2}^{(e)} \end{Bmatrix} = \begin{Bmatrix} P_{1}^{(e)} \\ P_{2}^{(e)} \end{Bmatrix}$$

$$q_i^{(e)} = $$ Axial displacement of element e at the the local node i

$$P_i^{(e)} = $$ Axial force of element e at the local node i

Goal : Derive [1] from [2]
The relationship can be expressed in form of,

$$ \mathbf{q_{2x1}^{(e)}} = \mathbf{T_{2x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} $$

Now consider the displacement vector of the local node 1, which is denoted by $$ \mathbf{d_1^{(e)}} $$

$$ \mathbf{d_1^{(e)}} = d_i^{(e)} \vec i + d_2^{(e)} \vec j $$

Now $$ q_1^{(e)} $$ equals the axial displacement of node i which is the orthogonal projection of displacement vector of $$ \mathbf{d_1^{(e)}} $$ of node 1 on the axis $$ \tilde{x} $$ of element (e)

Node 1:

$$ \Rightarrow q_1^{(e)} = \vec d_1^{(e)} \cdot \vec \tilde{i} $$

$$ q_1^{(e)} = (d_1^{(e)} \tilde{i} + d_2^{(e)} \tilde{j}) \cdot \vec \tilde{i} $$

$$ q_1^{(e)} = d_1^{(e)} (\vec i \cdot \vec \tilde{i}) + d_2^{(e)} (\vec j \cdot \vec \tilde{i}) $$

REMEMBER: from lecture Sept. 8th that

$$ l^{(e)} = \mathbf{\tilde{i}} \cdot \mathbf{i} = cos(\theta^{(e)}) $$

$$ m^{(e)} = \mathbf{\tilde{i}} \cdot \mathbf{j} = sin(\theta^{(e)}) $$

Therefore,

$$ q_1^{(e)} = l^{(e)} d_1^{(e)} + m^{(e)} d_2^{(e)} $$

$$ q_1^{(e)} = \left[ \begin{matrix} l^{(e)} & m^{(e)} \end{matrix} \right] \left[ \begin{matrix} d_1^{(e)} \\ d_2^{(e)} \end{matrix} \right] $$

Node 2:

$$ \Rightarrow q_2^{(e)} = \vec d_3^{(e)} \cdot \vec \tilde{j} $$

$$ q_2^{(e)} = (d_3^{(e)} \tilde{i} + d_4^{(e)} \tilde{j}) \cdot \vec \tilde{j} $$

$$ q_2^{(e)} = d_3^{(e)} (\vec i \cdot \vec \tilde{j}) + d_4^{(e)} (\vec j \cdot \vec \tilde{j}) $$

$$ q_2^{(e)} = m^{(e)} d_3^{(e)} +  l^{(e)} d_4^{(e)} $$

$$ q_2^{(e)} = \left[ \begin{matrix} l^{(e)} & m^{(e)} \end{matrix} \right] \left[ \begin{matrix} d_3^{(e)} \\ d_4^{(e)} \end{matrix} \right] $$

Now you can combine $$ q_1^{(e)} $$ and $$ q_2^{(e)} $$ into one matrix,

$$ \begin{Bmatrix} q_1^{(e)} \\ q_2^{(e)} \end{Bmatrix} = \left[ \begin{matrix} l^{(e)} & m^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \end{matrix} \right] \begin{Bmatrix} d_1^{(e)} \\ d_2^{(e)} \\ d_3^{(e)} \\ d_4^{(e)} \end{Bmatrix} $$

Class Notes: September 26, 2008
Equation [2]

$$ \mathbf{k^{(e)}} \mathbf{q^{(e)}} = \mathbf{P^{(e)}} $$

Fallowing the same argument above we get the relation

$$ \begin{Bmatrix} p_1^{(e)} \\ p_2^{(e)} \end{Bmatrix} = \left[ \begin{matrix} l^{(e)} & m^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \end{matrix} \right] \begin{Bmatrix} f_1^{(e)} \\ f_2^{(e)} \\ f_3^{(e)} \\ f_4^{(e)} \end{Bmatrix} $$

Equation [3]

$$ \mathbf{p_{2x1}^{(e)}} = \mathbf{T_{2x4}^{(e)}} \mathbf{f_{4x1}^{(e)}} $$

recalling element axial FD relation

Equation [4]

$$ \mathbf {k_{2x4}^{(e)}} \mathbf{q_{4x1}^{(e)}} =\mathbf{p_{2x1}^{(e)}}  $$

and Equation [5]

$$ \mathbf{q_{2x1}^{(e)}} = \mathbf{T_{2x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} $$

Substituting equations 3 and 5 into equation 4 yields to

$$ \mathbf {k_{2x4}^{(e)}} \mathbf{T_{2x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} =\mathbf{T_{2x4}^{(e)}} \mathbf{f_{4x1}^{(e)}}  $$

The Goal is to have $$ \mathbf{k^{(e)}} \mathbf{d^{(e)}} = \mathbf{f^{(e)}} $$ so move $$ \mathbf{T^{(e)}} $$ from R.H.S. to L.H.S. by premultiplying the equation by $$ \mathbf{{(T^{(e)})}^{(-1)}} $$        ( inverse $$ \mathbf{T^{(e)}} $$)

But unfortunately $$ \mathbf{T^{(e)}} $$ is a rectangular matrix [2x4], can't inverse $$ \mathbf{T^{(e)}} $$

to solve this problem the transpose property was used as follows

$$ \mathbf{{T_{4x2}^{(e)}}^{T}} \mathbf {k_{2x4}^{(e)}} \mathbf{T_{2x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} = \mathbf{f_{4x1}^{(e)}} $$

$$ \mathbf{{T_{4x2}^{(e)}}^{T}} \mathbf {k_{2x4}^{(e)}} \mathbf{T_{2x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} = \mathbf{f_{4x1}^{(e)}} $$

$$ \mathbf {k_{2x4}^{(e)}} \mathbf{d_{4x1}^{(e)}} = \mathbf{f_{4x1}^{(e)}}  $$

Homework

Verify that

$$ \mathbf {I_{4x4}} = \mathbf{{T_{4x2}^{(e)}}^{T}} \mathbf{T_{2x4}^{(e)}}   $$

$$ \mathbf {I_{4x4}} =  \left[ \begin{matrix} l^{(e)} & 0 \\ m^{(e)} & 0 \\ 0 & l^{(e)} \\ 0 & m^{(e)}  \end{matrix} \right]  \left[ \begin{matrix} l^{(e)} & m^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \end{matrix} \right]   $$

Applying the principle of virtual work will result reduction of global FD relation:

$$ \mathbf {k_{2x4}^{(e)}} = \mathbf{{T_{4x2}^{(e)}}^{T}}  \mathbf{\tilde{k}_{4x1}^{(e)}} \mathbf{T_{2x4}^{(e)}}   $$

Remark

why not solve as follows $$ \mathbf{d_{4x1}^{(e)}} = \mathbf{{(T^{(e)})}^{(-1)}}  \mathbf{f_{4x1}^{(e)}}  $$  ?

The answer is that the used methods to compute $$ \mathbf{{(K^{(e)})}^{(-1)}} $$ could not be used since K is singular. this will result a determinants value of zero and thus K is not inverted.

Recall

for an unconstrained structure system there are three possible body motion in 2-D (plane) two translation and one rotation.

Homework

Compute the eigenvalues of K and make observation about the number of zero eigenvalues.

Stiffness matrix K is defined as

$$\textbf{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}\,$$

with values from sample problem in homework-2 K becomes

$$\textbf{K}= \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}\,$$

Solving for the eigenvalues of K using MATLAB results in

Matlab: Computing eigenvalues
The following is a matlab script and output for the two-bar K example.

Zero Eigenvalues correspond to zero stored elastic energy, rigid body modes- (shape, eigenvictors)

Class Notes: September 29, 2008
Using the Global Force-Displacement Relation

Use the global force-displacement relationship matrix from lecture 10:

$$ \left[ \begin{array}{cc} k_{13}^{(1)} & k_{14}^{(1)} \\ k_{23}^{(1)} & k_{24}^{(1)} \\ (k_{33}^{(1)} + k_{11}^{(2)}) & (k_{34}^{(1)} + k_{12}^{(2)}) \\  (k_{43}^{(1)} + k_{21}^{(2)}) & (k_{44}^{(1)} + k_{22}^{(2)})  \\  k_{31}^{(2)} & k_{32}^{(2)} \\ k_{41}^{(2)} & k_{42}^{(2)}  \\ \end{array} \right] $$ $$\left[ \begin{array}{l} d_1 \\ d_2  \\ d_3  \\ d_4 \\ d_5 \\ d_6 \end{array} \right] = \left[ \begin{array}{l} F_1 \\ F_2  \\ F_3  \\ F_4 \\ F_5 \\ F_6 \end{array} \right]$$ (6x2) x (2x1)  =  (6x1)

Only do computations for rows 1,2,5,6 to get F1,F2,F5,F6. ....

Closing the loop between the Finite Element Method and statics: Virtual Displacement

For the two bar truss system:

Going from computing reactions to computing the displacements in statics "closes the loop".

by statics, reactions are known and therefor the member forces $$P_1^{(1)}$$ and $$P_2^{(1)}$$

Compute the axial displacement of the degrees of freedom (extension of bars):

$$q_2^{(1)} = \frac{P_1^{(1)}}{k^{(1)}} = \frac{P_2^{(1)}}{k^{(1)}}$$

$$ q_1^{(1)} = 0$$, fixed global node 1

$$q_1^{(2)} = \frac{-P_2^{(2)}}{k^{(2)}} = -AB$$

$$q_2^{(2)} = 0 $$, fixed global node 3

Question: How to get the displacement DOF's from the above results for node 2?

Class Notes: Wednesday, October 1, 2008


In order to find the coordinates of point D, one must first know the coordinates of points B and C. The magnitudes of AC and AB are solved for below.


 * $$AC=\frac{\vert P_2^{(1)}\vert}{k^{(1)}}=\frac{5.1243}{\frac{3}{4}}=6.8324$$


 * $$AB=\frac{\vert P_1^{(2)}\vert}{k^{(2)}}=\frac{6.2760}{5}=1.2552$$

Now that the distances of points B and C are known from the reference point, point A, we must solve for the coordinates of points B and C by constructing lines parallel to initial elements and measuring that distance along the line to find $$(x_B,y_B)$$ and $$(x_C,y_c)$$.

To do this we need the equations for the line AB and line AC. The figure below shows a line being constructed that runs through P and some arbitrary point Q at $$(x,y)$$. This line has the same slope as $$\vec \tilde{i}$$.



Constructing this line is as follows


 * $$\overrightarrow{P Q} = \vert PQ \vert \vec \tilde{i}=(PQ)[\cos\theta \vec i + \sin\theta \vec j]$$


 * $$=(x-x_p)\vec i + (y-y_p)\vec j$$

From this one can see that
 * $$(x-x_p)=(PQ)\cos\theta$$


 * $$(y-y_p)=(PQ)\sin\theta$$

This is sufficient to solve for the coordinates of points B and C, however to make one further useful extension, these two expressions can be combined to create the equation of a line.


 * $$\frac{(y-y_p)}{(x-x_p)}=\frac{\sin\theta}{\cos\theta}=\tan\theta \qquad \Rightarrow \qquad y-y_p = (\tan\theta)(x-x_p)$$

Using these expressions the coordinates for B and C are now to be solved for. Chosing A as the origin will make calculations convenient, so for these instances that is what was chosen.


 * $$(x_B,y_B)\Rightarrow(x_B - x_A, y_B - y_A)=((AB)\cos\theta^{(2)},(AB)\sin\theta^{(2)})$$
 * $$\quad =((1.2552)\cos 135^\circ,(1.2552)\sin 135^\circ)$$
 * $$(x_B,y_B) = (-0.8876,0.8876)$$


 * $$(x_C,y_C)\Rightarrow(x_C - x_A, y_C - y_A)=((AC)\cos\theta^{(1)},(AC)\sin\theta^{(1)})$$
 * $$\quad =((6.8324)\cos 30^\circ,(6.8324)\sin 30^\circ)$$
 * $$(x_C,y_C) = (5.9170,3.4162)$$

Once the equation of a line is known in the form shown earlier, one can easily derive the equation of a line perpendicular to the original, passing through a given point by adding a multiple of $$\frac {\pi}{2}$$ to the argument of the tangent term. The equation of the line perpendicular to the equation of the line derived above, passing through point P is given below.


 * $$ y - y_p = (\tan{\theta + \frac{\pi}{2}})(x-x_p)$$

The last unknown position is that of point D. This can either be calculated via the intersection of perpendicular lines passing through B and C, or via the results of the FEM solution. Choosing point A to be the origin once again for the same convenience, the coordinates of point D can be expressed as follows:


 * $$\overrightarrow{AD} = (x_D-x_A)\vec i +(y_D-y_a)\vec j = x_D \vec i + y_D \vec j$$

By definition $$\overrightarrow{AD}$$ is the displacement vector which is defined as


 * $$\overrightarrow{AD} = d_3 \vec i + d_4 \vec j$$

where $$d_3$$ and $$d_4$$ are the displacement DoFs previously solved for.


 * $$\therefore (x_D,y_D) = (d_3,d_4)$$


 * $$ \Rightarrow (x_D,y_D) = (4.3502,6.1254)$$

3-Bar Truss System


When numbering local nodes be sure to consider convenience in assembling the global stiffness matrix, $$\mathbf{\underline{K}}$$. Careful selection of local nodes can greatly ease the process and structure of the global stiffness matrix. The figure below shows the most convenient local node numbering.



Class Notes: Friday, October 3, 2008
There is a clarification/correction to the angle for element 3, $$\theta^{(3)}$$, as shown by the following image and equation...



$$ \theta^{(3)}=-90^{\circ}-45^{\circ}=-135^{\circ} $$

Now that the angle for element 3 is clear, we can examine the three-bar truss system with respect to the reaction forces and the external force P...



We know the following equations to be true for this three-bar truss...

$$\Sigma{F}_{x}=0\, $$

$$\Sigma{F}_{y}=0\, $$

$$\Sigma{M}_{A}=0\, $$

where $$\Sigma{M}_{A}=0$$ is trivial, due to the fact that all force vectors pass through this point.

But what about point B ~ $$\Sigma{M}_{B}\, $$?

We will use the following image and subsequent equations to explain the moment about point B. Note: Points A and B are unrelated to the original three-bar truss figure...



The following is true for all $$A'\,$$ on the line of action of $$\vec{F}$$... $$ \Sigma{\vec{M}}_{B}=\vec{BA}\times\vec{F}=\vec{BA'}\times\vec{F}$$ $$ \vec{BA'}=\vec{BA}+\vec{AA'}$$ $$ \Sigma{\vec{M}}_{B}=(\vec{BA}+\vec{AA'})\times\vec{F}=\vec{BA}\times\vec{F}+\vec{AA'}\times\vec{F}$$

where we know that $$\vec{AA'}\times\vec{F}=\vec{0}\, $$. Now we go back to the three-bar truss: Node $$A$$ is in equilibrium: $$ \sum_{i=0}^3\vec{F}=\vec{0}\, $$

And we know that $$ \sum_{i}\vec{M}_{B,i}=\sum_{i}\vec{BA'_{i}}\times\vec{F_{i}}\, $$

$$A'_{i}\, $$ is any point on the line of action of $$\vec{F_{i}}\, $$

$$ \sum_{i}\vec{M}_{B,i}=\sum_{i}\vec{BA'_{i}}\times\vec{F_{i}}=\vec{BA}\times\sum_{i}\vec{F_{i}}=\vec{0}\, $$, therefore proving $$ \Sigma{M}_{B}=\vec{0}\, $$.

Next, we move on to the global stiffness matrix for the three-bar truss. A general outline of the overlapping nature of the three individual element stiffness matrices is shown below...



Next, we fill this matrix...

HW: Complete the global stiffness matrix for the three-bar truss system.
Using the method described in HW2 and the values $$ k_{(1)}=1.2 $$, $$ k_{(2)}=0.8 $$ and $$ k_{(3)}=0.6 $$, we can find each individual element stiffness matrix...

$$ k^{(1)} = \left[ \begin{array}{llll} 0.9 & 0.52 & -0.9 & -0.52\\ 0.52 & 0.3 & -0.52 & -0.3\\ -0.9 &-0.52 & 0.9 & 0.52\\ -0.52 & -0.3 & 0.52 & 0.3 \end{array} \right] $$

$$ k^{(2)} = \left[ \begin{array}{llll} 0.6 & -0.346 & -0.6 & 0.346\\ -0.346 & 0.2 & 0.346 & -0.2\\ -0.6 & 0.346 & 0.6 & -0.346\\ 0.346 & -0.2 & -0.346 & 0.2 \end{array} \right] $$

$$ k^{(3)} = \left[ \begin{array}{llll} 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 \end{array} \right] $$

$$ \underline{K} = \left[ \begin{array}{llllllll} 0.9 & 0.52 & -0.9 & -0.52 & 0 & 0 & 0 & 0\\ 0.52 & 0.3 & -0.52 & -0.3 & 0 & 0 & 0 & 0\\ -0.9 & -0.52 & 1.8 & 0.474 & -0.6 & 0.346 & -0.3 & -0.3\\ -0.52 & -0.3 & 0.474 & 0.8 & 0.346 & -0.2 & -0.3 & -0.3\\ 0 & 0 & -0.6 & 0.346 & 0.6 & -0.346 & 0 & 0\\ 0 & 0 & 0.346 & -0.2 & -0.346 & 0.2 & 0 & 0\\ 0 & 0 & -0.3 & -0.3 & 0 & 0 & 0.3 & 0.3\\ 0 & 0 & -0.3 & -0.3 & 0 & 0 & 0.3 & 0.3\end{array} \right] $$

Matlab Homework - Truss Example
The following image is the modified output of the above code...

People who contributed to HW3
Jose Jayma 21:07, 5 October 2008 (UTC)

Eml4500.f08.FEABBQ.vitt 21:40, 5 October 2008 (UTC)

Stark 22:20, 5 October 2008 (UTC)

Esmail Hadman 03:27, 8 October 2008 (UTC)

William Kurth 14:39, 8 October 2008 (UTC)

Andy Koby 16:12, 8 October 2008 (UTC)