text:
User:Eas4200c.f08.carbon.orear/HW4/Matlab

 Comment: This page is identical to the previous one, except that the MATLAB section has been modified with additional work in the torsional constant calculation and centroid calculation. Here is a comparison between these two versions. Eas4200c.f08.carbon.orear 19:09, 7 November 2008 (UTC)

= Code =

The airfoil now has a centroid calculation missing from Homework 3. Instead of resubmitting HW 3, these changes are here in the new code.

AirfoilArea.m
= Debugging =

Circle
With a circular airfoil, the Number of Segments such that the computed Area of the airfoil is accurate to within 1% of the exact value is 23 segments at which the computed area was 3.1115 m3. This is 0.95% away from the actual area of a unit circle, 3.14195 m3. The Centroid of the 8 segment airfoil was found to be (1.0000, 4.2743e-018).

Code used to find ns_max

The area for a various number of segments is tabulated below.

Plotting the Area vs Number of segments up to the point of 1% accuracy we get the graph below.

At P0 = (0,0), Area computed equaled 3.1115

At P0 = (1,0), Area computed equaled 3.1115

At P0 = (1,-4), Area computed still equaled 3.1115

NACA Airfoil
With a circular airfoil, the Number of Segments such that the computed Area of the airfoil is accurate to within 1% of the previous value is 8 segments at which the computed area was 0.024616 m3. This is 0.77% away from the previous area of a unit circle, 0.024373 m3. At 10,000 segments, the computed area was 0.02569 m3. The Centroid of the 8 segment airfoil was found to be ( 0.2167, 0.0078).

Code used to find ns_max

The area for a various number of segments is tabulated below.

Plotting the Area vs Number of segments up to the point of 1% difference accuracy we get the graph below.

At P0 = (0,0), Area computed equaled 0.024616

At P0 = ( 0.2167, 0.0078) (Centroid), Area computed equaled 0.024616

At P0 = (.25,-1), Area computed still equaled 0.024616

=Torsion=

Single Cell
With a NACA 2415 airfoil, the Number of Segments such that the computed torsional constant (J) of the airfoil is accurate to within 1% of the previous value is 12 segments at which the computed area torsional constant (J) 4.7209 x 106 m4.

Code used to find ns_max

Multiple Cell
In the multiple cell case, using the code below, the torsional constant was found to be higher at torsional constant (J) 4.7518 x 106 m4, slightly higher then the single cell airfoil. The area summation of each cell was the same at 0.0251 m2


text:
Globulars in M31

M31, the Great Andomeda Galaxy

The first recorded observation of this object as a fuzzy patch was by the Persian astronomer, Abd al-Rahman al-Sufi. The nature of this "little cloud" has been the source of speculations throughout history. The answer was partially revealed shortly after the invention of the telescope which allowed Simon Marius to study it in 1612.

As a nebulousity it was cataloged in the 19th century and many other similar 'spiral nebulas' were discovered.

Photography institued a new era in the exploration of such objects and in 1917 Heber Curtis noticed that nova from the distant object were 10 magnitudes fainter making the object 500,000 light-years away. He hypothesized that such objects were "island universes" and were large copies of the stars, dust and gas that we called the milky way.

This remained controversial for a full decade, until Edwin Hubble, using the large palomar telescope, resolved some of its variable stars. As the timing of changes in the variable star was conclusively related to its brightness, he was able to conclusively show that M31, was indeed far beyond this galaxy. We now know it to be 2.2 million light-years away.

A globular cluster

A globular cluster is a very old and tightly bound group of stars. A globular cluster may have some million stars and therefore it shines very brightly. But on the other hand, globular clusters are found in the Halo around a galaxy and thus are much further away than many objects seen in the telescopes (galaxies being an obvious exception). The milky way galaxy has 158 known globular clusters, however many galaxies have far more. The Andromeda galaxy is suspected to have 500.

Observing globulars in M31
It is amazing to see any star groups 2 million light-years away, yet nevertheless globular clusters are a definite (but demanding observing target for a 10" reflecting telescope). At this distance, all sense of it being a cluster is lost and it shines as one of M31 brighter stars.  While some maps show many hundreds, many of which can be imaged with modest equiptment, the limit for a 10" scope is somewhere around 6-12.

The following sketch was carefully made of all field stars I saw during a clear night with a 10" dobsonian. Two globulars are shown, a third was clearly missed (having a brightness similar to other observed field stars)

Please feel free to attach any sketches of images which capture M31 globulars.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

title: Globulars in M31


text: File:102308 globM31.jpg

title: File:102308 globM31.jpg


text:
Talk:Breton/Breton pronunciation

IPA transliteration
Hello, would you have much difficulty appending full IPA transliteration, in addition to the simplified one, or instead of it? Most of IPA symbols are intuitive anyhow and moreover it could become even more helpful for the users with deeper knowledge of IPA. Thanks --Magairlin 13:45, 24 October 2008 (UTC) ==> OK, provided I can find how to type in IPA symbols. Any help welcome !--AC'hM 18:27, 26 October 2008 (UTC)

title: Talk:Breton/Breton pronunciation


title: User:Eml4500.f08.echo.sell.cm/HW4


text: File:Connectivity.JPG

title: File:Connectivity.JPG


text: File:LMM.JPG

title: File:LMM.JPG


title: File:Connectivity1.JPG


text: File:2 bar eig vectors.JPG

title: File:2 bar eig vectors.JPG


text: Introduction to US History/United We Stand:The History of the United States of America/Foreign Exploration


title: User:Eml4500.f08.team.guan/fixed101


text: File:IMG Katha.JPG

title: File:IMG Katha.JPG


title: File:2cell.jpg


text: File:Matrixa.jpeg

title: File:Matrixa.jpeg


text:
User talk:Uf.team.aero/HW5

Nov. 3 goes on HW6.

title: User talk:Uf.team.aero/HW5


title: File:Fig1101new.JPG


text:
User:Eml4500.f08.gravy.jad/NotesRedo3

 Comment: The differences between this page and the original page are the additions in the student contribution tables added to the main text, and an equation was corrected and added in a table. Here is a link to the original page HW3 report original And here is a link to the differences between the two versions. Comparison Eml4500.f08.gravy.mmm 23:24, 5 November 2008 (UTC)

FEA Notes and Homework 3

Force Displacement
Derivation of element Force Displacement with respect to global coordinate system.

$$\ k^{(e)}_{4x4} d^{(e)}_{4x1}=f^{(e)}_{4x1} \,$$ equation 1

The figure is labeled with the correct conventions to aid with the solving method. The four forces and displacements are collapsed to two forces and displacements by combining the x and y directions into a single vector.

Force Displacement matrices calculations
The force displacement can be written as

$$ \ K^{(e)} * \begin{bmatrix} 1 & -1 \\ -1 &1 \end{bmatrix}_{\hat{K}_{2x2}^{(e)}}* \begin{Bmatrix} q_1^{(e)} \\ q_2^{(e)} \end{Bmatrix}_{q_{2x1}^{(e)}} = \begin{Bmatrix} p_1^{(e)} \\ p_2^{(e)} \end{Bmatrix}_{p_{2x1}^{(e)}} \,$$ equation 2

$$\ q_i^{(e)} = \,$$ axial disp.of element at local node i

$$\ p_i^{(e)} = \,$$ axial disp.of element at local node i

The force displacement formula simplifies to the above and the unknowns can be solved for.

Goal: Derive Equation 1 from Equation 2

Want the relationship between $$\underline {q^{(e)}}$$ and $$\underline {d^{(e)}}$$ also $$\underline {p^{(e)}}$$ and $$\underline {f^{(e)}}$$

These relationships can be expressed in the form:

$$\ \underline {q^{(e)}_{2x1}}= \underline { T^{(e)}_{2x4}} * \underline {d^{(e)}_{4x1}} \,$$

To relate the q and d matrices and the p and f matrices a new matrix must be introduced. This matrix is the T matrix.

Consider the displacement vector of local node 1 denoted by

$$ \vec{d^{(e)}_{1}} $$

$$ \vec{d^{(e)}_{1}}= \vec{d^{(e)}_{1}}* \vec{i} + \vec{d^{(e)}_{2}} *\vec{j} $$

$$\ q_i^{(e)} = \,$$ axial displacement of node 1 is the orthogonal projection of the displacement.

Vector $$ \vec{d^{(e)}_{1}} $$ of node 1 on the axis $$\tilde{x}$$ of element $$e^{(1)}$$

$$\ = ( d_1^{(e)}*\vec{i} + d_2^{(e)}*\vec{j}) \cdot \vec{\tilde{i}} \,$$

$$\ = d_1^{(e)}* (\vec{i} \cdot \vec{\tilde{i}}) + d_2^{(e)}*(\vec{j} \cdot \vec{\tilde{i}} ) \,$$

$$\ (\vec{i} \cdot \vec{\tilde{i}} )= cos \theta^{(e)} = l^{(e)} \,$$

$$\ (\vec{j} \cdot \vec{\tilde{i}} )=sin \theta^{(e)} = m^{(e)} \,$$

$$\ q_1^{(e)}=l^{(e)}*d_1^{(e)} + m^{(e)}* d_2^{(e)} \,$$

$$\ = \begin{bmatrix} l^{(e)} & m^{(e)} \end{bmatrix}_{1x2} * \begin{Bmatrix} d_1^{(e)} \\ d_2^{(e)} \end{Bmatrix}_{2x1} \,$$

Similarly for node 2 : HW.
$$\ q_2^{(e)}=l^{(e)}*d_1^{(e)} + m^{(e)}* d_2^{(e)}= \begin{bmatrix} l^{(e)} & m^{(e)} \end{bmatrix}_{1x2} * \begin{Bmatrix} d_3^{(e)} \\ d_4^{(e)} \end{Bmatrix} \,$$

$$ \begin{Bmatrix}q_1^{(e)} \\ q_2^{(e)} \end{Bmatrix}_{\underline{q_{2x1}^{(e)}}} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0& 0\\0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix}_{\underline{T_{2x4}^{(e)}}}* \begin{Bmatrix} d_1^{(e)} \\ d_2^{(e)} \\d_3^{(e)} \\ d_4^{(e)} \end{Bmatrix}_{\underline{d_{4x1}^{(e)}}}$$

$$\ \underline{ q_{2x1}^{(e)} }= \underline{ T_{2x4}^{(e)} }* \underline{ d_{4x1}^{(e)} } \,$$

Similarly (same argument) :

$$\ \begin{Bmatrix} P_1^{(e)}\\P_2^{(e)} \end{Bmatrix}_{\underline{ P_{2x1}^{(e)} } } = \underline{ T_{2x4}^{(e)} } * \begin{Bmatrix} f_1^{(e)} \\ f_2^{(e)} \\ f_3^{(e)} \\ f_4^{(e)} \end{Bmatrix} _{\underline{ f_{4x1}^{(e)} } } \,$$

$$\underline{ P^{(e)} } = \underline{ T^{(e)}* f^{(e)} }$$

Recall element axial Force Displacement relation:

$$\ \underline{\hat{k}_{2x2}^{(e)}} * \underline{q_{2x1}^{(e)}}= \underline{P_{2x1}^{(e)}} \,$$

$$\underline{\hat{k}^{(e)}} *( \underline{ T^{(e)} } * \underline{ d^{(e)} })_{\underline{q^{(e)}}} = ( \underline{ T^{(e)} } * \underline{ f^{(e)} })_{\underline{P^{(e)}}}$$

Goal: Want to have $$\ \underline{k^{(e)}} * \underline{d^{(e)}}= \underline{f^{(e)}} \,$$

So "move " $$\underline{ T^{(e)}}$$ from the right hand side to the left hand side by multiplying the equation by $$ \underline{ T^{(e) -1}}$$ (inverse of $$\underline{ T^{(e)}} $$)

Unfortunately, $$\underline{ T^{(e)}}$$ is a rectangular matrix (2x4) so you can not invert it.

Ans: $$[\underline{ T_{4x2}^{(e)T} } * \underline{\hat{k}_{2x2}^{(e)}} * \underline{ T_{2x4}^{(e)}} ]_{is the matrix we are looking for} * d_{4x1}^{(e)}= f_{4x1}^{(e)}$$

$$\ \underline{k^{(e)} } * \underline{d^{(e)} }= \underline{f^{(e)} } \,$$

1) $$\ \underline{k^{(e)} } = \underline{T^{(e)T}} * \underline{\hat{k}^{(e)} }* \underline{T^{(e)} } \,$$

HW: verify above.
Verify the following equation:

$$\ \underline{k}^{(e)} = \underline{T}^{(e)T} * \underline{k}^{(e)} * \underline{T}^{(e)}$$

Where: $$\ \underline{T}^{(e)} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix}$$ and $$\ \underline{k}^{e)} = k^{(e)} * \begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix}$$

Since $$k^{(e)}$$ is a scaler, we can move it within the equation without changing the problem. So after moving the scaler and inserting in the matrices, the following equation is found:

$$\ \underline{k}^{(e)} = k^{(e)} * \begin{bmatrix} l^{(e)} & 0\\ m^{(e)} & 0\\ 0 & l^{(e)}\\ 0 & m^{(e)}\end{bmatrix} * \begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix} * \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix}$$

