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: K_6x6d_6x1 = F_6x1

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 w_6x1

w is known as the "weighting mat."

Note that the dot product is used in this equation.

Equation (1)↔ Equation (2)

Force Displacement Relation Proof

Proof that Kd=F

A) Equation(1) → Equation(2): Trivial B) Want to Show Equation(2) → Equation(1) There are six choices that can be used to show this. Choice 1 Select w at w1=1, w2=w3 = ... = w6=0 $\mathbf{w}^T_{1x6} = [1\; 0\; 0\; 0\; 0\; 0]_{1x6}$ $\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F}) = 1[\sum_{j=1}^6{K_{1j}d_{j}-F_{1}}] + 0[\sum_{j=1}^6{K_{2j}d_{j}-F_{2}}]+...+0[\sum_{j=1}^6{K_{6j}d_{j}-F_{6}}]$ $\Rightarrow \sum_{j=1}^{6}{K_{ij}d_j=F_1}$ (This is the 1st Equation) Choice 2 Select w at w1=0, w2=1, w3=w4=...=w6=0 $\mathbf{w}^T_{1x6} = [0\; 1\; 0\; 0\; 0\; 0]_{1x6}$ $\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F}) = 0[\sum_{j=1}^6{K_{1j}d_{j}-F_{1}}] + 1[\sum_{j=1}^6{K_{2j}d_{j}-F_{2}}]+0[\sum_{j=1}^6{K_{3j}d_{j}-F_{3}}]+...+0[\sum_{j=1}^6{K_{6j}d_{j}-F_{6}}]$ $\Rightarrow \sum_{j=1}^{6}{K_{ij}d_j=F_2}$ (This is the 1st Equation) Choice 3 Select w at w1=0, w2=0, w3=1,w4=w5 =w6=0 $\mathbf{w}^T_{1x6} = [0\; 0\; 1\; 0\; 0\; 0]_{1x6}$ $\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F}) = 0[\sum_{j=1}^6{K_{1j}d_{j}-F_{1}}] + 0[\sum_{j=1}^6{K_{2j}d_{j}-F_{2}}]+1[\sum_{j=1}^6{K_{3j}d_{j}-F_{3}}]+0[\sum_{j=1}^6{K_{4j}d_{j}-F_{4}}]...+0[\sum_{j=1}^6{K_{6j}d_{j}-F_{6}}]$ $\Rightarrow \sum_{j=1}^{6}{K_{ij}d_j=F_3}$ (This is the 3rd Equation) Choice 4 Select w at w1=w2=w3=w5=w6=0,w4=1 $\mathbf{w}^T_{1x6} = [0\; 0\; 0\; 1\; 0\; 0]_{1x6}$ $\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F}) = 0[\sum_{j=1}^6{K_{1j}d_{j}-F_{1}}] + ...+ 0[\sum_{j=1}^6{K_{3j}d_{j}-F_{3}}]+1[\sum_{j=1}^6{K_{4j}d_{j}-F_{4}}]+0[\sum_{j=1}^6{K_{5j}d_{j}-F_{5}}]+0[\sum_{j=1}^6{K_{6j}d_{j}-F_{6}}]$ $\Rightarrow \sum_{j=1}^{6}{K_{ij}d_j=F_4}$ (This is the 4th Equation) Choice 5 Select w at w1=w2=w3=w4=w6=0,w5=1 $\mathbf{w}^T_{1x6} = [0\; 0\; 0\; 0\; 1\; 0]_{1x6}$ $\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F}) = 0[\sum_{j=1}^6{K_{1j}d_{j}-F_{1}}] + ...+ 0[\sum_{j=1}^6{K_{4j}d_{j}-F_{4}}]+1[\sum_{j=1}^6{K_{5j}d_{j}-F_{5}}]+0[\sum_{j=1}^6{K_{6j}d_{j}-F_{6}}]$ $\Rightarrow \sum_{j=1}^{6}{K_{ij}d_j=F_5}$ (This is the 5th Equation) Choice 6 Select w at w1=...w5=0,w6=1 $\mathbf{w}^T_{1x6} = [0\; 0\; 0\; 0\; 0\; 1]_{1x6}$ $\mathbf{w}\bullet(\mathbf{Kd}-\mathbf{F}) = 0[\sum_{j=1}^6{K_{1j}d_{j}-F_{1}}] + ...+ 0[\sum_{j=1}^6{K_{5j}d_{j}-F_{5}}]+1[\sum_{j=1}^6{K_{6j}d_{j}-F_{6}}]$ $\Rightarrow \sum_{j=1}^{6}{K_{ij}d_j=F_6}$ (This is the 6th Equation) Hence Kd=F or Kd-F=0

