User:Eml4500.f08.RAMROD.A/HW5

 This version is very similar to the original but some sections have been rearranged to be placed in the correct order, as there was some confusion as to whose portions went where. Also, a parsing issue was resolved. I, (Eml4500.f08.RAMROD.A 01:11, 14 November 2008 (UTC),) have fixed the HW for final submission and the comparison to the original can be found here: [Comparsion]

Justification of Eliminated rows 1, 2, 5, and 6 to obtain K 2-bar truss:
Force Displacement Relation: K6x6d6x1 = F6x1

Moving F to the other side of the equation gives us Equation (1):

 Equation (1)

$$\mathbf{Kd}-\mathbf{F}=\mathbf{0}_{6x1} $$

Using the Principle of Virtual Work, we can arrive at Equation (2):

 Equation (2)

$$\mathbf{w}_{6x1}\bullet(\mathbf{Kd}-\mathbf{F})=\mathbf{0}_{1x1} $$
 * This is true for all w6x1
 * w is known as the "weighting mat."
 * Note that the dot product is used in this equation.

Equation (1)↔ Equation (2)

Using Principle of Virtual Work
Taking account for the boundary conditions of the 2 bar truss:
 * d1=d2=d5=d6=0

Weighting coefficient must be "kinematically admissible," ie cannot violate the boundary conditions
 * w1=w2=w5=w6=0

Weighting coefficients = virtual displacement


 * (calculus of variations)

From before:

$$\begin{Bmatrix} w_3\\ w_4\end{Bmatrix}\bullet(\mathbf{Kd}-\mathbf{F})=0$$

$$\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F})$$

$$\begin{Bmatrix} w_3\\ w_4 \end{Bmatrix}\bullet[\begin{bmatrix} 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{bmatrix} \begin{Bmatrix} d_3\\ d_4 \end{Bmatrix} - \begin{Bmatrix} F_1\\ F_2\\ F_3\\ F_4\\ F_5\\ F_6 \end{Bmatrix}] = 0$$
 * This is true for all $$\begin{Bmatrix}

w_3\\ w_4\end{Bmatrix}$$

$$\mathbf{\overline{K}} = \begin{bmatrix} K_{33} & K_{34} \\ K_{43} & K_{44} \end{bmatrix}$$

$$\mathbf{\overline{d}} = \begin{Bmatrix} d_3\\ d_4 \end{Bmatrix}$$

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

Comparing Equations (3) and (4)
Equation (3) : $$(kd-F) = 0$$

Equation (4) : $$w(kd-F) = 0$$

Equation (3) → Equation (4) is trivial

Equation (4) → Equation (3)
 * Since Equation (4) is valid for all w, select w=1


 * $$1(kd-F) = 0$$ Thus proving Equation (3) is true

Principal of Virtual Work Matrix
We are going to try to derive element stiffness matrix by Principle of Virtual Work.

$$K^{(e)}=\underline{T}^{(e)^\tau} \hat{\underline{K}}^{(e)}\underline{T}^{(e)}$$

Recalling the force displacement with axial degree of freedom $$\underline{q}^{(e)}$$

$$\hat{\underline{K}}^{(e)}_{2x2}\cdot \underline{q}^{(e)}_{2x1} = \underline{P}^{(e)}_{2x1}$$

$$\hat{\underline{K}}^{(e)}_{2x2}\cdot \underline{q}^{(e)}_{2x1} - \underline{P}^{(e)}_{2x1}= \mathbf{0}_{2x1}$$ Equation $$(1)$$

Principle of Virtual Work applied here: $$\hat{\mathbf{W}}\bullet \left(\hat{\underline{K}}^{(e)}_{2x2}\cdot \underline{q}^{(e)}_{2x1} - \underline{P}^{(e)}_{2x1}\right)= \mathbf{0}_{2x1}$$ Equation $$(2)$$

We showed equation $$(1) \Leftrightarrow (2)$$

Recall $$\mathbf{\underline{q}}^{(e)}= T^{(e)}\mathbf{\underline{d}}^{(e)}$$