Solving the first set of matrices gives:

$$\ \underline{k}^{(e)} = k^{(e)} * \begin{bmatrix} l^{(e)} & -l^{(e)}\\ m^{(e)} & -m^{(e)}\\ -l^{(e)} & l^{(e)}\\ -m^{(e)} & m^{(e)}\end{bmatrix} * \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix}$$

Solving the second set of matrices gives:

$$\ \underline{k}^{(e)} = k^{(e)}\begin{bmatrix} (l^{(e)})^2 & l^{(e)}m^{(e)} & -(l^{(e)})^2 & -(l^{(e)}m^{(e)})  \\  l^{(e)}m^{(e)} & (m^{(e)})^2 & -(l^{(e)}m^{(e)}) & -(m^{(e)})^2  \\ -(l^{(e)})^2 & -(l^{(e)}m^{(e)}) & (l^{(e)})^2 & l^{(e)}m^{(e)}\\ -(l^{(e)}m^{(e)}) & -(m^{(e)})^2 & l^{(e)}m^{(e)} & (m^{(e)})^2 \end{bmatrix} $$

Which is the definition of the stiffness matrix, therefore the equation is true.

Justification for 1) : PVW (later see 10-1 for fist approach of PVW Reduction of Global Force Displacement relation)

$$ \underline{K_{6x6}}* \underline{d_{6x1}}= \underline{F_{6x1}} which reduces to \underline{K_{2x2}}* \underline{d_{2x2}}= \underline{F_{2x2}}$$

Remark: Why not solve as follows

$$ \underline{d} = \underline{K^{-1}}* \underline{F} $$ ?

The reason is that K can't be unversed because it is singular

Singular matrix has a determinate of 0 so it is not invertible

Recall: need to compute $$\frac{1}{det(K)}$$ to find K-1

Why? For an unconstrained structure system there are three possible rigid body motion in 2-D (2 translation, and 1 rotation)

HW: Eigenvalues of the Global Stiffness Matrix K
For an unconstrained structure system, there are three possible rigid body motionsin 2-D (2 translation and 1 rotation).

Dynamic Evaluation Problem

 K v = λ Mv

$$\ \underline{k} = \begin{bmatrix} K_{11}^{(1)} & K_{12}^{(1)} & K_{13}^{(1)} & K_{14}^{(1)} & 0 & 0  \\ K_{21}^{(1)}& K_{22}^{(1)} & K_{23}^{(1)} & K_{24}^{(1)} & 0 & 0 \\ K_{31}^{(1)} & K_{32}^{(1)} & (K_{33}^{(1)} + K_{11}^{(2)}) & (K_{34}^{(1)} + K_{12}^{(2)}) & K_{13}^{(2)} & K_{14}^{(2)} \\ K_{41}^{(1)}& K_{42}^{(1)} & (K_{43}^{(1)} + K_{21}^{(2)}) & (K_{44}^{(1)} + K_{22}^{(2)}) & K_{23}^{(2)} & K_{24}^{(2)} \\ 0 & 0 & K_{31}^{(2)} & K_{32}^{(2)} & K_{33}^{(2)} & K_{34}^{(2)} \\ 0 & 0 & K_{41}^{(2)} & K_{42}^{(2)} & K_{43}^{(2)} & K_{44}^{(2)}  \end{bmatrix}_{6x6}\,$$

$$\ \ = \begin{bmatrix} \frac{9}{16}& \frac{3\sqrt{3}}{16} & -\frac{9}{16} & -\frac{3\sqrt{3}}{16} & 0 & 0 \\  \frac{3\sqrt{3}}{16} & \frac{3}{16} & -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & 0 & 0 \\ -\frac{9}{16} & -\frac{3\sqrt{3}}{16} & \frac{49}{16} & \frac{3\sqrt{3} - 40}{16} & -\frac{5}{2} & \frac{5}{2} \\ -\frac{3\sqrt{3}}{16} & -\frac{3}{16} & \frac{3\sqrt{3} - 40}{16} & \frac{43}{16} & \frac{5}{2} & -\frac{5}{2} \\  0 & 0 & -\frac{5}{2} & \frac{5}{2} & \frac{5}{2} & -\frac{5}{2} \\ 0 & 0 & \frac{5}{2} & -\frac{5}{2} & -\frac{5}{2} & \frac{5}{2}  \end{bmatrix}_{6X6}\,$$

MATLAB Code

Values of the stiffness matrix are assigned to matrix K >> K= [9/16 (3*sqrt(3))/16 -9/16 -(3*sqrt(3))/16 0 0; (3*sqrt(3))/16 3/16 -(3*sqrt(3))/16 -3/16 0 0; -9/16 -(3*sqrt(3))/16 49/16 ((3*sqrt(3)-40)/16) -5/2 5/2; -(3*sqrt(3))/16 -3/16 ((3*sqrt(3)-40)/16) 43/16 5/2 -5/2; 0 0 -5/2 5/2 5/2 -5/2; 0 0 5/2 -5/2 -5/2 5/2]

Values of K in decimal form 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

Matlab function to calculate eigenvalues >> d = eig(K)

List of the eigenvalues for global stiffness matrix K d = -0.0000  -0.0000    0.0000    0.0000    1.4705   10.0295

λ = 0, 0, 0, 0, 1.4705, 10.0295

It is important to note that four of the six eigenvalues are zero. These values represent zero stored elastic energy (rigid body at the modes).

//(end Remark)

HW: FD Rel. on p. 10-1
$$(\underline{k}_{6x2}) (\underline{d}_{2x1}) = \underline{F}_{6x1} $$

$$\ \begin{bmatrix} K_{13}^{(1)} & K_{14}^{(1)} \\ K_{23}^{(1)} & K_{24}^{(1)} \\ (K_{33}^{(1)} + K_{11}^{(2)}) & (K_{34}^{(1)} + K_{12}^{(2)}) \\ (K_{43}^{(1)} + K_{21}^{(2)}) & (K_{44}^{(1)} + K_{22}^{(2)}) \\ K_{31}^{(2)} & K_{32}^{(2)} \\ K_{41}^{(2)} & K_{42}^{(2)} \end{bmatrix}_{6x2}\,$$ $$\ \begin{bmatrix} d_3 \\ d_4 \end{bmatrix}_{2x1}\,$$ $$\ \ = \begin{bmatrix} F_1 \\ F_2 \\ F_3 \\ F_4 \\ F_5 \\ F_6 \end{bmatrix}_{6x1}\,$$

$$\ \begin{bmatrix} -\frac{9}{16} & -\frac{3\sqrt{3}}{16} \\ -\frac{3\sqrt{3}}{16} & -\frac{3}{16} \\ \frac{49}{16} & \frac{3\sqrt{3} - 40}{16}\\ \frac{3\sqrt{3} - 40}{16} & \frac{43}{16} \\ -\frac{5}{2} & \frac{5}{2} \\ \frac{5}{2} & -\frac{5}{2} \end{bmatrix}_{6X2}\,$$ $$\ \begin{bmatrix} 4.3527 \\ 6.1271 \end{bmatrix}_{2x1}\,$$ $$\ \ = \begin{bmatrix} F_1 \\ F_2 \\ F_3 \\ F_4 \\ F_3 \\ F_4 \end{bmatrix}_{6x1}\,$$

$$F_1 = (-\frac{9}{16})(4.3527) + (-\frac{3\sqrt{3}}{16})(6.1271) = -4.4382$$

$$F_2 = (-\frac{3\sqrt{3}}{16})(4.3527) + (-\frac{3}{16})(6.1271) = -2.5624$$

$$F_5 = (-\frac{5}{2})(4.3527) + (\frac{5}{2})(6.1271) = 4.436$$

$$F_5 = (\frac{5}{2})(4.3527) + (-\frac{5}{2})(6.1271) = -4.436$$

Note: F3 and F4 were calculated in previous section to be 0 and 7 respectively

Closing the loop between the FEM and statics: Virtual Displacement
There are two different ways on solving these types of problems: using FEM, or using statics.

Element 1:

Element 2:

These pictures are FBD's of 2 force bodies of element 1 & 2. By Statics reactions are known and therefore the member forces $$\ [P_{1}^{(1)}, P_{2}^{(2)}] \,$$ are equal and opposite in direction. We need to compute the axial displacement dofs ( amount of extension of bars):

$$\ q_{2}^{(1)}= P^{(1)}/k^{(1)}= \underline{AC} \,$$ . ( $$\ q_{1}^{(1)}=0 \,$$ because it's fixed at global node 1.

$$\ q_{1}^{(2)}= -P^{(2)}/k^{(2)}= -\underline{AC} \,$$ . ( $$\ q_{2}^{(2)}=0 \,$$ because it's fixed at global node 3.

Q : How to back out from above results of displacement of dofs of global node 2?

MATLAB code for truss displacement

The above picture is the output graph from MATLab after the above code was run. The solid line represents the deformed truss, while the dashed lines show the initial, or un-deformed, truss.

$$\ u_{y}=R*sin(\alpha) \simeq R*\alpha \,$$

$$\ u_{x}=R*(1-cos(\alpha)) \simeq 0 \,$$ $$\ (cos(\alpha) \simeq 1) \,$$

hence

Here we are assuming that alpha is very small. When we make that assumption sin(alpha) becomes alpha, and cos(alpha) becomes one. The picture shown above is exaggerated, and alpha is very small.

Infinitesimal Displacement (Related to Virtual Displacement)
We can solve for the displacement vectors AC and AB by first referencing all displacement from point A which will be (x,y)=(0,0) on the coordinate system. The actual displacement vector can be defined by the x and y displacement vectors in the form of d(3)(x-coordinates) and d(4)(y-coordinates). Points B and C are the axial displacement vectors in regards to element 1 and element 2.

$$\ AC= \frac{|P_{2}^{(1)}|}{k^{(1)}} = \frac{5.1243}{3/4} = 6.8324 \,$$

HW: Solving AC and AB
From previous lecture

Element Length

$$\ L^{(1)}=4 \,$$

$$\ L^{(2)}=2 \,$$

Young's Modules

$$\ E^{(1)}=3 \,$$

$$\ E^{(2)}=5 \,$$

Cross-sectional Area

$$\ A^{(1)}=1 \,$$

$$\ A^{(2)}=2 \,$$

Inclination angle

 * $$\ \theta^{(1)} = 30^o \,$$

 * $$\ \theta^{(2)} = -45^o \,$$

 * $$\ \theta^{(1)} = 30^o \,$$

Since $$ \ P_{(e)} = T_{(e)} * f_{(e)}\,$$

And the T matrix is

$$\ T^{(e)} = \begin{bmatrix} l^{(e)} & m^{(e)}& 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)}\end{bmatrix} \,$$

we can solve for P

The f matrix was solved in a previous lecture and is:

$$\ F= \begin{bmatrix} -4.4378 \\ -2.5622 \\ +4.4378 \\ +2.5622 \end{bmatrix}_{4X1}\,$$

l m and k were also solved in a previous lecture and are:

$$\ l^{(1)}=cos\theta^{(1)}= cos 30^o= \frac{\sqrt{3}}{2} \,$$

$$\ m^{(1)}=sin\theta^{(1)}= sin 30^o= \frac{1}{2} \,$$

$$\ k^{(1)}= \frac{E^{(1)}A^1}{L^1} \frac{\left(3\right ) \left(1\right )}{4}= \frac {3}{4}\,$$

$$\ l^{(2)}=cos\theta^{(2)}= cos -45^o= .707 \,$$

$$\ m^{(2)}=sin\theta^{(2)}= sin -45^o= -.707 \,$$

$$\ k^{(1)}= \frac{E^{(5)}A^1}{L^1} \frac{\left(2\right ) \left(1\right )}{2}= 5\,$$

$$T^{(1)}$$ and $$T^{(2)}$$ are solved as the following:

$$\ T^{(1)} = \begin{bmatrix} \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2}\end{bmatrix} \,$$

$$\ T^{(2)} = \begin{bmatrix} .707 & -.707 & 0 & 0 \\ 0 & 0 & .707 & -.707\end{bmatrix} \,$$

With T and f the P s were found as being,

$$P_{2}^{(1)}$$= $$P_{1}^{(2)}$$=5.1243

so

$$\ AC= \frac{|P_{2}^{(1)}|}{k^{(1)}} = \frac{5.1243}{3/4} = 6.8324 \,$$

$$\ AC= \frac{|P_{1}^{(2)}|}{k^{(2)}} = \frac{5.1243}{5} = 1.0249 \,$$

HW sove (x,y) for points B and C
Knowing that the angle to find x_b and y_b is 30 degrees, and to find x_c and y_c is 135 degrees, we can use simple trigonometry to solve for the values:

$$\ x_i = L*cos(\theta), yi = L*sin(\theta)\,$$

$$\ x_b = AC*cos(30)= 6.8324 * cos(30) = 5.9170, y_b = AC*sin(30) = 6.8324 * sin(30) = 3.4162\,$$

$$\ x_c = AB*cos(135) = 1.0249 * cos(135) = -.7247, y_c = AB*sin(135) = 1.0249 * sin(135) = .7247 \,$$

