User:Uf.team.wiki/HW6

Principle of Virtual Work (continuous) of dynamics of the elastic bar problem
Partial Differential Equation:

$$ \frac{\partial }{\partial x} \left [(EA) \frac{\partial u}{\partial x}\right ] + f = m \ddot u \qquad \qquad (1)$$

Equation (1) is equivalent to the multiple degree of freedom system (MDOF):

$$- \mathbf {Kd} + \mathbf F = \mathbf {M \ddot d} \Rightarrow $$

$$ \mathbf {M \ddot d} + \mathbf {Kd} = \mathbf F \qquad \qquad (2)$$

Derive (2) from (1):

$$ \int\limits_{x=0}^{x=L} W(x) \left \{ \frac{\partial }{\partial x}\left [(EA) \frac{\partial u}{\partial x}\right ] + f - m \ddot u \right \} dx = 0 \qquad \qquad (3) \qquad$$


 * NOTE: Equation (3) applies to all possible W(x), where W(x) is the weighting function


 * $$ (1) \Rightarrow (3)$$ : trivial
 * $$ (3) \Rightarrow (1)$$ : not trivial


 * (3) written as $$ \int W(x)g(x)dx = 0 $$ [for all W(x)]

Since (3) holds for all W(x), select W(x)=g(x), then (3) becomes:

$$ \int \underbrace{g^2}_{\ge 0} dx = 0 \quad \Rightarrow \quad g(x)=0$$

Integration by parts
Recall: Continuous PVW {Equation (3), above}


 * 1st term: $$ \qquad r(x) = (EA)\frac{\partial u}{\partial x} $$
 * 2nd term: $$ \qquad s(x) = W(x) $$

Integration by parts:

$$ \int\limits_{x=0}^{x=L} \underbrace{W(x)}_{s} \frac{\partial }{\partial x} \underbrace{\left [(EA) \frac{\partial u}{\partial x}\right ]}_{r} dx = 0 \left [ W (EA) \frac{\partial u}{\partial x}  \right ] - \int\limits_{0}^{L}\frac{dW}{dx}(EA)\frac{\partial u}{\partial x}dx $$

$$ \Rightarrow \underbrace{W(L)(EA)(L)\frac{\partial u}{\partial x}(L,t)}_{F=N(L,t)} - \underbrace{W(0)(EA)(0)\frac{\partial u}{\partial x}(0,t)}_{N(0,t)}- \int\limits_{0}^{L}\frac{dW}{dx}(EA)\frac{\partial u}{\partial x}dx $$

Model Problem


Looking at the figure above, it is shown that there are two Boundary Conditions
 * @ x=0, select W(x) such that W(0)=0 (ie: kinematically admissible)

Motivation: Apply the discrete Principle of Virtual Work:

$$ \underset{6 \times 1}{\mathbf W } \cdot \underbrace {\left (\underset{6 \times 2}{[\qquad] }\underset{2 \times 1}{\begin{Bmatrix} d_{3}\\ d_{4}\\ \end{Bmatrix}} - \underset{2 \times 1}{\mathbf F} \right )}_{6 \times 1} = \underset{1 \times 1}{0}$$

$$ \mathbf F^T = \begin{bmatrix} F_1 &  F_2  &  F_3  &  F_4  &  F_5  &  F_6 \\ \end{bmatrix} $$

(where F1, F2, F5, F6 are known reactions)

Since W can be selected arbitrarily, select W such that W1 = W2, W5 = W6 so to eliminate equations involving unknown reactions → eliminate rows 1,2,5,6:

$$ \underset{2 \times 2}{\overline \mathbf K} \underset{2 \times 1}{\overline \mathbf d} = \underset{2 \times 1}{\overline \mathbf F} $$

NOTE: Recall that:

$$\overline \mathbf W \cdot (\mathbf {Kd}- \mathbf F) = 0 \hookrightarrow for all \quad{\underset{6 \times 1}\mathbf W} \qquad (2)$$

Where: W_{3}\\ W_{4}\\ \end{Bmatrix} $$ d_{3}\\ d_{4}\\ \end{Bmatrix} $$ F_{3}\\ F_{4}\\ \end{Bmatrix} $$
 * $$\underset{2 \times 1}{\overline \mathbf W} \begin{Bmatrix}
 * $$\underset{2 \times 1}{\overline \mathbf d} \begin{Bmatrix}
 * $$\underset{2 \times 1}{\overline \mathbf F} \begin{Bmatrix}

