User:Eml4500.f08.qwiki.bishop/HW6

=Lecture 29 - Monday, November 3, 2008=

Principle of Virtual Work
The initial conditions for the above elastic bar are described as follows:
 * At t=0, prescribe $$u(x,t=0)=\bar{u}(x)$$
 * $$\bar{u}(x)$$ is a known function of displacement
 * Also at t=0, $$\frac{\partial u}{\partial x}(x,t=o)=\dot{u}(x,t=0)=\bar{v}(x)$$
 * $$\bar{v}(x)$$ is a known function of velocity

PVW (continued) of dynamics of the elastic bar from Lecture 28 includes the following Partial Differential Equation of Motion:

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

Knowing that the Discrete Equation of Motion is as follows: $$-\mathbf{k}*\mathbf{d}+\mathbf{F}=\mathbf{M}*\mathbf{\ddot{d}}$$

Noticing that there is a one to one correspondence between the two equations, for example $$\mathbf{M}*\mathbf{\ddot{d}}$$ is similar to $$m(x)\ddot{u}$$, the above Discrete EOM can be rewritten as:

Equation (2) $$\mathbf{M}*\mathbf{\ddot{d}}+\mathbf{k}*\mathbf{d}=\mathbf{F}$$

Equation (2) is a multiple degree of freedom equation (MDOF).

Single Degree of Freedom Example
A single degree of freedom (SDOF) situation is shown below. The equation of motion is: $$m\ddot{d}+kd=F$$

In this situation, 2 initial conditions are needed since it is a 2nd order derivative.

Deriving Equation (2) from Equation (1)
To derive Equation (2) from Equation (1) the following integral is used.

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

Where W(x) is the weighting function and the integral is for all possible W(x).

Going from Equation (1) to Equation (3) is a trivial task, but when going the opposite direction, it can be quite challenging. Thus the following steps must be taken to provide the proper derivation.

Steps to solve for Equation (1)

 * 1) Equation (3) can be rewritten as: $$\int {W(x)g(x)dx=0  }$$  for all w(x)
 * 2) Since Equation (3) holds for all W(x), select W(x)=g(x)
 * 3) Equation (3) then becomes: $$\int {g^{2}dx}=0$$
 * 4) Solving Step 3 equation yields g(x)=0

Two-Bar Truss and Rectangular Truss Zero Eigenvalues Plot Revisited
=Lecture 30 - Wednesday, November 5, 2008=

Integration by Parts
Knowing the steps of integration by parts is necessary to continue the Principal of Virtual Work. A basic review is shown below.

$$ \left(rs\right)^{'}=r^{'}s+rs^{'} $$

$$r'=\frac{dr}{dx}$$

$$s'=\frac{ds}{dx}$$

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

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

Recall Continuous Principal of Virtual Work

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

Using integration by parts, the following term is derived:

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

$$=W(L)(EA)(L)\frac{\partial u}{\partial x}(L,t)-W(0)(EA)(0)\frac{\partial u}{\partial x}(0,t)-\int_{0}^{L}\frac{\partial W}{\partial x}(EA)\frac{\partial u}{\partial x}dx$$

Now consider the following model:



At x=0, select W(x) such that W(0)=0. Doing so makes the system kinematically admissble.


 * Motivation: Discrete Principal of Virtual Work applied the equation from lecture 10:

$$\mathbf{W}_{6x4}([\mathbf{K}]_{6x2}\begin{Bmatrix} d_{1}\\ d_{2} \end{Bmatrix}_{2x1}-\mathbf{F}_{6x1})=0_{1x1}$$

For all $$\mathbf{W}$$

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

=0$$ so to eliminate equations involving unknown reactions. Doing so will eliminate rows 1, 2, 5, and 6.
 * For this statement, $$F_1$$, $$F_2$$, $$F_5$$, and $$F_6$$ are unknown reactions.
 * Because $$\mathbf{W}$$ can be selected arbitrarily, select $$\mathbf{W}$$ such that $$W_1=W_2=W_5=W_6
 * The previous equation now reduces to the following equation:

$$\mathbf{K}_{2x2}\mathbf{d}_{2x1}=\mathbf{F}_{2x1}$$

Note that $$\mathbf{W}(\mathbf{K}\mathbf{d}-\mathbf{F})=0$$ for all  $$\mathbf{W}$$

Back to Continuous PVW

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


 * Using continuous PVW the following term is derived:

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

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


 * These terms yield the following final equation:

$$\int_{0}^{L}W(m\ddot{u})dx+\int_{0}^{L}\frac{\partial W}{\partial x}(EA)\frac{\partial u}{\partial x}dx=W(L)F(t)+\int_{0}^{L}Wfdx$$

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

=Lecture 31 - Friday, November 7, 2008=

Continuous vs. Discrete Settings

 * To begin the process of finding the stiffness term, note the following figure:
 * Where n equals the number of nodes.


 * Assume displacement u(x) for $$x\epsilon [x_i,x_{i+1}]$$

Motivation for linear interpolation of u(x), 2-bar truss

 * Deformed shape is a straight line, i.e. there was an implicit assumption of linear interpolation of displacement between 2 nodes.
 * Consider the case where there are only axial displacements. (i.e., zero transverse displacements, see p.31-2)


 * 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)=N_i(x)d_i+N_{i+1}(x)d_{i+1}$$
 * Where Ni(x) and Ni+1(x)  are linear functions of x.
 * Homework:
 * $$N_{i}(x)=\frac{x-x_{i+1}}{x_{i}-x_{i+1}}$$
 * $$N_{i+1}(x)=\frac{x-x_{i}}{x_{i+1}-x_{i}}$$