(xb,yb) = (5.9170, 3.4162) (xc,yc) = (-0.7247, 0.7247)

- P is fixed, Q is "sliding" on the line

- the two lines are parallel $$\ \bar{PQ} =(PQ)*\vec {\tilde i} = (PQ)[cos\theta \vec {i} + sin\theta  \vec {j}] = (x-x_p)\vec i + (y-y_p)\vec j \,$$

$$\ \Rightarrow x-x_p=PQ cos\theta \,$$

$$\ \Rightarrow y-y_p=PQ sin\theta \,$$

$$\ \frac{y-y_p}{x-x_p}= tan\theta \,$$

$$\ (y-y_p)= (tan\theta)(x-x_p) \,$$

eq. for a line perpendicular to the line, passing P: $$\ (y-y_p)= \tan(\theta+\frac{\pi}{2})(x-x_p) \,$$

HW solve (x,y) for point D
To find x_d and y_d, the following equation is used:

$$\vec{AD} = x_d\vec{i}+ y_d\vec{j}$$

From HW2, we know that:

$$\vec{AD} = d_3\vec{i} + d_4\vec{j}$$

Also from HW2, we know that d3=4.352 and d4=6.1271. This means that:

$$\vec{AD} = 4.352\vec{i} + 6.1271\vec{j}$$

hence, (xd,yd) = (4.352, 6.1271).

since $$\ x_A=0 and y_A=0 \,$$ then $$\ \vec {AD} =(x_D-x_A)\vec i + (y_D - y_A)\vec j = x_D-\vec i + y_D\vec j \,$$