Similarly $$\hat{\underline{W}}^{(e)}_{2x1}= \bar{T}^{(e)}_{2x4}\mathbf{\underline{W}}^{(e)}_{4x1}$$

Virtual Displacement Terms
ŵ2x1 = virtual axial displacement, with corresponding q(e)2x1

w4x1 = virtual axial displacement in global coordinates system, with corresponding d(e)4x1

Equation (5)
$$ (T^{(e)}w)\cdot [\hat{K}^{(e)}(T^{(e)}d^{(e)})-P^{(e)}] = 0 $$   for all w4x1

We want to prove equation (5) using matrix identity equation (6):

Equation 6: $$ (\mathbf{AB})^T=\mathbf{B}^T\mathbf{A}^T $$

Homework Example to verify equation (6)
$$ A=\begin{bmatrix}1 & 2 & 3\\4 & 5 & 6 \end{bmatrix}, B=\begin{bmatrix}7 & 8 & 9\\1 & 2 & 3\\4 & 5 & 6\end{bmatrix} , B^T = \begin{bmatrix}7 & 1 & 4\\ 8 & 2 & 5\\ 9 & 3 & 6\end{bmatrix} , A^T=\begin{bmatrix}1 & 4\\ 2 & 5\\ 3 & 6\end{bmatrix} $$

$$ AB=\begin{bmatrix}1\cdot7+2\cdot1+3\cdot4 & 1\cdot8+2\cdot2+3\cdot5 & 1\cdot9+2\cdot3+3\cdot6 \\4\cdot7+5\cdot1+6\cdot4 & 4\cdot8+5\cdot2+6\cdot5 & 4\cdot9+5\cdot3+6\cdot6 \end{bmatrix} =\begin{bmatrix}21 & 27 & 33\\57 & 72 & 87\end{bmatrix} $$

$$(AB)^T = \begin{bmatrix}21 & 57\\27 & 72\\ 33 & 87\end{bmatrix} $$

$$ B^TA^T = \begin{bmatrix}21 & 57\\27 & 72\\ 33 & 87\end{bmatrix} $$

This proves Equation (6): $$ (\mathbf{AB})^T=\mathbf{B}^T\mathbf{A}^T $$

We can recall that $$ \underline{a}\cdot\underline{b} = \underline{a}^T\underline{b} $$, where a and b are nx1 and aT is 1xn

It follows that from Eq(1):

$$ (T^{(e)}w)^T[\hat{k}^{(e)}(T^{(e)}d^{(e)})-P^{(e)}] = 0 $$, for all w,

and using Eq (6) above, that:

$$ w^TT^{(e)T}[\hat{k}^{(e)}(T^{(e)}d^{(e)})-P^{(e)}] = 0 $$

and that:

$$ w\cdot[(T^{(e)T}\hat{k}^{(e)}T^{(e)})d^{(e)}-(T^{(e)T}P^{(e)})] = 0 $$, for all w4x1

Also:

$$ w\cdot[\hat{k}^{(e)}d^{(e)}-f^{(e)}] = 0 $$, for all w.

All the above proves:

$$ k^{(e)}_{4x4} = T^{(e)T}_{4x2}\hat{k}^{(e)}_{4x2}T^{(e)}_{2x4} $$

Principal of Virtual Work
Equation 2:

$$ \hat{w}\cdot(\hat{k}^{(e)}q^{(e)}-P^{(e)}) = 0 $$, for all w

Equation 3:

$$ q^{(e)}_{2x1}=T^{(e)}_{2x4}d^{(e)}_{4x1} $$

and finally equation 4:

$$ \hat{w} = T^{(e)}w $$

