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

Lecture Notes
Continuing where we left off from last homework...

Step 2
The second step is to draw the element pictures. In these pictures you should label the element dofs and the element forces in either global or local coordinates.





Step 3
The third step is to develop the global force-displacement relationship at the element level. As we saw earlier, for a specified element the force-displacement relationship is: $$\textbf{k}^{(e)}*\textbf{d}^{(e)}=\textbf{f}^{(e)}$$

Where k(e) is the element stiffness matrix, d(e) is the element displacement matrix and f(e) the element force matrix.

The axial stiffness of a specified element is determined by the equation:

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

Where E is the Young's modulus for the element, A is the cross-sectional area for the element and L is the length if the element.

Now we need to find the director cosines (l(e) and m(e)) for each element using the following equations:

$$l^{(e)}=\textbf{i}^{\sim}\cdot\textbf{i}=cos\theta^{(e)}$$

$$m^{(e)}=\textbf{j}^{\sim}\cdot\textbf{j}=sin\theta^{(e)}$$



Using the director cosines, axial stiffness and the following formula we can determine the element stiffness matrix.

$$\textbf{k}^{(e)}=k^{(e)}\left[\begin{array}{c c c c} (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{array}\right ]$$

For element 1:

$$\theta=30$$

$$l=cos\theta=\frac{\sqrt{3}}2$$

$$m=sin\theta=\frac{1}2$$

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

This gives us the stiffness matrix for element 1:

$$k^{(1)}=\left[\begin{array}{c c c c} \frac{9}{16} & \frac{3\sqrt{3}}{16} & -\frac{9}{16} & -\frac{3\sqrt{3}}{16}\\ \frac{3\sqrt{3}}{16} & \frac{3}{16} & -\frac{3\sqrt{3}}{16} & -\frac{3}{16}\\ -\frac{9}{16} & -\frac{3\sqrt{3}}{16} & \frac{9}{16} & \frac{3\sqrt{3}}{16}\\ -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & \frac{3\sqrt{3}}{16} & \frac{3}{16} \end{array}\right]$$

Notice that the stiffness matrix for element 1 is symmetric.

For element 2:

$$\theta=-45$$

$$l=cos(-45)=\frac{\sqrt{2}}2$$

$$m=sim(-45)=\frac{\sqrt{2}}2$$

$$k=\frac{EA}L=\frac{(5)(2)}2=5$$

This gives us the stiffness matrix for element 2:

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

Notice that the stiffness matrix for element 2 is also symmetric.

Now we have to assemble the global stiffness matrix. This is illustrated in the picture below where the shaded area represents the addition of overlapping parts of $$\textbf{k}^{(1)}$$ and $$\textbf{k}^{(2)}$$.



$$K_{33}=k_{33}^{(1)} + k_{11}^{(2)}$$

$$K_{34}=k_{34}^{(1)} + k_{12}^{(2)}$$

$$K_{43}=k_{43}^{(1)} + k_{21}^{(2)}$$

$$K_{44}=k_{44}^{(1)} + k_{22}^{(2)}$$

For this example we get:

$$\textbf{K}=\left[\begin{array}{c c c c c c} \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{49}{16} & \frac{-40+3\sqrt{3}}{16} & -\frac{5}2 & -\frac{5}2\\ -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & \frac{-40+3\sqrt{3}}{16} & \frac{43}{16} & -\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{array}\right]$$

Step 4
Now on the fourth step we eliminate known dofs to reduce the global force-displacement relationship.

Due to the applied constraints, $$d_1=d_2=d_5=d_6=0$$. This allows us to eliminate columns 1,2,5 and 6 due to the principal of virtual work. We also note that $$F_1=F_2=F_3=F_5=F_6=0$$. This allows us to eliminate rows 1,2,5 and 6 by the same principal. The resulting, simplified force-displacement relationship follows:

$$\left[\begin{array}{c c} K_{33} & K_{34}\\ K_{34} & K_{44} \end{array}\right] \begin{Bmatrix} d_3\\ d_4 \end{Bmatrix}=\begin{Bmatrix} F_3\\ F_4 \end{Bmatrix}$$

Let $$F_3=0$$ and $$F_4=7$$.Now plugging in numbers we get:

$$\left[\begin{array}{c c} \frac{49}{16} & \frac{-40+3\sqrt{3}}{16}\\ \frac{-40+3\sqrt{3}}{16} & \frac{43}{16} \end{array}\right] \begin{Bmatrix} d_3\\ d_4 \end{Bmatrix}= \begin{Bmatrix} 0\\ 7 \end{Bmatrix}$$

Now we can solve for the unknown displacement dofs by multiplying both sides by:

$$\textbf{K}^{-1}=\frac{1}{det\textbf{K}} \left[\begin{array}{c c} K_{44} & -K_{34}\\ -K_{43} & K_{33} \end{array}\right]$$.

Plugging in numbers we get:

$$\textbf{K}^{-1}=\frac{1}{3.4988} \left[\begin{array}{c c} \frac{49}{16} & \frac{40-3\sqrt{3}}{16}\\ \frac{40-3\sqrt{3}}{16} & \frac{43}{16} \end{array}\right]$$

It is important to note that: $$\textbf{K}^{-1}\neq\frac{1}{det\textbf{K}}\textbf{K}^T$$

Following through with this, we get:

$$\begin{Bmatrix} d_3\\ d_4\end{Bmatrix}= \textbf{K}^{-1}\begin{Bmatrix} F_3\\ F_4\end{Bmatrix}$$

Plugging in the numbers now gives:

$$\begin{Bmatrix} d_3\\ d_4\end{Bmatrix}= \left[\begin{array} {c c} .8753 & .6217\\ .627 & .7681 \end{array}\right] \begin{Bmatrix} 0\\ 7\end{Bmatrix}$$

Solving this system gives us:

$$d_3=4.352$$

$$d_4=6.127$$

Step 5
Now in step five we can compute reactions by one of two methods.

The first method is to use the element force-displacement relationship.

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

For element 1 we would have:

$$\left[\begin{array}{c c c c} \frac{9}{16} & \frac{3\sqrt{3}}{16} & -\frac{9}{16} & -\frac{3\sqrt{3}}{16}\\ \frac{3\sqrt{3}}{16} & \frac{3}{16} & -\frac{3\sqrt{3}}{16} & -\frac{3}{16}\\ -\frac{9}{16} & -\frac{3\sqrt{3}}{16} & \frac{9}{16} & \frac{3\sqrt{3}}{16}\\ -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & \frac{3\sqrt{3}}{16} & \frac{3}{16} \end{array}\right]\begin{Bmatrix} 0\\ 0\\ 4.352\\ 6.127\end{Bmatrix}= \begin{Bmatrix} f_1^{(1)}\\ f_2^{(1)}\\ f_3^{(1)}\\ f_4^{(1)}\end{Bmatrix}$$

Solving this system we get:

$$f_1^{(1)}=-4.438$$

$$f_2^{(1)}=-2.562$$

$$f_3^{(1)}=4.438$$

$$f_4^{(1)}=2.562$$

For element 2 we have:

$$\left[\begin{array}{c c c c} \frac{5}2 & -\frac{5}2 & \frac{5}2 & -\frac{5}2\\ -\frac{5}2 & \frac{5}2 & -\frac{5}2 & \frac{5}2\\ -\frac{5}2 & \frac{5}2 & -\frac{5}2 & \frac{5}2\\ \frac{5}2 & -\frac{5}2 & \frac{5}2 & -\frac{5}2 \end{array}\right]\begin{Bmatrix} 4.352\\ 6.127\\ 0\\ 0\end{Bmatrix}= \begin{Bmatrix} f_1^{(2)}\\ f_2^{(2)}\\ f_3^{(2)}\\ f_4^{(2)}\end{Bmatrix}$$

Solving this system we get:

$$f_1^{(2)}=-4.438$$

$$f_2^{(2)}=4.438$$

$$f_3^{(2)}=4.438$$

$$f_4^{(2)}=-4.438$$

Alternate Method
Another method to solve the truss problem is by using statics. Since the system is statically indeterminate we should use the Euler Cut method. In the Euler Cut method, the global free body diagram can be cut at any point and into and dimensions that make the problem easier to solve.



Proof that P1 = P2
We are asked to prove that P1 = P2. From previous calculations we have a full FD calculation for element 1:

$$\begin{bmatrix} -0.5625 & -0.32476 \\ -0.32476 & -0.1875 \\ 0.5625 & 0.32476 \\ 0.32476 & 0.1875 \end{bmatrix}\begin{Bmatrix} 4.352 \\ 6.1271 \end{Bmatrix} = \begin{Bmatrix} -4.4378 \\ -2.5622 \\ 4.4378 \\ 2.5622 \end{Bmatrix} = \begin{Bmatrix} f_1^{(1)} \\ f_2^{(1)} \\ f_3^{(1)} \\ f_4^{(1)} \end{Bmatrix}$$

The root sum squared method can be used for the proof required. Generally the root sum squared equation is:

$$P_1^{(e)} = \sqrt{\left(f_1^{(e)}\right)^2 + \left(f_2^{(e)}\right)^2}$$

We can then expand that to our condition and create a system of equations for P1 and P2.

$$P_1^{(1)} = \sqrt{\left(f_1^{(1)}\right)^2 + \left(f_2^{(1)}\right)^2}$$

$$P_2^{(1)} = \sqrt{\left(f_3^{(1)}\right)^2 + \left(f_3^{(1)}\right)^2}$$

From the FD matrix we have values for all of the local forces used in this system so P1 and P2 can be calculated directly and their equality proven.

$$P_1^{(1)} = \sqrt{(-4.4378)^2 + (-2.5622)^2} = 5.1243$$

$$P_2^{(1)} = \sqrt{(4.4378)^2 + (2.5622)^2} = 5.1243$$

Derivation of $$\vec{\tilde{i}}\cdot\vec{j}$$
The definition of $$\vec{\tilde{i}}$$ given in the lecture notes is: $$\vec{\tilde{i}}=\cos\theta^{(e)}\vec{i}+\sin\theta^{(e)}\vec{j}$$

Taking the dot product of $$\vec{\tilde{i}}$$ and $$\vec{j}$$ yields: $$\vec{\tilde{i}}\cdot\vec{j}=\cos\theta^{(e)}\vec{i}\cdot\vec{j}+\sin\theta^{(e)}\vec{j}\cdot\vec{j}$$

By definition, $$\vec{i}\cdot\vec{j}$$ is equal to 0, and $$\vec{j}\cdot\vec{j}$$ is equal to 1. Using these relationships reduces the above equation to: $$\vec{\tilde{i}}\cdot\vec{j}=\sin\theta^{(e)}$$

Previously, the direction cosine $$m^{\left(e\right)}$$ was defined as $$ \cos \left(\frac{\pi }{2}-\theta ^{(e)}\right)$$, which is equivalent to $$\sin\theta^{\left(e\right)}$$. Thus the derived expression for $$\vec{\tilde{i}}\cdot\vec{j}$$ is: $$\vec{\tilde{i}}\cdot\vec{j}=m^{(e)}$$

Derivation of Element Stiffness Matrices
$$ \bar{k}^{(e)}= \begin{bmatrix} k_{11}^{(e)} & k_{12}^{(e)} & k_{13}^{(e)} & k_{14}^{(e)} \\ k_{21}^{(e)} & k_{22}^{(e)} & k_{23}^{(e)} & k_{24}^{(e)} \\ k_{31}^{(e)} & k_{32}^{(e)} & k_{33}^{(e)} & k_{34}^{(e)} \\ k_{41}^{(e)} & k_{42}^{(e)} & k_{43}^{(e)} & k_{44}^{(e)} \end{bmatrix}; e = 1,2 $$

For an elastic bar element:

$$ \bar{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}; K^{(e)} = \frac{E^{(e)}A^{(e)}}{L^{(e)}} $$

$$ l^{(e)} = \frac{(x_j-x_i)}{L^{(e)}}; m^{(e)} = \frac{(y_j-y_i)}{L^{(e)}} $$

Where x and y are the global coordinates of the nodes i and j at either end of the bar element.

Using MATLAB
Global coordinates of nodes:

nodes = [x1, y1; x2, y2; x3, y3];

>> nodes = [0, 0; 4*cos(pi/6) 4*sin(pi/6); 4*cos(pi/6)+2*cos(pi/4) 4*sin(pi/6)-2*sin(pi/4)]

nodes =

0           0       3.4641            2       4.8783      0.58579

Directional forces l and m: >> l_1 = (3.4641-0)/4

l_1 =

0.86603

>> m_1 = (2-0)/4

m_1 =

0.5

>> l_2 = (4.8783-3.4641)/2

l_2 =

0.7071

>> m_2 = (0.5858-2)/2

m_2 =

-0.7071

Elemental Stiffness Matrices: >> k_1 = 3*1/4*[l_1^2 l_1*m_1 -l_1^2 -l_1*m_1; l_1*m_1 m_1^2 -l_1*m_1 -m_1^2; -l_1^2 -l_1*m_1 l_1^2 l_1*m_1; -l_1*m_1 -m_1^2 l_1*m_1 m_1^2]

k_1 =

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_2 = 5*2/2*[l_2^2 l_2*m_2 -l_2^2 -l_2*m_2; l_2*m_2 m_2^2 -l_2*m_2 -m_2^2; -l_2^2 -l_2*m_2 l_2^2 l_2*m_2; -l_2*m_2 -m_2^2 l_2*m_2 m_2^2]

k_2 =

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

Node Equation Verification
The homework question on page 11-4 of the class notes asked the class to verify the node equation for the following node:



Using the laws of statics ($$\sum $$Fx = 0 and $$\sum $$Fy = 0) we know that sin(30) * p(1)1 must equal sin(45) * p(2)2. We also know that P must equal cos(30) * p(1)1 + cos(45) * p(2)2. These two equations simplify as follows:

$$(1/2)p^1_1 - (\sqrt(2)/2)p^2_2 = 0$$

$$(\sqrt(3)/2)p^1_1 + (\sqrt(2)/2)p^2_2 - P = 0$$

The force P is a known value therefore this is a system of two equations with two unknowns ($$p^1_1, p^2_2)$$. Solving for these unknowns in terms of P gives us:

$$p^1_1 = P/(1/2+\sqrt(3)/2)$$

$$p^2_2 = P/[\sqrt(2)(1/2+\sqrt(3)/2)]$$

MATLAB Code
The code below represents the program that calculates the displacements and reactions of the nodes of the system described above. The code utilizes three functions which are below the main code. The professor provided the main code while the three functions were copied from code provided by the online resources of the textbook

Two Bar Truss System.m
% Two bar truss example 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);