Unknown reaction: $$ N(0,t) = (EA)(0) \frac {\partial u}{\partial x} (0,t) $$

Continuous PVW: $$ W(L)F(t) - \int\limits_{0}^{L} \frac {\partial W}{\partial x}(EA)\frac {\partial u}{\partial x}dx + \int W(x)[f-m \ddot u ] dx = 0 \qquad \qquad $$

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

$$ \Rightarrow \int\limits_{0}^{L} W (m \ddot u )dx + \int \limits_{0}^{L} \frac {\partial W}{\partial x}(EA)\frac {\partial u}{\partial x}dx = W(L)F(t) + \int \limits_{0}^{L} W f dx \qquad $$

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

Sample table for comparisons
When looking at the Principles of Virtual Work, a comparison can be made between the equations for the continuous setting and the discrete setting. A summary of this comparison is made below and shown in the table.

Looking at the stiffness term, and focusing on the general element i with the use of of a dotted area shown in Figure 2 below, represents an example of interpolation using an elastic bar.



Focusing in on a general element, i in the interpolation process, we may look at the figure below to aid in the calculation and assuming that the displacement u(x) is for $$x_i\leq x\leq x_{i+1}$$.



Figure 4 shown below is an example of a 2 - Bar truss example that serves as motivation when calculating the interpolation of u(x). Here there is an assumed linear interpolation between both nodes 1 and 2. Here ther deformed shape is a straight line, that was part of the implicit assumption made in the linear interpolation.



Continuation: Using the Principle of Virtual Work to obtain the Discrete Principal of Virtual Work
Rube Goldberg device

Walter Benjamin

- "honesty, imagination, ethics"

Going from continuous PVW to Discrete PVW
Lagrangian interpolation

motivation for form of $$\mathbf{N_{i}(x)}$$ and $$\mathbf{N_{i+1}(x)}$$

1. It is assumed that $$\mathbf{N_{i}(x)}$$ and $$\mathbf{N_{i+1}(x)}$$ are straight lines, and thus any linear combination of $$\mathbf{N_{i}(x)}$$ and $$\mathbf{N_{i+1}(x)}$$ is also linear, and in particular expression for $$\mathbf{u(x)}$$:

$$\mathbf{u(x)=N_i(x)d_i+N_{i+1}(x)d_{i+1}}$$

where $$N_i(x) and N_{i+1}(x)$$ are linear functions in x and equations are shown below.


 * $$\mathbf{N_{i}(x)={\alpha_i} + {\beta_i}x}$$


 * $$\mathbf{N_{i+1}(x)={\alpha_{i+1}} + {\beta_{i+1}}x}$$

where ($$\mathbf{\alpha_{i}, \beta_i, \alpha_{i+1}, \beta_{i+1}}$$) are all real numbers

Looking at the linear combination of $$\mathbf{N_{i}}$$ and $$\mathbf{N_{i+1}(x)}$$ :

$$\mathbf{N_{i+1}(x)}$$ : $$\mathbf{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} }$$

$$ =\mathbf{(\alpha_{i}d_{i} + \alpha_{i+1}d_{i+1}) +(\beta_{i}d_{i} + \beta_{i+1}d_{i+1})x} $$

Looking at the solution above, one can see that it is clearly a linear function of x

2) For the second case, we can recall equation for $$u(x)$$ (interpolation of $$u(x))$$

$$ \mathbf{u(x_{i})= N_{i}(x_{i})d_{i}+ N_{i+1}(x_{i})d_{i+1}}$$


 * where $$ \mathbf{u(x_{i})= d_{i}}$$ and $$\mathbf{ N_{i+1}(x_{i})d_{i+1} }= 0$$

Use of Interpolation
Looking at the interpolation result for u(x), the same interpolation is applied to form W(x) shown below for an approximation of the weighting coefficient:

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

Next, we are ready to construct the element stiffness matrix for any standard element i, from the elemental stiffness matrix:

$$\beta = \int_{ x }^{{x}_{ i+1}}  [{ N}_{ \acute{i}}{w}_{i}+ {N }_{\acute{i}+1 }{w}_{i+1}](EA) [{ N}_{\acute{i}}{d}_{i}+ {N }_{\acute{i}+1 }{d}_{i+1}]dx$$

Where : $${N}_{\acute{i}}=\frac{d{N}_{i}(x)}{dx}$$ and $${N}_{\acute{i}+1}=\frac{d{N}_{i+1}(x)}{dx}$$

From β:

$$\acute{w}(x)=[{ N}_{ \acute{i}}{w}_{i}+ {N }_{\acute{i}+1 }{w}_{i+1}]$$

$$\acute{u}(x)=[{ N}_{ \acute{i}}{d}_{i}+ {N }_{\acute{i}+1 }{d}_{i+1}]$$

It is important to note:

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

where N(x) can be approximated as:

$$\overline \underset{1\times 2}{\mathbf N(x)}=[{N}_{i}(x) {N}_{{i}+1}(x)]$$

Next, we can take the derivative of u with respect to x:

$$\frac{du}{dx} = {\begin{bmatrix}N'_i(x) & N'_{i+1}(x)\end{bmatrix}}\begin{Bmatrix}d_i \\ d_{i+1}\end{Bmatrix}$$

where B(x) can be approximated as:

$$\overline \underset{1\times 2}{\mathbf B(x)}=[{N}_{\acute{i}}(x) {N}_{{\acute{i}}+1}(x)]$$

Similarly the weighting coefficient, W(x) can be approximated as:

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

where taking the derivative with respect to x yields:

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

Recall: Use of element with degrees of freedom
In the figure below, both the global and local degrees of freedom are shown on each of the axial displacement vectors



From the above figure, one may truly appreciate the power of simplifying both notations for global and local degree of freedoms by using matricies to express the displacements and weighting coefficients

$$\begin{Bmatrix} d_i \\ d_{i+1} \end{Bmatrix} = \begin{Bmatrix} d_1^{(i)} \\ d_2^{(i)} \end{Bmatrix} = \underline\textbf{d}^{(i)} $$

$$\begin{Bmatrix} W_i \\ W_{i+1} \end{Bmatrix} = \begin{Bmatrix} W_1^{(i)} \\ W_2^{(i)} \end{Bmatrix} = \underline\textbf{W}^{(i)} $$

$$\beta =\int_{0}^{L}{\frac{dW}{dx}(EA)\frac{\delta u}{\delta x}dx}=\int_{x_{i}}^{x_{i+1}}{(\underline B \times \underline W^i)(EA)(\underline B \times \underline d^i)dx } $$

$$=w^i\cdot (\underline k ^i \underline d ^i)$$

The above equation for β is rearranged:

$$\beta = \underline w^i \cdot (\int \underline B^T (EA)\underline B dx)\underline d^i$$

The elemental stiffness matrix can be rewritten as seen below:

$$\underline k^i=\int_{x_{i}}^{X_{i+1}}{\underline B(x)^T (EA)\underline B(x) dx)\underline dx}$$

For this assignment EA is taken to be constant, then the equivalent value for K becomes:

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

Coordinate Transformation
The first step is to transform the variable coordinate system from $$x$$ to $$\tilde{x}$$

Therefore, the notation can be abbreviated as:

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

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

Using this relationship, the equation below for $$\textbf{k}^{(i)}:$$

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

Below is a figure of an element with varying areas, and moduli of elasticity and will be solved.



Associated with the figure above, the equations similarly for the area and the modulus of elasticities are shown below for the element i.

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

Incorporation of the Mean Value Theorm (MVT)


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


 * Shape functions $$N_{1}^{(i)},N_{2}^{(i)} $$ for the Basis

 Comparing Stiffness Matices 


 * $$k^{i}$$ to the matrix obtained from:

$$\frac{1}{2}(A_{1}+A_{2})$$ and $$\frac{1}{2}(E_{1}+E_{2})$$     (where $$E_{1} \neq E_{2} $$)

Provides:

$$\frac{(E_{1}+E_2)(A_{1}+A_{2})}{4L^{i}}\begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix}= k_{ave}^{(i)}$$

 Find:  k (i) - k (i)ave


 * Recall the Mean Value Theorem (MVT) related to a centroid:

$$MVT: \int_{x=a}^{x=b}{f(x)dx=f(\bar{x})[b-a]}$$ $$\bar{x} \epsilon [a,b] $$ where $$a\leq \bar{x}\leq b$$
 * For:

$$\int_{a}^{}{X dA }=\bar{X}\int_{a}^{}{dA}=\bar{X}A$$

$$\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 most cases


 * The average value of  f :

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


 * The average value of  g :

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

Analysis of the Electric Pylon
This section analyzes an Electric Pylon with:

Number of Elements = 91

Material = 300M Steel

Young's Modulus of Each Element = 200GPa

Cross Sectional Area of Each Element = 4$${cm}^{2}$$

Required MATLAB Functions
Before running the 2-Bar Truss System program, the following three MATLAB program functions must be in the same directory. They are obtained from the textbook website (See HW1 for textbook information). Matlab Truss Files

PlaneTrussElement.m
This function takes the Modulus of Elasticity (e), the Cross-Sectional Area (A) and the Nodal Coordinates (coord) of the element and generates the local stiffness matrix.

NodalSoln.m
This function takes the Global Stiffness Matrix (K), the Global Force vector (R), and two inputs that are the boundary conditions. The output is an array of the displacement and reaction forces.

PlaneTrussResults.m
This function inputs the Modulus of Elasticity (e), Cross-Sectional Area (A), Nodes of Each Element (coord), and the Displacement calculated in the previous function (disps). The output are the Strain, Stress, and Axial Force of each element.

Main MATLAB Function
The following MATLAB program (electric.m) is the main function to break down the 6-Bar Truss System.

Results
The following are the results to the Electric Pylon

The displacement matrix of the Electric Pylon is as follows (Units are in millimeters (mm)):

$$\begin{Bmatrix} d_{1} \\ d_{2} \\ d_{3} \\ d_{4} \\ d_{5} \\ d_{6} \\ d_{7} \\ d_{8} \\ d_{9} \\ d_{10} \\ d_{11} \\ d_{12} \\ d_{13} \\ d_{14} \\ d_{15} \\ d_{16} \\ d_{17} \\ d_{18} \\ d_{19} \\ d_{20} \\ d_{21} \\ d_{22} \\ d_{23} \\ d_{24} \\ d_{25} \\ d_{26} \\ d_{27} \\ d_{28} \\ d_{29} \\ d_{30} \\ d_{31} \\ d_{32} \\ d_{33} \\ d_{34} \\ d_{35} \\ d_{36} \\ d_{37} \\ d_{38} \\ d_{39} \\ d_{40} \\ d_{41} \\ d_{42} \\ d_{43} \\ d_{44} \\ d_{45} \\ d_{46} \\ d_{47} \\ d_{48} \\ d_{49} \\ d_{50} \\ d_{51} \\ d_{52} \\ d_{53} \\ d_{54} \\ d_{55} \\ d_{56} \\ d_{57} \\ d_{58} \\ d_{59} \\ d_{60} \\ d_{61} \\ d_{62} \\ d_{63} \\ d_{64} \\ d_{65} \\ d_{66} \\ d_{67} \\ d_{68} \\ d_{69} \\ d_{70} \\ d_{71} \\ d_{72} \\ d_{73} \\ d_{74} \\ d_{75} \\ d_{76} \\ d_{77} \\ d_{78} \\ d_{79} \\ d_{80} \\ d_{81} \\ d_{82} \\ d_{83} \\ d_{84} \\ d_{85} \\ d_{86} \\ d_{87} \\ d_{88} \\ d_{89} \\ d_{90} \\ d_{91} \\ d_{92} \\ \end{Bmatrix}$$ = $$\begin{Bmatrix} 0 \\            0  \\            0  \\            0  \\     -0.24614  \\      0.21895  \\     -0.15704  \\   -0.0027926  \\     -0.21101  \\      -0.3432  \\     -0.21576  \\      0.38648  \\     -0.21763  \\     -0.14083  \\     -0.21754  \\     -0.62552  \\      0.30459  \\      0.35724  \\      0.33289  \\      -0.6888  \\      0.84792  \\      0.74729  \\      0.87754  \\     -0.24461  \\      0.91043  \\      -1.1024  \\       1.5418  \\       1.2099  \\       1.5409  \\      0.23911  \\         1.48  \\     -0.75546  \\       1.4796  \\      -1.4891  \\        2.174  \\      0.72439  \\       1.9197  \\      -1.1672  \\       4.5898  \\       2.8773  \\       3.4349  \\        -2.63  \\       6.0714  \\       6.5191  \\       6.0714  \\       4.4011  \\       6.0714  \\        3.522  \\       6.0738  \\       2.8633  \\       6.0862  \\       2.0639  \\       6.1061  \\       1.2663  \\       6.0417  \\      0.11623  \\        5.919  \\      -1.1413  \\       5.7852  \\      -2.4771  \\        5.643  \\      -3.7848  \\       5.4056  \\      -6.3546  \\       4.8496  \\      -18.249  \\       6.4456  \\       4.9206  \\        6.754  \\      -8.1536  \\        6.649  \\       4.0415  \\       7.1973  \\      -5.2394  \\       6.7724  \\       3.5162  \\       6.7648  \\       2.4998  \\       6.7485  \\       1.6336  \\       6.7999  \\      0.71154  \\       6.9414  \\     -0.44367  \\       7.1235  \\      -1.8243  \\       7.3567  \\      -3.8516  \\         8.42  \\       3.5162  \\       10.425  \\      -3.8516  \\ \end{Bmatrix}$$