Problems of the Continuous Case (PDE's)
The new problems will contain varying parameters

Example of problem:

Elastic bar with varying cross section A(x), varying modulus E(x), subjective to varying axial load (distributed), a concentrated load and inertia force (dynamic).

note: load will be time dependent

Problem Setup


The following is a FBD of an infinitesimal section area of the elastic bar with varying properties.



The bottom force that is not labeled in the picture can be represented as the following equation;

$$\displaystyle A(x) \partial x \rho(x) \frac{\partial^2 u}{\partial t^2}$$

Within this equation $$ A(x) \partial x \rho(x)$$ is equal to $$\displaystyle m(x)$$ which has the units of $$ \frac{mass}{length}$$. Now we proceed to sum the forces in the x-direction and the results are as follows;

$$\displaystyle \sum{F_x}= -N(x,t) + N(x+dx,t) + f(x,t)dx - m(x)\ddot{u} = 0$$

where $$\ddot{u}$$ is equal to the acceleration as given in Newton's notation. Now we perform a Taylor Series expansion on this equation but we will be ignoring the higher order terms. This is performed below.

$$\displaystyle \sum{F_x}= \frac{\partial N}{\partial x}(x,t) + (Higher Order Terms ) + f(x,t) - m(x)\ddot{u}$$  $$\displaystyle f(x+dx) = f(x) + \frac{df}{dx}(x)dx + \frac{1}{2}\frac{d^2f}{dx^2} +.....$$  $$\displaystyle \frac{\partial N}{\partial x} + f = m\ddot{u}$$

The second equation is part of the expansion while the last equation is known as the equation of motion for the elastic bar which ignores the higher order terms of the expansion. The equation of motion is otherwise abbreviated as $$\displaystyle EOM$$. An expression for the force $$N(x,t)$$ can be obtained and the results are found below.

$$\displaystyle N(x,t) = A(x)\sigma(x) = A(x) E(x) \varepsilon (x,t) = A(x) E(x) \frac{\partial u}{\partial x}(x,t)$$

Putting this expression into the $$EOM$$ yields the partial differential equation of motion (PDE of Motion) for the elastic bar.

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

In order to solve this equation two boundary conditions along with two initial conditions are needed. The two initial conditions are initial displacement and initial velocity. For the boundary conditions, two cases need to be considered. They are depicted below.



For case a, the two boundary conditions are as follows;

$$\displaystyle u(0,t) = u(L,t) = 0$$

$$\displaystyle \frac{\partial u}{\partial x}(L,t) = \frac{F(t)}{A(L)E(L)}$$

For case two the boundary conditions are as follows;

$$\displaystyle u(0,t) = o$$

$$\displaystyle N(L,t) = F(t)$$

The initial conditions are specific to different problems

Debugged 2 Bar Truss System
Below is the code to the debugged two bar truss system. Since the modulus of elasticity and the areas for the given example were not the same the code had to be modified in order for it to run properly. Before the correction the code was running as if the all of the elements had the same modulus of elasticity and the same area. Below is the corrected code and the results from the corrected code. Functions PlaneTrussElemtens.m, NodalSoln.m and PlaneTrussResults.m were used in this code and the discription for these functions can be found in the functions section.

Results from the debugged two bar truss code

Comparison of Results
Once the code was corrected the same reaction results were obtained when comparing it to the statics methods that have been done in class. Below are the results to the statics method.

Finite Element Method:

Using this equation, we can now solve for the f(1) matrix and thus we find the values for the reactions: $$f^{(1)}_{1} = F_{1} = -4.4378$$ $$ f^{(1)}_{2} = F_{2} = -2.5622$$

Using this same procedure applied to element 2, we can solve for the values of the other two unknown reactions giving us: $$f^{(2)}_{3} = F_{5} = 4.4378$$ $$f^{(2)}_{4} = F_{6} = -4.4378$$

Statics

Using the Euler cut method, you find the following components from the resultant P forces.

$$P_x^1 = -4.4378$$ $$P_y^1=-2.5622$$ $$P_x^2=4.4378$$ $$P_y^2=-4.4378$$

Six Bar Truss Code
Below is the code for the six bar truss system found on page 226 in the the textbook. All six of the bars have the same modulus of elasticity (e) and the same area (A). The MATLAB functions of PlaneTrussElement.m, NodalSoln.m, and PlaneTrussResults.m are all found in the code. The descriptions and functions can be found below the results of the six bar truss system.

Results of the above code.

Functions
The functions PlaneTrussElement.m, NodalSoln.m, PlaneTrussResults.m were all called in the code for the 5 truss system. They are all necessary for the completion of the truss problem.

PlaneTrussElement.m

This function take the Young's Modulus (e), the Area of the element (A) and the coordinates of the element ends (coord) and generates the stiffness matrix for a plane truss element. The function first calculates the lengths of the the elements. Once those are calculated the stiffness matrix is calculated.

NodalSoln.m

The function takes the global coefficient matrix (K), the global right hand side vector (R), the list of degrees of freedom with specified values, and the specified values to determine the displacements and reactions at each node. The dof of the the system is first found by using the  command which finds the longest dimension of R. df then finds the difference between the dofs that have known values (a value of zero) and the dof that were found in the previous line. The displacements and the reactions are then calculated.

PlaneTrussResults.m

This function computes the plane truss element results. It takes in the modulus of elasticity (e), the area of the cross-section (A), the coordinates at the element ends (coord), and the displacements at the elemtent ends (disps) and calculates the axial strain (eps), stress(sigma) and force (force). The results are stored in the variable  and sent back to the main program.

Plot of Six Bar Truss System
Below is the code and the plot of the six bar truss system. The plot was made from the results that were found in the above code. The results demonstrate that only node 2 and node 5 had a displacement of (0.2131, 0.2500) and (-0.0061, 0.0122) respectively. The undeformed system is plotted with the dotted lines and the deformed system is plotted with solid lines. The deformed displacements have been scaled up by a factor of 1.2 in order to enhance in the plot.

Six Bar Truss Code for Elements with Varying E
Next, the six bar truss Matlab code used above, was modified to allow for different values of E for each element. This code uses all the same functions as above (i.e. PlaneTrussElement.m, PlaneTrussResults.m and NodalSoln.m). The modified code is shown below. The results from this situation, as well as the constant E case and the undeformed truss system are shown in the plot below. The varying E case is shown by the dotted nodes while the undeformed case is shown dashed. The element and node numbers are suppressed to avoid cluttering the diagram. The code used to create this plot is very similar to that shown above for the six bar truss with constant E, but repeated for the two deformation cases.

Three Bar Space Truss Example
Below is the code for the the 3D three bar truss system found in the textbook on page 230. There are four nodes in the system with each node containing 3 degrees of freedom. The code uses three additional functions for it to be run. SpaceTrussElement.m, NodalSoln.m, and NodalSoln.m were all used for the code. Below this code a description of each function can be found.

Functions
The functions SpaceTrussElement.m, NodalSoln.m, SpaceTrussResults.m were all called in the code for the 5 truss system. They are all necessary for the completion of the truss problem.

SpaceTrussElement.m This function creates the stiffness matrix of a space truss element. It takes in the modulus of elacticity (e), the truss area (A), and the corresponding coordinates for the truss. The function first calculates the length of each element and then calculates the stiffness matrix (k).

NodalSoln.m

The function takes the global coefficient matrix (K), the global right hand side vector (R), the list of degrees of freedom with specified values, and the specified values to determine the displacements and reactions at each node. The dof of the the system is first found by using the  command which finds the longest dimension of R. df then finds the difference between the dofs that have known values (a value of zero) and the dof that were found in the previous line. The displacements and the reactions are then calculated.

SpaceTrussResults.m

This function takes in the modulus of elasticity (e), the area of each truss (A), the coordinate of the nodes and the displacements (disps) at the element ends. It calculates the axial strain (eps), stress(sigma) and force (force). The results are stored in the variable  and sent back to the main program.

3 Bar Space Truss Results
Below are the results from the three bar space truss example.

Code for the lot of the 3 Bar Space Truss
Below is the code for the 3 Bar Space Truss.

Plot of 3 Bar Space Truss
Below is the plot of the 3 bar Space truss with the viewing mode of x = 1m, y = 1m, and z= 1m. The dotted lines represent the undeformed truss system and the solid lines represent the deformed truss system.

]

