User:Eml4500.f08.ateam.nobrega/Homework Report Six

The Principle of Virtual Work: A Continuation of the Dynamics of an Elastic Bar
For the initial conditions, at t = 0:


 * $$ U(x,t=0) = \overline{U}(x) $$       The known functional displacement
 * $$ \dot{U}(x,t=0) = \overline{V}(x) $$       The known functional velocity

To create an equation for a multiple degree of freedom system, start with


 * $$ \frac{d}{dz}[EA \frac{du}{dx}] + f = m \ddot{u} $$    Equation (1)

This becomes,


 * $$ -kd + f = m \ddot{d} $$

or


 * $$ m \ddot{d} + kd = f $$    Equation (2)

To derive Equation 2 from 1 move all the terms to one side and take the integral with respect to x. This produces:


 * $$ \int_{x=0}^{x=L} w(x)(\frac{d}{dx}[EA \frac{du}{dx}] + f - m \ddot{u}) \,dx = 0 $$    Equation (3)

This is true for all w(x), the weighting function. It should be noted that producing Equation 3 from Equation 1 is a trivial solution, but producing Equation 1 from Equation 3 is not. Equation 3 can be re-written as,


 * $$ \int_0^L w(x)g(x) \,dx = 0 $$

Here g(x) is substituted into Equation 3 for


 * $$ \frac{d}{dx}[EA \frac{du}{dx}] + f - m \ddot{u} $$

Because Equation 3 holds true for all w(x), select w(x) such that it is equal to g(x). Doing so allows Equation 3 to become:


 * $$ \int_0^L g^2(x) \,dx = 0 $$

for $$ g(x) \ge 0 $$

Integration by Parts
When integrating an argument with two functions dependent upon the same variable, you must use the method of integration by parts.


 * $$ \int r(x)s(x)dx

(rs)'=r's + rs'

r'=\frac{dr}{dx}

s'=\frac{ds}{dx}

\int (rs)' = \int r's + \int rs'

\int (rs)' = rs - \int rs' $$

(Recall continuation of Principle of Virtual Work Equation 3 pg 29-1 also featured above)

First term: $$ r(x)=EA\frac{du}{dx}

s(x)=w(x) $$

by integration by parts:


 * $$ \int_{x=0}^{x=L} w(x)(\frac{d}{dx}[EA \frac{du}{dx}])dx = WEA\frac{du}{dx} - \int_{0}^{L} \frac{dw}{dx}(EA)[\frac{du}{dx}]dx

$$

The negative sign in front of the integral on the left hand of the equation corresponds to the negative sign in the equation:

$$ F-kd = M\ddot{d} $$

Elastic Beam Problem
Boundary Conditions: At x=0, select w(x) such that w(0)=0 (i.e. Kinematically Admissible)

The motivation is the discrete PVW applied to the equation below (also found on pg 10-1)

$$ W([k] \begin{bmatrix}     d_3 \\      d_4 \end{bmatrix} -F) = 0 $$

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

Where F3 and F4 are known and F1, F2, F5 and F6 are unknown reactions. Since W can be selected arbitrarily, select W such that W1=W2=W5=W6=0 so to eliminate equations involving unknown reactions. Thus, eliminating rows 1,2,5 and 6.

$$ kd=F $$

Where k is a 2x2 matrix, d is a 2x1 matrix and F is a 2x1 matrix. Note: W(kd-F)=0 for all values of W.

Unknown Reaction:

$$ N(0,t)=EA(0)\frac{du}{dx}(0,t) $$

Continue PVW with functions W(L) and F(t)

$$ -\int_{0}^{L}\frac{dW}{dx}(EA)\frac{du}{dx}dx + \int_{0}^{L} W(x)[f-m\ddot{u}]dx = 0 $$

For all W(x) such that W(0)=0

$$ \int_{0}^{L}w(m\ddot{u})dx + \int_{0}{L}\frac{dW}{dx}EA\frac{du}{dx}dx = W(L)F(t) + \int_{0}^{L}W f dx $$

For all W(x) such that W(0)=0.



Stiffness Terms:





Assume displacement u(x) for xi <= x <= xi+1
 * (i.e. x belongs to [xi, xi+1])

Motivation for linear interpolation of u(x): 2 Bar Truss