=Lecture 32 - Monday, November 10, 2008=

Rube-Goldberg Devices
The topic of the Rube-Goldberg device was introduced in lecture to stimulate students' imaginations.

Rube Goldberg: accomplishing by complex means what seemingly could be done simply - Merriam-Webster's Online Dictionary

Goldburger To Go!
The above screen shot is of the [http://www.donpixel.com/play/en/060302124616/ Goldberger To Go! by PBS Kids] website. The website contains a machine which needs to be fixed in order to make the Rube-Goldberg device work properly. In the set up there are 13 different parts that can be changed to get the right result. The user can change whichever part of the device they want and then run it to see where the errors are located.

After fiddling with the parts, clicking on the Tips link, shows the user the scientific and engineering aspects of each portion. For instance, scrolling over the pendulum portion describes how gravity causes it to move and about the inertia the movement contains.

Lagrangian Interpretation
Motivation for form of Ni(x) and Ni+1(x)

1. Ni(x) and Ni+1(x) are linear- (represented by straight lines), thus any linear combination of Ni and Ni+1 is also linear, and in particular expressions for u(x):

Algebra $$\mathbf{u(x)=N_{i}(x)d_{i}+N_{i+1}(x)d_{i+1}}$$ $$\mathbf{N_{i}(x)=\alpha _{i}+\beta _{i}x}\;\Rightarrow \;with\; \alpha _{i} \;and \;\beta _{i} = \;real\; numbers$$ $$\mathbf{N_{i+1}(x)=\alpha _{i+1}+\beta _{i+1}x}\;\Rightarrow \;with\; \alpha _{i+1} \;and \;\beta _{i+1} = \;real\; numbers$$

Looking at the linear combination of Ni(x) and Ni+1(x) :

$$\mathbf{N_{i}(x)d_{i}+N_{i+1}(x)d_{i+1}=(\alpha _{i}+\beta _{i}x)d_{i}+(\alpha _{i+1}+\beta _{i+1}x)d_{i+1}}$$

Therefore:

$$\mathbf{N_{i}(x)d_{i}+N_{i+1}(x)d_{i+1}=(\alpha _{i}d_{i}+\alpha _{i+1}d_{i+1})+(\beta _{i}d_{i}+\beta _{i+1}d_{i+1})x}$$

Verification
2. Recall equation for u(x) (interpretation of u(x)) shown above. At x = xi the equation becomes:

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

Here, $$\mathbf{N_{i}( x_{i} )d_{i}}$$ becomes 1 and $$\mathbf{N_{i+1}(x_{i})d_{i+1}}$$ becomes 0. Thus, at x = xi, u = di which is what we set out to prove.

=Lecture 33 - Wednesday, November 12, 2008=

FEM via PVW (continued)
As we saw in Lecture ___, the assumed solution u(x) for a two-node element is a linear interpolation between the nodal unknowns, and is expressed as $$ u(x) = N_i(x) d_i + N_{i+1}(x) d_{i+1}$$. To obtain a solution for the linear displacement at $$x_{i+1} $$, we must evaluate $$ u(x) $$ for $$ u(x_{i+1})$$. Doing so yields:
 * $$ u(x_{i+1}) = N_i(x_{i+1}) d_i + N_{i+1}(x_{i+1}) d_{i+1} = d_{i+1} $$
 * $$ u(x_{i+1}) = \Bigl(\frac{x - x_{i+1}}{x_i - x_{i+1}}*d_i + \frac{x - x_i}{x_{i+1} - x_i} * d_{i+1}\Bigr) \Bigl\vert_{x = x_{i+1}} $$
 * $$ u(x_{i+1}) = \Bigl(\frac{x_{i+1} - x_{i+1}}{x_i - x_{i+1}}*d_i + \frac{x_{i+1} - x_i}{x_{i+1} - x_i} * d_{i+1}\Bigr) = \Bigl(\frac{0}{x_i - x_{i+1}}*d_i + \frac{x_{i+1} - x_i}{x_{i+1} - x_i} * d_{i+1}\Bigr)

= 0 * d_i + 1 * d_{i+1} $$ Therefore, one can see that $$ u(x_{i+1}) = d_{i+1} $$.

Since the weighting function w(x) is a linear function similar to $$ u(x) $$, one can linearly interpolate between $$ x_i \text{ and } x_{i+1} $$ to solve for $$w(x) \text{ : } $$
 * $$ w(x) = \Bigl(\frac{x - x_{i+1}}{x_i - x_{i+1}}*w_i + \frac{x - x_i}{x_{i+1} - x_i} * w_{i+1}\Bigr)

= N_i(x)w_i + N_{i+1}(x)w_{i+1} $$.

Element Stiffness Matrix for Element i
In practice, most finite element solutions are written in matrix notation to ease understanding of the system. For simple elements such as the one pictured above from Lecture ___, matrix notation is not necessary. However, for complicated elements it is nearly impossible to express element solutions without using matrix notation. In this section, we will show how to transform from the continuous stiffness equation $$ \beta = \int_{0}^{L} \frac{\delta w}{\delta x}(EA)\frac{\delta u}{\delta x}dx $$ to a discrete stiffness matrix $$ \mathbf{w}*(\mathbf{K} \mathbf{d}) $$. To begin, we must first note that
 * $$ \frac{\delta w}{\delta x} = [N_i^{'}w_i + N_{i+1}^{'}w_{i+1}] \text { and }

\frac{\delta u}{\delta x} = [N_i^{'}d_i + N_{i+1}^{'}d_{i+1}] \text{.}$$

$$\text{Therefore, } \beta = \int_{0}^{L} \frac{\delta w}{\delta x}(EA)\frac{\delta u}{\delta x}dx  \text{ can be expressed as } $$
 * $$ \beta = \int_{x_i}^{x_{i+1}} [N_i^{'}w_i + N_{i+1}^{'}w_{i+1}](EA)[N_i^{'}d_i + N_{i+1}^{'}d_{i+1}]dx \text{, where } N_i^{'} = \frac {\delta N_i(x)}{\delta x} \text{ and } N_{i+1}^{'} = \frac {\delta N_{i+1}(x)}{\delta x}.$$

Note that u(x) and $$ \frac{\delta u(x)}{\delta x} \quad $$can also be expressed in matrix form, where
 * $$ u(x) =

\lfloor N_i(x) \quad N_{i+1}(x)\rfloor \begin{bmatrix} d_i \\ d_{i+1} \\ \end{bmatrix} = \mathbf{N(x)}_{1x2} \begin{bmatrix} d_i \\ d_{i+1} \\ \end{bmatrix} \text{, and }$$


 * $$\frac{\delta u(x)}{\delta x} =

\lfloor N_i^{'}(x) \quad N_{i+1}^{'}(x)\rfloor \begin{bmatrix} d_i \\ d_{i+1} \\ \end{bmatrix} = \mathbf{B(x)}_{1x2} \begin{bmatrix} d_i \\ d_{i+1} \\ \end{bmatrix} $$

Similarly, w(x) and $$ \frac{\delta w(x)}{\delta x} \quad $$can also be expressed in matrix form, where
 * $$ w(x) = \mathbf{N(x)}_{1x2}

\begin{bmatrix} w_i \\ w_{i+1} \\ \end{bmatrix} \text{, and } $$


 * $$ \frac{\delta w(x)}{\delta x} = \mathbf{B(x)}_{1x2}

\begin{bmatrix} w_i \\ w_{i+1} \\ \end{bmatrix} $$

The element DOFs for element i are shown in the picture below:

The picture above shows that: Therefore, we may express equivalent DOFs in matrix format:
 * Node 1 of element i is equivalent to node i
 * Node 2 of element i is equivalent to node i+1
 * $$ \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)} $$