By definition: $$\ \vec{AD} = d_3 \vec i + d_4 \vec j \,$$ (these are displacement vectors of A from FEM

///end closing loop

3-Bar Truss Syst:

The FEA method of solving for a 3-Bar Truss system is done similarly to the 2-Bar Truss system except with an extra element. A 3-Bar Truss system cannot be solved though statics due to the fact that there are too many unknowns.

$$\ E^{(1)} =2, E^{(2)} =4, E^{(3)} =3 \,$$

$$\ A^{(1)} =3, A^{(2)} =1, A^{(3)} =2 \,$$

$$\ L^{(1)} =5, L^{(2)} =5, L^{(3)} =10 \,$$

$$\ P=30, \theta^{(1)}=30^\circ = -\theta^{(2)}, \theta^{(3)}=45^\circ \,$$

Local node numbering is equivalent to the convenience in assembly k.

Convenient local node numbering

$$ \sum{F_x} =0$$

$$ \sum{F_y} =0$$

$$ \sum{M_A} =0$$ trivial

The sum of the moments about A equals 0 so it is trivial. This means that there are three unknowns and only two equations. So the problem can not be solved using statics and must be solved using FEA.

Question: How about $$\sum{M_B}?$$ It still equals 0

$$ \sum{M_B} = \vec{BA'}\times \vec{F}$$

for A' on line of action of $$\vec{F}$$

$$\vec{BA'}= \vec{BA} + \vec{AA'}$$

$$ \sum{M_B}=(\vec{BA}+ \vec{AA'})\times \vec{F}$$ $$ = \vec{BA}\times F + \vec{AA'}\times \vec{F}= 0 $$

Back to 3-bar truss

Node A is in equilibrium : $$\sum_{i=0}{\vec{F_i}} =\vec{0}$$

$$ \sum_i{M_{B,i}}=\sum_i{\vec{BA'_i}}\times \vec{F_i}$$

A'i= any point on line of action of $$\vec{F_i}$$

$$ \sum_i{M_{B,i}}=\sum_i{\vec{BA}}\times \vec{F_i} \vec{BA} \times \sum_i{\vec{F_i}} = 0 $$

Converting from k matrices to K matrix

This figure is a representation of the k matrices forming the K matrix. The third element k matrix is split so that the ks corresponding to node 2 are added to the ks corresponding to node 2 for both element one and two.

Where $$K_{33}= k^{(1)}_{33}+ k^{(2)}_{11}+k^{(3)}_{11} $$ and

$$K_{34}= k^{(1)}_{34}+ k^{(2)}_{13}+k^{(3)}_{12} $$

HW: fill out the rest
$$ \begin{bmatrix} k^{(1)}_{11} & k^{(1)}_{12} & k^{(1)}_{13} & k^{(1)}_{14} & 0 & 0 & 0 & 0\\ k^{(1)}_{21} & k^{(1)}_{22} & k^{(1)}_{23} & k^{(1)}_{24} & 0 & 0 & 0 & 0\\ k^{(1)}_{31} & k^{(1)}_{32} & k^{(1)}_{33}+k^{(2)}_{11}+k^{(3)}_{11} & k^{(1)}_{34}+k^{(2)}_{12}+k^{(3)}_{12} & k^{(2)}_{13} & k^{(2)}_{14} & k^{(3)}_{13} & k^{(3)}_{14}\\ k^{(1)}_{41} & k^{(1)}_{42} & k^{(1)}_{43}+k^{(2)}_{21}+k^{(3)}_{21} & k^{(1)}_{44}+k^{(2)}_{22}+k^{(3)}_{22} & k^{(2)}_{23} & k^{(2)}_{24} & k^{(3)}_{23} & k^{(3)}_{24}\\ 0 & 0 & k^{(2)}_{31} & k^{(2)}_{32} & k^{(2)}_{33} & k^{(2)}_{34} & 0 & 0\\ 0 & 0 & k^{(2)}_{41} & k^{(2)}_{42} & k^{(2)}_{43} & k^{(2)}_{44} & 0 & 0\\ 0 & 0 & k^{(3)}_{31} & k^{(3)}_{32} & 0 & 0 & k^{(3)}_{33} & k^{(3)}_{34}\\ 0 & 0 & k^{(3)}_{41} & k^{(3)}_{42} & 0 & 0 & k^{(3)}_{43} & k^{(3)}_{44} \end{bmatrix} $$

Contributing Team Members
The following students contributed to this report:

Jessica Dyer Eml4500.f08.gravy.jad 01:45, 8 October 2008 (UTC)

Mikal McFarlane Eml4500.f08.gravy.mmm 15:37, 7 October 2008 (UTC)

Patrick Clayton Eml4500.f08.gravy.psc 22:49, 7 October 2008 (UTC)

Sean Sullivan Eml4500.f08.gravy.sms 01:01, 8 October 2008 (UTC)

Brian Eggeman Eml4500.f08.gravy.bje 01:43, 8 October 2008 (UTC)

Justin Tennant Eml4500.f08.gravy.JMT 01:47, 8 October 2008 (UTC)

title: User:Eml4500.f08.gravy.jad/NotesRedo3


title: File:Drawing1.svg


title: File:Drawing11.svg


title: File:FiveBarTruss1.jpg


title: File:Drawing2.svg


text: File:Spars ns max.jpg

title: File:Spars ns max.jpg


title: File:ThreeBardef1.jpg


text: Template:Inactive user

title: Template:Inactive user


text:
User:Eml4500.f08.team.guan/fixed926

The Transformation Matrix
Let's begin with the general equation transforming the forces in global coordinates into axial forces.

$$\begin{bmatrix}

P^{(e)}_{1}\\ P^{(e)}_{2} \end{bmatrix} = T^{(e)}

\begin{bmatrix}

f^{(e)}_{1}\\ f^{(e)}_{2}\\ f^{(e)}_{3}\\ f^{(e)}_{4} \end{bmatrix}$$ Here P is the axial force. Here, the transformation matrix must be a 2 by 4 matrix in order for the math to be correct. Recall the elemental axial force displacement relationship: $$\hat{k}^{(e)}_{2x2}q^{(e)}_{2x1} = p^{(e)}_{2x1}$$ This translates to $$\hat{k}^{(e)}(T^{(e)}d^{(e)}) = (T^{(e)}f^{(e)})$$ The goal here is get the elemental force displacement relationship, so, in order for that to happen, one must "move" the transformation matrix from the right hand side to the left hand side by multiplying by the inverse of the transformation matrix. However, it is unfortunate that the transformation matrix is a rectangular matrix (2 by 4), so it cannot be inverted. How can we obtain the elemental force displacement matrix then? The answer is the following equation. $$[T^{(e)T}_{4x2}\hat{k}^{(e)}_{2x2}(T^{(e)}_{2x4}]d^{(e)}_{4x1} = f^{(e)}_{4x1}$$

After deriving this, it is necessary to indicate that the principle of virtual work is responsible for the elimination of columns 1, 2, 5, and 6 in the stiffness matrix so that $$\displaystyle k_{6x6}d_{6x1}=F_{6x1}$$ becomes $$\displaystyle k_{2x2}d_{2x1}=F_{2x1}$$

Why not solve the global force displacement relationship by multiplying both sides by the inverse of the stiffness matrix K -1 ? This is because the inverse of the stiffness matrix could not be computed and there is a singularity (i.e. the determinant is zero and thus K is not invertible). Recall, you need to compute 1 over the determinant of the matrix to find the inverse. Why? For an unconstrained system, there are three possible rigid body motions in 2 dimensions (2 translation and 1 rotation). Thus, we calculate the eigenvalues of the stiffness matrix to give us information of the motions of this constrained system.

The Finite Element Method eigenvalues can be compared to a dynamics eigenvalue problem such as $$ Kv = \lambda MV$$ where K is a stiffness matrix, $$\lambda$$ is an eigenvalue, and M is a mass matrix. In both the dynamics approach and the Finite Element approach, a zero eigenvalue represents zero stored energy.

Because eigenvalues are used in a broad range of applications, the terminology associated with them varies across the mathematical disciplines, and at times even within a single discipline. Eigenvalues are often referred to as modes, or mode shapes. These modes refer to rigid body modes and correspond to elastic energy for the problems covered in this course of Finite Element Analysis.

title: User:Eml4500.f08.team.guan/fixed926


title: File:Eigvec2bar.jpg


text: File:Eigveccasea.jpg

title: File:Eigveccasea.jpg


text: File:Gaar curve.jpg

title: File:Gaar curve.jpg


text: File:HW5 a.jpeg

title: File:HW5 a.jpeg


text:
User:Eas4200c.f08.wiki.d/HW4

Idealized Single Cell Airfoil
function IdealizedSingleCell

a = 4; b = 2; r = b/2; t1 = .008; t2 = .01; t3 = .01;

y = linspace(0, -r, 25);

y(26:50) = linspace(-r, 0, 25);

y(51) = a; y(52) = 0;

z = sqrt(r^2 - y(1:25).^2);

z(26:50) = -sqrt(r^2 - y(26:50).^2);

z(51) = -r; z(52) = r;

plot(y,z)

A_bar = pi/2*(b/2)^2 + .5*b*a;

ThetaG = 1/(2*A_bar)*(pi*(b/2)/t1 + a/t2 + sqrt(a^2+b^2)/t3); % *q

J = 2*A_bar/ThetaG

EDU>> IdealizedSingleCell

J =

0.1001

Rectangular 2-Cell Airfoil
function Rectangle2Cell

a = .3; b = .6; c = .4;

t1 = .003; t2 = .005; t12 = .004;

A1_bar = a*c; A2_bar = b*c;

X = 2*a/t1 + c/t1 + c/t12;

Y = 2*b/t2 + c/t2 + c/t12;

Z = -c/t12;

alpha = (A2_bar*X - A1_bar*Z)/(A1_bar*Y - A2_bar*Z);

J = (4*A1_bar^2 + 4*alpha*A1_bar*A2_bar)/(X + Z*alpha)

EDU>> Rectangle2Cell

J =

8.5507e-004

Circle Validation
function Circle(c,P_0)

R = c/2; check = 1; j = 1;

while check == 1 ns = 2*j; y = linspace(0,c,ns); for i = 1:ns z_u(i) = sqrt(R^2-(y(i)-R)^2); z_l(i) = -sqrt(R^2-(y(i)-R)^2); end

plot(y,z_u,y,z_l)

A_bar = 0; y_bar = 0; z_bar = 0; s = 0; y_u = y;   y_l = y;

for i = 1:ns-1 PQ = [0 y(ns-i)-y(ns-i+1) z_u(ns-i)-z_u(ns-i+1)]; r = [0 y(ns-i+1)-P_0(2) z_u(ns-i+1)-P_0(3)]; a = cross(r,PQ)/2; A_bar = A_bar + a;       y_bar = y_bar + (y_u(i)+y_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); z_bar = z_bar + (z_u(i)+z_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end

for i = 1:ns-1 PQ = [0 y(i+1)-y(i) z_l(i+1)-z_l(i)]; r = [0 y(i)-P_0(2) z_l(i)-P_0(3)]; a = cross(r,PQ)/2; A_bar = A_bar + a;       y_bar = y_bar + (y_l(i)+y_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); z_bar = z_bar + (z_l(i)+z_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end NS(j) = ns; A_BAR(j) = A_bar(1); if abs(pi*R^2-A_bar(1))/(pi*R^2)<.01 check = 0; else j = j+1; end end

y_bar = y_bar/s; z_bar = z_bar/s;

A_bar = A_BAR(j)

subplot(2,1,1), plot(y_u,z_u,y_l,z_l,y_bar,z_bar,'x'), xlabel('y'), ylabel('z'), title('NACA 2415 Airfoil')

subplot(2,1,2), plot(NS,A_BAR), xlabel('ns'), ylabel('A_b_a_r'), title('ns vs. A_b_a_r')

EDU>> Circle(1,[0 0 0])

A_bar =

0.7779

NACA 2415 Single Cell Airfoil
function Airfoil(NACA,c,P_0)

m = floor(NACA/1000)/100;

p = (floor(NACA/100)-floor(NACA/1000)*10)/10;

ta = (NACA-floor(NACA/100)*100)/100;

ts = .002;

check = 1; j = 1;

while check ==1 ns = 5*j; y = linspace(0,1,ns); for i = 1:ns if y(i)<=p z_c(i) = (m/p^2)*(2*p*y(i)-y(i)^2); theta(i) = atan((m/p^2)*(2*p-2*y(i))); else z_c(i) = (m/(1-p)^2)*((1-2*p)+2*p*y(i)-y(i)^2); theta(i) = atan((m/(1-p)^2)*(2*p-2*y(i))); end z_t(i) = 5*ta*(.2969*y(i)^.5-.1260*y(i)-.3516*y(i)^2+.2843*y(i)^3-.1015*y(i)^4); y_u(i) = (y(i) - z_t(i)*sin(theta(i)))*c; z_u(i) = (z_c(i) + z_t(i)*cos(theta(i)))*c; y_l(i) = (y(i) + z_t(i)*sin(theta(i)))*c; z_l(i) = (z_c(i) - z_t(i)*cos(theta(i)))*c; end Area = 0; y_bar = 0; z_bar = 0; s = 0; for i = 1:ns-1 PQ = [0 y_u(ns-i)-y_u(ns-i+1) z_u(ns-i)-z_u(ns-i+1)]; r = [0 y_u(ns-i+1)-P_0(2) z_u(ns-i+1)-P_0(3)]; a = cross(r,PQ)/2; Area = Area + a;       y_bar = y_bar + (y_u(i)+y_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); z_bar = z_bar + (z_u(i)+z_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end for i = 1:ns-1 PQ = [0 y_l(i+1)-y_l(i) z_l(i+1)-z_l(i)]; r = [0 y_l(i)-P_0(2) z_l(i)-P_0(3)]; a = cross(r,PQ)/2; Area = Area + a;       y_bar = y_bar + (y_l(i)+y_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); z_bar = z_bar + (z_l(i)+z_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end NS(j) = ns; A_BAR(j) = Area(1); J(j) = 4*Area(1)^2*ts/s; if j ~= 1 if abs((J(j)-J(j-1))/J(j-1))*100<1 check = 0; end end j = j+1; end

y_bar = y_bar/s; z_bar = z_bar/s;

A_bar = A_BAR(j-1)

subplot(2,1,1), plot(y_u,z_u,y_l,z_l,y_bar,z_bar,'x'), xlabel('y'), ylabel('z'), title('NACA 2415 Airfoil')

subplot(2,1,2), plot(NS,J), xlabel('ns'), ylabel('Torsional Constant, J'), title('ns vs. J')

NACA 2415 3-Cell Airfoil
function Airfoil3Cell(NACA,c,P_0)

m = floor(NACA/1000)/100;

p = (floor(NACA/100)-floor(NACA/1000)*10)/10;

ta = (NACA-floor(NACA/100)*100)/100;

ts = .002; tp = .003;

check = 1; j = 1;

while check ==1 ns = 20*j; y = linspace(0,1,ns); for i = 1:ns if y(i)<=p z_c(i) = (m/p^2)*(2*p*y(i)-y(i)^2); theta(i) = atan((m/p^2)*(2*p-2*y(i))); else z_c(i) = (m/(1-p)^2)*((1-2*p)+2*p*y(i)-y(i)^2); theta(i) = atan((m/(1-p)^2)*(2*p-2*y(i))); end z_t(i) = 5*ta*(.2969*y(i)^.5-.1260*y(i)-.3516*y(i)^2+.2843*y(i)^3-.1015*y(i)^4); y_u(i) = (y(i) - z_t(i)*sin(theta(i)))*c; z_u(i) = (z_c(i) + z_t(i)*cos(theta(i)))*c; y_l(i) = (y(i) + z_t(i)*sin(theta(i)))*c; z_l(i) = (z_c(i) - z_t(i)*cos(theta(i)))*c; end for i = 1:ns if y_u(i) > c/4 P1UpperIndex = i;           break end end for i = 1:ns if y_l(i) > c/4 P1LowerIndex = i;           break end end for i = P1UpperIndex:ns if y_u(i) > 3*c/4 P2UpperIndex = i;           break end end for i = P1LowerIndex:ns if y_l(i) > 3*c/4; P2LowerIndex = i;           break end end %Cell 1 C1y = (y_u(P1UpperIndex)+y_l(P1LowerIndex))/2; for i = 1:P1UpperIndex C1y(i+1) = y_u(P1UpperIndex-i+1); end for i = 1:P1LowerIndex C1y = [C1y, y_l(i)]; end C1y = [C1y, (y_u(P1UpperIndex)+y_l(P1LowerIndex))/2]; C1z = (z_u(P1UpperIndex)+z_l(P1LowerIndex))/2; for i = 1:P1UpperIndex C1z(i+1) = z_u(P1UpperIndex-i+1); end for i = 1:P1LowerIndex C1z = [C1z, z_l(i)]; end C1z = [C1z, (z_u(P1UpperIndex)+z_l(P1LowerIndex))/2]; %Cell 2 for i = 1:P2UpperIndex-P1UpperIndex C21y(i) = y_u(P2UpperIndex-i+1); end C21y = [C21y, y_l(P1LowerIndex), y_u(P2UpperIndex); for i = 1:P2UpperIndex-P1UpperIndex C21z(i) = z_u(P2UpperIndex-i+1); end C21z = [C21z, z_l(P1LowerIndex), z_u(P2UpperIndex); for i=1:P2LowerIndex-P1LowerIndex C22y(i) = y_l(P1LowerIndex+i-1); end C22y = [C22y, y_u(P2UpperIndex), y_l(P1LowerIndex)] for i=1:P2LowerIndex-P1LowerIndex C22z(i) = z_l(P1LowerIndex+i-1); end C22z = [C22z, z_u(P2UpperIndex), z_l(P1LowerIndex)] %Cell 3 C3y = (y_u(P2UpperIndex)+y_l(P2LowerIndex))/2; for i = 1:ns-P2LowerIndex C3y(i+1) = y_l(P2LowerIndex+i-1) end for i = 1:ns-P2UpperIndex C3y'(i) = y_u(ns-i+1) end C3y = [C3y, C3y', (y_u(P2UpperIndex)+y_l(P2LowerIndex))/2]; C3z = (z_u(P2UpperIndex)+z_l(P2LowerIndex))/2; for i = 1:ns-P2LowerIndex C3z(i+1) = z_l(P2LowerIndex+i-1) end for i = 1:ns-P2UpperIndex C3z'(i) = z_u(ns-i+1) end C3z = [C3z, C3z', (z_u(P2UpperIndex)+z_l(P2LowerIndex))/2]; %Area & Torsional Constant Area1 = 0; Area21 = 0; Area22 = 0; Area3 = 0; for i = 1:size(C1y) PQ = [0 C1y(ns-i)-C1y(ns-i+1) C1z(ns-i)-C1z(ns-i+1)]; r = [0 C1y(ns-i+1)-(y_u(P1UpperIndex)+y_l(P1LowerIndex))/2 C1z(ns-i+1)-0]; a = cross(r,PQ)/2; Area1 = Area1 + a;   end for i = 1:size(C21y) PQ = [0 C21y(i+1)-C21y(i) C21z(i+1)-C21z(i)]; r = [0 C21y(i)-y_l(P1LowerIndex) C21z(i)-z_l(P1LowerIndex)]; a = cross(r,PQ)/2; Area21 = Area21 + a;   end for i = 1:size(C22y) PQ = [0 C22y(i+1)-C22y(i) C22z(i+1)-C22z(i)]; r = [0 C22y(i)-y_u(P2UpperIndex) C22z(i)-z_u(P2UpperIndex)]; a = cross(r,PQ)/2; Area22 = Area22 + a;   end for i = 1:size(C3y) PQ = [0 C3y(i+1)-C3y(i) C3z(i+1)-C3z(i)]; r = [0 C3y(i)-(y_u(P2UpperIndex)+y_l(P2LowerIndex))/2 C3z(i)-0]; a = cross(r,PQ)/2; Area3 = Area3 + a;   end NS(j) = ns; A_BAR(j) = Area(1); if j ~= 1 if abs((A_BAR(j)-A_BAR(j-1))/A_BAR(j-1))*100<.01 check = 0; end end j = j+1; end

y_bar = y_bar/s; z_bar = z_bar/s;

A_bar = A_BAR(j-1)

P1y = [y_u(P1UpperIndex) y_l(P1LowerIndex)];

P1z = [z_u(P1UpperIndex) z_l(P1LowerIndex)];

P2y = [y_u(P2UpperIndex) y_l(P2LowerIndex)];

P2z = [z_u(P2UpperIndex) z_l(P2LowerIndex)];

subplot(2,1,1), plot(y_u,z_u,y_l,z_l,y_bar,z_bar,'x',P1y,P1z,P2y,P2z), xlabel('y'), ylabel('z'), title('NACA 2415 Airfoil')

subplot(2,1,2), plot(NS,A_BAR), xlabel('ns'), ylabel('A_b_a_r'), title('ns vs. A_b_a_r')

title: User:Eas4200c.f08.wiki.d/HW4


text: File:Single-graph.png

title: File:Single-graph.png


text:
User:Shengtaoliu/ENES-100/project 1/Design

Problem Statement

 * Figure out how many ohm resistor for voltmeter.
 * Get the answer in the class.
 * Test the converted ammeter.

The Design Process

 * Use R=V/I get, get the voltmeter's resistor is 10K and shunt resistor is 0.3 ohm.
 * Shunt resistor should be parallel with voltmeter to reduce the voltmeter's resistor.
 * The resistor in the circuit series with converted ammeter is 10K.
 * Here is the diagram I draw.

The initial design

 * Replaced 1M voltmeter's resistor with 1 ohm resistor.
 * Tested the converted ammeter, the needle flip to another side of voltmeter.
 * Because the voltmeter's resistor is too small, and the voltmeter is in series connection.

Operations

 * solder 1k resistor off.

 * Here is my youtube video test converted ammeter
 * Turn on the 3amp power supply.

Maintainability, reliability, and safety

 * The resistor may get burn up and melt when test a couple times.
 * Keep alert. When smell something burning check around.
 * This resistor get burn and melt.

Next Steps

 * Meter is broken, get a new meter.
 * Test converted ammeter in the circuit.
 * Connect ammeter to the power supply.
 * Do tutorial for this project.


text: File:Idealized Single Cell Airfoil.jpg

title: File:Idealized Single Cell Airfoil.jpg


text:
User:Eml4500.f08.delta 6.krueger/HW

Verify the element stiffness matrix. Start with equation 1: $$k^{(e)}={\tilde{T}^{(e)T}}{\tilde{k}^{(e)}}{\tilde{T}^{(e)}}$$ The substitute in with: equation 2 '''$$R^{(e)T}=\begin{bmatrix} l^{(e)} & -m^{(e)}\\ m^{(e)} & l^{(e)} \end{bmatrix}$$''' and equation 3 '''$$R^{(e)}=\begin{bmatrix} l^{(e)} & m^{(e)}\\ -m^{(e)} & l^{(e)} \end{bmatrix}$$''' and equation 4 '''$${\tilde{T}^{(e)T}}=Diag\begin{bmatrix} R^{(e)T} & R^{(e)T}\end{bmatrix}$$''' and equation 5 '''$${\tilde{T}^{(e)}}=Diag\begin{bmatrix} R^{(e)} & R^{(e)}\end{bmatrix}$$''' also '''$${\tilde{k}^{(e)}}=k^{(e)}\begin{bmatrix} 1 & -1\\ -1 & 1 \end{bmatrix}$$''' From previous work $$l^{(1)}=cos(30)=.866025$$ and $$m^{(1)}=sin(30)=.5$$. so '''$$R^{(e)T}=\begin{bmatrix} .866025 & -.5\\ .5 & .866025 \end{bmatrix}$$''' and '''$$R^{(e)}=\begin{bmatrix} .866025 & .5\\ -.5 & .866025 \end{bmatrix}$$''' and '''$$k^{(e)}=\begin{bmatrix} .5625 & 3248 & -.5625 & -.3248\\ .3248 & .1875 & -.3248 & -.1875\\ -.5625 & -.3248 & .5625 & .3248\\ -.3248 & -.1875 & .3248 &.1875 \end{bmatrix}$$''' Now the known matrices are substituted into equations 2, 3, 4 and 5. Then equations 2, 3, 4 and 5 are substituted into equation 1 to verify that: '''$$k^{(e)}=\begin{bmatrix} .5625 & 3248 & -.5625 & -.3248\\ .3248 & .1875 & -.3248 & -.1875\\ -.5625 & -.3248 & .5625 & .3248\\ -.3248 & -.1875 & .3248 &.1875 \end{bmatrix}$$''' Prove equation 1 $${\tilde{f}^{(e)}}={\tilde{T}^{(e)}}f^{(e)}$$ known formulas: equation 2 $${\tilde{f}^{(e)}}={\tilde{k}^{(e)}}{\tilde{d}^{(e)}}$$ equation 3 $${\tilde{d}^{(e)}}={\tilde{T}^{(e)}}d^{(e)}$$ equation 4, where I is the identity matrix $${\tilde{k}^{(e)}}=k^{(e)}I$$ equation 5 $$f^{(e)}=k^{(e)}d^{(e)}$$ Substituting equation 2 into equation 1: $${\tilde{f}^{(e)}}={\tilde{k}^{(e)}}{\tilde{T}^{(e)}}d^{(e)}$$ The substituting equation 3 and equation 4 into equation 1: $${\tilde{f}^{(e)}}=k^{(e)}I{\tilde{T}^{(e)}}d^{(e)}$$ Finally, equation 5 is substituted into equation 1 and the Identity matrix is dropped. This gives the wanted results. $${\tilde{f}^{(e)}}={\tilde{T}^{(e)}}f^{(e)}$$

title: User:Eml4500.f08.delta 6.krueger/HW


text:
User:Eml4500.f08.wiki1.schaet/hw4

MatLab Homework
Part of this homework was to do the example of the 5-bar Truss found on page 25 of the book in Matlab and then to do the 3-bar example established previously in Matlab. Below is the Matlab code for the 5-bar truss with extensive comments.

 MATLAB Code

The image of the deformed beams versus the undeformed beams are as follows:

The Matlab code for the 3-bar truss example with extensive comments are as follows:  MATLAB Code

The image for this result is as follows:

The following are the three M-files given online:  PlaneTrussElement

 NodalSoln

 PlaneTrussResults

title: User:Eml4500.f08.wiki1.schaet/hw4


title: File:Equilnode2.jpg


title: File:Node21.jpg


text: File:Hexagon.GIF

title: File:Hexagon.GIF


text:
User:Eml4500.f08.team.allen./hw4

Eigenvalues and Eigenvectors
Zero Eigenvalues correspond to rigid body motion and zero stored elastic energy.

$$ \mathbf{k}\mathbf{u}\mathbf{i} = \mathbf{0}\mathbf{u}\mathbf{i} $$

The linear combination of $$ \left\{\mathbf{u}_{i}, i=1,...,4 \right\}$$ The equation which defines $$ \mathbf{w} $$ is $$ \sum_{i=1}^{4}{\alpha i\mathbf{u}_{i}}=:\mathbf{w} $$

Homework
The following Matlab code was used to solve for the eigenvalues of the Global Stiffness Matrix of the two truss problem:

Where the eigenvector is given as $$\lambda = \begin{Bmatrix} 0\\ 0\\ 0\\ 0\\ 1.4706\\ 10.0294 \end{Bmatrix} $$

title: User:Eml4500.f08.team.allen./hw4


title: File:Mtg19ex1.JPG


text:
Talk:Theory of relativity/Special relativity

Exercise
According to the classical theory, the kinetic energy of a body can be written as: $$\frac{1}{2}mv^2$$, where $$m$$ is the mass, and $$v$$ is the velocity of the body. According to the special theory of relativity, the kinetic energy of a body can be written as $$(\gamma-1)mc^2$$, where $$\gamma=\frac{1}{\sqrt{1-\frac{v^2}{c^2}} }$$ is called the Lorentz factor. Prove that the classical theory would be correct, according to the special theory of relativity, if the speed of light would be infinite. In other words, prove that: $$ \lim_{c \to \infty}(\gamma-1)mc^2 = \frac{1}{2}mv^2$$

[Solution] We have to prove that

 * $$\lim_{c \to \infty}(\gamma-1)mc^2 = \frac{1}{2}mv^2$$

By extending we get

 * $$ =\lim_{c \to \infty} \frac{(\gamma+1)(\gamma-1)mc^2} {(\gamma+1)} $$

Siplifying gives

 * $$ =\lim_{c \to \infty} \frac{(\gamma^2-1^2)mc^2} {\gamma+1} $$

And further

 * $$ =\lim_{c \to \infty} \frac{(\gamma)^2mc^2-(1)^2mc^2} {\gamma+1} $$

By substituting $$\gamma$$ we get

 * $$=\lim_{c \to \infty} \frac{\mathbf{\left(\frac{1}{\sqrt{1-\frac{v^2}{c^2}}}\right)}^2mc^2-mc^2} {\gamma+1} $$

Simplification gives

 * $$=\lim_{c \to \infty}\frac{\frac{mc^2} {\mathbf{\left (\sqrt{1-\frac{v^2}{c^2}} \right)}^2}-mc^2}{\gamma+1} $$

And so we get

 * $$=\lim_{c \to \infty}\frac{\frac{mc^2} {\mathbf{\left (1-\frac{v^2}{c^2} \right)}}-mc^2}{\gamma+1} $$

Which gives

 * $$=\lim_{c \to \infty}\frac{\frac{mc^2c^2}{c^2-v^2}-mc^2}{\gamma+1}$$

Further

 * $$=\lim_{c \to \infty}\frac{\frac{mc^4}{c^2-v^2}\frac{-(mc^2)(c^2-v^2)}{c^2-v^2}}{\gamma+1}$$

And naturally

 * $$=\lim_{c \to \infty}\frac{\frac{mc^4-mc^4+mv^2c^2}{c^2-v^2}}{\gamma+1}$$

This gives

 * $$=\lim_{c \to \infty}\frac{\frac{mv^2c^2}{c^2(1-\frac{v^2}{c^2})}}{\gamma+1}$$

And so

 * $$=\lim_{c \to \infty}\frac{mv^2}{(1-\frac{v^2}{c^2})(\gamma+1)}$$

By substituting $$\gamma$$ again we get

 * $$=\lim_{c \to \infty}\frac{mv^2}{(1-\frac{v^2}{c^2})(\sqrt{1-\frac{v^2}{c^2}}+1)}$$

By inserting $$\infty$$ we get

 * $$=\frac{mv^2}{(1-\frac{v^2}{\infty^2})(\sqrt{1-\frac{v^2}{\infty^2}}+1)}$$

The definition of infinity gives

 * $$=\frac{mv^2}{(1-0)(\sqrt{1-0}+1)}$$

Simplification gives

 * $$=\frac{1}{2}mv^2$$

$$Q.E.D.$$

--Ofey 18:58, 24 October 2008 (UTC)]

Put all relativity theory resources in a single resource
I think it would be better, if the user comes to a resource called "Theory of Relativity", where there are sections about "General relativity", "Special relativity" and other subsections. Currently, there is no main page that links to all relativity theory resources. But I have not enough experience here to do it.--Sae1962 (discuss • contribs) 10:03, 30 June 2015 (UTC)


title: File:Oct15FEABBQ.jpg


text: File:NACA 2415 Single Cell Airfoil.jpg

title: File:NACA 2415 Single Cell Airfoil.jpg


text:
User:Eml4500.f08.team.foskey.ckf/hw4b

 HW Redo: A better figure for Node 2 was added. Further explanation of the math and equations were provided.

For a comparison between this page, and the original, click here

Eml4500.f08.team.foskey.ckf 02:01, 3 November 2008 (UTC)

Justifying the Global Stiffness Matrix
The assembly of the elemental stiffness matrices, k(e) into a global stiffness matrix, K, can be justified by reviewing the two bar truss problem.

If we look back at the Two bar truss problem, and recall the element force-displacement relation $$ k^{(e)}d^{(e)} = f^{(e)}$$, we can use the Euler Cut Principle to give the equilibrium of the node relationship.

Then, if the free body diagrams of elements 1 and 2 are considered, and the global degrees of freedom for global node 2 are compared to the degrees of freedom for elements 1 and 2, the equilibrium of the node can be confirmed.

Looking at the equilibrium of node 2 for the two bar truss pictured above, the addition of k values is because the forces in the x and y direction sum up to equal zero. Using a statics approach and calculating the forces in each direction gives:

$$ \Sigma F_x = 0= -f^{(1)}_3 - f^{(2)}_1=0$$     (1) $$ \Sigma F_y = 0= P-f^{(1)}_4 - f^{(2)}_2=0$$     (2)

When these two equations are rewritten they become:

$$f^{(1)}_{3} + f^{(2)}_{1}= 0$$

$$f^{(1)}_{4} + f^{(2)}_{2}= 0$$

The next step is to involve the force displacement relationship ( $${k^{(e)}} {d^{(e)}} = {f^{(e)}} $$ ). Replacing f with k*d will get the middle of the 6 by 6 global stiffness matrix.

The element forces can then be written in terms of k and d.

$$f_{3}^{(1)} = \begin{bmatrix} k_{31}^{(1)}d_{1}^{(1)} + k_{32}^{(1)}d_{2}^{(1)} + k_{33}^{(1)}d_{3}^{(1)} + k_{34}^{(1)}d_{4}^{(1)} \end{bmatrix}$$

$$f_{1}^{(2)} = \begin{bmatrix} k_{11}^{(2)}d_{1}^{(2)} + k_{12}^{(2)}d_{2}^{(2)} + k_{13}^{(2)}d_{3}^{(2)} + k_{14}^{(2)}d_{4}^{(2)} \end{bmatrix}$$

Then, the equation for the sum of the forces in the x-direction becomes:

$$ \begin{bmatrix} k_{31}^{(1)}d_{1}^{(1)} + k_{32}^{(1)}d_{2}^{(1)} + k_{33}^{(1)}d_{3}^{(1)} + k_{34}^{(1)}d_{4}^{(1)} \end{bmatrix} + \begin{bmatrix} k_{11}^{(2)}d_{1}^{(2)} + k_{12}^{(2)}d_{2}^{(2)} + k_{13}^{(2)}d_{3}^{(2)} + k_{14}^{(2)}d_{4}^{(2)} \end{bmatrix}=0$$

The sum of the forces in the y-direction can be found in a similar manner.

Thus, statics and equilibrium evaluation of global node 2 has verified the middle of the global stiffness matrix.

title: User:Eml4500.f08.team.foskey.ckf/hw4b


text: File:Hw5 airfoil1.png

title: File:Hw5 airfoil1.png


text: File:Hw5 airfoil2.png

title: File:Hw5 airfoil2.png


title: File:3cellairfoil.png


title: File:Eqnode2.jpg


text: File:CircleValidation.jpg

title: File:CircleValidation.jpg


text:
User:Eml4500.f08.team.Bengtson/teamHw3

Lecture 14
Goal: Transform from a 4x2 matrix to a 2x2. The transformation relationship

\begin{bmatrix} \boldsymbol{P}_1^{(e)}\\ \boldsymbol{P}_2^{(e)}\\ \end{bmatrix} =\boldsymbol{T}^{(e)} \begin{bmatrix} \boldsymbol{f}_1^{(e)}\\ \boldsymbol{f}_2^{(e)}\\ \boldsymbol{f}_3^{(e)}\\ \boldsymbol{f}_4^{(e)}\\ \end{bmatrix} $$

There is a more general form

$$\mathbf{P}^{(e)}_i=\mathbf{T}^{(e)}\mathbf{f}^{(e)}_i$$, where

$$\mathbf{P}^{(e)}_{i}$$ of element (e) axial force at local node $$i$$

$$\mathbf{T}^{(e)}$$ transforms force components to an axial force

$$\mathbf{f}^{(e)}_{i}$$ is the $$i$$-th force component of element (e)

Goal: to find the relationship between

$$\mathbf{q}^{(e)}_{(2x2)}$$ and $$\mathbf{d}^{(e)}_{(4x2)}$$, (axial displacements of each element node and their respective displacement component)

$$\mathbf{P}^{(e)}_{(2x2)}$$ and $$\mathbf{f}^{(e)}_{(4x2)}$$, (axial forces at each element node and their respective force components)

Let us consider the force transformation relationship: $$\mathbf{P}^{(e)}_{(2x2)}=\mathbf{T}^{(e)}_{(2x4)}\mathbf{f}^{(e)}_{(4x2)}$$

FD relationship is shown by:

$$\mathbf{\hat{k}}^{(e)}_{(2x2)}\mathbf{q}^{(e)}_{(2x2)}=\mathbf{P}_{(2x2)}^{(e)}$$

By introducing a transformation matrix into the FD relationship: $$\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}\mathbf{d}^{(e)}=\mathbf{T}^{(e)}\mathbf{f}^{(e)}$$

The ultimate goal is to express the above equation in the form $$\mathbf{k}^{(e)}\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$$. If $$\mathbf{T}^{(e)}$$ were an nxn matrix, then one could simply "move" $$\mathbf{T}^{(e)}$$ from the right hand side to the left hand side by multiplying the entire equation by the inverse $$\mathbf{T^{(e)}}^{-1}$$. Unfortunately, $$\mathbf{T}^{(e)}$$ is a rectangular 2x4 matrix, which means it is non-invertible.

For a more thorough understanding of the Finite Element Method, it is wise to derive the element force displacement with respect to the global coordinate system.

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

The element force displacement matrix can be written as

$$\mathbf{k^{(e)}}\begin{bmatrix} 1 & -1 \\ -1 & 1\\ \end{bmatrix}\begin{pmatrix} q_{1}^{(e)}\\ q_{2}^{(e)}\\ \end{pmatrix}=\begin{pmatrix} P_{1}^{(e)}\\ P_{2}^{(e)}\\ \end{pmatrix} $$

$$q^{(e)}_{i}$$=axial displacementof element e at local node $$i$$

$$P^{(e)}_{i}$$=axial force of element e at local node $$i$$

We want to find equation (1) from equation (2).

We want to find the relationship between:

$$q^{(e)}_{2x1}$$ and $$d^{(e)}_{4x1}$$

$$P^{(e)}_{2x1}$$ and $$f^{(e)}_{4x1}$$

The relationships can be expressed as:

$$q^{(e)}_{2x1}=T^{(e)}_{2x4}d^{(e)}_{4x1}$$

Consider the displacement vector of local node i, denoted by $$d^{(e)}_{i}$$:

$$\vec{d^{(e)}_{[i]}}=d^{(e)}_{1}\vec{i}+d^{(e)}_{2}\vec{j}$$

$$q^{(e)}_{i}$$=axial displacement of node [1] is the orthogonal projection of the displacement vector $$\vec{d_{[1]}^{(e)}}$$onto the $$\tilde{x}$$ axis of element e

Therefore, for node [1], $$q_1^{(e)}$$ can be derived using the following steps.

$$q_{1}^{(e)}=d_{[1]}^{(e)}*\vec{\tilde{i}}$$ $$q_{1}^{(e)}=(d_{1}^{(e)}\vec{i}+d_{2}^{(e)}\vec{j})\vec{\tilde{i}}$$ $$q_{1}^{(e)}=d_{1}^{(e)}(\vec{i}*\vec{\tilde{i}})+d_{2}^{(e)}(\vec{j}*\vec{\tilde{i}})$$ $$q_{1}^{(e)}=l^{(e)}d_{1}^{(e)}+m^{(e)}d_{2}^{(e)}$$ $$q_{1}^{(e)}=\begin{bmatrix} l^{(e)} & m^{(e)}\\ \end{bmatrix}\begin{pmatrix} d_{1}^{(e)}\\ d_{2}^{(e)}\\ \end{pmatrix}$$

The values for $$q_2^{(e)}$$ and $$q_2^{(e)}$$ can now be substituted into the following matrix:

$$\begin{pmatrix} q_{1}^{(e)}\\ q_{2}^{(e)}\\ \end{pmatrix}=\begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)}\\ \end{bmatrix}\begin{pmatrix} d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{4}^{(e)}\\ \end{pmatrix}$$ We cannot solve for the inverse of ($$\mathbf{T}^{(e)}$$ becayse it is not an nxn matrix, we must instead make use of the transpose:

$$\mathbf{T}^{(e)T}_{(4x2)}\hat{k}^{(e)}_{(2x2)}\mathbf{T}^{(e)}_{(2x4)}\mathbf{d}^{(e)}_{(4x1)}=\mathbf{f}^{(e)}_{(4x1)}$$

Note that $$\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}_{(4x2)}\mathbf{\hat{k}}^{(e)}_{(2x2)}\mathbf{T}^{(e)}_{(2x4)}$$

Therefore, $$\mathbf{k}^{(e)}\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$$

Homework: Verify that

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

title: User:Eml4500.f08.team.Bengtson/teamHw3


text:
User:Eml4500.f08.a-team.melvin/Lecture 22

Justification of assembly of $$ \underbrace{ element \ stiffness \ matrix }_{\underline{k}^{(e)} } $$ into  $$ \underbrace{ global \ stiffness \ matrix }_{\underline{K}} $$

e = 1,...,n (# of elements)

Consider the example of a two-bar truss.

Recall element FD relation: k(e) d(e) = f(e)

p.11-3: Euler cut principle, method 2 (equilibrium of global node 2)

p.4-2, 4-3: FBDs of element 1 and element 2, element dofs. d(e)

p.8-3: For node 2, identify global dofs to element dofs for both element 1 and 2.

$$ \sum Fx = 0 $$ = -f3(1) - f1(2) = 0

$$ \sum Fy = 0 $$ = P - f4(1) - f2(2) = 0

Next, use element FD relation k(e) d(e) = f(e)


text: File:3-Bar Truss1024.jpg

title: File:3-Bar Truss1024.jpg


text: File:5 bar truss1024.jpg

title: File:5 bar truss1024.jpg


text:
User:Eml4500.f08.delta 6.guzman/MATLAB Three Bar Truss Example with Plot

Three Bar Truss Example (MATLAB)
It is important to understand how to solve a bar truss system employing a finite element MATLAB code in order to compare the resulting matrices with those obtained from hand calculations.

Take this truss example with the following:

Data

MATLAB Code (Analysis of the system)

MATLAB Functions

In order to run this MATLAB code the following functions are needed:

Note: this MATLAB code was obtained from the MATLAB.zip Chapter 9 folder given by the Student Companion Site Bhatti: Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations website http://bcs.wiley.com/he-bcs/Books?action=resource&bcsId=3105&itemId=0471648078&resourceId=7629

Note: this MATLAB code was obtained from the MATLAB.zip Common folder given by the Student Companion Site Bhatti: Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations website http://bcs.wiley.com/he-bcs/Books?action=resource&bcsId=3105&itemId=0471648078&resourceId=7629

Note: this MATLAB code was obtained from the MATLAB.zip Chapter 9 folder given by the Student Companion Site Bhatti: Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations website http://bcs.wiley.com/he-bcs/Books?action=resource&bcsId=3105&itemId=0471648078&resourceId=7629

Results MATLAB Code (Plot of the system)


title: File:Lalala1.jpg


text:
User:Eml4500.f08.bike.mcdonald/homework report 5

 See my comments below. Please don't remove these comments in case you amend your work; just add a comment in these comment boxes to say what you did. Eml4500.f08 12:56, 10 November 2008 (UTC)

Principal of Virtual Work
To justify the elimination of rows 1, 2, 5, and 6 in order to obtain K as a 2x2 matrix in a 2-bar truss, the FD relation

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

Can be rearranged to be:

$$ \mathbf{Kd} - \mathbf{F} = \mathbf{0}$$ (1)

Which is designated to be equation 1.

Utilizing the principle of virtual work, equation 1 becomes:

$$ \mathbf{W}\cdot (\mathbf{Kd} - \mathbf{F}) = \mathbf{0}$$ (2)

Where W is a 6x1 matrix and Kd - F is a 6x1, and 0 is a scalar number. This equation is designated as equation 2. This equation is true for all values of W and W is called the weighting matrix.

Considering this, Equation 2 being true implies Equation one is also true, and vice versa.

The proof going from Equation 1 to 2 is trivial. However, going from two to one is not trivial and to do this we must choose different values for W

1st Choice:

Choose W so that

$$W_{1} = 1, W_{2} = W_{3} = W_{4} = W_{5} = W_{6} = 0$$

Therefore

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

and

$$ \mathbf{W} \cdot ( \mathbf{Kd-F} ) = 1[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]  $$

Thus, the result becomes

$$\sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j= \mathbf{F}_1 $$

Which is designated to be equation (a).

2nd Choice:

Choose W such that W2=1, W1=W3=W4=W5=W6=0 so that:

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

Therefore:

$$ \mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 1[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]  $$

Thus, the result becomes

$$\sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j= \mathbf{F}_2 $$

Which is designated to be equation (b)

3rd Choice: C hoose W such that W3=1, W1=W2=W4=W5=W6=0 so that

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

$$ \mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 1 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]  $$

Thus, the result becomes  $$\sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j= \mathbf{F}_3 $$

Which is designated to be equation (c)

Choice 4:

Choose W such that W4=1, W1=W2=W3=W5=W6=0 so that

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

$$ \mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 1 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]  $$

Thus, the result becomes  $$\sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j= \mathbf{F}_4 $$

This is designated to be equation (d)

Choice 5:

Choose W such that W5=1, W1=W2=W3=W4=W6=0 so that

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

$$ \mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 1 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 0 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]  $$

Thus, the result becomes  $$\sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j= \mathbf{F}_5 $$

This is designated to be equation (e)

Choice 6:

Choose W such that W6=1, W1=W2=W3=W4=W5=0 so that

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

$$ \mathbf{W} \cdot ( \mathbf{Kd-F} ) = 0[ \sum^6_{j=1} \mathbf{K}_{1j} \mathbf{d}_j- \mathbf{F}_1] + 0[ \sum^6_{j=1} \mathbf{K}_{2j} \mathbf{d}_j- \mathbf{F}_2] + 0 [ \sum^6_{j=1} \mathbf{K}_{3j} \mathbf{d}_j- \mathbf{F}_3] + 0 [ \sum^6_{j=1} \mathbf{K}_{4j} \mathbf{d}_j- \mathbf{F}_4] + 0 [ \sum^6_{j=1} \mathbf{K}_{5j} \mathbf{d}_j- \mathbf{F}_5] + 1 [ \sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j- \mathbf{F}_6]  $$

so that the result becomes  $$\sum^6_{j=1} \mathbf{K}_{6j} \mathbf{d}_j= \mathbf{F}_6 $$

This equation is designated to be equation (f)

Thus, when the equations (a) through (f) are combined, they become $$\mathbf{Kd}= \mathbf{F} $$ which is equation (1) from before.

The weighting coefficients must be kinetically admissible which means they cannot violate the boundary conditions so W1 = W2 =  W5=  W6=0  where the weighting coefficients represent the virtual displacement.

Now,

$$\mathbf{W} \cdot (\mathbf{Kd}-\mathbf{F})$$

Where Kd is a 6x2 matrix minus a 2x1 vector.

Thus

$$= \begin{Bmatrix} W_3 \\ W_4 \end{Bmatrix} (\mathbf{\overline{K} \overline{d}- \overline{F} })$$ (3)

Which is designated as Equation 3.

and is true for all values of $$\begin{Bmatrix} W_3 \\ W_4 \end{Bmatrix} $$

and the K matrix is reduced to $$\mathbf{\overline{K}}=\begin{bmatrix} K_{33} && K_{34} \\ K_{43} && K_{44}\end{bmatrix} $$

as well as:

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

$$\overline{F } = \begin{Bmatrix}F_{3} \\ F_{4} \end{Bmatrix} $$

being the reduced F and d matrices.

Analysis of The Axial Force Displacement Relation Using The Principle of Virtual Work.
Recall: $$ \mathbf{q}^{(i)}=\mathbf{Td}^{(i)} $$

An example of the use of a weighting factor is shown in the computation of our course grade. This computation is shown below.

Course Grade = $$\alpha _{o} \cdot (HW Grade)+\sum_{i=1}^{3}{\alpha _{i}}\cdot (Exam_{i})$$

The weighting factors used in the above computation are $$\alpha _{o}$$ and $$\alpha _{i}$$

Deriving the following relation, $$ \mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)} $$

and also recalling the following Force Displacement relation with axial degrees of freedom,

$$\mathbf{\hat{k}}^{e}_{2\times 2}\mathbf{q}^{(e)}_{2\times 1}= \mathbf{p}^{(e)}_{2\times 1} $$

By subtracting the force matrix $$\mathbf{p}$$, the relation is now written as follows.

$$\mathbf{\hat{k}}^{(e)}\mathbf{q}^{(e)}- \mathbf{p}^{(e)}=0 \qquad $$ Equation (1)

Now, the Principle of Virtual Work can be applied, and it is shown below.

$$ \mathbf{W}\cdot (\mathbf{\hat{k}}^{e}\mathbf{q}^{(e)}-\mathbf{p}^{(e)})=0 \qquad $$   Equation (2)

This equation holds true for all $$\mathbf{W}$$.

Principle of Virtual Work Conversion to Force Displacement Equation
$$\hat{\textbf{w}}_{2x1}$$ = virtual axial displacement corresponding to $$\textbf{q}^{(e)}$$. $$\textbf{w}_{4x1}$$ = virtual displacement in global coordinate system, corresponding to $$\textbf{d}^{(e)}$$.

Substituting Equations (3) and (4) into Equation (2) yields the following Equation: $$\left(\textbf{T}^{(e)}\textbf{W} \right)\cdot \left[\mathbf{\hat{k}}^{(e)}\left(\mathbf{T}^{(e)}\mathbf{\bar{d}}^{(e)}-\mathbf{p}^{(e)} \right) \right]=0$$ for all $$\textbf{w}_{4x1} \; \; \; \; \; \; \;(5)$$

In order to solve Equation (5), a few matrix operations must be recalled:

$$\left(\textbf{AB} \right)^T=\textbf{B}^T\textbf{A}^T \; \; \; \; \; \; \; (6)$$

$$\textbf{a}_{nx1}\cdot\mathbf{b}_{nx1} = \textbf{a}^t_{1xn}\mathbf{b}_{nx1} \; \; \; \; \; \; \; (6)$$

To solve Equation (5), Equation (6) and (7) must be applied:

The dot product must be replaced with the transpose as shown below from Equation (7). $$\left(\textbf{T}^{(e)}\textbf{W} \right)^T\left[\mathbf{\hat{k}}^{(e)}\left(\mathbf{T}^{(e)}\mathbf{\bar{d}}^{(e)}-\mathbf{p}^{(e)} \right) \right]=0$$ for all $$\textbf{w}_{4x1}$$

The matrix transpose is replaced with the product of individual transposes from Equation (6). $$\textbf{W}^{T}\textbf{T}^{(e)T}\left[\mathbf{\hat{k}}^{(e)}\left(\mathbf{T}^{(e)}\mathbf{\bar{d}}^{(e)}-\mathbf{p}^{(e)} \right) \right]=0$$ for all $$\textbf{w}_{4x1}$$

After distributing the transpose matrix T the following equation is the result. $$\textbf{W}^{T}\cdot\left[\left(\textbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}\right)\mathbf{\bar{d}}^{(e)}-\left(\textbf{T}^{(e)T}\mathbf{p}^{(e)} \right) \right]=0$$ for all $$\textbf{w}_{4x1}$$

$$\mathbf{k}^{(e)}$$ replaced $$\textbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$$, and $$\mathbf{f}^{(e)}$$ replaced $$\textbf{T}^{(e)T}\mathbf{p}^{(e)}$$.

The resulting equation is: $$\mathbf{w}\cdot\left[\mathbf{k}^{(e)}\mathbf{d}^{(e)}-\mathbf{f}^{(e)}\right]=0$$ for all $$\mathbf{w}$$

After dividing both sides by w, and moving f(e) to the other side yielded: $$\mathbf{k}^{(e)}\mathbf{d}^{(e)}=\mathbf{f}^{(e)}$$

The previous methods have all been for discrete, non-continuous cases (matrices), and the new analytic methods will be for continuous cases (PDEs).

Motivational Model Problem: An elastic bar with varying A(x) and E(x) is subject to a varying axial load (distributed), concentrated load, and inertia force (dynamic). Both the axial and concentrated loads are time dependent. The figure below depicts the problem.

Equation of Motion, Constitutive Relation, PDE of Motion, and Boundary Condition Derivation
The free body diagram of section dx can be seen below in Figure 1.

Summing the forces in the x direction yields

$$\sum F_x = 0 = -N(x,t) + N(x+dx,t) + f(x,t)dx - m(x){\ddot{u}}dx $$

$$ = \frac{\partial N(x,t)}{\partial x}dx + h.o.t + f(x,t)dx - m(x){\ddot{u}}dx $$ ⇨ $$ Equation \left(1\right)$$

The dx terms are canceled out of this equation, as well as the high order terms(h.o.t). Let this equation be call equation (1).

Why are the high order terms neglected?

Recall Taylor series expansion: $$ f(x+dx) = f(x) + \frac{df(x)}{dx}dx + \frac{1}{2}\frac{d^{2}f(x)}{dx^{2}}dx^{2} + ... $$ where $$\frac{1}{2}\frac{d^{2}f(x)}{dx^{2}}dx^{2}$$ is considered a h.o.t and is therefore neglected.

From equation (1) the equation of motion(EOM) of the system can be derived. Let the EOM be called equation (2).

$$ \frac{\partial N}{\partial x} + f = m{\ddot{u}} $$ ⇨ $$ Equation \left(2\right)$$

Equation (3), the constitutive relation, can then be determined.

$$ N\left(x,t\right) = A(x)\sigma(x,t) $$ ⇨ $$ Equation \left(3\right)$$ where $$ \sigma\left(x,t\right) = E(x)\epsilon(x,t) $$ and $$ \epsilon(x,t) = \frac{\partial u(x,t)}{\partial x} $$

Substituting equation (3) into equation (2) yields the partial differential equation of motion, equation (4) shown below.

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

There are 2 initial conditions(displacement and velocity), and therefore 2 boundary conditions(2nd order derivative with respect to x).

For Figure 2 the boundary conditions are $$ u\left(0,t\right) = 0 = u(L,t) $$

For Figure 3 the boundary condition at 0 is also $$ u\left(0,t\right) = 0 $$, but the boundary condition at L can be derived as follows

$$ N\left(L,t\right) = F(t) $$ where $$ N = A\left(L\right)\sigma(L,t)$$ and $$\sigma\left(L,t\right) = E(L)\epsilon(L,t)$$ and $$ \epsilon(L,t) = \frac{\partial u(L,t)}{\partial x} $$

The boundary condition at L can therefore be simplified to

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

Axial FD Relation
When an additional displacement DOF is added to the system, the Axial Force-Displacement Relation does not change. This is because a member can only have two-reaction forces, and those forces are constrained to be along the direction of the member. Because deformation in a member is caused by these forces, the only axial displacements in a space truss member are found at the bar ends in the direction of the member.

Thus, the Axial Force-Displacement Relation is as shown below, where P represents the axial forces at each node in the member and q represents the axial displacement of each node.

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

$$ \begin{Bmatrix}P_{1}^{(e)} \\ P_{2}^{(e)}\end{Bmatrix} = k^{(e)}\begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}\begin{Bmatrix}q_{1}^{(e)}\\ q_{2}^{(e)}\end{Bmatrix} $$

