User:Eml4500.f08.echo.mott/ECHO.HW2

This page contains material required for EML4500 Homework Report No.2. The previous version of this page was submitted to the class.

Statically Determinate Problem
Here is an example of a statically determinate problem. It can be solved using the traditional methods learned in statics.

2 Force body example:




 * Drawing and Labeling the Global Picture:




 * Breaking it down into two elements:



$$\sum F_x = 0; R_1cos( \theta_1) - R_2cos( \theta_2) + Pcos( \theta_P) = 0$$

$$\sum F_y = 0; R_1sin( \theta_1) - R_2sin( \theta_2) + Psin( \theta_P) = 0$$

$$\sum M_1 = 0; R_2sin( \theta_2)*(l_2cos( \theta_2) + l_1cos( \theta_1)) + Psin( \theta_P)*(l_1cos( \theta_1)) = 0$$

Three equations and three unknowns, therefore it is statically determinate.

Statically Indeterminate
Some problems cannot be solved using the 3 equations from statics. That is, if the problem contains more than 3 unknowns it cannot be solved using traditional statics methods.

The example on the left can be solved because each member is a 2-force member. Hence, there are only 3 unknowns. The example on the right is a statically indeterminate structure because there are 6 unknowns.



The Element Stiffness Matrix
After the diagrams have been drawn for each element in the structure, the element stiffness matrix should be computed for each element. The stiffness matrix k looks like this:

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

The scalar $$k^{(e)}$$ can be found from the properties of the structure where E is the Elastic Modulus, A is the cross-sectional area of the element, and L is the length of the element.

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

$$l^{(e)}$$ is the x-length of the element in reference to global coordinates and $$m^{(e)}$$ is the y-length of the element from a global perspective. They can be thought of as the length of the element projected onto global coordinates from local coordinates.

The x-length can be found by using: $$l_{s}= \frac{x_{2}-x_{1}}{L}$$

The y-length of the element can be found by using: $$l_{s}= \frac{y_{2}-y_{1}}{L}$$

With the inclined angle referenced to the horizontal, the horizontal and vertical lengths can be found by:

$${l_s}^x = {cos \theta}$$

$${l_s}^y = sin \theta$$



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

$$m^{(e)}= \mathbf{\tilde{i}} \cdot \mathbf{j}=cos( \frac{ \pi}{2} - \theta^{(e)})=sin( \theta^{(e)})$$

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

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

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

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

After reviewing the stiffness matrix two observations can be made:

1) Of all the numbers in the matrix, only three need to be computed. The other coefficients have the same absolute values.
 * $$l_s^{}m_s $$
 * $${l_s}^2$$
 * $${m_s}^{2}$$

2) Matrix k is symmetric. That is $$k_{i j}^{1}=k_{j i}^{1}$$.

The Global Stiffness Matrix
Consider the Global Force-Displacement Relationship for a Truss with 2 elastic bars. There are 6 degrees of freedom in the global picture.

Written in compact notation we have the following:

$$\begin{bmatrix}K_{ij}^{}\end{bmatrix}_{6x6} \left\{d_{ij}\right\}_{6x1} = \left\{F_{i}\right\}_{6x1}$$

In the more general form $${nxn}$$ :

$$\sum_{j=1}^{6} K_{ij}d_{j}= F_{i}$$; i = 1, 2, ..., 6

$$K_{nxn} = \begin{bmatrix}K_{ij}^{}\end{bmatrix}_{nxn}$$ = global stiffness matrix

$$d_{nx1} = \left\{d_{j}^{}\right\}_{nx1}$$ = global displacement matrix

$$F_{nx1} = \left\{F_{i}^{}\right\}_{nx1}$$ = global force matrix

Now recall the Element Force-Displacement relationship:

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

$$k_{4x4}^{(e)} = \begin{bmatrix}k_{ij}^{(e)}\end{bmatrix}_{4x4}$$ = element stiffness matrix

$$d_{4x1}^{(e)} = \left\{d_{j}^{(e)}\right\}_{4x1}$$ = element displacement matrix

$$f_{4x1}^{(e)} = \left\{f_{i}^{(e)}\right\}_{4x1}$$ = element force matrix

Use an assembly process to go from element matrices(stiffness, displacement, force) to global matrices.


 * First, Identify the correspondence between the element displacement degrees of freedom's and the global displacement degrees of freedom's.

Global Level: $$\left\{d_1, d_2, \dots, d_6\right\}$$

Element level:
 * $$\left\{d_1^{(1)}, d_2^{(1)}, d_3^{(1)}, d_4^{(1)}\right\}$$
 * $$\left\{d_1^{(2)}, d_2^{(2)}, d_3^{(2)}, d_4^{(2)}\right\}$$

The relationship between global-local(element) dof's are as follows:

$$\text{Global Node 1} = \begin{cases} d_1 = d_1^{(1)} \\ d_2 = d_2^{(1)}\end{cases}$$

$$\text{Global Node 2} = \begin{cases} d_3 = d_3^{(1)} = d_1^{(2)} \\ d_4 = d_4^{(1)} = d_2^{(2)}\end{cases}$$

$$\text{Global Node 3} = \begin{cases} d_5 = d_3^{(2)} \\ d_6 = d_4^{(2)}\end{cases}$$

Conceptual step of the assembly:



Taking the global stiffness matrix and substituting in the appropriate values:

$$\mathbf{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}= \begin{bmatrix} K_{11} & K_{12} & K_{13} & K_{14} & K_{15} & K_{16}\\K_{21} & K_{22} & K_{23} & K_{24} & K_{25} & K_{26}\\K_{31} & K_{32} & K_{33} & K_{34} & K_{35} & K_{36}\\K_{41} & K_{42} & K_{43} & K_{44} & K_{45} & K_{46}\\K_{51} & K_{52} & K_{53} & K_{54} & K_{55} & K_{56}\\K_{61} & K_{62} & K_{63} & K_{64} & K_{65} & K_{66}\end{bmatrix}$$

Plugging in values,

$$K_{11}=k_{11}^{(1)}=\frac{9}{16}=0.5625$$

$$K_{12}=K_{21}=k_{12}^{(1)}=\frac{3\sqrt{3}}{16}=0.3248$$

$$K_{13}=K_{31}=k_{13}^{(1)}=-k_{11}^{(1)}=-0.5625$$

$$K_{14}=K_{23}=K_{32}=K_{41}=k_{14}^{(1)}=-k_{12}^{(1)}=-0.3248$$

$$K_{22}=k_{22}^{(1)}=\frac{3}{16}=0.1875$$

$$K_{24}=K_{42}=k_{24}^{(1)}=-k_{22}^{(1)}=-0.1875$$

$$K_{33}=k_{33}^{(1)}+k_{11}^{(2)}=	\begin{matrix} \frac{9}{16} \end{matrix}+\begin{matrix} \frac{5}{2} \end{matrix}=3.0625$$

$$K_{34}=k_{34}^{(1)}+k_{12}^{(2)}=	\begin{matrix} \frac{3\sqrt{3}}{16} \end{matrix}+\begin{matrix} \frac{-5}{2} \end{matrix}=-2.1752$$

$$K_{35}=K_{53}=k_{13}^{(2)}=-k_{11}^{(2)}=\frac{-5}{2}=-2.5$$

$$K_{36}=K_{45}=K_{54}=K_{63}=k_{14}^{(2)}=-k_{12}^{2}=2.5$$

$$K_{43}=K_{34}=k_{43}^{(1)}+k_{21}^{(2)}=k_{34}^{(1)}+k_{12}^{(2)}=-2.1752$$

$$K_{44}=k_{44}^{(1)}+k_{22}^{(2)}=\frac{3}{16}+\frac{5}{2}=2.6875$$

$$K_{46}=K_{64}=k_{24}^{(2)}=\frac{-5}{2}=-2.5$$

$$K_{55}=k_{33}^{(2)}=\frac{5}{2}=2.5$$

$$K_{56}=K_{65}=k_{34}^{(2)}=\frac{-5}{2}=-2.5$$

$$K_{66}=k_{44}^{(2)}=\frac{5}{2}=2.5$$

$$K_{15}^{}=K_{16}=K_{25}=K_{26}=K_{51}=K_{52}=K_{61}=K_{62}=0$$