By substituting equivalent values for $$\frac{\delta u(x)}{\delta x} \quad$$ and $$ \frac{\delta w(x)}{\delta x} \quad $$ into $$ \beta = \int_{x_i}^{x_{i+1}} \frac{\delta w}{\delta x}(EA) \frac{\delta u}{\delta x}dx \quad $$, the continuous stiffness equation reduces to
 * $$ \beta = \int_{x_i}^{x_{i+1}} (\mathbf{B} \mathbf{w}^{(i)})(EA)(\mathbf{B} \mathbf{d}^{(i)})dx $$

(Move this and/or expain) Goal: $$ \beta = \mathbf{w}^{(i)} * (\mathbf{k}^{(i)} \mathbf{d}^{(i)}) $$

Since $$(EA) $$, $$(\mathbf{B} \mathbf{w}^{(i)})  $$ , and $$(\mathbf{B} \mathbf{d}^{(i)}) $$ are scalars, we may rearrange the terms in any fashion. Doing so allows us to express the continuous stiffness equation as


 * $$ \beta = \int_{x_i}^{x_{i+1}} (EA)(\mathbf{B} \mathbf{w}^{(i)})(\mathbf{B} \mathbf{d}^{(i)})dx $$

Recall that by the definition of a transpose:
 * $$ (\mathbf{B} \mathbf{w}^{(i)})^T = (\mathbf{w}^{(i)T} \mathbf{B}^T) = (\mathbf{w}^{(i)} \mathbf{B}^T)$$
 * $$ (\mathbf{B} \mathbf{d}^{(i)})^T = (\mathbf{d}^{(i)T} \mathbf{B}^T) = (\mathbf{d}^{(i)} \mathbf{B}^T)$$