Below is the same plot but in plane view down the x axis, plane view down the y axis, and plane view down the z axis.

=3-D Truss Relationships=

This process is very similar to that of the 2-D case. Therefore we will be taking the same steps to derive the Force-Displacement relationships. The elemental axial force displacement relationship is identical to the 2-D case and is shown below.

$$\displaystyle (k^e)(q^e)=p^e$$

Since the 3-D case will have three displacements at each node, the element displacement matrix will be a six by one (6x1) as will the element force matrix. These matrices are shown below.

 $$\displaystyle d^e = \begin{Bmatrix}d_1^e\\ d_2^e\\ d_3^e\\ d_4^e\\ d_5^e\\ d_6^e \end{Bmatrix}$$ <p style="text-align:center;">$$\displaystyle f^e = \begin{Bmatrix}f_1^e\\ f_2^e\\ f_3^e\\ f_4^e\\ f_5^e\\ f_6^e \end{Bmatrix}$$

For the 3-D case the two by six (2x6) transformation matrix $$T^e$$ that relates the axial degrees of freedom to the global degrees of freedom is shown below.

<p style="text-align:center;"> $$\displaystyle T^e = \begin{bmatrix} l^e & m^e & n^e & 0 & 0 & 0 \\ 0 & 0 & 0 & l^e & m^e & n^e \end{bmatrix}$$