Using Principle of Virtual Work

Taking account for the boundary conditions of the 2 bar truss:

d₁=d₂=d₅=d₆=0

Weighting coefficient must be "kinematically admissible," ie cannot violate the boundary conditions

w₁=w₂=w₅=w₆=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}= \bold{0}_{2x1}$ Equation $(1)$

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

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

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

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

Virtual Displacement Terms

ŵ_2x1 = virtual axial displacement, with corresponding q^(e)_2x1

w_4x1 = 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 w_4x1

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

Equation 6: $(\bold{AB})^T=\bold{B}^T\bold{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): $(\bold{AB})^T=\bold{B}^T\bold{A}^T$

We can recall that $\underline{a}\cdot\underline{b} = \underline{a}^T\underline{b}$, where a and b are nx1 and a^T 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 w_4x1

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

Composite Literature
An example of a material with a varying modulus of elasticity is that of reinforced walls. The rebar embedded in the wall has a different modulus than the concrete so as a whole the modulus varies throughout the wall.

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.

% 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(i), A(i), ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
format short g
results

Results from the debugged two bar truss code

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

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⁽¹⁾ 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 Example

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.

% Six bar truss example
e = 200*10^3; A = 0.001*1000^2; P = 20000.; 
alpha = pi/6; 
nodes = 1000*[0, 0; 4, 0; 0, 3; 4, 3; 2, 2]; 
dof=2*length(nodes);
conn=[1,2; 2,5; 5,3; 2,4; 1,5; 5,4];
lmm = [1, 2, 3, 4; 3, 4, 9, 10; 9, 10, 5, 6;
    3, 4, 7, 8; 1, 2, 9, 10; 9, 10, 7, 8];
elems=size(lmm,1);
K=zeros(dof); R = zeros(dof,1);
debc = [1, 2, 5, 6, 7, 8];
ebcVals = zeros(length(debc),1);
 
%load vector
R = zeros(dof,1); R(3) = P*sin(alpha); R(4) = P*cos(alpha); 
 
% Assemble global stiffness matrix
K=zeros(dof);
for i=1:elems
    lm=lmm(i,:);
    con=conn(i,:);
    k=PlaneTrussElement(e, A, 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

Results of the above code.

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

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.

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

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 length 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.

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);

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 results and sent back to the main program.

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];

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.

%Plot Six beam truss system
%
clear;
close;
 
%model with 2-D beam elements
dof = 2; %dof per node: axial disp x, 2= disp y
%obtain the coordinatates of all nodes
%
n_node = 5;             %number of nodes
n_elem = 6;            %number of elements
total_dof = 2 * n_node; %total dof of system
 
position(:, 1) = [0; 0];
position(:, 2) = [4; 0];
position(:, 3) = [0; 3];
position(:, 4) = [4; 3];
position(:, 5) = [2; 2];
% print the node coord.
for i = 1 : n_node;
    x(i) = position (1,i);
    y(i) = position (2,i);
end
 
node_connect (1, 1) = 1;    %element 1
node_connect (2, 1) = 2;
 
node_connect (1, 2) = 2;    %element 2
node_connect (2, 2) = 5;
 
node_connect (1, 3) = 3;    %element 3
node_connect (2, 3) = 5;
 
node_connect (1, 4) = 2;    %element 4
node_connect (2, 4) = 4;
 
node_connect (1, 5) = 1;    %element 5
node_connect (2, 5) = 5;
 
node_connect (1, 6) = 5;    %element 6
node_connect (2, 6) = 4;
 
% connect all beam elements by connectivity array
for i = 1 : n_elem;
    node_1 = node_connect(1,i);
    node_2 = node_connect(2,i);
    xx = [x(node_1),x(node_2)];
    yy = [y(node_1),y(node_2)];
    axis([ -1 7 -1 5])
    plot(xx,yy,'--')
    hold on
end
 
text(x(node_connect(1,1)),y(node_connect(1,1))-.2,'Global Node 1', 'HorizontalAlignment', 'center')
text(x(node_connect(1,2)),y(node_connect(1,2))-.2,'Global Node 2', 'HorizontalAlignment', 'center')
text(x(node_connect(1,3))-.22,y(node_connect(1,3))+.2,'Global Node 3', 'HorizontalAlignment', 'center')
text(x(node_connect(2,4))+.6,y(node_connect(2,4))+.8,'Global Node 4', 'HorizontalAlignment', 'center')
text(x(node_connect(1,6))-1.27,y(node_connect(1,6)),'Global Node 5', 'HorizontalAlignment', 'center')
 