Due to the fact that $$\mathbf{w}^{(i)} $$ and $$ \mathbf{d}^{(i)} $$ are independent from x, we may remove them from the integral:
 * $$ \beta = \mathbf{w}^{(i)} * \Bigl (\int_{x_i}^{x_{i+1}} \mathbf{B}^T(EA) \mathbf{B} dx \Bigr) \mathbf{d}^{(i)} $$

From the above equation, one may easily see that the stiffness equation $$ \beta = \mathbf{w}^{(i)} * \Bigl(\int_{x_i}^{x_{i+1}} \mathbf{B}^T(EA) \mathbf{B} dx \Bigr) \mathbf{d}^{(i)} $$ reduces to $$ \beta = \mathbf{w}^{(i)} * (\mathbf{k}^{(i)} \mathbf{d}^{(i)})$$. Therefore, we have solved for the stiffness matrix $$ \mathbf{k}^{(i)} $$:
 * $$ \mathbf{k}^{(i)} = \int_{x_i}^{x_{i+1}} \mathbf{B}^T(EA) \mathbf{B} dx $$

It is also known that
 * $$ B(x) = \lfloor N_i^{'}(x) \quad N_{i+1}^{'}(x)\rfloor = \Bigl \lfloor \frac{1}{x_i - x_{i+1}} \quad \frac{1}{x_{i+1} - x_i} \Bigr \rfloor = \Bigl \lfloor \frac{1}{-L^{(i)}} \quad \frac{1}{L^{(i)}} \Bigr \rfloor \quad$$, where the length of element i equals $$\quad L^{(i)} = x_{i+1} - x_i $$