This transformation matrix can be used to relate the axial and global degrees of freedom and forces as shown below.

<p style="text-align:center;">$$\displaystyle \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_1\\ d_2\\ d_3\\ d_4\\ d_5\\ d_6\end{Bmatrix}=\begin{Bmatrix}q_1\\ q_2\end{Bmatrix}$$

It can also be used to relate the global forces to the local axial forces. This is shown below.

<p style="text-align:center;">$$\displaystyle \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_1\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6\end{Bmatrix}=\begin{Bmatrix}p_1\\ p_2\end{Bmatrix}$$

Using the Principle of Virtual Work, the global FD relationship is given by:

<p style="text-align:center;">$$\displaystyle K*(d)=F$$

Also, we can represent the elemental matrix $$k^e$$ in terms of the axial matrix $$\hat{k}$$ and the transformation matrix, $$\displaystyle T$$. This is shown below:

<p style="text-align:center;">$$k^e = \hat{k^e}(T)(T^T)$$

Lastly, the relationship between the element forces in global coordinates, $$\displaystyle f^e$$, and the axial forces, $$\displaystyle p^e$$, can be found using the transformation matrix, $$\displaystyle T$$.

<p style="text-align:center;">$$\displaystyle f^e(T)=p^e$$

Or more detailed:

<p style="text-align:center;">$$\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_1\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6\end{Bmatrix}=\begin{Bmatrix}p_1\\ p_2\end{Bmatrix}$$

Using the above derived equation for $$\displaystyle k^e$$, we can give the actual values for its elements. Multiplying out the right side of the equation gives us the following:

<p style="text-align:center;">$$k^e = \hat{k^e}(T)(T^T)$$

<p style="text-align:center;">$$k^e = \begin{bmatrix} 1 &-1 & 1\\ -1 & 1 & -1\\ 1&  -1&1 \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} l^e & 0\\ m^e & 0\\ n^e & 0\\ 0 & l^e\\ 0 & m^e\\ 0 & n^e \end{bmatrix}$$

This is the elemental stiffness matrix, $$\displaystyle k^e$$ for element e.



Contributing Members
Eml4500.f08.RAMROD.A 21:11, 7 November 2008 (UTC)

EML4500.F08.RAMROD.B 19:54, 7 November 2008 (UTC)

EML4500.f08.RAMROD.E 18:05, 6 November 2008 (UTC)

Eml4500.f08.ramrod.D 01:59, 7 November 2008 (UTC)

Eml4500.f08.ramrod.c 19:11, 7 November 2008 (UTC)

Eml4500.f08.RAMROD.F 20:36, 7 November 2008 (UTC)