The reactions at the legs of the Electric Pylons are as follows (Units are in Newtons (N)):

$$\begin{Bmatrix} R_1 \\ R_2 \\ R_3 \\ R_4 \\ \end{Bmatrix}$$ = $$\begin{Bmatrix} 183.78 \\     -501.65 \\      -183.78 \\       1501.7 \\ \end{Bmatrix}$$

In the results array, the left column is the Axial Strain, the middle column is the Axial Stress (MPa), and the right column is the Axial Force (N) for each element.

Looking at the results array, the element with the highest tensile and compressive stress are as follows

$$\sigma_{81}=10.981MPa$$

$$\sigma_{55}=-10.679MPa$$

Plot
The following MATLAB function is used to plot the Electric Pylon. It is a modified version of the function provided here.There are comments embedded in the function. The elements of highest compression and tension stresses are shown by arrows.

The following plot shows the undeformed and deformed shape of the Electric Pylon. The deformed displacement is magnified by a factor of 1000 to get a better look.



Statics Method
The Electric Pylon is Statically Determinate. Unlike previous 2-D problems, taking the moments about a point does not yield a 0 = 0. The first two equations are the sum of the forces in the x and y directions

$$\sum{F_{x}}=R_{1}+R_{3}=0$$

$$\sum{F_{y}}=R_{2}+R_{4}-1000=0$$

The third equations is the sum of the moments taken about Node 1. The convention is that positive are moments pointing into the pylon. The distances are calculated from the MATLAB plot.

$$\sum{M_{1}}=1000(6-1.488)+R_{4}(4.493-1.488)=0$$

From this:

$$\R_{4}=1501.5$$

$$\R_{2}=-501.5$$

Using the nodal locations of Element 1 for relations of the X and Y component of the force, it is calculated that:

$$\R_{1}=197.9$$

And since:

$$\R_{1}=-\R_{3}$$

$$\R_{3}=-197.9$$

The solution of the Statics Method only differ by a small margin versus the FEA Method.

Contributing Members
The following Team Wiki members contributed to this report.

Renee Hood --Eml4500.f08.wiki.hood 00:35, 18 November 2008 (UTC)

Rodney Dagulo --Eml4500.f08.wiki.dagulo.r 19:45, 20 November 2008 (UTC)

Gonzalo Barcia --Eml4500.f08.wiki.barcia 10:45, 21 November 2008 (UTC)

Thomas Kehoe --Eml4500.f08.wiki.Kehoe 12:12, 21 November 2008 (UTC)

Cortland Glowacki --Eml4500.f08.wiki.cort 12:57, 21 November 2008 (UTC)

Ricardo Lopez --Eml4500.f08.wiki.Lopez 4:12, 21 November 2008 (UTC)