User:Eml4500.f08.qwiki.bishop/HW2

 Comment: The pictures that were previously used in this homework submission have been changed. The previous images ( and ) were not originally authored by Team Qwiki, however I did not intentionally plagiarize Team Echo's work. While working on Homework 2, both Team Qwiki and Team Echo created two images, with the same name: and . Due to the nature of MediaWiki, if two different teams create an image with the same name, only the latest version is displayed when called upon in an article. Therefore, even though Team Qwiki uploaded two of our own images (Element1.jpg and Element2.jpg) during homework 2, Team Echo's images were displayed because their upload was the most recent. In the future, we will make sure to upload files that have a unique name. The correct images that were authored by Team Qwiki are included in this updated version of Homework 2.

A comparison between the current version of Homework 2 and our prior submission can be found here. Eric Gatch 

=Lecture 6 - Monday, September 8, 2008=

Continued from Lecture 5
2. Element Picture

3. Global FD at element level

Element Stiffness Matrix
Element Stiffness Matrix:

$$\mathbf{k^{(e)}}=\begin{bmatrix} \mathbf{k}_{11} & \mathbf{k}_{12} & \mathbf{k}_{13} & \mathbf{k}_{14}\\ \mathbf{k}_{21} & \mathbf{k}_{22} & \mathbf{k}_{23} & \mathbf{k}_{24}\\ \mathbf{k}_{31} & \mathbf{k}_{32} & \mathbf{k}_{33} & \mathbf{k}_{34}\\ \mathbf{k}_{41} & \mathbf{k}_{42} & \mathbf{k}_{43} & \mathbf{k}_{44}\\ \end{bmatrix} $$

Knowing the above element stiffness matrix, is the result of:

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

Where (e) is the element number.

$$k^{e}=\frac{E^{e}A^{e}}{L^{e}}$$ This is the axial stiffness of bar element "e"

$$l^{e},m^{e}$$ are director cosines of $$\tilde{x}$$ axis (from local node 1 to local node 2) with respect to the global (x,y) coordinate system.

$$l^{e}=\vec{\tilde{i}}\cdot \vec{i}=cos\theta^{e}$$

$$m^{e}=\vec{\tilde{i}}\cdot \vec{j}=sin\theta^{e}$$

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

Distributing by unit vector "i" it is easily seen that the equation reduces to $$l^{e}$$ because a unit vector dotted with itself equals 1, whereas a unit vector dotted with its perpendicular counter part equals 0.

Note: the director cosines are the components of $$\vec{\tilde{i}}$$

$$\vec{\tilde{i}}=cos\theta ^{e}\vec{i}+sin\theta ^{e}\vec{j}$$

=Lecture 7 - Wednesday, September 10, 2008= To continue with modeling the 2-bar truss system, element 1 is further examined to find the stiffness matrix k (1). Element 1 is located at an angle $$\theta ^{(1)}={\color{red}30^{\circ}}$$

Therefore inputting this angle into the following equations of the director cosines yields:

$$l^{(1)}=cos\theta ^{(1)}=cos30^{\circ} =\frac{\sqrt{3}}{2}=0.866$$

$$m^{(1)}=sin\theta ^{(1)}=sin30^{\circ} =\frac{1}{2}=0.50$$

Next, to find the stiffness coefficient the below equation is used.

$$k^{(1)}=\frac{E^{(1)} A^{(1)}}{L^{(1)}}=\frac{3*1}{4}=0.75$$

An example of finding out what each element of the stiffness matrix is shown below.

$$\mathbf{k_{11}}=k^{(1)}*(l^{(1)})^{2}=\frac{3}{4}*(\frac{\sqrt{3}}{2})^{2}=\frac{9}{16}=0.5625 $$

To compute element 2, the same steps are taken. Element 2 is at an angle $$\theta ^{(2)}={\color{red}45^{\circ}}$$ Again inputting this angle into the following equations of the director cosines yields:

$$l^{(2)}=cos\theta ^{(2)}=cos45^{\circ} =\frac{\sqrt{2}}{2}=0.707$$

$$m^{(2)}=sin\theta ^{(2)}=sin45^{\circ} =\frac{\sqrt{2}}{2}=0.707$$

Next, to find the stiffness coefficient the below equation is used.

$$k^{(2)}=\frac{E^{(2)} A^{(2)}}{L^{(2)}}=\frac{5*2}{2}=5$$

An example of finding out what each element of the stiffness matrix is shown below.

$$\mathbf{k_{11}}=k^{(2)}*(l^{(2)})^{2}=5*(\frac{\sqrt{2}}{2})^{2}=\frac{5}{2}=2.5 $$

=Lecture 8 - Friday, September 12, 2008=
 * Global Force-Displacement relation:
 * $$\mathbf{K_{(nxn)}d_{(nx1)} = F_{(nx1)}}$$
 * For this example, n=6