The deformed shape of the truss is a straight line as well because of the implicit assumption of linear interpolation of displacement between 2 nodes.

Consider the case where there are only axial displacements (i.e. 0 translational displacement).

Question: Express u(x) in terms of di=u(xi) and di+1=u(xi+1) as a linear function in x (i.e. linear interpolation)

u(x) = Ni(x)di + Ni+1(x)di+1

where Ni(x) and Ni+1(x) are linear functions of x.

N i+1 (x)= (x - xi) * (xi+1 - xi)-1



Honesty, Imagination, Ethics
The following quotes are excerpts from Dr. L. Vu-Quoc main homework page for the Finite Element Analysis: EML 4500 course at the University of Floida.

Continuation PVW to Discrete PVW
Langrangian Interpolation

Motivation for the form of $$ N_{i}(x) $$ and $$ N_{i+1}(x) $$

1. $$ N_{i}(x) $$ and $$ N_{i+1}(x) $$ are linear (straight lines). Thus any linear combination of $$ N_{i} $$ and $$ N_{i+1} $$ is also linear and in particular the expression for $$ u(x) $$ on p.31-3.

$$ N_{i}(x) = \alpha_{i} + \beta_{i}(x) $$, where $$ \alpha_{i} $$ and $$ \beta_{i} $$, are real numbers

$$ N_{i+1}(x) = \alpha_{i+1} + \beta_{i+1}(x) $$, where $$ \alpha_{i+1} $$ and $$ \beta_{i+1} $$, are real numbers

$$ N_{i}d_{i}+N_{i+1}d_{i+1} = (\alpha_{i}+\beta_{i}x)d_{i}+(\alpha_{i+1}+\beta_{i+1}x)d_{i+1} $$, is clearly a linear function in x

2. Recall the equation for $$ u(x) $$

$$ u(x_i) = N_{i}(x_i)d_i+N_{i+1}(x_i)d_{i+1} $$

Where $$ N_{i}(x_{i})d_i = 1 $$ and $$ N_{i+1}(x_{i})d_{i+1} = 0 $$

Continueing the FEM via PVW
From a previous lecture on meeting 31 about the linear interpolation of U(x), W(x) can be shown by applying the same principle.

i.e.

$$ W(x)=N_{i}(x)W_{i}+N_{i+1}(x)W_{i+1} $$

Now the element stiffness matrix for element I can be written as:

$$ \beta = \int_{x_{i}}^{x_{i+1}}{[N_{i}'W_{i}+N_{i+1}'W_{i+1}](EA)} $$

with

$$[N_{i}'W_{i}+N_{i+1}'W_{i+1}]=W'(x)$$

$$[N_i'd_i+N_{i+1}'d_{i+1}]=U'(x)$$

$$ N_i':=\frac{dN_i(x)}{dx} $$

Likewise for $$N_{i+1}'$$

NOTE

$$ u(x)=\begin{bmatrix} N_i(x) & N_{i+1}(x) \end{bmatrix} \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} =\mathbf{N(x)_{(1x1)}} \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} $$

$$ \frac{du(x)}{dx}=\begin{bmatrix} N_i'(x) & N_{i+1}'(x) \end{bmatrix} \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} =\mathbf{B(x)} \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} $$

Similarly W(x) can be shown as:

$$ W(x)=\mathbf{N(x)}\begin{Bmatrix} W_i\\ W_{i+1} \end{Bmatrix} $$

$$ \frac{dW(x)}{dx}=\mathbf{B(x)}\begin{Bmatrix} W_i\\ W_{i+1} \end{Bmatrix} $$

Now, let us recall the element degrees of freedom as discussed in a previous lecture



therefore

$$ \begin{Bmatrix} d_i\\ d_{i+1}

\end{Bmatrix} = \begin{Bmatrix} d_1^{(i)}\\ d_2^{(i)} \end{Bmatrix} = \mathbf{d^{(i)}} $$

$$ \begin{Bmatrix} W_i\\ W_{i+1}

\end{Bmatrix} = \begin{Bmatrix} W_1^{(i)}\\ W_2^{(i)} \end{Bmatrix} = \mathbf{W^{(i)}} $$

The element stiffness matrix can now be written as