Element DOFs and Element Forces in Global xyz Coordinates
The following image shows the proper setup of the element displacement DOFs and element forces for a space truss element.

Since each element in a space truss system has 6 displacement DOFs and forces, the element displacement DOF and force matrices are of size 6 x 1. These matrices are shown below.

$$ \mathbf{d}^{(e)} = \begin{Bmatrix}d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{4}^{(e)}\\ d_{5}^{(e)}\\ d_{6}^{(e)}\end{Bmatrix} \qquad \mathbf{f}^{(e)} = \begin{Bmatrix}f_{1}^{(e)}\\ f_{2}^{(e)}\\ f_{3}^{(e)}\\ f_{4}^{(e)}\\ f_{5}^{(e)}\\ f_{6}^{(e)}\end{Bmatrix} $$

Transformation Matrix
For plane truss members, the transformation matrix is as shown below, where $$l^{(e)}$$ and $$m^{(e)}$$ are the director cosines for a 2D system.

$$ \mathbf{T}^{(e)}_{2 x 4} = \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix} $$

For space truss members, the transformation matrix is slightly different, due to the fact that an additional coordinate, z, is added to the system. The director cosine corresponding to the z-coordinate is represented by $$n^{(e)}$$. The transformation matrix for a space truss system is shown below.