text (x(node_connect(2,1))*1.1-2,y(node_connect(2,1))/2 -.2,'Element 1', 'HorizontalAlignment', 'center')
text (x(node_connect(2,1))-1.6, y(node_connect(2,1))+1 ,'Element 2', 'HorizontalAlignment', 'center')
text( x(node_connect(2,3))/2 -1, y(node_connect(2,3))/2 ,'Element 3', 'HorizontalAlignment', 'center')
text( x(node_connect(2,4))/2 + 1, y(node_connect(2,4))*1.1 ,'Element 4', 'HorizontalAlignment', 'center')
text( x(node_connect(1,4))*1.5-1.2, y(node_connect(2,3))*.85+.65 ,'Element 5', 'HorizontalAlignment', 'center')
hold on
 
 
%Displacements were multiplied by a factor of 1.2 in order for them to show up
%in the plot
 
position_disp(:, 1) = [0; 0];
position_disp(:, 2) = [5.2133 * 1.2; .25 * 1.2];
position_disp(:, 3) = [0; 3 * 1.2];
position_disp(:, 4) = [4 * 1.2; 3 * 1.2];
position_disp(:, 5) = [1.5439 * 1.2; 2.0122 * 1.2];
 
 
% print the node coord.
for i = 1 : n_node;0
    x(i) = position_disp (1,i);
    y(i) = position_disp (2,i);
end
 
node_connect_disp (1, 1) = 1;    %element 1
node_connect_disp (2, 1) = 2;
 
node_connect_disp (1, 2) = 2;    %element 2
node_connect_disp (2, 2) = 5;
 
node_connect_disp (1, 3) = 3;    %element 3
node_connect_disp (2, 3) = 5;
 
node_connect_disp (1, 4) = 2;    %element 4
node_connect_disp (2, 4) = 4;
 
node_connect_disp (1, 5) = 1;    %element 5
node_connect_disp (2, 5) = 5;
 
node_connect_disp (1, 6) = 5;    %element 6
node_connect_disp (2, 6) = 4;
 
% connect all beam elements by connectivity array
for i = 1 : n_elem;
    node_1 = node_connect_disp(1,i);
    node_2 = node_connect_disp(2,i);
    xx = [x(node_1),x(node_2)];
    yy = [y(node_1),y(node_2)];
    axis([ -1 7 -1 5])
    plot(xx,yy,'-')
    hold on
end
 
 
title ('Five Bar Truss System')
xlabel('x')
ylabel('y')

File:Six bar truss.JPG

Undeformed and deformed six-bar truss system

I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

This is a candidate to be copied to the Wikimedia Commons. Appropriately licensed media are more accessible to other Wikimedia projects if placed on Commons. Any user may perform this move—please see Moving images to the Commons for more information. Please remove this tag after this image has been copied to Commons. Copy to Commons via CommonsHelper

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

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.

% Three bar space truss example
a1 = 200; a2 = 600; e = 200000; P = 20000;
nodes = 1000*[.96, 1.92, 0; -1.44, 1.44, 0; 0, 0, 0; 0, 0, 2]; 
dof=3*length(nodes);
conn=[1,4; 2,4; 3,4];
lmm = [1, 2, 3, 10, 11, 12; 
    4, 5, 6, 10, 11, 12; 
    7, 8, 9, 10, 11, 12];
debc = [1:9];
ebcVals = zeros(length(debc),1);
 
%load vector
R = zeros(dof,1); R(11) = -P; 
 
% Assemble global stiffness matrix
K=zeros(dof);
for i=1:2
    lm=lmm(i,:);
    con=conn(i,:);
    k=SpaceTrussElement(e, a1, nodes(con,:));
    K(lm, lm) = K(lm, lm) + k;
end
lm=lmm(3,:);
con=conn(3,:);
k=SPaceTrussElement(e, a2, nodes(con,:));
K(lm, lm) = K(lm, lm) + k
 
% Nodal solution and reactions
[d, reactions] = NodalSoln(K, R, debc, ebcVals)
results=[];
for i=1:2
    results = [results; SpaceTrussResults(e, a1, ...
            nodes(conn(i,:),:), d(lmm(i,:)))];
end
results = [results; SpaceTrussResults(e, a2, ...
        nodes(conn(3,:),:), d(lmm(3,:)))];
format short g
results

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).