$$ \beta =\int_{x_{i}}^{x_{i+1}}{{BW}^{(i)}_{(1x1)}(EA)_{(1x1)}{Bd}^{(i)}_{(1x1)}}dx = {W}^{(i)}\cdot({k}^{(i)}{d}^{(i)}) $$

$$ \beta =\int_{x_i}^{X_{i+1}}(EA)({BW^{(i)}})_{(1x1)}\cdot({Bd^{(i)}}){(1x1)}dx

$$

$$ ({BW^{(i)}})^T({Bd^{(i)}}) $$

where

$$ {BW^{(i)T}} = {W^{(i)T}B^T} = {W^{(i)} \cdot B^T} $$

$$ \beta= \mathbf{W^{(i)}}\cdot\int \mathbf{B^T}EA\mathbf{B}dx\mathbf{d^{(i)}} $$

$$ \mathbf{k^{(i)}}_{(2x2)}=\int_{x_i}^{x^{i+1}}{\mathbf{B(x)^T}}_{(2x1)}(EA)_{1x1}\mathbf{B(x)}_{(1x2)}dx $$

If we consider $$ EA $$ to constant the matrix $$ {k^{(i)}} $$ can be written as

$$\textbf{k}^{(i)} = \frac{EA}{L^{(i)}}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$$

The next step is the transformation of variable coordinate system from:

$${x}$$ to $$ \tilde{x} $$

Given

$$ d\tilde{x}=dx $$

$$ \tilde{x}:=x-x_i $$

Now $$ {k^{(i)}} $$ can be written as

$$ \mathbf{k^{(i)}}=\int_{\tilde{x}=0}^{\tilde{x}=L^{(i)}}{\mathbf{B^T(\tilde{x})}}EA(\tilde{x})\mathbf{B(\tilde{x}})d(\tilde{x}) $$

The following problem outlined by the image below can be solved using the transformed coordinate system.

Given:



With the following functions for the modulus elasticity and cross sectional area $$ {k^{(i)}} $$ matrix can be written as

$$A(\tilde{x}) = N_1^{(i)}(\tilde{x})A_1 + N_2^{(i)}(\tilde{x})A_2$$

$$E(\tilde{x}) = N_1^{(i)}(\tilde{x})E_1 + N_2^{(i)}(\tilde{x})E_2$$

$$ \mathbf{k^{(i)}}=\int_{\tilde{x}=0}^{\tilde{x}=L^{(i)}}{\mathbf{B^T(\tilde{x})}}{[N_1^{(i)}(\tilde{x})A_1 + N_2^{(i)}(\tilde{x})A_2]}{[N_1^{(i)}(\tilde{x})E_1 + N_2^{(i)}(\tilde{x})E_2]}(\tilde{x})\mathbf{B(\tilde{x}})d(\tilde{x}) $$

$$N_{1}^{(i)}(\tilde{x})=1-\frac{\tilde{x}}{L^{(i)}}= \begin{cases} 1 & \text{ at } \tilde{x}=0 \\ 0 & \text{ at } \tilde{x}=L^{(i)} \end{cases}$$

$$N_{2}^{(i)}(\tilde{x})=\frac{\tilde{x}}{L^{(i)}}= \begin{cases} 0 & \text{ at } \tilde{x}=0 \\ 1 & \text{ at } \tilde{x}=L^{(i)} \end{cases} $$

Homework:

Book p. 159 ; Set $$ E_{1}=E_{2}=E_{3}$$, let $$A(\tilde{x})$$ be linear as previously defined, and obtain $$\mathbf{k}^{(i)}$$ and compare to expression given in the book

Where:

$$A(\tilde{x})=N_{1}^{(i)}(\tilde{x})A_{1}+N_{2}^{(i)}(\tilde{x})A_{2}$$

and

$$E(\tilde{x})=N_{1}^{(i)}(\tilde{x})E_{1}+N_{2}^{(i)}(\tilde{x})E_{2}$$

Compare general $$\mathbf{k}^{(i)}$$ from before to $$\mathbf{k}_{avg}^{(i)}$$ where,

$$ \frac {(E_{1}+E_{2})}{2}\frac{(A_{1}+A_{2})}{2}\frac{1}{L^{(i)}} \begin{bmatrix} 1 &-1 \\ -1 & 1 \end{bmatrix}=\mathbf{k}_{avg}^{(i)}$$