$$\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}\begin{pmatrix} d_1\\ d_2\\ d_3\\ d_4\\ d_5\\ d_6\\ \end{pmatrix}=\begin{pmatrix} F_1\\ F_2\\ F_3\\ F_4\\ F_5\\ F_6\\ \end{pmatrix}$$


 * In compact notation, this matrix can be written in the form:

$$\begin{bmatrix} k_{ij}\\ \end{bmatrix}_{(6x6)}\begin{pmatrix} d_j\\ \end{pmatrix}_{(6x1)}=\begin{pmatrix} F_i\\ \end{pmatrix}_{(6x1)}$$


 * (More generally, nxn)

$$\sum_{j=1}^{6}{K_{ij}d_j}=F_i, i=1,...,6$$

The above global matrices can be define as follows:


 * $$K_{(nxn)}=[K_{ij}]_{(nxn)}$$ = Global Stiffness Matrix
 * $$d_{(nx1)}=(d_j)_{(nxn)}$$ = Global Displacement Matrix
 * $$F_{(nx1)}=(f_i)_{(nxn)}$$ = Global Force Matrix

When dealing with a single element, the matrices can be written as:


 * $$K_{(4x4)}^{(e)}=[K_{ij}^{(e)}]_{(4x4)}$$ = Element Stiffness Matrix
 * $$d_{(4x1)}^{(e)}=(d_j^{(e)})_{(4x4)}$$ = Element Displacement Matrix
 * $$F_{(4x1)}^{(e)}=(f_i^{(e)})_{(4x4)}$$ = Element Force Matrix

An assembly process is used to go form an element matrix to a global matrix.

=Lecture 9 - Monday, September 15, 2008= $$\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{pmatrix} d_1\\ d_2\\ d_3\\ d_4\\ d_5\\ d_6\\ \end{pmatrix}=\begin{pmatrix} F_1\\ F_2\\ F_3\\ F_4\\ F_5\\ F_6\\ \end{pmatrix}$$
 * Global Force/Displacement Relationship
 * (Global stiffness matrix, K) X (global displacement, d) = (global force, F)
 * $$\mathbf{K_{(6x6)}d_{(6x1)} = F_{(6x1)}}$$
 * The global stiffness matrix, K, is found by combining the element stiffness matrices $$k^{(1)}$$ & $$k^{(2)}$$ as seen below.
 * NOTE: There is an overlap for displacements $$d_1$$ & $$d_1$$, which results from Element 1 and Element 2 sharing a node; and their k values are then summed together.

=Lecture 10 - Wednesday, September 17, 2008=

Elimination of Known Degree's of Freedom
If d1 = d2 = d5 = d6 = 0, then rows 1, 2, 5 and 6 can be eliminated from the global stiffness matrix. This is possible because when the stiffness matrix is multiplied by a displacement matrix containing those zero values, the respective values for the force components will equal zero. So. ..

$$	\mathbf{D} = \begin{Bmatrix} d_1=0\\d_2=0\\d_3\\d_4\\d_5=0\\d_6=0 \end{Bmatrix}$$

Respective Global Stiffness Matrix
Now. ..

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

Then. . . by the Principle of Virtual Work (PVW), we delete the Global Stiffness Matrix rows corresponding to those rows that equal zero in the Global Stiffness Matrix.

Therefore, the Global Stiffness Matrix, K, reduces to

$$	\mathbf{K} = \begin{vmatrix} k_{33} & k_{34}\\k_{43} & k_{44}\end{vmatrix}$$

Resulting Overall Force Displacement Relationship
$$\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}=\begin{Bmatrix} 0\\P\end{Bmatrix}$$

Now we go back to the 2 bar truss problem, for which the solution is:

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

Computing the Reactions
There are two different methods for finding the reactions:

Method 1

Use the Element Force Displacement Relationship

For this case, the element stiffness matrix is known and the element displacements are known. Some of the element forces are known.

Element 1

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



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

Element 2

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



$$d^{(2)}=\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}$$

=Lecture 11 - Friday, September 19, 2008=

Element 1 Equilibrium
$$\mathbf{f}^{(1)}=\mathbf{k}^{(1)}\mathbf{d}^{(1)}= \begin{bmatrix} 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\\ \end{bmatrix} \begin{bmatrix} 0\\ 0\\ 4.352\\ 6.127\\ \end{bmatrix}= \begin{bmatrix} 0.5625 & 0.3248\\ 0.3248 & 0.1875\\ \end{bmatrix} \begin{bmatrix} 4.325\\ 6.127\\ \end{bmatrix}= \begin{bmatrix} 4.438\\ 2.562\\ \end{bmatrix} $$ $$\mathbf{f}^{(1)}= \begin{bmatrix} f^{(1)}_1\\ f^{(1)}_2\\ f^{(1)}_3\\ f^{(1)}_4\\ \end{bmatrix}= \begin{bmatrix} -4.438\\ -2.562\\ 4.438\\ 2.562\\ \end{bmatrix} $$ Note that:
 * Recall that for Element 1, the local force relationship is given by:
 * To be correct, we must include both the internal and reaction forces of element 1. Since the internal and reaction forces are equal and opposite, the local force matrix for element 1 is given by:
 * $$ f^{(1)}_3 $$,and $$ f^{(1)}_4 $$ are the internal forces of Element 1
 * $$ f^{(1)}_1 $$ and $$ f^{(1)}_2 $$ are reaction forces of Element 1