function k = SpaceTrussElement(e, A, coord)
% k = SpaceTrussElement(e, A, coord)
% Generates stiffness matrix of a space 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); z1=coord(1,3);
x2=coord(2,1); y2=coord(2,2); z2=coord(2,3);
L=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L; ns=(z2-z1)/L;
k = e*A/L*[ls^2, ls*ms, ls*ns, -ls^2, -(ls*ms), -(ls*ns);
    ls*ms, ms^2, ms*ns, -(ls*ms), -ms^2, -(ms*ns);
    ls*ns, ms*ns, ns^2, -(ls*ns), -(ms*ns), -ns^2;
    -ls^2, -(ls*ms), -(ls*ns), ls^2, ls*ms, ls*ns;
    -(ls*ms), -ms^2, -(ms*ns), ls*ms, ms^2, ms*ns;
    -(ls*ns), -(ms*ns), -ns^2, ls*ns, ms*ns, ns^2];

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 length 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.

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);

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 results and sent back to the main program.

function results = SpaceTrussResults(e, A, coord, disps)
% results = SpaceTrussResults(e, A, coord, disps)
% Compute space truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
 
x1=coord(1,1); y1=coord(1,2); z1=coord(1,3);
x2=coord(2,1); y2=coord(2,2); z2=coord(2,3);
L=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L; ns=(z2-z1)/L;
T=[ls,ms,ns,0,0,0; 0,0,0,ls,ms,ns];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;
results=[eps, sigma, force];

3 Bar Space Truss Results

Below are the results from the three bar space truss example.

K =
 
  1.0e+004 *
 
  Columns 1 through 7 
 
    0.1460    0.2919   -0.3041         0         0         0         0
    0.2919    0.5839   -0.6082         0         0         0         0
   -0.3041   -0.6082    0.6335         0         0         0         0
         0         0         0    0.3567   -0.3567    0.4954         0
         0         0         0   -0.3567    0.3567   -0.4954         0
         0         0         0    0.4954   -0.4954    0.6880         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
         0         0         0         0         0         0         0
   -0.1460   -0.2919    0.3041   -0.3567    0.3567   -0.4954         0
   -0.2919   -0.5839    0.6082    0.3567   -0.3567    0.4954         0
    0.3041    0.6082   -0.6335   -0.4954    0.4954   -0.6880         0
 
  Columns 8 through 12 
 
         0         0   -0.1460   -0.2919    0.3041
         0         0   -0.2919   -0.5839    0.6082
         0         0    0.3041    0.6082   -0.6335
         0         0   -0.3567    0.3567   -0.4954
         0         0    0.3567   -0.3567    0.4954
         0         0   -0.4954    0.4954   -0.6880
         0         0         0         0         0
         0         0         0         0         0
         0    6.0000         0         0   -6.0000
         0         0    0.5026   -0.0647    0.1913
         0         0   -0.0647    0.9405   -1.1036
         0   -6.0000    0.1913   -1.1036    7.3216
 
 
d =
 
         0
         0
         0
         0
         0
         0
         0
         0
         0
   -0.1871
   -2.5920
   -0.3858
 
 
reactions =
 
  1.0e+004 *
 
    0.6667
    1.3333
   -1.3889
   -0.6667
    0.6667
   -0.9259
         0
         0
    2.3148
 
 
results =
 
   0.00050936       101.87        20375
   0.00033036       66.072        13214
   -0.0001929       -38.58       -23148
 
>>

Code for the lot of the 3 Bar Space Truss

Below is the code for the 3 Bar Space Truss.

%Plot three beam truss system
%
clear;
close;
 
%model with 2-D beam elements
dof = 3; %dof per node: axial disp x, 2= disp y
%obtain the coordinatates of all nodes
%
n_node = 4;             %number of nodes
n_elem = 3;            %number of elements
total_dof = 3 * n_node; %total dof of system
 
position(:, 1) = [-.96; 1.92; 0 ];
position(:, 2) = [-1.44; 1.44; 0 ];
position(:, 3) = [0; 0; 0 ];
position(:, 4) = [0; 0; 2]
% print the node coord.
for i = 1 : n_node;
    x(i) = position (1,i);
    y(i) = position (2,i);
    z(i) = position (3,i);
end
 
node_connect (1, 1) = 1;    %element 1
node_connect (2, 1) = 4;
 
node_connect (1, 2) = 2;    %element 2
node_connect (2, 2) = 4;
 
node_connect (1, 3) = 3;    %element 3
node_connect (2, 3) = 4;
 
