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

=Homework 5=

Justification of Reduction of Stiffness Matrix for 2-bar Truss Problem
From the general force displacement relationship, the following relation can be derived:


 * $$\mathbf{kd}= \mathbf{F}$$ => $$\mathbf{kd - F}= \mathbf{0}\;\;(1)$$

Using this relation, the Principle of Virtual Work can be expressed as follows:


 * $$\mathbf{W}\cdot\mathbf{(kd - F)} = \mathbf{0}$$ for all $$\mathbf{W}\;\;(2)$$

Where $$\mathbf{W}\;$$ refers to the weighting matrix. It is useful to show the equivalency of the above two expressions. While is is trivial to show how equation 2 comes from equation 1, it is more challenging to show how equation 1 comes from equation 2. In order to prove that this is the case, consider the following weighting coefficients:


 * $$W_{1} = 1\;$$, $$\;W_{2}=...=W_{6} = 0$$

Which can be expressed in transposed form as:


 * $$\mathbf{W^{T}} = \begin{bmatrix} 1&0&0&0&0&0\end{bmatrix}$$

Applying this choice for $$\mathbf{W}$$ to equation 2 yields the following result:


 * $$\mathbf{W}\cdot\mathbf{(kd - F)} = [1\cdot\sum_{j=1}^{6} (k_{1j}d_{j} - F_{1}) + 0\cdot\sum_{j=1}^{6} (k_{2j}d_{j} - F_{2}) + ... + 0\cdot\sum_{j=1}^{6} (k_{6j}d_{j} - F_{6})] = 0$$

It is easily shown that the result of the above expression is:


 * $$\sum_{j=1}^{6} k_{1j}d_{j} = F_{1}$$

Consider another set of weighting coefficients and resulting transposed form:


 * $$W_{1} = 0\;$$, $$W_{2} = 0\;$$ , $$\;W_{3}=...=W_{6} = 0$$


 * $$\mathbf{W^{T}} = \begin{bmatrix} 0&1&0&0&0&0\end{bmatrix}$$

As before, applying this matrix to equation 2 yields a result of $$\sum_{j=1}^{6} k_{2j}d_{j} - F_{2} = 0$$. Thus, it is apparent that for any choice of weighting matrix $$\mathbf{W}$$, the result is equivalent to $$\mathbf{kd - F}= \mathbf{0}$$, which is equation 1.

Using this method, the rest of the force-displacement relations can be found.

Choice 3:


 * $$W_1=0\;, W_2=0, W_3=1, W_4=0, W_5=0, W_6=0$$


 * $$\mathbf{W^T}=\begin{bmatrix} 0&0&1&0&0&0 \end{bmatrix}$$


 * $$\sum_{j=1}^{6} k_{3j}d_{j} - F_{3} = 0$$


 * $$\sum_{j=1}^{6} k_{3j}d_{j} = F_{3}$$

Choice 4:


 * $$W_1=0\;, W_2=0, W_3=0, W_4=1, W_5=0, W_6=0$$


 * $$\mathbf{W^T}=\begin{bmatrix} 0&0&0&1&0&0 \end{bmatrix}$$


 * $$\sum_{j=1}^{6} k_{4j}d_{j} - F_{4} = 0$$


 * $$\sum_{j=1}^{6} k_{4j}d_{j} = F_{4}$$

Choice 5:


 * $$W_1=0\;, W_2=0, W_3=0, W_4=0, W_5=1, W_6=0$$


 * $$\mathbf{W^T}=\begin{bmatrix} 0&0&0&0&1&0 \end{bmatrix}$$


 * $$\sum_{j=1}^{6} k_{5j}d_{j} - F_{5} = 0$$


 * $$\sum_{j=1}^{6} k_{5j}d_{j} = F_{5}$$

Choice 6:


 * $$W_1=0\;, W_2=0, W_3=0, W_4=0, W_5=0, W_6=1$$


 * $$\mathbf{W^T}=\begin{bmatrix} 0&0&0&0&0&1 \end{bmatrix}$$


 * $$\sum_{j=1}^{6} k_{6j}d_{j} - F_{6} = 0$$


 * $$\sum_{j=1}^{6} k_{6j}d_{j} = F_{6}$$