Remember: Recall the mean value Theorem (MVT) and its relation to the centroid:

Mean Value Theorem: for $$\bar{x}\epsilon [a,b]$$,     $$a\leq \bar{x}\leq b$$,

$$ \int_{x=a}^{x=b}{f(x)dx}=f(\tilde{x})[b-a] $$

Where the centroid is defined as:

$$\int_{A}^{}{xdA} = \bar{x} \int_{A}^{}{dA} = \bar{x}A$$

Therefore:

$$\int_{x=a}^{x=b}{f(x)g(x)dx} = f(\bar{x})g(\bar{x})[b-a]$$ where  $$ a\leq \bar{x}\leq b$$

In general:

$$f(\bar{x})\neq \frac{1}{b-a} \int_{a}^{b}{f(x)dx}$$

$$g(\bar{x})\neq \frac{1}{b-a} \int_{a}^{b}{g(x)dx} $$

The Electric Pylon
Code was provided that creates a two dimensional truss system of an electric pylon. This code was to be scaled to 60 meters and loaded with a downward 1000 N force on the right most tip. The stresses in the element were to be calculated and the largest compressive and tensile stresses identified in the truss. The mass of each element needed to be taken into account with the mass of each element being redistributed to its end most points. The lowest eigenvectors were then plotted as well as the lowest periods of vibration.

Functions Needed to Load and Solve the Truss System
Length.m

The Length function calculates the element length based on the elements coordinates. function L = Length(coord) x_1 = coord(1,1); x_2 = coord(2,1); y_1 = coord(1,2); y_2 = coord(2,2); L = sqrt((x_2 - x_1)^2 + (y_2 - y_1)^2); %calculates element length using coordinates

ReduceM.m

This function takes the mass matrix and reduces it. function Mr = ReduceM(Mass, F, debc, ebcVals) %calculates the reduced matrix dof = length(F); df = setdiff(1:dof, debc); Mr = Mass(df, df);

ReduceS.m

This function reduces the stiffness matrix. function Kr = ReduceS(K, F, debc, ebcVals) %Reduces the stiffness matrix dof = length(F); df = setdiff(1:dof, debc); Kr = K(df,df);

PlaneTrussElement.m

This function calculates the length of the bar and the director cosines as well as the elemental stiffness matrix before returning it to the call function. 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

This function uses the global stiffness matrix to calculate the number of global degrees of freedom. It calculates the global displacements and the reactions seen in the bars. function [d, reactions] = NodalSoln(K, R, deg, num) dof = length(R); diff = setdiff(1:dof, deg); Ksub = K(diff, diff); Rsub = R(diff) - K(diff, deg)*num; value = Ksub\Rsub; d = zeros(dof,1); d(deg) = num; d(diff) = value; reactions = K(deg,:)*d - R(deg);

PlaneTrussResults.m

This function creates the transform matrix T after calculating the length of the bars and the director cosines. This transform matrix turns the global displacements into axial displacements so it can compute the axial forces, stress, and strain acting within the bar. 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];

Results
Highest Tensile Stress : 9.05 MPa

Element with Highest Tensile Stress: 81

Highest Compressive Stress: -8.70 MPa

Element with Highest Compressive Stress: 55

Lowest three vibrational periods :    0.5465 secs, 0.1265 secs , 0.1159 secs

Deformed Truss








Problem Statically Determinant?
This problem is not statically determinate. Each element has two degrees of freedom for each node, or a total of four degrees of freedom per element. With only three static equations, the problem is statically indeterminate. The two force member method does not apply as the forces cannot be assumed to act only in the direction of member.

Eigenvalue Plot Deformation










Eigenvalue Plot Deformation








Contributing Team Members
Eml4500.f08.ateam.nobrega 01:19, 21 November 2008 (UTC)

Eml4500.f08.ateam.boggs.t 01:37, 21 November 2008 (UTC)

Eml4500.f08.ateam.carr 02:50, 21 November 2008 (UTC)

Eml4500.f08.ateam.shah 18:21, 21 November 2008 (UTC)

Eml4500.f08.ateam.mcnally 18:51, 21 November 2008 (UTC)

Eml4500.f08.ateam.rivero 19:31, 21 November 2008 (UTC)