$$ \mathbf{T}^{(e)}_{2 x 6} = \begin{bmatrix} l^{(e)} & m^{(e)} & n^{(e)} & 0 & 0 & 0 \\ 0 & 0 & 0 & l^{(e)} & m^{(e)} & n^{(e)} \end{bmatrix} $$

Relation between Axial DOFs and Element DOFs in Global Coordinates
The relation between the axial displacement DOFs and element DOFs using the transformation matrix is shown below.

$$ \mathbf{q}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{d}_{6\times 1}^{(e)} $$

$$ \begin{Bmatrix} q_{1}^{(e)}\\ q_{2}^{(e)} \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} d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{4}^{(e)}\\ d_{5}^{(e)}\\ d_{6}^{(e)} \end{Bmatrix} $$

Relation between Axial Forces and Element Forces in Global Coordinates
The relation between the axial forces and element forces using the transformation matrix is shown below.

$$ \mathbf{P}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)} $$

$$ \begin{Bmatrix} P_{1}^{(e)}\\ P_{2}^{(e)} \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} f_{1}^{(e)}\\ f_{2}^{(e)}\\ f_{3}^{(e)}\\ f_{4}^{(e)}\\ f_{5}^{(e)}\\ f_{6}^{(e)} \end{Bmatrix} $$