%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

PlaneTrussElement.m
function k = PlaneTrussElement(e, A, coord) % PlaneTrussElement(e, A, coord) % Generates stiffness matrix for a plane truss element % e = modulus of elasticity % A = area of cross-section % coord = coordinates at the element ends

x1=coord(1,1); y1=coord(1,2); x2=coord(2,1); y2=coord(2,2); L=sqrt((x2-x1)^2+(y2-y1)^2); ls=(x2-x1)/L; ms=(y2-y1)/L; k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms; ls*ms, ms^2,-ls*ms,-ms^2; -ls^2,-ls*ms,ls^2,ls*ms; -ls*ms,-ms^2,ls*ms,ms^2];

NodalSoln.m
function [d, rf] = NodalSoln(K, R, debc, ebcVals) % [nd, rf] = NodalSoln(K, R, debc, ebcVals) % Computes nodal solution and reactions % K = global coefficient matrix % R = global right hand side vector % debc = list of degrees of freedom with specified values % ebcVals = specified values dof = length(R); df = setdiff(1:dof, debc); Kf = K(df, df); Rf = R(df) - K(df, debc)*ebcVals; dfVals = Kf\Rf; d = zeros(dof,1); d(debc) = ebcVals; d(df) = dfVals; rf = K(debc,:)*d - R(debc);