% connect all beam elements by connectivity array
for i = 1 : n_elem;
    node_1 = node_connect(1,i);
    node_2 = node_connect(2,i);
    xx = [x(node_1),x(node_2)];
    yy = [y(node_1),y(node_2)];
    zz = [z(node_1),z(node_2)];
    axis([ -3 3 -3 3 -3 3])
    plot3(xx,yy,zz,'--')
    view([1,1,1])
    hold on
end
 
text(x(node_connect(1,1)),y(node_connect(1,1)),z(node_connect(1,1))-.5,'Global Node 1', 'HorizontalAlignment', 'center')
text(x(node_connect(1,2))-1,y(node_connect(1,2)),z(node_connect(1,2))+.3,'Global Node 2', 'HorizontalAlignment', 'center')
text(x(node_connect(1,3)),y(node_connect(1,3)),z(node_connect(1,3))-.5,'Global Node 3', 'HorizontalAlignment', 'center')
text(x(node_connect(2,3)),y(node_connect(2,3)),z(node_connect(2,3))+.5,'Global Node 4', 'HorizontalAlignment', 'center')
 
%text (x(node_connect(1,1))/2,y(node_connect(1,1))/2 +1,'Element 1', 'HorizontalAlignment', 'center')
%text (x(node_connect(2,2))/2 -1, y(node_connect(2,2))/2 -.5,'Element 2', 'HorizontalAlignment', 'center')
%text( x(node_connect(2,3))/2, y(node_connect(2,3))*.75 ,'Element 3', 'HorizontalAlignment', 'center')
 
 
 
%Scaled by 1/10
 
position_disp(:, 1) = [-.96; 1.92; 0 ];
position_disp(:, 2) = [-1.44; 1.44; 0 ];
position_disp(:, 3) = [0; 0; 0 ];
position_disp(:, 4) = [-0.1871; -2.5920; 2 + -0.3858]
 
% print the node coord.
for i = 1 : n_node;0
    x(i) = position_disp (1,i);
    y(i) = position_disp (2,i);
    z(i) = position_disp (3,i);
end
 
node_connect_disp (1, 1) = 1;    %element 1
node_connect_disp (2, 1) = 4;
 
node_connect_disp (1, 2) = 2;    %element 2
node_connect_disp (2, 2) = 4;
 
node_connect_disp (1, 3) = 3;    %element 3
node_connect_disp (2, 3) = 4;
 
% connect all beam elements by connectivity array
for i = 1 : n_elem;
    node_1 = node_connect_disp(1,i);
    node_2 = node_connect_disp(2,i);
    xx = [x(node_1),x(node_2)];
    yy = [y(node_1),y(node_2)];
    zz = [z(node_1),z(node_2)];
    axis([ -3 3 -3 3 -3 3])
    view([1,1,1])
    plot3(xx,yy,zz,'-')
    hold on
end
 
 
title ('Three Bar Space Truss System')
xlabel('x')
ylabel('y')
zlabel('z')

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.

File:3barspace111.JPG

Undeformed and deformed three bar space truss with a viewing mofe of [1 1 1]

I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

This is a candidate to be copied to the Wikimedia Commons. Appropriately licensed media are more accessible to other Wikimedia projects if placed on Commons. Any user may perform this move—please see Moving images to the Commons for more information. Please remove this tag after this image has been copied to Commons. Copy to Commons via CommonsHelper

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.

File:3barspacex.JPG

Undeformed and deformed three bar space truss with a viewing down plane x

I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

This is a candidate to be copied to the Wikimedia Commons. Appropriately licensed media are more accessible to other Wikimedia projects if placed on Commons. Any user may perform this move—please see Moving images to the Commons for more information. Please remove this tag after this image has been copied to Commons. Copy to Commons via CommonsHelper

File:3barspacey.JPG

Undeformed and deformed three bar space truss with a viewing down plane y

I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

This is a candidate to be copied to the Wikimedia Commons. Appropriately licensed media are more accessible to other Wikimedia projects if placed on Commons. Any user may perform this move—please see Moving images to the Commons for more information. Please remove this tag after this image has been copied to Commons. Copy to Commons via CommonsHelper

File:3barspacez.JPG

Undeformed and deformed three bar space truss with a viewing down plane z

I, the copyright holder of this work, hereby release my work worldwide into the public domain. Where this is not legally possible, I grant anyone the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

This is a candidate to be copied to the Wikimedia Commons. Appropriately licensed media are more accessible to other Wikimedia projects if placed on Commons. Any user may perform this move—please see Moving images to the Commons for more information. Please remove this tag after this image has been copied to Commons. Copy to Commons via CommonsHelper

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

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

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

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

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

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

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

$\displaystyle f^e(T)=p^e$

Or more detailed:

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

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

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