Derivations Using the Principal of Virtual Work
In this section, various relations for space truss elements are derived using the Principal of Virtual Work.

Global FD Relation
The element axial force-displacement relation is shown below.

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

Relations for the element axial displacement DOFs and element axial forces in terms of the transformation matrix and the element displacement DOFs and forces in global coordinates are shown below.

$$ \mathbf{q}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{d}_{6\times 1}^{(e)} $$

$$ \mathbf{P}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)} $$

Substituting these relations into the element axial force-displacement relation, the following relation is developed.

$$ \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)} = \mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6}\mathbf{d}_{6\times 1}^{(e)} $$

Applying the Principal of Virtual Work, the above relation can be written as shown below, greatly simplifying the problem.

$$ \mathbf{f}_{6\times 1}^{(e)} = [\mathbf{T}^{(e)^{T}}_{2\times 6}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6}]\mathbf{d}_{6\times 1}^{(e)} $$

 Error: No, you cannot "solve" the equation further above to obtain the above equation; you need to use the PVW. Eml4500.f08 11:06, 10 November 2008 (UTC)  The above derivation has been modified to explain how the PVW is used. Eml4500.f08.bike.mcdonald 21:55, 11 November 2008 (UTC)

Multiplying out the term inside the brackets, it can be shown that it is equal to the element stiffness matrix, as shown below.

$$ \mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} = \begin{bmatrix} l^{(e)}&0 \\ m^{(e)}&0 \\ n^{(e)}&0 \\ 0&l^{(e)} \\ 0&m^{(e)} \\ 0&n^{(e)} \end{bmatrix}k^{(e)}\begin{bmatrix} 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} $$

$$ \mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} = k^{(e)}\begin{bmatrix} l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} & -l^{(e)^{2}} & -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} \\ m^{(e)}l^{(e)}& m^{(e)^{2}} & m^{(e)}n^{(e)} & -m^{(e)}l^{(e)} & -m^{(e)^{2}} & -m^{(e)}n^{(e)} \\ n^{(e)}l^{(e)}& m^{(e)}n^{(e)} & n^{(e)^{2}} & -n^{(e)}l^{(e)} & -m^{(e)}n^{(e)} & -n^{(e)^{2}} \\ -l^{(e)^{2}}& -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} & l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} \\ -m^{(e)}l^{(e)}& -m^{(e)^{2}} & -m^{(e)}n^{(e)} & m^{(e)}l^{(e)} & m^{(e)^{2}} & m^{(e)}n^{(e)} \\ -n^{(e)}l^{(e)}& -m^{(e)}n^{(e)} & -n^{(e)^{2}} & n^{(e)}l^{(e)} & m^{(e)}n^{(e)} & n^{(e)^{2}} \end{bmatrix} $$

$$ \mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} = \mathbf{k}^{(e)}_{6\times 6} $$

Substituting this relation into the above relation for the element force matrix in global coordinates gives the desired result, the element force-displacement relation, proved using the Principle of Virtual Work.

$$ \mathbf{f}_{6\times 1}^{(e)} = \mathbf{k}^{(e)}_{6\times 6}\mathbf{d}_{6\times 1}^{(e)} $$

Element Stiffness Matrix in Global Coordinates in Terms of Transformation Matrix and Element Stiffness Matrix from Axial FD Relation
The following relation shows the element stiffness matrix of a space truss element in terms of the transformation matrix, $$T^{(e)}$$, and the element stiffness matrix from the axial force-displacement relation, $$\hat{k}^{(e)}$$. This relation was derived in the above derivation of the global force-displacement relation.

$$ \mathbf{k}^{(e)}_{6\times 6} = \mathbf{T}^{(e)^{T}}_{6\times 2}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}_{2\times 6} $$

Element Forces in Global Coordinates in Terms of Transformation Matrix and Element Axial Forces
The following relation was derived and shown in the above derivation of the space truss element force-displacement relation.

$$ \mathbf{P}_{2\times 1}^{(e)} = \mathbf{T}^{(e)}_{2\times 6}\mathbf{f}_{6\times 1}^{(e)} $$

Multiplying both sides of this relation by $$\mathbf{T}^{(e)^{T}}_{6\times 2}$$ gives the relation shown below for the element forces in global coordinates in terms of the transformation matrix and the element axial forces.

$$ \mathbf{f}_{6\times 1}^{(e)} = \mathbf{P}_{2\times 1}^{(e)}\mathbf{T}^{(e)^{T}}_{6\times 2} $$

Element Stiffness Matrix in Global Coordinates in Terms of Director Cosines
The following relation shows the element stiffness matrix of a space truss element in terms of the director cosines, $$l^{(e)}$$, $$m^{(e)}$$, and $$n^{(e)}$$. This relation was derived in the above derivation of the global force-displacement relation.