Plane Truss Results.m
function results = PlaneTrussResults(e, A, coord, disps) % results = PlaneTrussResults(e, A, coord, disps) % Compute plane truss element results % e = modulus of elasticity % A = Area of cross-section % coord = coordinates at the element ends % disps = displacements at element ends % The output quantities are eps = axial strain % sigma = axial stress and force = axial force.

x1=coord(1,1); y1=coord(1,2); x2=coord(2,1); y2=coord(2,2); L=sqrt((x2-x1)^2+(y2-y1)^2); ls=(x2-x1)/L; ms=(y2-y1)/L; T=[ls,ms,0,0; 0,0,ls,ms]; d = T*disps; eps= (d(2)-d(1))/L; sigma = e.*eps; force = sigma.*A; results=[eps, sigma, force];

The output to the command window when the program is executed is shown below. The matrices k_local and k are the local and elemental stiffness matrices. The global stiffness matrix is K. The displacement matrix is d. The x and y reactions of global nodes 1 and 3 are shown in the reactions matrix. The results matrix is composed of a 2 by 1 and two 2 by 2 matrices concatenated together. The first column in results represents the strain for elements 1 and 2. The second and third columns represent the stress matrix of the elements where the diagonals are the correct values. The fourth and fifth columns represent the axial force matrix of the elements where, again, the diagonals are the correct values.