$$\mathbf{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 Unknown Displacements
Next, the elimination of known degree's of freedom can create the following matrices.

$$ d_1^{} = d_2^{} = d_5^{} = d_6^{}=0 $$

$$	\mathbf{K}\begin{Bmatrix} d_1=0\\d_2=0\\d_3\\d_4\\d_5=0\\d_6=0 \end{Bmatrix}=\begin{vmatrix} k_{13} & k_{14} \\ k_{23} & k_{24}\\k_{33} & k_{34}\\k_{43} & k_{44}\\k_{53} & k_{54}\\k_{63} & k_{64}\end{vmatrix}\begin{Bmatrix} d_3\\d_4\end{Bmatrix}=\mathbf{F}$$

Applying fixed boundary conditions allows us to delete the corresponding columns in the global stiffness matrix K

By the Principal of Virtual Work (PVW) we're able to delete the corresponding rows (in this case - 1, 2, 5, 6). Corresponding rows of F are also deleted.

The resulting FD relation is: $$ \begin{bmatrix} k_{33} & k_{34} \\ k_{43} & k_{44} \end{bmatrix} \begin{Bmatrix} d_3\\d_4\end{Bmatrix}=\begin{Bmatrix} F_3\\F_4\end{Bmatrix}$$

$$\mathbf{F}=\begin{Bmatrix} F_3\\F_4\end{Bmatrix}=\begin{Bmatrix} 0\\P\end{Bmatrix}$$

And P =>

What is the inverse of K without using a calculator?

Taking the determinant of K:

$$det \mathbf{K}=k_{33}k_{44}-k_{k34}k_{43}$$

$$\mathbf{K}^{-1}=\frac{1}{det\mathbf{K}}\begin{bmatrix} k_{44} & -k_{34} \\ -k_{43} & k_{33} \end{bmatrix}$$

Verifying:

$$\mathbf{K}\mathbf{K}^{-1}=\mathbf{K}^{-1}\mathbf{K}=\mathbf{I}=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$$

This is essentially:

$$\frac{1}{det(\mathbf{K})}\begin{bmatrix} det(\mathbf{K}) & 0 \\ 0 & det(\mathbf{K}) \end{bmatrix}=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$$

Note: $$\mathbf{K}^T=\begin{bmatrix} k_{33} & k_{43} \\ k_{34} & k_{44} \end{bmatrix}$$

Therefore: $$\mathbf{K}^{-1}\ne \ \frac{1}{det(\mathbf{K})}\mathbf{K}^T$$

Going back to the 2 bar truss problem...

$$\begin{Bmatrix} d_3 \\ d_4 \end{Bmatrix}=\mathbf{K}^{-1}\begin{Bmatrix} 0 \\ P \end{Bmatrix}=\begin{Bmatrix} 4.352 \\ 6.1271 \end{Bmatrix}$$

Solving for Reactions
Now we must compute the reactions:


 * Method 1:

Use element FD relations.

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

Taking a look at element 1 again and applying the knowns:



$$d^{(1)}=\begin{Bmatrix} 0 \\ 0 \\ 4.352 \\ 6.1271 \end{Bmatrix}$$

Taking a look at element 2:



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

$$d^{(1)}=\begin{Bmatrix} 4.352 \\ 6.1271 \\ 0 \\ 0 \end{Bmatrix}=\begin{Bmatrix} d_1^{(2)} \\ d_2^{(2)} \\ d_3^{(2)} \\ d_4^{(2)} \end{Bmatrix}$$

Code

 * The 2-bar truss example outlined above can be solved using MATLAB. The code is shown here:

clear all;

e = [3 5]; A = [1 2]; P = 7; L=[4 2]; alpha = pi/3; beta = pi/4;

nodes = [0, 0]; L(1)*cos(pi/2-alpha), L(1)*sin(pi/2-alpha); L(1)*cos(pi/2-alpha)+L(2)*sin(beta),L(1)*sin(pi/2-alpha)-L(2)*cos(beta)];

dof = 2*length(nodes);

conn = [1,2; 2,3];

lmm = [1, 2, 3, 4; 3, 4, 5, 6];

elems = size(lmm,1);

K=zeros(dof); R = zeros(dof,1); debc = [1, 2, 5, 6];

ebcVals = zeros(length(debc),1);


 * The load vector

R = zeros(dof,1); R(4) = P;


 * Assembly of 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


 * The 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

Results

 * The code shown above generates the following results:

k_local = 0.75       -0.75        -0.75         0.75

k = 0.5625     0.32476      -0.5625     -0.32476      0.32476       0.1875     -0.32476      -0.1875      -0.5625     -0.32476       0.5625      0.32476     -0.32476      -0.1875      0.32476       0.1875

k_local = 5   -5           -5     5

k = 2.5        -2.5         -2.5          2.5         -2.5          2.5          2.5         -2.5         -2.5          2.5          2.5         -2.5          2.5         -2.5         -2.5          2.5

K = 0.5625     0.32476      -0.5625     -0.32476            0            0      0.32476       0.1875     -0.32476      -0.1875            0            0      -0.5625     -0.32476       3.0625      -2.1752         -2.5          2.5     -0.32476      -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

R = 0    0     0     7     0     0

d = 0           0        4.352       6.1271            0            0

reactions = -4.4378     -2.5622       4.4378      -4.4378

results = 1.7081      5.1244       8.5406       5.1244       17.081       0.6276       1.8828        3.138       1.8828        6.276


 * The previous code was taken from the EML4500 course website. A text version of the data can be found                                                                                                             here

Contributing Members HW2
Chris Mott 18:40, 26 September 2008 (UTC)

Charles 18:41, 26 September 2008 (UTC)

Bradley LaCroix Eml4500.f08.echo.lacroix 18:54, 26 September 2008 (UTC)

C. Sell 18:58, 26 September 2008 (UTC)

Brandon Sell Eml4500.f08.echo.sell.br 19:10, 26 September 2008 (UTC)

Jason Bruce 18:40, 26 September 2008 (UTC)