Returning to the specific 2-bar truss problem, accounting for the boundary conditions allows the removal of the following global displacements:


 * $$d_{1}=d_{2}=d_{5}=d_{6}=0\;$$

The weighting coefficients must be "kinematically admissable" - this means that it must not violate the boundary conditions. Thus, the following is a result of the above expression:


 * $$W_{1}=W_{2}=W_{5}=W_{6}=0\;$$

This leads to the conclusion that the weighting coefficients are essentially virtual displacements. As a result of the previous 2 relations, equation 2 can now be rewritten as:


 * $$\begin{Bmatrix}W_{3}\\W_{4}\end{Bmatrix}\cdot\mathbf{(\bar{k}\bar{d} - \bar{F})} = 0$$ for all $$\begin{Bmatrix}W_{3}\\W_{4}\end{Bmatrix}$$

Where $$\mathbf{\bar{k}}\;$$, $$\mathbf{\bar{d}}\;$$, and $$\mathbf{\bar{F}}$$ are defined as follows:


 * $$\mathbf{\bar{k}} = \begin{bmatrix}k_{33}&k_{34}\\k_{43}&k_{44}\end{bmatrix}$$


 * $$\mathbf{\bar{d}} = \begin{Bmatrix}d_{3}\\d_{4}\end{Bmatrix}$$


 * $$\mathbf{\bar{F}} = \begin{Bmatrix}F_{3}\\F_{4}\end{Bmatrix}$$

An example of weighting is the calculation of a student's final grade in this course, which is given as follows: $$course\;grade = \alpha_{0} + HW Grade + \Sigma(\alpha Exam Grade)$$

Prove Kd=F is equivalent to W(kd-F)=0
$$\mathbf{W}\cdot(\mathbf{K}\mathbf{d}-\mathbf{F})=\begin{Bmatrix} w_3\\ w_4\\ \end{Bmatrix}\cdot(\mathbf{\bar{K}}\mathbf{\bar{d}}-\mathbf{\bar{F}})=0 $$ for all $$ \begin{Bmatrix} w_3\\ w_4\\ \end{Bmatrix}$$

Remember and apply our previous derivations for K, d, and F:

k_{33} & k_{34}\\ k_{43} & k_{44}\\ \end{bmatrix}$$ d_3\\ d_4\\ \end{bmatrix}$$ F_3\\ F_4\\ \end{bmatrix}$$
 * $$\mathbf{K}=\begin{bmatrix}
 * $$\mathbf{d}=\begin{bmatrix}
 * $$\mathbf{F}=\begin{bmatrix}

$$\begin{Bmatrix} w_3\\ w_4\\ \end{Bmatrix}\cdot(\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})=0$$

This matrix can be converted by the dot product and simple matrix math to a single algebraic expression:

$$w_3(k_{33}d_3+k_{34}d_4-F_3)+w_4(k_{43}d_3+k_{44}d_4-F_4)=0$$

This equation will have multiple solutions depending on $$w_3$$ and $$w_4$$. Since by definition the relationships apply for any $$w_3$$ and $$w_4$$ then any arbitrary value can be applied to solve the expression created above. Picking different sets of values will create a system of equations that, when written in matrix form, are identical to $$\mathbf{W}\cdot(\mathbf{K}\mathbf{d}-\mathbf{F})=0$$ and other values will also prove identity.

If we pick 0 for our values of w:

$$w_3=0, w_4=0$$