Transformation of Element Coordinate System
=Lecture 34 - Friday, November 14, 2008=

by: Eric Gatch

from p. 31-4

$$\mathbf{N}_{1}^{(i)}(\tilde{x}) = HW 6$$

and

$$\mathbf{N}_{2}^{(i)}(\tilde{x}) = \frac{\mathbf\tilde{x}}{{L}^{(i)}} - \begin{cases} 0 \; at \; \tilde{x}=0\\ 1 \; at \; \tilde{x}={L}^{(i)} \end{cases}$$

Shaper Functions (basis) $${N}_{1}^{(i)} and {N}_{2}^{(i)} $$

Homework 6 - Comparing the Stiffness Matrix from the notes to that from the Textbook

Part A

Note: When 'The Text' is mentioned in Part A refer to Page 159

Set $${E}_{1} = {E}_{2} = {E}$$

Let $$ A(\tilde{x})$$ be linear as in p. 33-5 so that $$\frac{{A}_{1} + {A}_{2}}{2}$$ is the average area.

Obtain $$\mathbf{\underline{k}^{(i)}}$$ from the previous homework problem and compare it to the stiffness matrix found in the book.

Previously we found that

$$\mathbf{\underline{k}^{(i)}}(\tilde{x}) = $$

And from 'The Text' we see that

$$\frac{E}{{L}^{(i)}}\frac{{A}_{1} + {A}_{2}}{2}\begin{bmatrix} 1 & -1 \\ -1 & 1\end{bmatrix} = \mathbf{\underline{{k}^{(i)}}}$$

Part B

Next compare the general $$\underline{{k}^{(i)}} $$ found above to the stiffness matrix obtained using $$ \frac{1}{2}{({A}_{1} + {A}_{2})} $$ and  $$ \frac{1}{2}{({E}_{1} +{E}_{2})} $$

Note: $$ {E}_{1} \neq {E}_{2} $$

This gives. ..

$$\underline{{k}_{ave}^{(i)}} = \frac{({E}_{1} + {E}_{2})({A}_{1} + {A}_{2})}{4*{L}^{(i)}}\begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}$$

Next. . . Find $$\underline{{k}^{(i)}} - \underline{{k}_{ave}^{(i)}}$$

Recall the Mean Value Theorem (MVT) and it's relation to the centroid:

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

for $$ a <= x <= b $$



Centroid: $$\int_{A} x * dA = \tilde{x}\int_{A} dA = \tilde{x} * A$$

so. ..

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

but. ..

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

Modify the Two Bar Truss Code to Accommodate the General Stiffness Matrix $$ \underline{k}^{(i)} $$ on page 33-5



Note: the only matlab function to change is the one that computes the element stiffness matrix.

=MATLAB Electric Pylon Problem=

Results

 * The highest tensile stress occurs at element number 77
 * The highest compressive stress occurs at element 41
 * The problem is not statically indeterminate because there numerous bars contributing to the reaction force, and there are not enough equations to solve for the unknowns.
 * Results (Global Node Displacements, Reactions, Axial Strain, Axial Stress, Axial Force)

=Contributing Members=
 * Mike Berry - Eml4500.f08.qwiki.berry 21:54, 17 November 2008 (UTC)
 * Andrea Booher - EML4500.f08.qwiki.booher 00:28, 19 November 2008 (UTC)
 * Chris Bishop - Eml4500.f08.qwiki.bishop 01:33, 19 November 2008 (UTC)
 * David Nobles - Eml4500.f08.qwiki.nobles 16:14, 20 November 2008 (UTC)
 * Eric Gatch - Eml4500.f08.qwiki.gatch 17:57, 21 November 2008 (UTC)
 * Lee Harris - Eml4500.f08.qwiki.harris 19:59, 21 November 2008 (UTC)