$$ \mathbf{k}_{6\times 6}^{(e)} = k^{(e)}\begin{bmatrix} l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} & -l^{(e)^{2}} & -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} \\ m^{(e)}l^{(e)}& m^{(e)^{2}} & m^{(e)}n^{(e)} & -m^{(e)}l^{(e)} & -m^{(e)^{2}} & -m^{(e)}n^{(e)} \\ n^{(e)}l^{(e)}& m^{(e)}n^{(e)} & n^{(e)^{2}} & -n^{(e)}l^{(e)} & -m^{(e)}n^{(e)} & -n^{(e)^{2}} \\ -l^{(e)^{2}}& -m^{(e)}l^{(e)} & -n^{(e)}l^{(e)} & l^{(e)^{2}} & m^{(e)}l^{(e)} & n^{(e)}l^{(e)} \\ -m^{(e)}l^{(e)}& -m^{(e)^{2}} & -m^{(e)}n^{(e)} & m^{(e)}l^{(e)} & m^{(e)^{2}} & m^{(e)}n^{(e)} \\ -n^{(e)}l^{(e)}& -m^{(e)}n^{(e)} & -n^{(e)^{2}} & n^{(e)}l^{(e)} & m^{(e)}n^{(e)} & n^{(e)^{2}} \end{bmatrix} $$

Rerunning the MATLAB Code for 2-Bar Truss System and Interpreting Array
The setup of the 2-bar truss system to be analyzed is shown below. This system was previously analyzed and verified using statics in Homework Report 2, however there was a bug in the MATLAB code that needed to be corrected. The steps taken to correct the bug will be shown in this section, and the new results that come from the modified MATLAB code will be presented.

The properties of the 2-bar truss system to be analyzed is shown in the table below.

Modified MATLAB Code
When the MATLAB code for solving the 2-bar truss system using the Finite Element Method was first analyzed in Homework Report 2, there was a bug that showed up when calculating the  array. This bug was due to the fact that the MATLAB code passed the entire  and   arrays into the   function, instead of passing only the proper array index in.

The section of code with the error is shown below prior to fixing the bug.

In order to fix the bug, the variables  and   were replaced with   and , which allowed for the values at each index in these arrays to be properly used when calculating the results. This is because the function  is located inside a   loop from the value   to the variable , which in this case has a value of 2.

The modified section of code is shown below.

The overall modified MATLAB code used to solve the 2-bar truss system is shown below.

Interpreting Array
The MATLAB command window output after running the MATLAB code shown in the above section is shown below.

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

The last value shown in the MATLAB command window output is the  array. The first row of the array is associated with Element 1, and the second row is associated with Element 2. The first column gives the value of axial strain in the corresponding elements, the second column gives the value of axial stress in the corresponding elements, and the third column gives the value of axial force in the corresponding elements.

A table showing the results of the Finite Element Analysis on the 2-bar truss system is given below.

In order to verify that the results obtained from the MATLAB code shown above, the calculation of axial force, axial stress, and axial strain using statics is shown below for both element 1 and 2 of the 2-bar truss system.

The axial forces of both elements were obtained by summing forces in both the x and y directions, as shown below.

$$ \sum{F_{x}}= F_{1}\cos{30^\circ} - F_{2}\cos{45^\circ} $$

$$ \sum{F_{y}}= F_{1}\sin{30^\circ} + F_{2}\sin{45^\circ} $$

$$ F_{1} = 5.124 $$

$$ F_{2} = 6.277 $$

Element 1
Knowing the axial force in element 1, the axial stress and axial strain can be calculated as shown below.

$$ F_{1} = 5.124 $$

$$ \sigma_{1} = \frac{F_{1}}{A_{1}} = \frac{5.124}{1} = 5.124 $$

$$ \epsilon_{1} = \frac{\sigma_{1}}{E_{1}} = \frac{5.124}{3} = 1.708 $$

Comparing these values with the values obtained from the Finite Element Method solution above, it can be seen that both methods are in complete agreement.

Element 2
Knowing the axial force in element 2, the axial stress and axial strain can be calculated as shown below.

$$ F_{2} = 6.277 $$

$$ \sigma_{2} = \frac{F_{2}}{A_{2}} = \frac{6.277}{2} = 3.139 $$

$$ \epsilon_{2} = \frac{\sigma_{2}}{E_{2}} = \frac{3.139}{5} = 0.628 $$

Comparing these values with the values obtained from the Finite Element Method solution above, it can be seen that both methods are in complete agreement.

Finite Element Analysis of Various 6-Bar Truss Systems
The setup of the 6-bar truss system as given on page 226 of FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti is shown in the figure below.

6-Bar Truss System With Same Young's Modulus and Cross-Sectional Area Values
The properties of the 6-bar truss system with the same values for Young's Modulus and cross-sectional area are shown in the table below.

MATLAB Code
The following MATLAB code was provided as an online resource from FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti.

In order for this code to run properly, it needs three additional MATLAB functions. These functions are listed and shown in detail below.

Plot of Undeformed and Deformed Structure
The plot of the undeformed and deformed shapes of the 6-bar truss system with the same values for Young's Modulus and cross-sectional area is shown below with a multiplication factor on the deformation of 3000.

The source code added to the original MATLAB code to produce the above plot is given below.

6-Bar Truss System Where Young's Modulus Values Are Not Equal
The properties of the same 6-bar truss system analyzed above but different and varying values for Young's Modulus are shown in the table below.

MATLAB Code
The following MATLAB code was provided as an online resource from FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti.

In order for this code to run properly, it needs three additional MATLAB functions. These functions are listed and shown in detail below.

Plot of Undeformed and Deformed Structure
The plot of the undeformed and deformed shapes of the 6-bar truss system with the same values for Young's Modulus and cross-sectional area is shown below with a multiplication factor on the deformation of 3000.

The source code added to the original MATLAB code to produce the above plot is given below.

Plot Comparing Deformed Structures of Both 6-Bar Truss Systems
The following plot shows a comparison of the deformed structures of the two 6-bar truss systems analyzed previously using the Finite Element Method. The blue dotted line represents the undeformed structure of the 6-bar truss system, the solid red line represents the deformed structure of the 6-bar truss system with the same values of Young's Modulus, and the solid green line represents the deformed structure of the 6-bar truss system with different values of Young's Modulus.

The MATLAB code used to produce the above plot is shown below.

Finite Element Analysis of 3-Bar Space Truss
The setup of the 3-bar space truss system as given on page 230 of FUNDAMENTAL Finite Element Analysis and Applications With Mathematica and MATLAB Computations, written by M. Asghar Bhatti is shown in the figure below.

The element and node properties of the 3-bar space truss system are shown in the following two tables.

In order to properly create MATLAB code to solve the 3-bar space truss system using the Finite Element Method, the  and   arrays must be populated. The proper assembly of the  array is shown below.

$$ \texttt{conn}=\left[ \begin{array}{cc} 1 & 4 \\ 2 & 4 \\ 3 & 4 \end{array} \right] $$

The proper assembly of the  array is shown below.

$$ \texttt{lmm}=\left[ \begin{array}{cccccc} 1 & 2 & 3 & 10 & 11 & 12 \\ 4 & 5 & 6 & 10 & 11 & 12 \\ 7 & 8 & 9 & 10 & 11 & 12 \end{array} \right] $$

MATLAB Code
The MATLAB code, shown below, and used to solve the 3-bar space truss system was provided as an online resource by FUNDAMENTAL Finite Element Analysis and Applications with Mathematica and MATLAB Computations, written by M. Asghar Bhatti.

The following additional functions were called in the above code and were necessary to properly solve the 3-bar truss system using MATLAB.

MATLAB Results
The output in the MATLAB command window after running the above code is shown below.

K =

1.0e+004 *

Columns 1 through 11

0.1460   0.2919   -0.3041         0         0         0         0         0         0   -0.1460   -0.2919    0.2919    0.5839   -0.6082         0         0         0         0         0         0   -0.2919   -0.5839   -0.3041   -0.6082    0.6335         0         0         0         0         0         0    0.3041    0.6082         0         0         0    0.3567   -0.3567    0.4954         0         0         0   -0.3567    0.3567         0         0         0   -0.3567    0.3567   -0.4954         0         0         0    0.3567   -0.3567         0         0         0    0.4954   -0.4954    0.6880         0         0         0   -0.4954    0.4954         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0    6.0000         0         0   -0.1460   -0.2919    0.3041   -0.3567    0.3567   -0.4954         0         0         0    0.5026   -0.0647   -0.2919   -0.5839    0.6082    0.3567   -0.3567    0.4954         0         0         0   -0.0647    0.9405    0.3041    0.6082   -0.6335   -0.4954    0.4954   -0.6880         0         0   -6.0000    0.1913   -1.1036

Column 12

0.3041   0.6082   -0.6335   -0.4954    0.4954   -0.6880         0         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

The computed reactions of the 3-bar space truss system are shown below.

$$ \mathbf{R}= \left[ \begin{array}{c} R_{1} \\ R_{2} \\ R_{3} \\ R_{4} \\ R_{5} \\ R_{6} \\ R_{7} \\ R_{8} \\ R_{9} \\ \end{array} \right] = \left[ \begin{array}{c} 6667 \\ 13333 \\ -13889 \\ -6667 \\ 6667 \\ -9259 \\ 0 \\ 0 \\ 23148 \\ \end{array} \right] $$

Determination of Statically Determinate or Indeterminate System
The 3-bar space truss system analyzed above is a statically determinate system. This is because each of the three members of the system have an axial force directed along the member, which gives the system 3 unknown forces. Because the space truss system has 3 displacement DOFs at each node, the summation of forces can be performed in the x-direction, y-direction, and z-direction. This gives 3 unique equations. Because there are 3 equations and 3 unknowns, this system can be solved using only the methods of statics.

In order to verify that the results obtained from the MATLAB code shown above, the calculation of axial force, axial stress, and axial strain using statics is shown below for the 3-bar space truss system.

Three unique equations were obtained by summing forces in the x, y, and z directions, as shown below.

$$ \sum{F_{x}}= R_{1}l^{(1)} + R_{2}l^{(2)} + R_{3}l^{(3)} = 0 $$

$$ \sum{F_{y}}= R_{1}m^{(1)} + R_{2}m^{(2)} + R_{3}m^{(3)} = -P $$

$$ \sum{F_{z}}= R_{1}n^{(1)} + R_{2}n^{(2)} + R_{3}n^{(3)} = 0 $$

$$ R_{1} = 20375 $$

$$ R_{2} = 13214 $$

$$ R_{3} = -23149 $$

Element 1
Knowing the axial force in element 1, the axial stress and axial strain can be calculated as shown below.

$$ R_{1} = 20375 $$

$$ \sigma_{1} = \frac{R_{1}}{A_{1}} = \frac{20375}{200} = 101.875 $$

$$ \epsilon_{1} = \frac{\sigma_{1}}{E_{1}} = \frac{101.875}{200\times10^{3}} = 5.094\times10^{-4} $$

Element 2
Knowing the axial force in element 2, the axial stress and axial strain can be calculated as shown below.

$$ R_{2} = 13214 $$

$$ \sigma_{2} = \frac{R_{2}}{A_{2}} = \frac{13214}{200} = 66.07 $$

$$ \epsilon_{2} = \frac{\sigma_{2}}{E_{2}} = \frac{66.07}{200\times10^{3}} = 3.3035\times10^{-4} $$

Element 3
Knowing the axial force in element 3, the axial stress and axial strain can be calculated as shown below.

$$ R_{3} = -23149 $$

$$ \sigma_{3} = \frac{R_{3}}{A_{3}} = \frac{-23149}{600} = -38.58 $$

$$ \epsilon_{3} = \frac{\sigma_{3}}{E_{3}} = \frac{-38.58}{200\times10^{3}} = -1.929\times10^{-4} $$

Comparing these values with the values obtained from the Finite Element Method solution above, it can be seen that both methods are in complete agreement.

Plots of Undeformed and Deformed Structure using MATLAB
The following sections show plots of the undeformed and deformed shapes of the 3-bar space truss system described above. These plots were created using the following MATLAB code, which was added to the end of the MATLAB code used to solve the 3-bar space truss system using the Finite Element Method.

Perspective Views
A perspective view of the undeformed and deformed shapes of the 3-bar space truss analyzed using the Finite Element Method is shown below. This perspective view is from the coordinates: x = -2 m, y = -2 m, z = 3 m

In order to obtain a more clear image of the deformation of the 3-bar space truss system, a perspective view from the coordinates of x = 1 m, y = 1 m, and z = 1 m is shown below.

Contributing Team Members
Andrew McDonald - Eml4500.f08.bike.mcdonald 21:25, 2 November 2008 (UTC) Garrett Pataky - Eml4500.f08.bike.pataky 20:25, 3 November 2008 (UTC) Sam Bernal - Eml4500.f08.bike.bernal 03:32, 7 November 2008 (UTC) Bobby Sweeting - Eml4500.f08.bike.sweeting 22:10, 2 November 2008 (UTC) Shawn Gravois - Eml4500.f08.bike.gravois 06:09, 5 November 2008 (UTC) Eric Viale - Eml4500.f08.bike.viale 21:44, 6 November 2008 (UTC)


title: File:Meet14-1ex.JPG


title: File:Area2.jpeg


title: File:Area1.jpeg


text: File:Hw4 poisson.jpg

title: File:Hw4 poisson.jpg


text:
File:Simple airfoil.JPG

Eas4200c.f08.ZYX.O'Mahoney 20:27, 24 October 2008 (UTC)

title: File:Simple airfoil.JPG


text: File:Newtrans.jpg

title: File:Newtrans.jpg


text: File:Three-bar Truss Structure (undeformed and deformed).GIF