$$(0)(k_{33}d_3+k_{34}d_4-F_3+(0)(k_{43}d_3+k_{44}d_4-F_4)=0$$

$$0=0$$

If we pick 0 for one value of w and 1 for the other:

$$w_3=0, w_4=1$$

$$(0)(k_{33}d_3+k_{34}d_4-F_3+(1)(k_{43}d_3+k_{44}d_4-F_4)=0$$

$$k_{43}d_3+k_{44}d_4-F_4=0$$

$$F_4=k_{43}d_3 + k_{44}d_4$$

If we pick 1 for our values of w:

$$w_3=1, w_4=1$$

$$(1)(k_{33}d_3+k_{34}d_4-F_3+(1)(k_{43}d_3+k_{44}d_4-F_4)=0$$

$$k_{33}d_3+k_{34}d_4+k_{43}d_3+k_{44}d_4=F_3+F_4$$

$$\begin{bmatrix} k_{33}d_3 + k_{34}d_4\\ k_{43}d_3 + k_{44}d_4\\ \end{bmatrix}=\begin{bmatrix} F_3\\ F_4\\ \end{bmatrix}$$

These three examples each bring us back to our proof that:

$$\mathbf{Kd}=\mathbf{F}$$

Derivation of $$\mathbf{k^{(e)}} = \mathbf{T^{(e)^{T}}}\mathbf{\hat{k}^{(e)}}\mathbf{T^{(e)}}$$
Recall the force-displacement relationship adjusted for axial degrees of freedom $$\mathbf{q^{(e)}}$$:


 * $$\mathbf{\hat{k}^{(e)}}\mathbf{q^{(e)}} = \mathbf{P^{(e)}}$$

This relationship can be rearranged to the following:


 * $$\mathbf{\hat{k}^{(e)}}\mathbf{q^{(e)}} - \mathbf{P^{(e)}} = \mathbf{0}\;\;(1)$$

Applying the Principle of Virtual Work results in the following:


 * $$\mathbf{\hat{W}}\cdot(\mathbf{\hat{k}^{(e)}}\mathbf{q^{(e)}} - \mathbf{P^{(e)}}) = \mathbf{0}$$ for all $$\mathbf{\hat{W}}$$

It has previously been shown that equations 1 and 2 are equivalent. Recalling the following two relations:


 * $$\mathbf{q^{(e)}} = \mathbf{T^{(e)}}\mathbf{d^{(e)}}\;\;(3)$$


 * $$\mathbf{\hat{W}^{(e)}} = \mathbf{T^{(e)}}\mathbf{W^{(e)}}\;\;(4)$$

$$\mathbf{\hat{W}}$$ refers to virtual axis displacements, corresponding to $$\mathbf{q^{(e)}}$$.

$$\mathbf{W}$$ refers to virtual displacements in global coordinate system, corresponding to $$\mathbf{d^{(e)}}$$.

Placing equations 3 and 4 into equation 2 yields the following:


 * $$(\mathbf{T^{(e)}}\mathbf{W})\cdot[\mathbf{\hat{k}^{(e)}}(\mathbf{T^{(e)}}\mathbf{d^{(e)}}) - \mathbf{P^{(e)}}] = \mathbf{0}$$ for all $$\mathbf{W}\;\;(5)$$

Recall the following two vector identities:


 * $$(\mathbf{A}\cdot\mathbf{B})^{T} = \mathbf{B}^{T}\mathbf{A}^{T}\;\;(6)$$

Let matrix A and matrix B be as follows:

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

Then matrix $$AB$$ and matrix $$AB^T$$ are as follows:

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

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

The matrix $$A^T, B^T and A^TB^T$$ are as follows:

$$A^T = \begin{bmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{bmatrix}$$

$$B^T = \begin{bmatrix} 7 & 1 & 4 \\ 8 & 2 & 5 \\ 9 & 3 & 6 \end{bmatrix}$$

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

Therefore, $$A^TB^T = (AB)^T$$


 * $$\mathbf{a}\cdot\mathbf{b} = \mathbf{a}^{T}\mathbf{b}\;\;(7)$$

Applying equations 7 and 8 into equation 5 yields the following derivation:

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

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

$$\mathbf{W}^{T}[\left(\mathbf{T^{(e)}}^{T}\mathbf{\hat{k}^{(e)}}\mathbf{T^{(e)}}\right)\mathbf{d^{(e)}} - \left\{\mathbf{T^{(e)}}^{T}\mathbf{P^{(e)}}\right\}] = \mathbf{0}$$

In the last expression, the term inside the parentheses is equivalent to $$\mathbf{k^{(e)}}$$, and the term inside the curly brackets is equivalent to $$\mathbf{f^{(e)}}$$. Making these substitutions yields the following result:

$$\mathbf{W^{T}}[\mathbf{k^{(e)}}\mathbf{d^{(e)}} - \mathbf{f^{(e)}}] = \mathbf{0}$$ for all $$\mathbf{W}$$.

Because the previous equation has been proven for all values of $$W$$ we can choose $$W$$ to be one and hence the equation becomes:

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

Continuous Case - Motivational Model Problem


Consider the elastic bar shown in Figure 1. The bar has a varying cross-section and Young's modulus, and is subject to a time-dependent loading condition consisting of a distributed, varying axial load as well as a concentrated load, signified by $$f(x,t)\;$$ and $$P(x)\;$$ respectively. A section of the bar is shown in Figure 2.

The product $$A(x)\rho (x)\;$$ can be expressed as $$m(x)\cdot dx\;$$, where $$m(x)\;$$ is the mass per unit length. The term $$\frac{\partial{^{2}u}}{\partial{x^{2}}}$$ refers to the acceleration of the deforming bar, and can be represented as $$\ddot{u}$$.

The sum of forces on the section can be written as follows:


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


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

"h.o.t" refers to the higher order terms.

Recall the general definition of the Taylor Series expansion, $$f(x+dx) = f(x) + \frac{\partial{f(x)}}{\partial x}dx + \frac{1}{2}\frac{\partial{^{2}x}}{\partial{x}^{2}} + ...$$, where $$\frac{1}{2}\frac{\partial{^{2}x}}{\partial{x}^{2}} + ...$$ are the higher order terms and are typically neglected. Dividing out the $$dx\;$$ terms and applying this concept to equation 1 yields equation 2, the equation of motion for the elastic bar:


 * $$\frac{\partial{N}}{\partial x} + f = m\ddot{u}\;\;(2)$$

The constitutive relation for the bar can be given as equation 3:


 * $$N(x,t) = A(x)\sigma (x,t)\;\;(3)$$

The $$\sigma (x,t)\;$$ term can be broken down as follows:


 * $$\sigma (x,t)= E(x)\epsilon (x,t)\;$$


 * $$\epsilon (x,t)= \frac{\partial{u}}{\partial{x}}(x,t)\;$$

The partial differential equation of motion for the elastic bar can be given as equation 4, which is the result of substituting equation 3 and resulting breakdown of terms into equation 2:


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

In order to manipulate this equation, boundary and initial conditions are required. Since the equation is of 2nd order in both x and t, two boundary and two initial conditions are necessary. The initial conditions consist of an initial displacement and velocity, given as $$[init. vel., init. disp]\;$$. The boundary conditions depend on the type of system being analyzed. Figure 3 shows a bar of length L with both ends constrained. The boundary conditions for this system are:


 * $$u(0,t) = 0 = u(L,t)\;$$

Figure 4 again shows a bar of length L, but with one end free and subject to force F. The first boundary condition, as with the constrained case, is:


 * $$u(0,t) = 0\;\;(bc\;1)$$

The second boundary condition can be derived using the following line of thought, similar to that used to derive equation 4:


 * $$N(L,t) = F(t)\;$$


 * $$N(L,t) = A(L)\sigma (L,t)\;$$


 * $$\sigma (L,t) = E(L)\epsilon (L,t)\;$$


 * $$\epsilon (L,t) = \frac{\partial u}{\partial x} (L,t)$$

Using the above relations, the second boundary condition for this case can be written as follows:


 * $$\frac{\partial u(L,t)}{\partial x} = \frac{F(t)}{A(L)E(L)}\;\;(bc\;2)$$

Composite Materials Research
Composite materials are materials made from at least two separate materials. The materials used in a composite material are of differing physical/chemical properties so that the final composite material has several material properties. Typical composites include carbon fiber, plywood, fiberglass, concrete, and some types of Kevlar reinforcement. Composites exist as a matrix element and a reinforcement element, both of which are important factors in composite creation. The reinforcement material combines with the matrix materials to create a finite combination of two or more materials, which creates the special composite material.

Composites are arranged so that any kind of loading is applied parallel to their fiber arrangement axis. They are not isotropic materials, so any loading that is applied to perpendicular to the fiber arrangement axis could result in failure of the composite material. If any of the material of the composite becomes damaged or separated from the matrix, imminent failure can occur.

An excellent source for composite material information is on Wikipedia at the following link: http://en.wikipedia.org/wiki/Composite_materials

MATLAB Code
The first part of the required MATLAB code debugs part of the two-bar truss system from Homework 2. The code is presented below for the corrected system representation and is called "TwoBarTrussSystemDebugged.m". It uses the same three functions (PlaneTrussElement.m, NodalSoln.m, PlaneTrussResults.m) defined in Homework 2. The only difference between "TwoBarTrussSystem.m" and "TwoBarTrussSystemDebugged.m" is Line 40 of code reproduced below the program.

The first part of the program defines the specifications of each bar and node including modulus of elasticity (e), cross-sectional area (A), and loading values (P). The "nodes" matrix defines the absolute locations of the three nodes. The "conn" matrix defines the global node numbers at the ends of each bar or element. Each row represents an element. As shown in the code, there are 2 sets of 2 numbers for the 2 nodes corresponding to each element. The "lmm" matrix represents the global degrees of freedom (dof) for each element. Each node has 2 dof's; therefore, each element has 4 corresponding dof's. Degrees of freedom are numbered from 1 to 6 where the first 2 correspond to the x and y direction of Global Node 1, the second 2 correspond to the x and y direction of Global Node 2, and the third 2 correspond to the x and y direction of Global Node 3. From all this data, the global stiffness matrix (K) is created. The displacements of the nodes are calculated by isolating the relevant parts of the global stiffness matrix, inverting them, and multiplying by the load array (R). Finally, the reactions (reactions) at each node are calculated by multiplying the displacement array (d) by the global stiffness matrix.

TwoBarTrussSystemDebugged.m
Line 40

Before

After

Command Window Output

>> TwoBarTrussSystemDebugged

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       5.1244       0.6276        3.138        6.276

The displacement array (d) shows that Global Node 2 was displaced because it was the only unfixed node. The output also displays the reactions at the fixed nodes in the "reactions" matrix. The results matrix (results) shows the strain, stress, and axial force in each element. Homework 2, the results matrix changed from a 3x5 to a 3x3 matrix, showing that the bug was removed. Below are the tabular results of the output above. These values match the values computed using the statics method from Homework 2.

The next part of the MATLAB code constructs and loads a truss system composed of six bars. The code below represents Example 4.1 (SixBarTrussEx.m) on page 226 of the textbook. The code for this program is structured similarly to the code for the two-bar truss system. However, there are six elements and five nodes instead of two elements and three nodes.

There are two cases for this example. The original system (SixBarTrussEx.m), Case 1, is the code provided by textbook. Case 2 (SixBarTrussExModified.m) is the same system with a different modulus of elasticity for each element.

SixBarTrussEx.m (Case 1)
Command Window Output

>> SixBarTrussEx

K =

Columns 1 through 5

85355       35355       -50000            0            0        35355        35355            0            0            0       -50000            0        85355       -35355            0            0            0       -35355  1.0202e+005            0 0           0            0            0        71554            0            0            0            0       -35777            0            0            0            0            0            0            0            0       -66667            0       -35355       -35355       -35355        35355       -71554       -35355       -35355        35355       -35355        35777

Columns 6 through 10

0           0            0       -35355       -35355            0            0            0       -35355       -35355            0            0            0       -35355        35355            0            0       -66667        35355       -35355       -35777            0            0       -71554        35777        17889            0            0        35777       -17889            0        71554        35777       -71554       -35777            0        35777        84555       -35777       -17889        35777       -71554       -35777  2.1382e+005            0 -17889      -35777       -17889            0  1.0649e+005

R =

0           0        10000        17321            0            0            0            0            0            0

d =

0           0      0.21311      0.24998            0            0            0            0   -0.0060971     0.012242

reactions =

-10873     -217.27       874.27      -437.13      -1.7279       -16666

results =

5.3276e-005      10.655        10655 -4.6334e-006    -0.92669      -926.69 -4.8873e-006    -0.97746      -977.46 -8.3326e-005     -16.665       -16665 1.5363e-006     0.30727       307.27 -9.659e-009  -0.0019318      -1.9318

The displacement array (d) shows that Global Nodes 2 and 5 were displaced because these were the only unfixed nodes. The output also displays the reactions at the fixed nodes in the "reactions" matrix. The results matrix (results) show the strain, stress, and axial force in each element. Below are the tabular results of the output above.

To visualize the loaded and unloaded system, the following code graphs the system in its loaded and unloaded state. In order to be able to realize the direction of node movement, the unfixed node displacements have been magnified by the displacement magnification factor (m) of 3000. As shown, the original values from the displacement array (d(3), d(4), d(9), d(10)) above were used below to get the right proportions of the nodal displacements.

SixBarTrussDisplay.m
Command Window Output

>> SixBarTrussDisplay



The dashed red line represents the undeformed system shape. The solid blue line represents the deformed system shape. The abbreviations "GN" and "E" stand for Global Node and Element, respectively. Note only nodes 2 and 5 are displaced because they are unfixed.

The code below is the code for Case 2. Below the program are the line changes to the program for Case 1.

SixBarTrussExModified.m (Case 2)
Lines 2, 22, and 31

Before

After

Command Window Output

>> SixBarTrussExModified

K =

Columns 1 through 5

76391       38891       -37500            0            0        38891        38891            0            0            0       -37500            0        69320       -31820            0            0            0       -31820        98486            0            0            0            0            0        71554            0            0            0            0       -35777            0            0            0            0            0            0            0            0       -66667            0       -38891       -38891       -31820        31820       -71554       -38891       -38891        31820       -31820        35777

Columns 6 through 10

0           0            0       -38891       -38891            0            0            0       -38891       -38891            0            0            0       -31820        31820            0            0       -66667        31820       -31820       -35777            0            0       -71554        35777        17889            0            0        35777       -17889            0        89443        44721       -89443       -44721            0        44721        89027       -44721       -22361        35777       -89443       -44721  2.3171e+005        16015 -17889      -44721       -22361        16015  1.1096e+005

R =

0           0        10000        17321            0            0            0            0            0            0

d =

0           0      0.26485      0.26083            0            0            0            0   0.00063864    -0.001246

reactions =

-9908.3      23.619      -90.274       45.137      -1.4008       -17389

results =

6.6213e-005      9.9319       9931.9 5.347e-007    0.096245       96.245 5.0465e-007     0.10093       100.93 -8.6943e-005     -17.389       -17389 -1.5183e-007   -0.033402      -33.402 -6.2645e-009  -0.0015661      -1.5661

The displacement array (d) shows that Global Nodes 2 and 5 were displaced because these were the only unfixed nodes. The output also displays the reactions at the fixed nodes in the "reactions" matrix. The results matrix (results) show the strain, stress, and axial force in each element. Below are the tabular results of the output above.

To visualize the loaded and unloaded system, the following code graphs the system in its unloaded state and loaded states for Cases 1 and 2. For this display, the magnification factor was changed to 10000. Also, the labels have been removed for the global nodes and elements. The displacement array for Case 1 was changed from "d" to "d1". For Case 2, a new array was created called "d2" using the values in the displacement array from the output above.

SixBarTrussDisplayModified.m
Command Window Output

>> SixBarTrussDisplayModified



The dashed red line represents the undeformed system shape. The solid blue line represents the deformed system shape for Case 1. The dotted black line represents the deformed system shape for Case 2.

The next part of the MATLAB code constructs and loads a three dimensional truss system composed of three bars. The code below represents Example 4.2 (ThreeBarSpaceTrussEx.m) on page 230 of the textbook. The code for this program is structured similarly to the code for the two previous systems. However, this system is a space system where displacement and reactions are measured in three dimensions.

The format of the array "conn" doesn't change. On the other hand, the number of columns in "nodes" and "lmm" increase because the matrices are accommodating for movement in the Z-direction.

This program utilizes two new functions: "SpaceTrussElement.m" and "SpaceTrussResults.m". These are displayed below the program.

SpaceTrussResults.m
Command Window Output

>> ThreeBarSpaceTrussEx

K =

Columns 1 through 6

1459.7      2919.3      -3040.9            0            0            0       2919.3       5838.6      -6081.9            0            0            0      -3040.9      -6081.9       6335.3            0            0            0            0            0            0       3566.7      -3566.7       4953.8            0            0            0      -3566.7       3566.7      -4953.8            0            0            0       4953.8      -4953.8       6880.3            0            0            0            0            0            0            0            0            0            0            0            0            0            0            0            0            0            0      -1459.7      -2919.3       3040.9      -3566.7       3566.7      -4953.8      -2919.3      -5838.6       6081.9       3566.7      -3566.7       4953.8       3040.9       6081.9      -6335.3      -4953.8       4953.8      -6880.3

Columns 7 through 12

0           0            0      -1459.7      -2919.3       3040.9            0            0            0      -2919.3      -5838.6       6081.9            0            0            0       3040.9       6081.9      -6335.3            0            0            0      -3566.7       3566.7      -4953.8            0            0            0       3566.7      -3566.7       4953.8            0            0            0      -4953.8       4953.8      -6880.3            0            0            0            0            0            0            0            0            0            0            0            0            0            0        60000            0            0       -60000            0            0            0       5026.4      -647.45       1912.9            0            0            0      -647.45       9405.4       -11036            0            0       -60000       1912.9       -11036        73216

d =

0           0            0            0            0            0            0            0            0     -0.18705       -2.592      -0.3858

reactions =

6666.7       13333       -13889      -6666.7       6666.7      -9259.3            0            0        23148

results =

0.00050936      101.87        20375   0.00033036       66.072        13214   -0.0001929       -38.58       -23148

The displacement array (d) shows that Global Node 4 was displaced because it was the only unfixed node. The output also displays the reactions at the fixed nodes in the "reactions" matrix. The results matrix (results) show the strain, stress, and axial force in each element. Below are the tabular results of the output above.

This problem is statically determinate. Because the system occupies a three-dimensional space, there exists three equations for equilibrium in the x, y, and z directions. There are also three unknown variables which are the axial forces reacting at Global Node 4 to accommodate the loading force (P).

The relationship between axial forces and directional forces is governed by the following equations.

The place holder "n" represents the nth Global Node or Element. The three equilibrium equations are shown below using the directional forces derived above.

When the equations for the directional forces are substituted into the above equations, the new equations only have the three axial forces as unknowns. The coefficients of these axial forces make up the matrix "F" in the following code. The "L" array defines the lengths of the three elements. The "R" array is the loading at Global Node 4 by force P in the negative Y-direction. The program below simultaneously solves the system of equations to determine the axial forces in the elements.

ThreeBarSpaceTrussStatics.m
Command Window Output

>> ThreeBarSpaceTrussStatics

L =

2.9339      2.8543            2

R =

0      20000           0

F =

0.32721     -0.5045            0      0.65441       0.5045            0     -0.68168     -0.70069           -1

axial_forces =

20375       13214       -23148

The axial forces produced using the statics method matches exactly to the values obtained in the table above for the FEA method. This proves that the system is statically determinate.

To visualize the loaded and unloaded system, the following code graphs the system in its loaded and unloaded state. In order to be able to realize the direction of node movement, the unfixed node displacements have been magnified by the displacement magnification factor (m) of 300. As shown, the original values from the displacement array (d(10), d(11), d(12)) above were used below to get the right proportions of the nodal displacements.

ThreeBarSpaceTrussDisplay.m
Command Window Output

>> ThreeBarSpaceTrussDisplay



The dashed red line represents the undeformed system shape. The solid blue line represents the deformed system shape. The abbreviations "GN" and "E" stand for Global Node and Element, respectively. Note only node 4 is displaced because it is unfixed. This view is a perspective view as if the observer was looking at the system from the global location (-2,-2,3).

The next three view are views along each axis (X, Y, Z), respectively. Line 27 from "ThreeBarTrussDisplay.m" was modified as shown below for each plot. The only other changes to the code were locations of the "GN" and "E" labels which are trivial.

Line 27

Before

After View along X-axis:

View along Y-axis:

View along Z-axis:







Contributing Team Members
Aaron Fisher - Eml4500.f08.Lulz.fisher 21:00, 7 November 2008 (UTC)

Eric Layton - Eml4500.f08.Lulz.Layton.Eric 06:39, 7 November 2008 (UTC)

Sam Miorelli - Eml4500.f08.lulz.abcd 20:48, 7 November 2008 (UTC)

Benjamin Mitchell - Eml4500.f08.Lulz.mitchell.bm 16:59, 7 November 2008 (UTC)

John Saxon - Eml4500.f08.Lulz.js 05:28, 6 November 2008 (UTC)

Andrew Strack - Eml4500.f08.lulz.strack 05:06, 29 October 2008 (UTC)