>> TwoBarTrussSystem

k_local =

0.7500  -0.7500   -0.7500    0.7500

k =

0.5625   0.3248   -0.5625   -0.3248    0.3248    0.1875   -0.3248   -0.1875   -0.5625   -0.3248    0.5625    0.3248   -0.3248   -0.1875    0.3248    0.1875

k_local =

5   -5    -5     5

k =

2.5000  -2.5000   -2.5000    2.5000   -2.5000    2.5000    2.5000   -2.5000   -2.5000    2.5000    2.5000   -2.5000    2.5000   -2.5000   -2.5000    2.5000

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

R =

0    0     0     7     0     0

d =

0        0    4.3520    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

Contributing Team Members
Aaron Fisher - Eml4500.f08.Lulz.fisher 18:12, 26 September 2008 (UTC)

Eric Layton - Eml4500.f08.Lulz.Layton.Eric 17:55, 26 September 2008 (UTC)

Sam Miorelli - Eml4500.f08.lulz.abcd 18:01, 26 September 2008 (UTC)

Benjamin Mitchell - Eml4500.f08.Lulz.mitchell.bm 17:58, 26 September 2008 (UTC)

John Saxon - Eml4500.f08.Lulz.js 16:27, 26 September 2008 (UTC)

Andrew Strack - Eml4500.f08.lulz.strack 18:50, 26 September 2008 (UTC)