User:Eml4500.f08.qwiki.bishop/HW5

 Comment:The change made for this new revision is that there are now labels on the plots. A comparison between the newly updated version of HW5 and our original submission can be found here

Eric GatchEml4500.f08.qwiki.gatch 18:23, 14 November 2008 (UTC)

=Lecture 24 - Monday, October 20, 2008=

=Lecture 25 - Wednesday, October 22, 2008= Test Day, No Lecture Notes

=Lecture 26 - Monday, October 27, 2008=

2-Bar Truss MATLAB Code
The 2-Bar Truss System shown below was previously analyzed via MATLAB and verified by statics calculations in Homework 2. However, Dr. Vu-Quoc pointed out during lecture that there was a bug in the author's code pertaining to the treatment of each element's modulus and cross-sectional area. The steps needed to fix the problem will be shown in this section, along with an explanation of the results. The following properties are given in the MATLAB code.

Debugging the 2-Bar Truss Code
The two-bar truss MATLAB code given by the author of Fundamental Finite Element Analysis and Applications does not take into account the differences in modulus and cross-sectional area when utilizing the  function to create a   array. One can see that the code below will cause  to use the same values of   and   for each element 1 and 2, even though they should not be equivalent.

i=1:elems results = [results; PlaneTrussResults(e, A , ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end

A simple fix to this problem is to replace  and   with   and. Doing so will yield a  array that is calculated from the true values of modulus and cross-sectional area for element 1 and 2.

i=1:elems results = [results; PlaneTrussResults(e(i), A(i) , ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end

The overall MATLAB code that correctly calculates the  array for a 2-bar truss system shown below:

% Two bar truss example provided by Wiley Pushishing % This code calculates axial stress and displacement of each element in a 2-bar truss system % Note: Code was modified on Nov. 4, 2008 to fix bug in calculating results 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(i), A(i), ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g results

After running the corrected code in MATLAB, the following results are obtained:

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 = Columns 1 through 5 0.5625     0.32476      -0.5625     -0.32476            0      0.32476       0.1875     -0.32476      -0.1875            0      -0.5625     -0.32476       3.0625      -2.1752         -2.5     -0.32476      -0.1875      -2.1752       2.6875          2.5            0            0         -2.5          2.5          2.5            0            0          2.5         -2.5         -2.5   Column 6 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       5.1244       0.6276        3.138        6.276

Analyzing MATLAB Results
The last section of code displays the  array for element one and two. Row 1 represents the results for element 1 and row 2 represents element 2. Each row provides values of axial strain, stress and force of its respective element. The values are summarized in the table below.

To verify that the above solutions are correct, it is wise to compare them using statics. Element 1 Thus, the MATLAB solution for element 1 is consistent with static solutions of material properties. Element 2 The MATLAB solutions for element 2 also agree with material property equations relating stress, modulus, strain, and force. Equilibrium of Element 1 and 2 To verify equilibrium, take the sum of the forces in the x- and y-direction $$F_x = (5.1244)(0.866) - (6.276)(0.707) = 4.437 - 4.437 = 0 $$ $$ F_y = 7 - (5.1244)(0.5) - (6.276)(0.707) = 0 $$
 * $$ \sigma_1 = \epsilon_1 E_1 = 1.7081 * 3 = 5.1244$$
 * $$ F_1 = \sigma_1 A_1 = 5.1244 * 1 = 5.1244 $$
 * $$ \sigma_2 = \epsilon_2 E_2 = 0.6276 * 5 = 3.138$$
 * $$ F_2 = \sigma_2 A_2 = 3.138 * 2 = 6.276 $$
 * $$ F_x = F_1 \cos 30^{\circ} - F_2 \cos -45^{\circ} $$
 * $$ F_y = P - F_1 \sin 30^{\circ} - F_2 \sin -45^{\circ} $$

Principle of Virtual Work
To help us understand the Principle of Virtual Work, consider the following situation: Note that $$\alpha_0$$ and $$\alpha_i$$ are weighting factors. By choosing different weighting factors $$(\alpha_0, \alpha_i)$$, etc. you can "filter" out influences. For example, if you let $$\alpha_0 = 0$$, the total grade at the end of the semester will only be based on the sum of each exam grade multiplied by its respective weighting factor, $$\sum^N_{i=1}{\alpha_i G_{Exam,i}}$$. Similarly, when we discuss the Principle of Virtual Work, we will make use of a weighting factor $$\mathbf{w}$$.
 * Suppose your course grade, GTotal, is based on the equation $$G_{Total}=\alpha_0G_{Homework} + \sum^N_{i=1}{\alpha_iG_{Exam, i}}$$

Our overall goal is to derive $$\mathbf{k}_{4x4}^{(e)} = \mathbf{T}_{4x2}^{(e)T} \hat\mathbf{k}_{2x2}^{(e)} \mathbf{T}_{2x4}^{(e)}$$ via the Principle of Virtual Work Recall: $$\hat\mathbf{k}_{2x2}^{(e)} \mathbf{q}_{2x2}^{(e)} = \mathbf{p}_{2x1}^{(e)}$$ $$ \Rightarrow \hat\mathbf{k}^{(e)}_{2x2} \mathbf{q}^{(e)}_{2x1} - \mathbf{p}^{(e)}_{2x1} = \mathbf{0}_{2x1} \qquad \rightarrow \qquad (1) $$ By the Principle of Virtual Work (PVW): $$\mathbf{w}_{2x1} * ( \hat\mathbf{k}^{(e)} \mathbf{q}^{(e)} - \mathbf{p}^{(e)} )_{2x1} = 0_{1x1} \quad $$ for all $$\; \mathbf{w}_{2x1}\qquad \rightarrow \qquad (2) $$ $$\therefore \;$$ We showed that $$ (1) \longleftrightarrow (2)$$ Recall: $$\mathbf{q}{(e)} = \mathbf{T}^{(e)} \mathbf{d}^{(e)} \qquad \rightarrow \qquad (3)$$ Similarly, $$ \; \hat\mathbf{w}_{2x1} = \mathbf{T}_{2x4}^{(e)} \mathbf{w}_{4x1}^{(e)} \qquad \rightarrow \qquad (4)$$
 * $$\mathbf{q}^{(e)} \;$$ is the real displacement matrix defined in the axial coordinate system
 * $$\mathbf{d}^{(e)} \;$$ is the real displacement matrix defined in the global coordinate system
 * $$\hat\mathbf{w}_{2x1} \;$$ is the virtual displacement matrix defined in the axial coordinate system
 * $$\mathbf{w}_{4x1} \;$$ is the virtual displacement matrix defined in the global coordinate system

=Lecture 27 - Wednesday, October 29, 2008=

Read Chapter 2: 2.1, 2.2, 2.6, 2.7

Continuing The Principal of Virtual Work
$$\mathbf{\hat{w}}_{(2x1)}$$ is the virtual displacement, corresponds to $$\mathbf{q}_{(2x1)}^{(e)}$$

$$\mathbf{w}_{(4x1)}$$ is the virtual displacement in the global coordinate system, corresponds to $$\mathbf{d}_{(4x1)}^{(e)}$$

Now, replace Equations (3) and (4) from the previous lectures into Equation (2). Doing so yields:

Equation (5) $$(\mathbf{T}^{(e)}\mathbf{W})[\mathbf{\hat{k}}^{(e)}(\mathbf{T}^{(e)}\mathbf{d}^{(e)})-\mathbf{p}^{(e)}]=0$$

for all $$\mathbf{w}_{(4x1)}$$

Recall the following relation:

Equation (6) $$(\mathbf{A}_{(pxq)}\mathbf{B}_{(qxr)})^{T}=\mathbf{B}^{T}\mathbf{A}^{T}$$

Also recall this following equation:

Equation (7) $$\mathbf{a}_{(nx1)}\mathbf{b}_{(nx1)}=\mathbf{a}_{(1xn)}^{T}\mathbf{b}_{(nx1)}$$

Now apply (7) and (6) into equation (5). Doing so yields the following:


 * $$(\mathbf{T}^{(e)}\mathbf{w})^T[\mathbf{\hat{k}}^{(e)}(\mathbf{T}^{(e)}\mathbf{d}^{(e)})-\mathbf{p}^{(e)}]=0$$

for all $$\mathbf{w}_{(4x1)}$$


 * $$\mathbf{W}^{T}\mathbf{T}^{(e)T}[\mathbf{\hat{k}}^{(e)}(\mathbf{T}^{(e)}\mathbf{d}^{(e)})-\mathbf{p}^{(e)}]=0$$


 * $$\mathbf{W}[(\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)})\mathbf{d}^{(e)}-(\mathbf{T}^{(e)T}\mathbf{p}^{(e)})]=0$$

for all $$\mathbf{w}_{(4x1)}$$


 * $$\mathbf{K}_{(4x4)}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$$


 * $$\mathbf{f}_{(4x1)}^{(e)}=\mathbf{T}^{(e)T}\mathbf{p}^{(e)}$$


 * $$\mathbf{W}[\mathbf{K}^{(e)}\mathbf{d}^{(e)}-\mathbf{f}^{(e)}]=0$$

Continuous Case (PDE's)
So far, only discrete cases have been dealt with (non continuous, matrices). Now, an elastic bar with varying properties with respect to x will be considered. The varying properties are:


 * Cross-sectional area, A(x)
 * Modulus of Elasticity, E(x)
 * Axial (distributed) and concentrated loads, f(x,t)
 * Mass density, $$\rho$$(x)

Note that the axial and concentrated loads are also time dependent.

The following figure illustrates the elastic bar to be examined with the varying properties.



=Lecture 28 - Friday, October 31, 2008=

Free Body Diagram
The following diagram is of an infinitely small element of the composite shown above.

Note: Time is involved because it is a dynamic problem. Also m(x) is the mass per unit length.

Now, using the useful equation $$\Sigma F_{x}=0$$ yields:

$$\Sigma F_{x}=-N(x,t)+N(x+dx,t)+f(x,t)dx-m(x)\ddot{u}=0$$

Expanding $$-N(x,t)+N(x+dx,t)$$ to $$-N(x,t)+N(x,t)+N(dx,t)=N(dx,t)$$

Equation 1

$$0=\frac{\partial N}{\partial x}(x,t)dx+h.o.t.+f(x,t)dx-m(x)\ddot{u}dx$$

Where H.O.T. stands for higher order terms which can be neglected by recalling the Taylor series expansion:

$$f(x+dx)=f(x)+\frac{df(x)}{dx}dx+1/2\frac{d^{2}f(x)}{dx^{2}}dx^{2}+...$$

Therefore Equation 1 can be re-written as Equation 2- The Equation of Motion for an Elastic Bar

$$\frac{\partial N}{\partial x}+f=m\ddot{u}$$

The N force can also be broken down to the following constitutive relationship, naming it Equation 3.



Replacing Equation 3 in Equation 2 yields Equation 4- The Partial Differential Equation of Motion

$$\frac{\partial }{\partial x}[A(x)E(x)\frac{\partial u}{\partial x}]+f(x,t)=m(x)\ddot{u}$$

Boundary Conditions
If there is a 2nd order partial differential equation, 2 conditions are needed. One constant arises after the first integral, and the second comes after the second integral.


 * Need 2 boundary conditions (2nd order derivative with respect to x)
 * Need 2 initial conditions (2nd derivative with respect to t) - initial displacement, initial velocity

Examples


The boundary conditions for the above fixed bar are u(0,t)= 0 = u(L,t).



Left Side Boundary: u(0,t)=0

Right Side Boundary: N(L,t)=F(t) which yields  $$A(L)E(L)\frac{\partial u}{\partial x}(L,t)= F(t)$$

Rearranging results in: $$\frac{\partial u}{\partial x}(L,t)= \frac{F(t)}{A(L)E(L)}$$

=Additional Assignment: Space Trusses=

FD Relation Derivations

 * The concepts of FEA for a plane (2-D) truss can be applied to a space (3-D) truss.
 * The axial FD relation remains the same for each element.
 * $$\begin{Bmatrix}

F_{1}\\ F_{2} \end{Bmatrix}=\begin{bmatrix} k & -k\\ -k & k \end{bmatrix}\begin{Bmatrix} d_{1}\\ d_{2} \end{Bmatrix}$$
 * The elemental degrees of freedom and elemental forces are transformed to global xyz coordinates by a elemental stiffness matrix that now includes a factor "n".
 * $$k^{(e)}\begin{bmatrix}

l^{(e)2} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} & -l^{(e)2} & -l^{(e)}m^{(e)} & -l^{(e)}n^{(e)}\\ l^{(e)}m^{(e)}& m^{(e)2} & n^{(e)}m^{(e)} & -l^{(e)}m^{(e)} & -m^{(e)2} & -n^{(e)}m^{(e)}\\ l^{(e)}n^{(e)} & m^{(e)}n^{(e)} & n^{(e)2} & -l^{(e)}n^{(e)}& -m^{(e)}n^{(e)} & -n^{(e)2}\\ -l^{(e)2} & -l^{(e)}m^{(e)} & -l^{(e)}n^{(e)} & l^{(e)2}& l^{(e)}m^{(e)} & l^{(e)}n^{(e)}\\ -l^{(e)}m^{(e)} & -m^{(e)2} & -n^{(e)}m^{(e)} & l^{(e)}m^{(e)}& m^{(e)2} & n^{(e)}m^{(e)}\\ -l^{(e)}n^{(e)} & - m^{(e)}n^{(e)}& -n^{(e)2} & l^{(e)}n^{(e)}& n^{(e)}m^{(e)} & n^{(e)2} \end{bmatrix}\begin{Bmatrix} d^{(e)}_{1}\\ d^{(e)}_{2}\\ d^{(e)}_{3}\\ d^{(e)}_{4}\\ d^{(e)}_{5}\\ d^{(e)}_{6} \end{Bmatrix}=\begin{Bmatrix} f^{(e)}_{1}\\ f^{(e)}_{2}\\ f^{(e)}_{3}\\ f^{(e)}_{4}\\ f^{(e)}_{5}\\ f^{(e)}_{6} \end{Bmatrix}\Rightarrow \mathbf{K^{(e)}}\mathbf{d^{(e)}}=\mathbf{f^{(e)}}$$
 * The derivation for the Transformation matrix, T (2x6), is similar to that of a plane truss and is as follows:
 * $$q^{(e)}_{1}=(\vec{i}*\tilde{i})d^{(e)}_{1}+(\vec{j}*\tilde{i})d^{(e)}_{2}+(\vec{k}*\tilde{i})d^{(e)}_{3}=l^{(e)}d^{(e)}_{1}+m^{(e)}d^{(e)}_{2}+n^{(e)}d^{(e)}_{3}$$


 * Which applied to the other axial displacement and then put into matrix form is (below is also the relation for element forces:
 * $$\begin{Bmatrix}

q^{(e)}_{1}\\ q^{(e)}_{2} \end{Bmatrix}=\begin{bmatrix} l^{(e)} & m^{(e)} & n^{(e)} & 0 & 0 & 0\\ 0& 0 & 0 & l^{(e)} & m^{(e)} & n^{(e)} \end{bmatrix}\begin{Bmatrix} d^{(e)}_{1}\\ d^{(e)}_{2}\\ d^{(e)}_{3}\\ d^{(e)}_{4}\\ d^{(e)}_{5}\\ d^{(e)}_{6} \end{Bmatrix}=\mathbf{T^{(e)}}\begin{Bmatrix} d^{(e)}_{1}\\ d^{(e)}_{2}\\ d^{(e)}_{3}\\ d^{(e)}_{4}\\ d^{(e)}_{5}\\ d^{(e)}_{6} \end{Bmatrix}$$
 * $$\begin{Bmatrix}

p^{(e)}_{1}\\ p^{(e)}_{2} \end{Bmatrix}=\begin{bmatrix} l^{(e)} & m^{(e)} & n^{(e)} & 0 & 0 & 0\\ 0& 0 & 0 & l^{(e)} & m^{(e)} & n^{(e)} \end{bmatrix}\begin{Bmatrix} f^{(e)}_{1}\\ f^{(e)}_{2}\\ f^{(e)}_{3}\\ f^{(e)}_{4}\\ f^{(e)}_{5}\\ f^{(e)}_{6} \end{Bmatrix}=\mathbf{T^{(e)}}\begin{Bmatrix} f^{(e)}_{1}\\ f^{(e)}_{2}\\ f^{(e)}_{3}\\ f^{(e)}_{4}\\ f^{(e)}_{5}\\ f^{(e)}_{6} \end{Bmatrix}$$ \tilde{f}^{(e)}_{1}\\ \tilde{f}^{(e)}_{2}\\ \tilde{f}^{(e)}_{3} \end{Bmatrix}=\begin{bmatrix} l^{(e)} & m^{(e)} & n^{(e)} \\ -m_{(e)}&l^{(e)} \\ n^{(e)} \end{bmatrix}\begin{Bmatrix} {f}^{(e)}_{1}\\ {f}^{e}_{2}\\ \end{Bmatrix}=\textbf{R}^{(e)}\begin{Bmatrix} {f}^{(e)}_{1}\\ {f}^{(e)}_{2}\\ \end{Bmatrix}$$ \tilde{f}^{(e)}_{1}\\ \tilde{f}^{(e)}_{2}\\ \tilde{f}^{(e)}_{3}\\ \tilde{f}^{(e)}_{4} \end{Bmatrix}=\begin{bmatrix} \mathbf{R}^{(e)}_{2x2} & \mathbf{0}_{2x2}\\ \mathbf{0}_{2x2} & \mathbf{R}^{(e)}_{2x2} \end{bmatrix}\begin{Bmatrix} {f}^{(e)}_{1}\\ {f}^{(e)}_{2}\\ {f}^{(e)}_{3}\\ {f}^{(e)}_{4} \end{Bmatrix}\Rightarrow \tilde{\mathbf{f}}^{(e)}_{4x1}=\tilde{\mathbf{T}}^{(e)}_{4x4}\mathbf{f}^{(e)}_{4x1}$$
 * $$\mathbf{\tilde{d}^{(e)}}=\mathbf{\tilde{T}^{(e)}d^{(e)}}$$
 * Similarly, $$\mathbf{\tilde{f}^{(e)}}=\mathbf{\tilde{T}^{(e)}f^{(e)}}$$
 * $$\begin{Bmatrix}
 * $$\begin{Bmatrix}
 * The set of R and zero matrices can be represented by the T matrix.
 * Applying the Principal of Virtual Work (PVW):
 * $$\mathbf{W}_{12x1}(\mathbf{k}_{12x12}\mathbf{d}_{12x1}-\mathbf{F}_{12x1})=0$$, for all W's.


 * The boundary conditions of this problem shows that:
 * $$d_1=d_2=d_3=d_4=d_5=d_6=d_7=d_8=d_9=0$$


 * Since, the weighting coefficients must be kinematically admissible, they cannot violate the boundary conditions, therefore:
 * $$W_1=W_2=W_3=W_4=W_5=W_6=W_7=W_8=W_9=0$$


 * In plugging these values for d and W into the force-displacement equation, the equation reduces to:
 * $$\bar{\mathbf{K}}_{3x3}\bar{\mathbf{d}}_{3x1}-\bar{\mathbf{F}}_{3x1}=\begin{bmatrix}

K_{10,10} & K_{10,11} & K_{10,12}\\ K_{11,10}& K_{11,11} &K_{10,12} \\ K_{12,12}& K_{12,11} & K_{10,12} \end{bmatrix}\begin{Bmatrix} d_{10}\\ d_{11}\\ d_{12} \end{Bmatrix}-\begin{Bmatrix} F_{10}\\ F_{11}\\ F_{12} \end{Bmatrix}=0$$
 * And typically the forces, F10, F11, F12 are known.
 * If we look at the force-displacement relationship, we can apply the T matrix conversion between global and axial coordiantes.
 * $$\tilde{\mathbf{k}}^{(e)}\tilde{\mathbf{d}}^{(e)}=\tilde{\mathbf{f}}^{(e)}\Rightarrow \tilde{\mathbf{k}}^{(e)}\tilde{\mathbf{T}}^{(e)}\mathbf{d}^{(e)}=\tilde{\mathbf{T}}^{(e)}\mathbf{f}^{(e)}$$
 * But if the T matrix can be inverted, then:
 * $$[\tilde{\mathbf{T}}^{(e) -1}\tilde{\mathbf{k}}^{(e)}\tilde{\mathbf{T}}^{(e)}]\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$$

Results for Matlab Code

 * The global stiffnes matrix, K (12x12)


 * The global displacement matrix, d (12x1)


 * The global reactions, f (9x1), excluding the known forces from given force P


 * The results
 * The conn and lmm arrays
 * This problem is statically indeterminate because there are 3 bars, and therefore too many unknowns for the amount of equations available to solve for them.
 * Here are the plots in plane view down the x axis, down the y axis, down the z axis, and from the point (-2,-2,3).

=Six Bar Truss Homework Problem=

Equivalent Modulus and Cross-Sectional Area for all Elements
Case 1

Plot for Case 1
The graph below shows the undeformed and deformed shaped of the 6-bar truss system for case 1:

Case 2 Plot
=Contributing Team Members= The following students contributed to this report


 * Andrea Booher EML4500.f08.qwiki.booher 17:27, 3 November 2008 (UTC)
 * Mike Berry Eml4500.f08.qwiki.berry 12:32, 4 November 2008 (UTC)
 * Chris Bishop Eml4500.f08.qwiki.bishop 19:56, 4 November 2008 (UTC)
 * David Nobles Eml4500.f08.qwiki.nobles 15:12, 7 November 2008 (UTC)
 * Lee Harris Eml4500.f08.qwiki.harris 01:14, 13 November 2008 (UTC)
 * Eric Gatch Eml4500.f08.qwiki.gatch 21:55, 7 November 2008 (UTC)