Observation: Element 1 in Equilibrium We know that Element 1 is in equilibrium because the sum of forces in the x- and y- directions equals zero and the force moment about any point in the element also equals zero. $$\sum F_x = f_1^{(1)} + f_3^{(1)} = -4.438 + 4.438 =0$$ $$\sum F_y = f_2^{(1)} + f_4^{(1)} = -2.562 + 2.562 =0$$ $$\sum M_{any pt} = \sum M_{node1} = -L^{(1)}sin(30)f_3^{(1)} + L^{(1)}cos(30)f_4^{(1)}= -4cos(30)(2.5622) + 4sin(30)4.4378 = -8.867 + 8.867 = 0$$

Element 2 Equilibrium
$$\mathbf{f}^{(2)}=\mathbf{k}^{(2)}\mathbf{d}^{(2)}= \begin{bmatrix} 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\\ \end{bmatrix} \begin{bmatrix} 4.352\\ 6.127\\ 0\\ 0\\ \end{bmatrix}= \begin{bmatrix} 2.5 & -2.5\\ -2.5 & 2.5\\ \end{bmatrix} \begin{bmatrix} 4.352\\ 6.127\\ \end{bmatrix}= \begin{bmatrix} f_1^{(2)}\\ f_2^{(2)}\\ \end{bmatrix}= \begin{bmatrix} -4.438\\ 4.438\\ \end{bmatrix} $$
 * Recall for Element 1 and 2 that the x- and y- nodal displacements are equal at their intersecting node. Simply put, $$ d_3^{(1)} = d_1^{(2)} $$ and $$ d_2^{(2)} = d_4^{(1)}$$. As a result, we can express the force-displacement relationship of Element 2 as:

$$\mathbf{f}^{(2)}= \begin{bmatrix} f^{(2)}_1\\ f^{(2)}_2\\ f^{(2)}_3\\ f^{(2)}_4\\ \end{bmatrix}= \begin{bmatrix} -4.438\\ 4.438\\ 4.438\\ -4.438\\ \end{bmatrix} $$
 * To be correct, we must include both the internal and reaction forces of element 2. Since the internal and reaction forces are equal and opposite, the local force matrix for Element 2 is given by:
 * $$ f^{(2)}_1 $$,and $$ f^{(2)}_2 $$ are the internal forces of Element 2
 * $$ f^{(2)}_3 $$ and $$ f^{(2)}_4 $$ are reaction forces of Element 2

Observation: Element 2 in Equilibrium We know that Element 1 is in equilibrium because the sum of forces in the x- and y- directions equals zero and the force moment about any point in the element also equals zero. $$\sum F_x = f_1^{(2)} + f_3^{(2)} = -4.438 + 4.438 =0$$ $$\sum F_y = f_2^{(2)} + f_2^{(2)} = 4.438 - 4.438 =0$$ $$\sum M_{node2} = -L^{(2)}sin(45)f_2^{(2)} + L^{(2)}cos(45)f_1^{(2)}= -2sin(45)(-4.4378) + -2sin(45)(4.4378) = 6.276 - 6.276 = 0$$

Analysis Methods
In general, there are two separate methods that can be used to solve for a 2-bar truss system using statics:

Method 1: FBD Analysis
In this analysis method:
 * Isolate the supports
 * Assign the applied force to only one FBD
 * Examine FBD 1 and FBD 2 Seaparately

Method 2: Nodal Finite Element Analysis
In this analysis method:
 * Isolate the supports AND the nodes
 * Assign the applied force to its respective node
 * Examine Element 1 and Element 2
 * Examine the node independently

=HW 2 Assignment: Matlab Code=

The following is an M-file created in MATLAB.

function two_bar_truss 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [eps, sigma, force] = PlaneTrussResults(e, A, coord, disps) % [eps, sigma, force] = 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;

The following are results from the above M-file.

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

=Contributing Team Members= The following students contributed to this report


 * Chris Bishop Eml4500.f08.qwiki.bishop 19:19, 25 September 2008 (UTC)
 * Michael Berry Eml4500.f08.qwiki.berry 19:41, 25 September 2008 (UTC)
 * Andrea Booher Eml4500.f08.qwiki.booher 20:21, 25 September 2008 (UTC)
 * David Nobles Eml4500.f08.qwiki.nobles 00:35, 26 September 2008 (UTC)
 * Eric Gatch Eml4500.f08.qwiki.gatch 00:01, 26 September 2008 (UTC)
 * Lee HarrisEml4500.f08.qwiki.harris 18:59, 26 September 2008 (UTC)