text:
User:Eml4500.f08.delta 6.krueger/Matlab sections 1-5

Section 1: Accessing Matlab
Once logged into the system being used, Matlab can be entered and exited using the matlab and quit system commands. When using certain systems, such as Unix systems and MS Windows, it is suggested that the user keep both the Matlab and local editor active.

Section 2: Entering Matrices
Matlab is a matrix manipulator. Matrices can be 1x1 matrices, scalars, matrices with one column or one row, vectors, or any variation of columns and rows. The inputs can be integers, variables or complex numbers.

Inputting matrices into Matlab can be done in many different ways. They can be input with an explicit list of elements, built in statements or functions, diskfile from local editor or from external data files or applications.

The following are examples of matrix inputs.

Complex Numbers can be put into matrices in the following ways.

If an external data file is used. The file can be loaded into Matlab using the load function. For example, if the file is named data.ext. The file can be uploaded into Matlab using the function load data.ext.

There are many built in function in Matlab that can be used to generate matrices. This section will look at three common matrix generators, rand, magic and hilb.

rand(n) generates a matrix with random entries between 0 and 1 into a nxn matrix. rand(m,n) generates a matrix with rando0m entries between 0 and 1 into an mxn matrix Below is an example of the rand function.

magic(n) will create an integral nxn matrix which is a magic square. This means that the rows, columns and diagonals have common sums. The following is an example of the magic(n) function.

hilb(n) will create a nxn hilbert matrix (m and n denote positive integers).

The following is an example of the hilb(n) function.

Matrix and vector entries can be referenced using A(3,2). This pulls the entry from the 3rd row and the 2nd column from matrix A.

The following is an example of how to reference an entry within a matrix.

Section 3: Matrix Operation, Array Operations
Section three references the basic matrix operators. These operators are addition,+, subtraction,-, Multiplication,*, power,^, transpose,', and left and right division,/ and \. These operators should be known by anyone with a programming background.

While using the multiplication, power and division operators, a period can be used that will make the operators operate entry-wise.

The following is an example of entry-wise operations.

Section 4: Statements, Expressions and Variables; Saving a Session
Matlab is an expression language. It evaluates and interprets expressions entered and produces matrices. Expressions are composed of operators, functions and variable names.

A Statement is normally terminated with the return button. Several statemnet can be place on the same line using commas or semicolons following the statements. Semicolons suppress the printing of the statement.

A session can be saved be using the save command and reloaded using the load command. The following are examples of the save and load commands. EDU>> save

Saving to: matlab.mat

EDU>> load

Loading from: matlab.mat

Section 5: Matrix Building Functions
This section references convenient matrix building functions.

The eye function give an identity matrix of a given matrix. The following is an example of the eye functiion. This example shows the identity matrix of matrix A. EDU>> eye(size(A))

ans =

1    0

0    1

The zeros function gives a matrix of zeros. The following is an example of the zero function.

EDU>> zeros(3,3)

ans =

0    0     0     0     0     0     0     0     0

The diag function creates or extracts diagonals. The following is an example of the diag function.

EDU>> diag(5)

ans =

5

Using two diag functions together a diagonal of a matrix can be retrieved. The following is an example of the diag(diag(A)) function that will retrieve the diagonal of the A matrix.

EDU>> diag(diag(A))

ans =

1.0000 + 5.0000i       0 0            4.0000 + 8.0000i

Matrices can also be built from blocks. The following is an example of a matrix built from matrix A, two zeros matices and one eye matrix. This function builds a 5x5 matrix.

EDU>> B=[A, zeros(3,2); zeros(2,3), eye(2)]

B =

1    0     0     0     0     0     1     0     0     0     0     0     1     0     0     0     0     0     1     0     0     0     0     0     1


text:
User:Eas4200c.f08.WIKI.E/hw3

Friday Octopber 8th
K) 1) Multicelled Section

For cells where $$i = 1, ... , n_{cell}$$

$$\tau = 2\sum_{i=1}^{N_c}{q_i \bar{A_i}} $$

$$q_i \;$$ = Shear flow in cell i.

$$\bar{A_i}$$ = "average" area in cell i.

Define: $$T_i = 2q_i\bar{A_i}$$ torque generated by one cell.

$$\tau = 2\sum_{i=1}^{N_c}{q_i \bar{A_i}} $$

$$T_i = \sum_{i=1}^{N_c}{T_i} $$

2) Shape of airfoil is "rigid" in the plane $$(X,Y)$$ however it can warp out of the plane.

$$\theta = \theta _1 = \theta _2= ...=\theta _{N_c}$$

The twist angle equation applied to cell i.

$$\theta_i =\frac{1}{2 G_i \bar{A_i}}\oint_{}^{}{\frac{q_i}{t_i}ds}$$

$$G_i\;$$ = Shear Modulus of cell i

$$t_i(s)\;$$ = thickness of walls on cell i. This is a linear curve coordinate along the cell wall.

''' Pb 1.1: Rectangular single cell section with modification. '''

Now more general with a single cell section:

$$t_i=0.08 \; m$$

$$t_2= t_3=0.01 \; m$$

$$\bar{A}=\frac{1}{2} \pi (\frac{b}{2})^2 +\frac{1}{2}ba= 5.5708$$

Shear flow comes from the torque equation: $$T=2q\bar{A}$$

$$q= \frac{T}{2\bar{A}}$$

Twist angle: $$\theta_i =\frac{1}{2 G_i \bar{A_i}}\sum_{j=1}^{3}{\frac{q_jl_j}{t_j}}$$ $$j=1,2,3$$ index for the segment number.

There is only one cell for $$\theta_i\;$$ so it reduces to $$\theta\;$$ in the equation above. You would normally have two indexes for multicelled section $$(i,j)\;$$.

$$\theta =\frac{1}{2 G \bar{A}}\sum_{j=1}^{3}{\frac{q_jl_j}{t_j}}$$

Note: $$\int$$ is and elongated S, standing for continuous summation. $$\sum{}$$ is a discrete sum.

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

title: User:Eas4200c.f08.WIKI.E/hw3


text:
User:Eml4500.f08.delta 6.krueger/Notes Monday September 15th

Global stiffness matrix K.

 * {| style="text-align:center; border-collapse:collapse;" cellpadding="4"

 * rowspan="6" | K  =
 * rowspan="6" | $$\left.\begin{align} \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \end{align} \right($$
 * title="stiffness matrix from element one" style="border-left: 1px solid black; border-top: 1px solid black;" | $$k_{11}^{(1)}$$
 * title="stiffness matrix from element one" style="border-top: 1px solid black;" | $$k_{12}^{(1)}$$
 * title="stiffness matrix from element one" style="border-top: 1px solid black;" | $$k_{13}^{(1)}$$
 * title="stiffness matrix from element one" style="border-top: 1px solid black; border-right: 1px solid black;" | $$k_{11}^{(1)}$$
 * 0
 * 0
 * rowspan="6" | $$\left)\begin{align} \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \end{align} \right.$$
 * title="stiffness matrix from element one" style="border-left: 1px solid black;" | $$k_{21}^{(1)}$$
 * title="stiffness matrix from element one" | $$k_{22}^{(1)}$$
 * title="stiffness matrix from element one" | $$k_{23}^{(1)}$$
 * title="stiffness matrix from element one" style="border-right: 1px solid black;" | $$k_{24}^{(1)}$$
 * 0
 * 0
 * title="stiffness matrix from element one" style="border-left: 1px solid black;" | $$k_{31}^{(1)}$$
 * title="stiffness matrix from element one" style="border-right: 1px solid black;" | $$k_{32}^{(1)}$$
 * title="stiffness matrix from element one and element two" style="border-top: 1px solid black;" | $$k_{33}^{(1)} + k_{11}^{(2)}$$
 * title="stiffness matrix from element one and element two" style="border-right: 1px solid black; border-top: 1px solid black;" | $$k_{34}^{(1)} + k_{12}^{(2)}$$
 * title="stiffness matrix from element two" style="border-top: 1px solid black;" | $$k_{13}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black; border-top: 1px solid black;" | $$k_{14}^{(2)}$$
 * title="stiffness matrix from element one" style="border-left: 1px solid black; border-bottom: 1px solid black;" | $$k_{41}^{(1)}$$
 * title="stiffness matrix from element one" style="border-right: 1px solid black; border-bottom: 1px solid black;" | $$k_{41}^{(1)}$$
 * title="stiffness matrix from element one and element two" style="border-bottom: 1px solid black;" | $$k_{43}^{(1)} + k_{21}^{(2)}$$
 * title="stiffness matrix from element one and element two" style="border-right: 1px solid black; border-bottom: 1px solid black;" | $$k_{44}^{(1)} + k_{22}^{(2)}$$
 * title="stiffness matrix from element two" | $$k_{23}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black;" | $$k_{24}^{(2)}$$
 * 0
 * style="border-right: 1px solid black;" | 0
 * title="stiffness matrix from element two" | $$k_{31}^{(2)}$$
 * title="stiffness matrix from element two" | $$k_{32}^{(2)}$$
 * title="stiffness matrix from element two" | $$k_{33}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black;" | $$k_{34}^{(2)}$$
 * 0
 * style="border-right: 1px solid black;" | 0
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{41}^{(2)}$$
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{42}^{(2)}$$
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{43}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black; border-bottom: 1px solid black;" | $$k_{44}^{(2)}$$
 * }
 * 0
 * style="border-right: 1px solid black;" | 0
 * title="stiffness matrix from element two" | $$k_{31}^{(2)}$$
 * title="stiffness matrix from element two" | $$k_{32}^{(2)}$$
 * title="stiffness matrix from element two" | $$k_{33}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black;" | $$k_{34}^{(2)}$$
 * 0
 * style="border-right: 1px solid black;" | 0
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{41}^{(2)}$$
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{42}^{(2)}$$
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{43}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black; border-bottom: 1px solid black;" | $$k_{44}^{(2)}$$
 * }
 * title="stiffness matrix from element two" style="border-bottom: 1px solid black;" | $$k_{43}^{(2)}$$
 * title="stiffness matrix from element two" style="border-right: 1px solid black; border-bottom: 1px solid black;" | $$k_{44}^{(2)}$$
 * }
 * }

The box, within the matrix on the top left, is the stiffness matrix from element one. The box to the lower right is the stiffness matrix from element two.

Where capital k (K) is the global stiffness coefficient and lowercase k (k) is the element stiffness coefficient. Below is the relationship between the global stiffness coefficients and the element stiffness coefficients for the global stiffness matrix above.

K11 = k11(1) ; K12 = k12(1) ; K13 = k13(1) ; K14 = k14(1) ; K15 = 0 ; K16 = 0

K21 = k21(1) ; K22 = k22(1) ; K23 = k23(1) ; K24 = k24(1) ; K25 = k25(1) ; K26 = k26(1)

K31 = k31(1) ; K32 = k32(1) ; K33 = k33(1) + k11(2) ; K34 = k34(1) + k12(2); K13 = k13(2) ; K14 = k14(2)

K41 = k41(1) ; K42 = k42(1) ; K43 = k43(1) + k21(2) ; K44 = k44(1) + k22(2); K23 = k23(2) ; K24 = k24(2)

K51 = 0 ; K52 = 0 ; K53 = k31(2) ; K54 = k32(2) ; K55 = k33(2) ; K56 = k34(2)

K61 = 0 ; K62 = 0 ; K63 = k41(2) ; K64 = k42(2) ; K65 = k43(2) ; K66 = k44(2)

Below are selected sample calculations of the global stiffness matrix.

$$K_{11}=\frac{9}{16}$$

$$K_{12}=\frac{3\sqrt{3}}{16}$$

$$K_{33}=K_{33}^{\left(1\right)}+K_{11}^{\left(2\right)}=\frac 9{16}+\frac 52=\frac{49}{16}=3.0625$$

$$K_{34}=K_{33}^{\left(1\right)}+K_{11}^{\left(2\right)}=\frac{3\sqrt 3}{16}+\left(-\frac 52\right)\approx -2.1752$$

$$K_{44}=K_{44}^{\left(1\right)}+K_{22}^{\left(2\right)}=\frac 3{16}+\frac 52=2.6875$$

4) Elimination of known degrees of freedom (dofs).
This Process reduces the global force displacement relationship by eliminating the first and last two columns of the global stiffness matrix.

d1 = d2 = d5 = d6 = 0

The displacements above are zero. Therefore, columns 1, 2, 5 and 6 of the global stiffness matrix go to zero. This reduces the matrix to a 6x2.


title: User:Eas4200c.f08.ZYX.reger/Week5,6


text:
User:EAS4200c.f08.Blue.E/Lecture 16

Continuing the discussion on multicell torsion:

The two important equations for single cell torsion are

$$T=2*q*A$$

and

$$\theta=\frac{1}{2(G)(A)}\oint\frac{q}{t}dl$$

Which is the torque and angle of twist respectively of the cell.

when moving to a multicell view in which cell $$i=1,....,n_{cell}$$

Torque on the multicell system now becomes equal to the sum of all the individual cell torques.

$$T=2*\sum{q_i*A_i}$$

in which

$$q_i$$ is the shear flow in cell i

$$A_i$$ is the "average" area in cell i

so therefore $$T_i=2*q_i*A_i$$

so the torque on the for the entire multicell system would then be the sum of all the $$T_i$$

$$T=\sum_{n=1}^{n_c}T_i$$

the same approach is used for the angle of twist $$\theta$$

in which

$$\theta= \theta_1,.....,\theta_i$$

Therefore

$$\theta_i=\frac{1}{2*G_i*A_i}\oint\frac{q_i}{t_i}dl$$

in which

$$G_i$$ is the shear modulus of cell i

and

$$t_i$$ is the thickness of cell i

so therefore the angl$$e$$ of twist for the entire multicell system would then be the sum of all the $$\theta_i$$

$$\theta=\sum_{n=1}^{n_c}\theta_i$$

Example

$$t_1=.008m$$

$$t_2=t_3=.01m$$

$$a=4m$$

$$b=2m$$

therefore the area will equal the sum of the semicircle and triangle which make up the overall airfoil shape.

$$A=\frac{\pi*(\frac{b}{a})^2}{2}+\frac{1}{2}*b*a=5.5708m $$

and

$$T=2*q*A$$

so therefore

$$q=\frac{T}{2*A}$$

in which T is the variable

so the twist angle will equal

$$\theta=\frac{1}{2*G*A}\sum_{j=1}^{3}\frac{q_j*l_j}{t_j}$$

in which l is the length of the path.

title: User:EAS4200c.f08.Blue.E/Lecture 16


text:
User:Eml4500.f08.delta 6.krueger/Notes Monday September 29th

Use global force displacement relation. $$ \begin{bmatrix} K_{13} & K_{14}\\ K_{23} & K_{24}\\ K_{33} & K_{34}\\ K_{43} & K_{44}\\ K_{53} & K_{54}\\ K_{63} & K_{64} \end{bmatrix} \begin{bmatrix} d_{3}\\ d_{4} \end{bmatrix} = \textbf{F} $$ This matrix is a 6x2 multiplied by a 2x1 matrix and equals a 6x1 matrix. The forces $$F_{3}$$ and $$F_{4}$$are the known applied loads. Therefore, there is no need to compute them. The computations that need to be done are for rows 1,2,5 and 6. This will give The reactions $$F_{1}$$, $$F_{2}$$, $$F_{5}$$ and $$F_{6}$$. Closing the loop between The finite element method (FEM) and statics: Virtual Displacement Sloving the two bar truss system uses a process shown in the figure below. Where side 1 is the finite element method and side two is the statics approach. The statics approach is used to confirm the finite element method and is known as closing the loop.  The process of determinig the displacements using statics was discussed above in "Determination of Element FD (wrt global Coordinates). In the statics process the reactions are known. Therefore, the member forces $$P^{(1)}_{1}$$ and $$P^{(2)}_{2}$$ are known. To compute the axial displacement dof's ($$q$$) or the amount of extension of the bars the member forces and the element stiffiness's need to be used. The following formulas show the relationship between the member forces, element stiffness's and the axial displacements. $$ P=Kq$$ so, $$ q=P/K$$ $$q^{(1)}_{2}= \frac{P^{(1)}_{1}}{K^{(1)}}= \frac{P^{(1)}_{2}}{K^{(1)}}= AC$$ $$q^{(2)}_{1}= \frac{-P^{(2)}_{2}}{K^{(2)}}= \frac{-P^{(2)}_{2}}{K^{(2)}}= -AB$$ The axial displacements dof's at node 1 and 3 are equal to zero.  $$q^{(1)}_{2}=q^{(2)}_{2}= 0$$  Question: How does one back out from the displacement dof's of node 2?  Where: A = node 2, AC is the extension of 1A and AB is the extension of 3A. First an extensions must be drawn perpendicular to AC and AB. Then a line can be drawn between the point D, where the two extensions meet, and node 1. A line is also drawn between node 3 and point D. These two lines are the dotted line and denote the deformed shape of the two-bar truss system. The vector R is the displacement of A or node 2. It is found by taking the square root of the sum of the squared magnitudes of AC and AB. $$R= (AC^{2} + AB^{2})^{\frac{1}{2}}$$ Review of the small angle approximation.

Where A is the angle between the two bars. $$U_{y}= Rsin(A)$$ $$U_{x}= R(1-cos(A))$$ If the angle is small then $$sin(A)$$ is approximately equal to A and $$cos(A)$$ is approximately equal to 1. The above equations become: $$U_{y}= R*A$$ $$U_{x}= 0$$


title: File:FEAHW3.JPG


text:
Enterprise Resource Planning/Project Management

Programme
This is a module part of the Open Source ERP/Executive Diploma course conducted by a private university.

Subject
ERP Project Management - From Inception To Achievement

Objectives
The objectives of this subject are:
 * To know the project schedule areas in a typical ERP Project
 * To go through the challenges and tasks of the phases
 * To be able to interface on the right page with the respective clients in each phase

Learning Outcomes
On completion of this subject the students shall be able to:
 * Walk a target audience through the various phases of a project
 * Analyse and discuss a best approach according to the environment
 * Address user concerns and suggest remedies or workarounds
 * Confident in acting as an ERP consultant that knows the right rules for the right occasion

Sypnosis
This subject goes through an ERP Project cycle, its phases and timelines after the succesful acceptance phase after the Pre-Sales. Each phase has its importance and dependency on the client's business and user environment. Its success is sought by knowing what experts know and how they dealt with the most common scenarios.

Implementation
(15 hrs)
 * Level of ERP Performance
 * Company Wide Implementation
 * Uneven ERP Phases

Training and Education
(15 hrs)
 * Education Plan For the Client Users
 * Case Studies
 * Training Materials

Project Scheduling
(15 hrs)
 * User Priority & Part Timers
 * Team Planning

Sales & Operations Planning
(15 hrs)

Basic ERP
(15 hrs)

Supply Chain Integration
(15 hrs)

Production ShopFloor
(15 hrs)

Future ERP
(15 hrs)

Download Resource

 * Fastrack Methodology


text:
User:Eml4500.f08.RAMROD.A/HW3

=HW 3 Ramrod=

Derivation of elemental force displacement with respect to the global system.
$$ K^{(e)}_{4x4} \cdot d^{(e)}_{4x1} = F^{(e)}_{4x1}$$

Statics dictates that this can be defined as a two-force member. Therefore, it can be drawn like this.

The force displacement relationship is as follows for P and q.

$$K^{(e)}_{2x2} \cdot q^{(e)}_{2x1} = P^{(e)}_{2x1}$$

$$K^{(e)}\begin{vmatrix} 1 & -1\\ -1 & 1 \end{vmatrix} \begin{Bmatrix} q^{(e)}_1\\ q^{(e)}_2 \end{Bmatrix}=\begin{Bmatrix} P^{(e)}_1\\ P^{(e)}_2 \end{Bmatrix}$$

$$q^{(e)}_i$$ = the axial displacement of element e at local node i.

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

Goal : Derive equation $$ K^{(e)}_{4x4} \cdot d^{(e)}_{4x1} = F^{(e)}_{4x1}$$ from $$K^{(e)}_{2x2} \cdot q^{(e)}_{2x1} = P^{(e)}_{2x1}$$ . Find the relationship between matrix $$q^{(e)}$$ and $$d^{(e)}$$ and between matrix $$P^{(e)}$$ and $$F^{(e)}$$.

These relations can be expressed in the form:

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

== Consider the displacement vector of local node 1, denoted by $$d^{(e)}_1 $$ == center text

$$ d^{(e)}_1 =  d^{(e)}_1\vec{i} \cdot d^{(e)}_2\vec{j}$$

$$q^{(e)}_1 $$ = the axial displacement of node 1 is the orthogonal projection of the displacement vector of node one. $$d^{(e)}_1 $$ on the zxis of the bar $$\tilde{x}$$ of element e.

$$q^{(e)}_1=\vec{d}^{(e)}_1\cdot \vec{\tilde{i}}$$

$$q^{(e)}_1=(d^{(e)}_1\vec{i}+d^{(e)}_2\vec{j})\cdot \vec{\tilde{i}}$$

$$q^{(e)}_1=d^{(e)}_1(\vec{i}\cdot \vec{\tilde{i}})+d^{(e)}_2(\vec{j}\cdot \vec{\tilde{i}})$$

$$q^{(e)}_1=l^{(e)}d^{(e)}_1+ m^{(e)}d^{(e)}_2$$

$$q^{(e)}_1=\begin{vmatrix} l^{(e)} & m^{(e)} \end{vmatrix}\begin{Bmatrix} d^{(e)}_1 \\d^{(e)}_2

\end{Bmatrix}$$

Consider the displacement vector of local node 2, denoted by $$d^{(e)}_2$$

$$ d^{(e)}_2 =  d^{(e)}_3\vec{i} \cdot d^{(e)}_4\vec{j}$$

$$q^{(e)}_2 $$ = the axial displacement of node 2 is the orthogonal projection of the displacement vector of node two. $$d^{(e)}_2 $$ on the zxis of the bar $$\tilde{x}$$ of element e.

$$q^{(e)}_2=\vec{d}^{(e)}_2\cdot \vec{\tilde{i}}$$

$$q^{(e)}_2=(d^{(e)}_3\vec{i}+d^{(e)}_4\vec{j})\cdot \vec{\tilde{i}}$$

$$q^{(e)}_2=d^{(e)}_3(\vec{i}\cdot \vec{\tilde{i}})+d^{(e)}_4(\vec{j}\cdot \vec{\tilde{i}})$$

$$q^{(e)}_2=l^{(e)}d^{(e)}_3+ m^{(e)}d^{(e)}_4$$

$$q^{(e)}_2=\begin{vmatrix} l^{(e)} & m^{(e)} \end{vmatrix}\begin{Bmatrix} d^{(e)}_3 \\d^{(e)}_4

\end{Bmatrix}$$

The result of the nodal analysis is combined to form this matrix.

$$\begin{Bmatrix} q^{(e)}_1\\ q^{(e)}_2 \end{Bmatrix}=\begin{vmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0& 0 & -l^{(e)} & -m^{(e)} \end{vmatrix}\begin{Bmatrix} d^{(e)}_1\\ d^{(e)}_2 \\ d^{(e)}_3 \\ d^{(e)}_4

\end{Bmatrix} $$

The Transformation Matrix
In order to relate the axial coordinate system's forces, f(e), to the element coordinate system's axial forces, p(e), it is necessary to employ a transformation matrix, T(e) such that:

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

Looking at this equation, we can tell that T(e) must be a 2X4 matrix in order for the dimensions to work out.

Knowing this, and using the element F-D relationship,

$$\hat{k}^{(e)}*q^{(e)} = p^{(e)}$$

we can substitute in the Transformation Matrix relationship above to get an equation relating the cartesion forces and displacements via the transformation matrix and \hat{k}. The resulting equation is shown below.

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

Rearranging this equation requires moving the T(e) term to the left side to allow us to solve for the internal forces in f(e). Since T(e) is rectangular, the transpose, rather than the inverse, is multiplied by each side resulting in:

$$[T^{(e)T}*\hat{k}^{(e)}*(T^{(e)}]*d^{(e)} = f^{(e)}$$

In order for this to be true, it must be proven that,

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

This proof is given below as the assigned homework.

Now that we have proven the derivation of k^e, we can continue to solving for the global displacement degrees of freedom.

Using the principle of virtual work (i.e. d1 = d2 = d5 = d6 = 0) we can reduce the global equation to two unknown DOF's, d3 and d4.

The next homework assignment was to solve for the eigenvalues of the global stiffness matrix, K.

2 Bar Truss System

Closing the loop refers to the two upward pointing arrows on the right side of the figure, making your way back to computing the reaction forces.

Using Statics
Taking a look back at the statics method, the reaction forces known and unknown therefore the member forces are P1(1), P2(2)
 * Note: Both members are in tension

How to Compute the axial displacement of dofs (degrees of freedom) or the amount of extension of the bars:
$$ q_2^{(1)} = \frac {P_1^{(1)}}{d^(1)} = \frac {P_2^{(1)}}{k^{(1)}} = AC $$
 * P1(1) = P2(1)
 * q1(1) = 0, Fixed global node 1

$$ q_1^{(2)} = \frac {-P_2^{(2)}}{k^{(2)}} = -AB $$

 * q2(2) = 0, Fixed global node 3

Finding The Displacements

Referring back to the picture of the truss system, it can be seen that the displacement of point $$\displaystyle D$$ is simply the displacement of the global node $$\displaystyle 2$$ at point $$\displaystyle A$$. In other words,

$$\displaystyle Y_{D} = d_{4} $$ and

$$\displaystyle X_{D} = d_{3} $$

We can now proceed to calculate the length of lines $$\displaystyle AC $$ and $$\displaystyle AB $$ be computed which as we already know are defined to be the axial forces. In other words it is the amount the bars extend or elongate due to the applied force. As previously stated the formula to find line $$\displaystyle AC $$ is as follows,

$$\displaystyle AC = \frac{\left|P_{2}^1 \right|}{k^1}$$

where

$$\displaystyle \left|P_{2}^1 \right| = 5.1243 $$ $$\displaystyle k^1 = \frac{3}{4}$$

Therefore,

$$\displaystyle AC = 6.8324$$

Remember that $$\displaystyle \left|P_{2}^1 \right|$$ is the internal force $$\displaystyle 2$$ at element $$\displaystyle 1 $$. Remember also that in order to obtain the value of $$\displaystyle \left|P_{2}^1 \right|$$, it is necessary to take the Root Sum Squared, RSS for short, of the two forces making acting at that node. The formula for the RSS of line AC extending from element 1 is,

$$\displaystyle \sqrt{(f_3)^2 + (f_4)^2}$$

The same can be done for $$\displaystyle AB $$ with the axial force from element two.

$$\displaystyle AB = \frac{\left|P_{1}^2 \right|}{k^2}$$ $$\displaystyle \left|P_{1}^2 \right| = 6.276 $$

$$\displaystyle k^2 = 5$$

Therefore

$$\displaystyle AB = 1.2552$$

From this point it is now possible to find the coordinates of points $$\displaystyle B $$ and $$\displaystyle C $$. Simple trigonometry will be used to carry out these calculations. First we will find the coordinates for point C, $$\displaystyle (x_C,y_C) $$. If you recall from previous lessons

$$\displaystyle \theta^1 = 30^{\circ} $$

and we already know the length of $$\displaystyle AC $$ so therefore the coordinates of $$\displaystyle C $$ can be expressed as

$$\displaystyle x_C = AC \cos (30^{\circ}) = 6.8324 *\frac{\sqrt{3}}{2} = 5.917$$

$$\displaystyle y_C = AC \sin (30^{\circ}) = 6.8324 *\frac{1}{2} = 3.4162$$

The coordinates of point $$\displaystyle B $$ can be found in the same way. the angle used is the angle the line $$\displaystyle AB $$ makes with the negative x-axis.

$$\displaystyle x_B = AB \cos (45^{\circ}) = 1.2552 *\frac{\sqrt{2}}{2} = -.88756 $$

$$\displaystyle y_B = AB \sin (30^{\circ}) = 1.2552*\frac{\sqrt{2}}{2} = .886756$$

We next need to find the coordinates of the point $$\displaystyle D $$. There are two unknowns at point $$\displaystyle D $$, the x and the y coordinate $$\displaystyle (x_D,y_D) $$. Since there are two unknowns, we need to find two equations in which these two variables are being utilized. The equation for the lines $$\displaystyle AB $$ and $$\displaystyle AC $$ are the two equations we will utliize. In order to do so we will work with the following figure.

The line $$\displaystyle PQ $$ can be represented in the following way;

$$\displaystyle \vec{PQ} = (PQ) \tilde{i} = PQ [\cos(\theta) \vec{i} + \sin(\theta) \vec{j}] $$

The cosine and the sine of the angle can then be expressed as;

$$\displaystyle \vec{PQ} = (x-x_p) \vec{i} + (y-y_p) \vec{j} $$

Dividing the previous two equations will yield the following;

$$\displaystyle \frac{(y-y_p)}{(x-x_p)} = \tan(\theta) $$

Multiplying through by $$\displaystyle (x-x_p)$$ and the equation of the line $$\displaystyle \vec{PQ} $$ is given by the following;

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

This equation for line $$\displaystyle \vec{PQ} $$ is equivalent to the equation for line $$\displaystyle \vec{AC} $$. The equation for a line that is perpendicular to line $$\displaystyle \vec{AC} $$ and passing throught the point $$\displaystyle A $$ is given by the following;

$$\displaystyle y-y_p= \tan(\theta +\frac{\pi}{2})(x-x_p) $$

This equation is equivalent to the equation for line $$\displaystyle \vec{AC} $$. When the coordinates for lines $$\displaystyle \vec{BD} $$ and $$\displaystyle \vec{CD} $$ are put into the above equation then it is possible to solve for the coordinates of point $$\displaystyle D $$.

$$\displaystyle y_D-y_B= \tan(\theta +\frac{\pi}{2})(x_D-x_B) $$

 $$\displaystyle y_D-y_C= \tan(\theta +\frac{\pi}{2})(x_D-x_C)$$

 $$\displaystyle y_D-.886756= \tan(\frac{\pi}{4})(x_D+.886756) $$

 $$\displaystyle y_D-3.4162= \tan (\frac{2\pi}{3}) - (x_D-5.917)$$

Solving for $$\displaystyle y_D $$ for each equation will allow you to then solve for $$\displaystyle x_D$$. Plug this value back into one of the equations will allow you to get the value of $$\displaystyle y_D$$. When done, the coordinate of point $$\displaystyle D $$ becomes;

$$\displaystyle D = (4.352, 6.1271) $$

Using the Finite Element Method (FEM), we have already found the displacement of point $$\displaystyle D $$. In other words;

$$\displaystyle \vec{AD} = (x_D-x_A) \vec{i} + (y_D-y_A) \vec{j} $$

In this equation the variables $$\displaystyle x_A $$ and $$\displaystyle y_A $$ are zero and the only variables that remain are the coordinates of point $$\displaystyle D $$. By definition;

 $$\displaystyle x_D = d_3 $$ and $$\displaystyle y_D = d_4 $$

And therefore,
$$\displaystyle \vec{AD} = d_3 \vec{i} + d_4 \vec{j} $$

where the coefficients are the global displacements of point $$\displaystyle A $$ to point $$\displaystyle D $$.

Three Truss Bar System
We will now look at the process of analyzing a three bar truss system as shown below.

It is important to note the helpfulness in numbering the local nodes of each element. Doing this in a certain manner make the assembly of the global stiffness matrix more convenient. The following numbering is the best and most convenient system.

The 3 Bar Truss and the Moment Equation
A look at the 3 Bar Truss Setup

From Statics, We know that there are typically three useful equations:

$$\displaystyle \sum{Fx}=0 $$

$$\displaystyle \sum{Fy}=0 $$

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

Statically Indeterminate
For this truss setup, even the moment equation renders no useful information because all reactions pass through the same point, A, and the setup is statically indeterminate

Now, we ask the question: what if we tried to sum the moments about some arbitrary point, B?

We will need to look at another picture and vector math to find out answers:

Can we get Information Summing Moments about B?

A look at the vector algebra:

$$\displaystyle \sum{\overline{M_B}}= \overline{BA} \times \overline{F} $$

Thus, $$\displaystyle \sum{\overline{M_B}}= \overline{BA'} \times \overline{F} $$ for all A' on line of action of F (see image right)

From Vector Algebra: $$\displaystyle \overline{BA'} = \overline{BA}+ \overline{AA'} $$

$$\displaystyle \overline{M_B} = (\overline{BA}+ \overline{AA'}) \times \overline{F} $$

$$\displaystyle \overline{M_B} = \overline{BA} \times \overline{F} + \overline{AA'} \times \overline{F} $$,

With The latter cross product in the previous equation equal to zero, we are left with: $$\displaystyle \overline{M_B} = \overline{BA} \times \overline{F} + 0 $$

More Generally:

$$\displaystyle \sum_{i}^{}{B_i}=\sum_{i}^{}{\overline{BA_i} \times \overline}$$ Where Ai is any point on the line of action of F.

Which Indicates: $$\displaystyle \sum{\overline{M_B}} = \overline{BA} \times \sum{\overline{F}} $$

We are already aware that: $$\displaystyle \sum{\overline{F}} = 0 $$

We can conclude that: $$\displaystyle \sum{\overline{M_B}} = \overline{BA} \times 0 = 0 $$

So, even if we take the moment from any arbitrary point, it will ultimately not provide any helpful information to solve the problem. The problem truly is statically indeterminate.

Three Truss Stiffness Matrix
Now we can solve for overall K by combining K(1), K(2) and K(3)

Matrix K $$\displaystyle K = \begin{bmatrix} k_{33}^{(1)}+k_{11}^{(2)}+k_{11}^{(3)} & k_{34}^{(1)}+k_{12}^{(2)}+k_{12}^{(3)}\\k_{43}^{(1)}+k_{21}^{(2)}+k_{21}^{(3)} & k_{44}^{(1)}+k_{22}^{(2)}+k_{22}^{(3)}\\k_{31}^{(2)}+k_{31}^{(3)} & k_{32}^{(2)}+k_{32}^{(3)}\\k_{41}^{(2)}+k_{41}^{(3)} & k_{42}^{(2)}+k_{42}^{(3)}\\k_{51}^{(3)} & k_{52}^{(3)}\\k_{61}^{(3)} & k_{62}^{(3)} \end{bmatrix} $$  =   $$\displaystyle K = \begin{bmatrix} k_{33} & k_{34}\\k_{43} & k_{44}\\k_{53} & k_{54}\\k_{63} & k_{64}\\k_{73} & k_{74}\\k_{83} & k_{84} \end{bmatrix} $$

MATLAB Code to Plot Undeformed and Deformed Two Bar Truss System
The code below plots the two bar truss system that was solved in class with the deformed system superimposed on top. The plot of the system is below the code.

% Matlab Script %********************************************************************* % filename: two_bar_truss.m  % % PURPOSE: %  Plot the two bar truss system discussed in class and the deformed %  system with matlab % % AUTHOR: %   Daniel Olivero %   Team Ramrod %   EML 4500 FEAD, Fall 2008. % % Modified on: vql, Fri, 03 Oct 2003, 08:37:46 EDT % Created on : Sun, 05 Oct 2008, 16:22:00 EDT % % DEPENDENCIES: %  call: % % REMARKS: % %*********************************************************************  %Plot Two beam truss system % clear; close;

%model with 2-D beam elements dof = 2; %dof per node: axial disp x, 2= disp y %obtain the coordinatates of all nodes % n_node = 3;             %number of nodes n_elem = 2;           %number of elements total_dof = 2 * n_node; %total dof of system

position(:, 1) = [0; 0]; position(:, 2) = [3.46; 2]; position(:, 3) = [4.9; .59];

% print the node coord. for i = 1 : n_node; x(i) = position (1,i); y(i) = position (2,i); end

node_connect (1, 1) = 1;   %element 1 node_connect (2, 1) = 2;

node_connect (1, 2) = 2;   %element 2 node_connect (2, 2) = 3;

% connect all beam elements by connectivity array for i = 1 : n_elem; node_1 = node_connect(1,i); node_2 = node_connect(2,i); xx = [x(node_1),x(node_2)]; yy = [y(node_1),y(node_2)]; axis([ -2 5 -2 5]) plot(xx,yy,'--') hold on end

text(x(node_connect(1,1)),y(node_connect(1,1)),'Global Node 1', 'HorizontalAlignment', 'center') text(x(node_connect(2,1)),y(node_connect(2,1)),'Global Node 2', 'HorizontalAlignment', 'center') text(x(node_connect(2,2)),y(node_connect(2,2)),'Global Node 3', 'HorizontalAlignment', 'center') text(x(node_connect(2,1))/2,y(node_connect(2,1))/2,'Element 1', 'HorizontalAlignment', 'center') text( x(node_connect(2,1)) + (x(node_connect(2,2))-x(node_connect(2,1)))/2, y(node_connect(2,1)) +(y(node_connect(2,2))-y(node_connect(2,1)))/2,'Element 2', 'HorizontalAlignment', 'center')

hold on

position_disp(:,1) = [0;0]; position_disp(:,2) = [7.812; 8.1271]; position_disp(:,3) = [4.9; .59];

% print the node coord. for i = 1 : n_node; x(i) = position_disp (1,i); y(i) = position_disp (2,i); end

node_connect_disp (1, 1) = 1;   %element 1 node_connect_disp (2, 1) = 2;

node_connect_disp (1, 2) = 2;   %element 2 node_connect_disp (2, 2) = 3;

% connect all beam elements by connectivity array for i = 1 : n_elem; node_1 = node_connect_disp(1,i); node_2 = node_connect_disp(2,i); xx = [x(node_1),x(node_2)]; yy = [y(node_1),y(node_2)]; axis([ -2 10 -2 10]) plot(xx,yy,'-') hold on end

text(x(node_connect_disp(2,1)),y(node_connect_disp(2,1)),'Global Node 2', 'HorizontalAlignment', 'center') text(x(node_connect_disp(2,1))/2,y(node_connect_disp(2,1))/2,'Element 1', 'HorizontalAlignment', 'center') text( x(node_connect_disp(2,1)) + (x(node_connect_disp(2,2))-x(node_connect_disp(2,1)))/2, y(node_connect_disp(2,1)) +(y(node_connect_disp(2,2))-y(node_connect_disp(2,1)))/2,'Element 2', 'HorizontalAlignment', 'center')

title ('Two Barr Truss System') xlabel('x') ylabel('y')

Contributing Team Members
Christopher Pergola EML4500.F08.RAMROD.B 19:32, 8 October 2008 (UTC)

EML4500.f08.RAMROD.E 19:19, 7 October 2008 (UTC)

Eml4500.f08.ramrod.D 19:24, 7 October 2008 (UTC)

Eml4500.f08.ramrod.c 03:14, 8 October 2008 (UTC

Chris Fontana Eml4500.f08.RAMROD.F 13:34, 8 October 2008 (UTC)

Eml4500.f08.RAMROD.A 16:24, 8 October 2008 (UTC)

title: User:Eml4500.f08.RAMROD.A/HW3


title: File:Displacement13.JPG


text: File:3Bar Truss.JPG

title: File:3Bar Truss.JPG


text: File:Skew Triangle.JPG

title: File:Skew Triangle.JPG


text: File:3Bar Elem.JPG

title: File:3Bar Elem.JPG


text: File:Skew Triangle2.JPG

title: File:Skew Triangle2.JPG


text:
User:Eas4200c.f08.blue.a/Lecture 16

Continuing the discussion on multicell torsion:

The two important equations for single cell torsion are

$$T=2*q*A$$

and

$$\theta=\frac{1}{2(G)(A)}\oint\frac{q}{t}dl$$

Which is the torque and angle of twist respectively of the cell.

when moving to a multicell view in which cell $$i=1,....,n_{cell}$$

Torque on the multicell system now becomes equal to the sum of all the individual cell torques.

$$T=2*\sum{q_i*A_i}$$

in which

$$q_i$$ is the shear flow in cell i

$$A_i$$ is the "average" area in cell i

so therefore $$T_i=2*q_i*A_i$$

so the torque on the for the entire multicell system would then be the sum of all the $$T_i$$

$$T=\sum_{n=1}^{n_c}T_i$$

the same approach is used for the angle of twist $$\theta$$

in which

$$\theta= \theta_1,.....,\theta_i$$

Therefore

$$\theta_i=\frac{1}{2*G_i*A_i}\oint\frac{q_i}{t_i}dl$$

in which

$$G_i$$ is the shear modulus of cell i

and

$$t_i$$ is the thickness of cell i

so therefore the angl$$e$$ of twist for the entire multicell system would then be the sum of all the $$\theta_i$$

$$\theta=\sum_{n=1}^{n_c}\theta_i$$

Example

$$t_1=.008m$$

$$t_2=t_3=.01m$$

$$a=4m$$

$$b=2m$$

therefore the area will equal the sum of the semicircle and triangle which make up the overall airfoil shape.

$$A=\frac{\pi*(\frac{b}{a})^2}{2}+\frac{1}{2}*b*a=5.5708m $$

and

$$T=2*q*A$$

so therefore

$$q=\frac{T}{2*A}$$

in which T is the variable

so the twist angle will equal

$$\theta=\frac{1}{2*G*A}\sum_{j=1}^{3}\frac{q_j*l_j}{t_j}$$

in which l is the length of the path.

title: User:Eas4200c.f08.blue.a/Lecture 16


text:
User:Eml4500.f08.wiki1.ambrosio/hw1

Section 1 How to Access MATLAB
After logging onto a computer with MATLAB installed, run the system command matlab or click on the desktop icon. To exit MATLAB, enter quit or exit in the MATLAB command window, or click on the File menu and select Exit MATLAB.

When the Editor in MATLAB is opened, it is convenient to have the Editor and Command windows opened separately. The specifics of the Editor are covered later on in this tutorial.

Section 2 Entering Matrices
MATLAB is based around operating with matrices. In fact, MATLAB is an abbreviation for Matrix Laboratory. The program can interpret these matrices as either scalars or vectors.

Matrices can be created by entering an explicit list of elements, generating built in statements and functions, creating a disk file with a local editor, or loading external data files or applications.

The following are all examples for creating matrix A, which is a 3x3 matrix. A = [1 2 3; 4 5 6; 7 8 9]

A = [1,2,3;4,5,6;7,8,9]

A = [ 1 2 3	4 5 6	7 8 9]

Notice that the elements in the rows of the matrix can be separated by commas or blank spaces.

The following is an example of how to enter a number in exponential form.

x = 2.34e-9

Notice when the number in entered, there are no spaces in between the e.

MATLAB is capable of doing computations with complex numbers. The following are examples of how to create complex matrices.

A = [1 2;3 4] + i*[5 6;7 8]

A = [1+5i 2+6i;3+7i 4+8i]

Notice when the complex numbers were listed inside the matrix, no blank spaces were used. It should also be noted that either i or j can be used as the imaginary unit. Also, if i or j are assigned as variables to represent a new value, they will not retain their value as an imaginary number. For example if the user enters the command

i = 5

then i will no longer be equal to 'sqrt(-1)', it will be equal to 5.

If a large matrix must be entered, creating an ASCII file is a good option because it can make it easier to correct errors. If an ASCII file is created and named data.ext, simply type load data.ext into the MATLAB command line and press enter. This command will read the file and save it as the variable data. MATLAB also has built in functions such as rand, magic, and hilb. If rand(n) is typed into the command line, MATLAB will create a n x n matrix with randomly generated numbers distributed uniformly between 0 and 1. If rand(m,n) is typed into the command line, MATLAB will create a m x n matrix with the same criteria as before. If magic(n) is typed into the command line, MATLAB will create an integral n x n matrix which is a magic square. This means all the rows, columns, and diagonals have a common sum. If hilb(n) is typed into the command line, MATLAB will create a n x n Hilbert matrix.

Individual entries in a matrix can be easily referenced. For example,

A(2,3)

selects the value in the second row, third column of matrix A. Also

X(3)

Would select the third coordinate of vector X.

Section 3 Matrix and Array Operations
The following matrix operations are available in MATLAB:

+	addition -	Subtraction ^	power ‘	conjugate transpose \	left division /	right division
 * multiplication

Since a scalar value entered into MATLAB is stored as a 1 x 1 matrix, these operations also apply to those values. If MATLAB is told to perform an operation that is not possible for a specific matrix, an error message will appear on the screen.

Matrix division can sometimes be confusing in MATLAB. The difference in right and left division is shown below.

x = A\b is the solution of A*x = b

x = b/A is the solution of x*A = b (another way of writing this could be x = (A’\b’)’ )

It is also important to note that addition and subtraction operate entry wise, but multiplication, power, and left and right division operate as matrix operations. These matrix operations can operate as entry wise operations by placing a period in front of them. An example is

[1,2,3,4].*[1,2,3,4]	or 	[1,2,3,4].^2

Both of these operations will give [1,4,9,16] as the solution.

Section 4 Statements, Expressions, and Variables and Saving Sessions
MATLAB uses expression language. This means when values are typed in, they are entered in the form of variable = expression or just as expression. For example,

x = 5*4

The evaluation of this expression is created as a matrix and displayed on the screen, such as

x = 20

If a variable and equal sign are not used, the result is automatically saved as the variable ans.

The user can also decide to continue the statement to the next line by using three or more periods (…) followed by enter, or several statements can be placed on a single line by separating them with commas (,) or semicolons (;).

If the variable is rather large, it may be useful to suppress the statement. By placing a semicolon (;) at the end of a statement, the assignment is carried out, but the result will not be printed on the screen.

It is also important to note that MATLAB is not case sensitive. For example, alpha is not the same as Alpha.

In order to see a list of variables that are currently being used in the workspace, type who (or whos) in the command line. If a variable needs to be cleared from the workspace, enter clear variablename'. If clear is entered into the command line, all nonpermanent variables will be cleared.

If the unit round off needs to be changed, type in eps (epsilon) to have MATLAB use a round off of 10^-16. This can be useful in specifying tolerances for convergence of iterative processes.

Sometimes a computer will enter a loop that does not stop or takes too long to calculate. These computations can be stopped in MATLAB by pressing CTRL-c.

When MATLAB is closed, all of the variables in the workspace are deleted. If these variables need to be saved, enter save in the command line before closing the program. These variables are written to a disk file named matlab.mat. When MATLAB is started again, enter the command load and the variables will be restored to the workspace.

Section 5 Matrix Building Functions
The following are convenient matrix building functions:

eye		identity matrix zeros 		matrix of zeros ones		matrix of ones diag		create or extract diagnols triu		upper triangular part of a matrix tril		lower triangular part of a matrix rand		randomly generated matrix hilb		Hilbert matrix magic		magic square toeplitz	see help toeplitz

For example, zeros(m,n) creates a m x n matix of zeros, and zeros(n) produces a n x n matrix of zeros. If A has already been defined as a matrix, then zeros(size(A)) will create a matrix of zeros that is the same size as matrix A.

If x has already been defined as a vector, diag(x) creates a diagonal matrix with x down the diagonal. If A has already been defined as a square matrix, then diag(A) creates a vector consisting of the diagonal of A.

Matrices can also be built from blocks consisting of other matrices. For example, if A is defined as a 3 x 3 matrix, then B = [A, zeros(3,2); zeros(2,3), eye(2)] will be 5 x 5 matrix.

title: User:Eml4500.f08.wiki1.ambrosio/hw1


text:
User:Eml4500.f08.wiki1.ambrosio/hw2

Notes 9/8
HW1 left off with the description of elements as springs and explained the 6 steps of the FEM process.

In order to continue on with the FEM process and develop the element stiffness matrix, an understanding of director cosines must exits. The director cosines of the $$\tilde{x}$$ axis with respect to the global ($$ \ x $$,$$ \ y $$) coordinates are derived below from the figure below.

\ l^{(e)} =  \cos \theta^{(e)} $$

\ m^{(e)} =  \cos (90^{o} - \theta^{(e)}) $$

\ m^{(e)} =  \sin \theta^{(e)} $$

A further proof is presented below to show the relationship of the director cosines.

 * $$ \bar{\tilde{i}} = \cos \theta^{(e)} \bar{i} + \sin \theta^{(e)} \bar{j}  $$

\begin{array}{lcl} \begin{align} \bar{\tilde{i}} \cdot \bar{i} & = (\cos \theta^{(e)} \bar{i} + \sin \theta^{(e)} \bar{j}) \cdot \bar{i} \\ & = \cos \theta^{(e)} \bar{i} \cdot \bar{i} + \sin \theta^{(e)} \bar{j} \cdot \bar{i} \\ & = \cos \theta^{(e)} (1) + \sin \theta^{(e)} (0)\\ & = \cos \theta^{(e)} \end{align} \\ l^{(e)} =  \cos \theta^{(e)} \end{array} $$

\begin{array}{lcl} \begin{align} \bar{\tilde{i}} \cdot \bar{j} & = (\cos \theta^{(e)} \bar{i} + \sin \theta^{(e)} \bar{j}) \cdot \bar{j} \\ & = \cos \theta^{(e)} \bar{i} \cdot \bar{j} + \sin \theta^{(e)} \bar{j} \cdot \bar{j} \\ & = \cos \theta^{(e)} (0) + \sin \theta^{(e)} (1)\\ & = \sin \theta^{(e)} \end{align} \\ m^{(e)} =  \sin \theta^{(e)} \end{array} $$

In order to describe the global force displacement at the element level, the k(e) matrix must first be determined. The k(e) matrix for a 4x4 matrix is defined below.

\mathbf{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} $$

Where

 * $$ k^{(e)} =

\frac{E^{(e)} A^{(e)}} {L^{(e)}} $$

And

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

m-file
% Two bar truss example clear all; e = [3 5]; A = [1 2]; P = 7; L=[4 2]; alpha = pi/3; beta = pi/4;

nodes = [0, 0; L(1)*cos(pi/2-alpha), L(1)*sin(pi/2-alpha); L(1)*cos(pi/2-alpha)+L(2)*sin(beta),L(1)*sin(pi/2-alpha)-L(2)*cos(beta)];

dof=2*length(nodes);

conn=[1,2; 2,3]; lmm = [1, 2, 3, 4; 3, 4, 5, 6]; elems=size(lmm,1); K=zeros(dof); R = zeros(dof,1); debc = [1, 2, 5, 6]; ebcVals = zeros(length(debc),1);

%load vector R = zeros(dof,1); R(4) = P;

% Assemble global stiffness matrix K=zeros(dof); for i=1:elems lm=lmm(i,:); con=conn(i,:); k_local=e(i)*A(i)/L(i)*[1 -1; -1 1] k=PlaneTrussElement(e(i), A(i), nodes(con,:)) K(lm, lm) = K(lm, lm) + k; end K R % Nodal solution and reactions [d, reactions] = NodalSoln(K, R, debc, ebcVals) results=[]; for i=1:elems results = [results; PlaneTrussResults(e, A, ...           nodes(conn(i,:),:), d(lmm(i,:)))]; end format short g

Results
k_local =

0.75       -0.75        -0.75         0.75

k =

0.5625     0.32476      -0.5625     -0.32476      0.32476       0.1875     -0.32476      -0.1875      -0.5625     -0.32476       0.5625      0.32476     -0.32476      -0.1875      0.32476       0.1875

k_local =

5   -5    -5     5

k =

2.5        -2.5         -2.5          2.5         -2.5          2.5          2.5         -2.5         -2.5          2.5          2.5         -2.5          2.5         -2.5         -2.5          2.5

K =

0.5625     0.32476      -0.5625     -0.32476            0            0      0.32476       0.1875     -0.32476      -0.1875            0            0      -0.5625     -0.32476       3.0625      -2.1752         -2.5          2.5     -0.32476      -0.1875      -2.1752       2.6875          2.5         -2.5            0            0         -2.5          2.5          2.5         -2.5            0            0          2.5         -2.5         -2.5          2.5

R =

0    0     0     7     0     0

d =

0           0        4.352       6.1271            0            0

reactions =

-4.4378     -2.5622       4.4378      -4.4378

results =

1.7081      5.1244       8.5406       5.1244       17.081       0.6276       1.8828        3.138       1.8828        6.276

title: User:Eml4500.f08.wiki1.ambrosio/hw2


title: File:3cellsection.JPG


title: File:Airfoil4.JPG


text:
User:Navy blue84

I want to learn HTML and some other stuff in my spare time. So I decided Wikiversity would be the best place to do it!!!

title: User:Navy blue84


text:
User:Eml4500.f08.gravy.bje/Notes3

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

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

title: User:Eml4500.f08.gravy.bje/Notes3


title: File:L11elem1.jpg


title: File:L11elem2.jpg


title: File:Diagram31.jpg


text: File:Compare Beam.png

title: File:Compare Beam.png


title: File:Infinitesimaldisp1.JPG


text:
Portal talk:European History/Medieval History

Is anyone still planning on developing this? If not, I'd like to merge this with Topic:European History. The Jade Knight (d'viser) 19:26, 7 October 2008 (UTC)


text:
User:EAS4200C.F08.WIKI.A/Homework3

Curved Panels
Below is an image of a curved panel in the (x,y,z) coordinates.

Note: Sometimes (y,z) coordinates are used for the axes in plane of the cross section of the panel instead of (x,y) coordinates.

Below is an image of another curved panel with thickness t. A force will be applied to the panel which can be seen in the next image below.

Below is a cross sectional view of the same curved panel with force F acting on it at an angle above the horizontal (y-axis). The cross section goes from point A to point B. It's height is b and it's width is a.

$$d\overrightarrow{F} = q\overrightarrow{dl} = q(dl_yj + dl_zk)$$

$$d\overrightarrow{F} = q(dl\cos\theta \overrightarrow{j} + dl\sin \theta \overrightarrow{k})$$

$$\displaystyle\ dl\cos\theta = dy$$

$$\displaystyle\ dl\sin\theta = dz$$

$$\displaystyle\ q = \tau t $$ = shear flow

Here's a closer look at the panel where the force is acting. You can see it's acting on and parallel with section dl.

Closed Thin-Walled Cross Section
$$\overrightarrow{T} = T\overrightarrow{i} $$

Below is the cross section of a closed thin-walled area. There is a force acting on section dl.

$$\overrightarrow{r}$$ = the distance from the origin to section dl

$$\rho$$ = the distance from the origin to the force that is perpendicular to the force vector.

$$ d\overrightarrow{T} = \overrightarrow{r}xd\overrightarrow{F} $$

$$ dT = \rho dF $$

Take a closer look where the force acts on on section dl.

$$dA = \frac{1}{2}\rho dl $$

$$T = \oint_{}^{}{dT} = q\oint_{}^{}{\rho dl} $$

 * $$ \rho dl = 2dA $$

$$ T = 2q\int_{\overline{A}}^{}{dA}$$

$$T = 2q\overline{A}$$

$$\overline{A}$$ = Average Area

Below is an image of the average area of the cross section of a curved panel:

=Open Thin Walled Cross-Sections=

The torque of the depicted thin walled structure is

 $$\displaystyle T = 2q\bar{A}$$

This is also equivalent to the location of the resultant force acting parallel to the line connecting the two ends of the object.

=Uniform Circular Cylinder: Non-Warping=

We will now take a look at a uniform circular cylinder cross-section. This cross-section behaves as a rigid disk.

We will look at the infinitesimal area $$\displaystyle dA $$ at a distance $$\displaystyle r $$ away from the origin. The torque here can be described by the following equation;

$$\displaystyle T = \int\int_{A} r\tau dA $$

It must be noted that;

$$\displaystyle \tau = G\gamma $$, $$\displaystyle \tau dA = F $$

where $$\displaystyle G $$ is the strain modulus and $$\displaystyle \gamma $$ is the shear strain. Together they comprise Hooks Law which is used to describe the deformation of the disk. The variable $$\displaystyle F $$ is the force and when multiplied by the displacement produces the moment around the origin.The shear strain is described by the following equation;

 $$\displaystyle \gamma = r \frac{d\alpha}{dx}$$

where $$\displaystyle \frac{d\alpha}{dx} = \theta $$ which is defined as the rate of twist. Since $$\displaystyle \tau $$ is proportional to $$\displaystyle \gamma $$ which is also proportional to $$\displaystyle r $$, the shear stress $$\displaystyle \tau $$ gets smaller in the material closer to the core because the distance away from the core $$\displaystyle r $$ gets smaller. Therefore it is more efficient and cost effective to remove material to make the cylinder lighter. The limits of the material will be tested but if calculations are done correctly then there shouldn't be anything to worry about.

The torque can be represented by the following;

$$\displaystyle T = \int _A r G (r\theta) dA = G\theta\int _Ar^2dA $$

where $$\displaystyle dA $$ is a function of $$\displaystyle dy $$ and $$\displaystyle dz $$.

$$\displaystyle dA = dy * dz = r(dr)(dw)$$

The strain modulus and rate of twist is not a function of $$\displaystyle z $$ or $$\displaystyle y $$ and is why they come out of the integral. The resulting integral is defined as the second polar moment of inertia.

$$\displaystyle J = \int _Ar^2dA $$

=Hollow Circular Cross-Section: Thin Walled=

The analysis of a hollow circular cross section with thin walls (ie. t << a), is conducted as follows:

$$r_{i} = a$$ (inner radius)  $$r_{o} = b$$ (outer radius)

The equation for the polar moment is done by superimposing 'empty' circle on top of one solid circle, and the resulting equation is as follows:

 $$J = \frac{1}{2}\pi (b^{4} - a^{4})$$

This can be expanded as:

 $$J = \frac{1}{2}\pi (b - a)(b + a)(b^{2} + a^{2}$$

If an average radius of the cross section is then defined as:

 $$\bar{r} = \frac{a + b}{2}$$

And it is noted that, since t is very small,

$$b \cong a$$

and

$$ b \cong \bar{r}$$

and

$$ a \cong \bar{r}$$

Therefore,

$$ a^{2} \cong \bar{r}^{2}$$

$$ b^{2} \cong \bar{r}^{2}$$

and noting that $$b - a = t$$

J can then be expressed as:

 $$J = 2\pi t\bar{r}^{3}$$

J can then be expressed in terms of the average radius by noting that:

$$J = (2\pi ^{-\frac{1}{2}}t)(\pi \bar{r}^{2})^{3/2}$$

With $$\pi \bar{r}^{2} = \bar{A}$$ it's easily seen that J is proportional to $$\bar{A}^{3/2}$$, with $$(2\pi ^{-\frac{1}{2}}t)$$ being a proportionality factor.

=Solid Circular Cross-Section vs Hollow Thin Walled Cross-Section=

In order to compare these two cross-sections, we will be using the figure from the text Mechanics of Aircraft Structures by C.T. Sun. The figure is depicted below with all dimensions labeled on the figure.

We will first compute the area for each cross-section.

$$\displaystyle  A_a = \pi r^2 = \pi (1cm)^2 = 3.1416 cm^2 $$

$$\displaystyle A_b = \pi r_o^2 - \pi r_i^2 = \pi (5.1cm)^2 - \pi (5cm)^2= 3.173cm^2$$

As the computtion shows, the areas are similar despite the size difference. We will next compute the polar moment of inertia for each cross section.

$$\displaystyle J_a = \frac{\pi r^4}{2} = \frac{\pi (1cm)^2}{2} = 1.5707cm^4 $$

$$\displaystyle J_b = \frac{\pi (r_0^4 - r_i^4)}{2} = \frac{\pi ((5.1cm)^4-(5cm)^4)}{2} = 80.927cm^4$$

Comparing the two Polar moment of inertia shows that the thin walled cross-section is a better torsional member than the solid cross-section. Since its polar moment of inertia (measure of its resistance to torsion) is much larger, it is a more desirable cross-section to use. The polar moment of inertia for the solid cross-section is only $$\displaystyle 2\% $$ of the polar moment of inertia for the hollow thin walled cross-section.

$$\displaystyle \frac{J_a}{J_b} = \frac{1.5707cm^4}{80.927cm^4} = .0194 $$

If the thickness to equal $$\displaystyle .02 r_i^{(c)}$$, where $$\displaystyle c $$ represents a new hollow thin walled cross-section, what would this inner radious $$\displaystyle r_i^{(c)}$$ equal such that the polar moment of inertia for figure (a) equals that of the new figure (c)?

Matlab Project Discription
For HW 3, the team was asked to analyze a NACA 4-digit airfoil series. The discription of the problem can be seen here. A diagram showing the shape of the airfoil in this series can be found on this website. The group is asked to break up the airfoil into a number (ns) of discrete sections in order to analyze the geometry of the airfoil.

In order to complete this problem the student must first understand how to discretize an object into triangular segments in order to find the area of the object. A diagram of how this is done is provided below.

Using this, the following equations can be used in order to find the area of the shaded region.

$$d\vec{T} = \vec{r} \times d\vec{F} = q\vec{r} \times \vec{PQ}$$

Here, $$d\vec{F} = \vec{PQ}$$ and $$q\vec{r} \times \vec{PQ} = (2dA)\vec{i}$$

With

$$||\vec{PQ}|| = d\mathit{l} $$

Now that the basic foundation of the analysis has been discussed, the MATLAB code along with its execution is displayed below.

{| class="toccolours collapsible collapsed" style="width:100%" ! MATLAB Code

Airfoil Function
function Airfoil(m,p,ta,c,P_0)

m = m/100; p = p/10; ta = ta/100; 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

Area = 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;   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;   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

A_bar = A_BAR(j-1)

subplot(2,1,1), plot(y_u,z_u,y_l,z_l), 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')

NACA 2415 Airfoil in MATLAB command window
EDU>> Airfoil(2,4,15,.5,[0 0 0])

A_bar =

0.0253

EDU>> Airfoil(2,4,15,.5,[0 0.19 0.01])

A_bar =

0.0255

EDU>> Airfoil(2,4,15,.5,[0 0.25 -1])

A_bar =

0.0256

 * }

Torsion
Torsion of uniform, non circular bars leads to warping of the cross section, which is defined as an axial displacement along the x-axis (ie, along the bars length due to the conventions used in class) of a point on the deformed, or rotated, cross section. A schematic representation is shown below.

$$u_{y}$$ = the y-component of the displacement vector $$\vec{PP^{'}}$$

and

$$u_{z}$$ = the z-component of the displacement vector $$\vec{PP^{'}}$$

It should be noted (and easily seen) from the figure that $$\vec{PP^{'}}$$ is not on a circle centered at the origin with radius $$\vec{OP^{'}}$$ and $$\vec{OP^{'}}$$ in fact has a larger magnitude than $$\vec{OP}$$

This is due to the line $$\vec{PP^{'}}$$ having a vertical and horizontal displacement from $$\vec{PP^{'}}$$

However, since the angles associated with twist in torsion are likely to be relatively small, one can use small angle approximations and assume that the magnitude of $$\vec{OP} = \vec{OP^{'}} = R$$

$$u_{z} = Rsin\alpha \cong R\alpha$$

$$u_{y} = R(1 - cos\alpha )$$

With $$\alpha = 0$$, this makes $$u_{y} = 0$$ and $$u_{z} = \alpha R$$

The rate of twist angle, $$\displaystyle \theta $$, is defined as

 * $$\displaystyle \theta = {\alpha \over x} $$

Where $$\displaystyle \alpha $$ is the angle of rotation. The horizontal displacement (displacement in the y-direction) is

 * (1) $$\displaystyle u_y = -\theta x z $$

The vertical displacement of a point of the cross section (displacement in the z-direction) can be determined from

 * $$\displaystyle u_z = +(PP') \cos \beta = +\alpha \zeta_p $$

Where $$\displaystyle (OP)\alpha $$ can be substituted in place of $$\displaystyle PP' $$ and $$\displaystyle \theta_x $$ can be substituted in place of $$\displaystyle \alpha $$, thus resulting in the following

 * (2) $$\displaystyle u_z = +\theta x y $$

The displacement due to warping along the x-axis is defined as

 * (3) $$\displaystyle u_x = \theta \psi (y,z) $$

Equations 1, 2, and 3 make up the kinematic assumptions.

=Roadmap for Torsional Analysis=

The following list lays out a general step by step roadmap of equations to apply for the torsional analysis of an aircraft wing:

 * A. Kinematic assumptions (equations 1, 2, and 3 above)
 * B. Strain-displacement relation given by the two equations:
 * $$\displaystyle \gamma_{yz} = {\partial u_x \over \partial y} + {\partial u_y \over \partial x}$$
 * $$\displaystyle \gamma_{zx} = {\partial u_x \over \partial z} + {\partial u_z \over \partial x}$$
 * C. Equilibrium equations for stresses:
 * $$\displaystyle {\partial \sigma_{yy} \over \partial y} + {\partial \tau_{zy} \over \partial z} + {\partial \tau_{xy} \over \partial x} = 0 $$
 * $$\displaystyle {\partial \tau_{yz} \over \partial y} + {\partial \sigma_{zz} \over \partial z} + {\partial \tau_{xz} \over \partial x} = 0 $$
 * $$\displaystyle {\partial \tau_{yx} \over \partial y} + {\partial \tau_{zx} \over \partial z} + {\partial \sigma_{xx} \over \partial x} = 0 $$
 * But, since $$\displaystyle \tau_{yx} $$ and $$\displaystyle \tau_{zx} $$ are independent of x in the displacement field, the equations reduce to the following:
 * $$\displaystyle {\partial \tau_{yx} \over \partial y} + {\partial \tau_{zx} \over \partial z} = 0 $$
 * D. Prandt stress funtion $$\displaystyle \phi (y,z) $$:
 * $$\displaystyle \tau_{yx} = {\partial \phi \over \partial z} $$
 * $$\displaystyle \tau_{zx} = -{\partial \phi \over \partial y} $$
 * E. Strain compatibility equation:
 * $$\displaystyle {\partial \gamma_{zx} \over \partial y} - {\partial \gamma_{yx} \over \partial z} = 2 \theta $$
 * F. Equation for $$\displaystyle \phi $$
 * $$\displaystyle {\partial^2 \phi \over \partial y^2} - {\partial^2 \phi \over \partial z^2} = -2G \theta $$
 * G. Boundary conditions for $$\displaystyle \phi $$:
 * The traction free boundary condition $$\displaystyle t_x = 0 $$ on a lateral surface yields the following:
 * $$\displaystyle {d \phi \over ds} = 0 $$
 * thus, $$\displaystyle \phi $$ is a constant. But, this constant becomes arbitrary since it is a solid section with only one contour boundary, thus
 * $$\displaystyle \phi = 0 $$
 * H. Resultant torque:
 * $$\displaystyle T = 2 \iint_A \phi dA $$
 * I. $$\displaystyle T = G J \theta $$ where
 * $$\displaystyle J = -{4 \over \triangledown^2 \phi} \iint_A \phi dA $$
 * $$\displaystyle T = 2q \bar A $$
 * J. Twist angle
 * $$\displaystyle \theta = {1 \over 2G \bar A} \oint {q \over t}ds $$
 * where $$\displaystyle s $$ is the curvilinear coordinate along a thin wall

 * K. 1) Multicelled Section

 * For cells where $$i = 1, ..., n_{cell}$$

 * [[Image:3cellsection.JPG|500px]]

 * $$\tau = 2\sum_{i=1}^{N_c}{q_i \bar{A_i}} $$

 * $$q_i \;$$ = Shear flow in cell i.

 * $$\bar{A_i}$$ = "average" area in cell i.

 * Define: $$T_i = 2q_i\bar{A_i}$$ torque generated by one cell.

 * $$\tau = 2\sum_{i=1}^{N_c}{q_i \bar{A_i}} $$

 * $$T_i = \sum_{i=1}^{N_c}{T_i} $$

 * 2) Shape of airfoil is "rigid" in the plane $$(X,Y)$$ however it can warp out of the plane.

 * $$\theta = \theta _1 = \theta _2= ...=\theta _{N_c}$$

The twist angle equation applied to cell i.

$$\theta_i =\frac{1}{2 G_i \bar{A_i}}\oint_{}^{}{\frac{q_i}{t_i}ds}$$

 * $$G_i\;$$ = Shear Modulus of cell i

 * $$t_i(s)\;$$ = thickness of walls on cell i. This is a linear curve coordinate along the cell wall.

''' Pb 1.1: Rectangular single cell section with modification. '''

Now more general with a single cell section:

$$t_i=0.08 \; m$$

$$t_2= t_3=0.01 \; m$$

$$\bar{A}=\frac{1}{2} \pi (\frac{b}{2})^2 +\frac{1}{2}ba= 5.5708$$

Shear flow comes from the torque equation: $$T=2q\bar{A}$$

$$q= \frac{T}{2\bar{A}}$$

Twist angle: $$\theta_i =\frac{1}{2 G_i \bar{A_i}}\sum_{j=1}^{3}{\frac{q_jl_j}{t_j}}$$ $$j=1,2,3$$ index for the segment number.

There is only one cell for $$\theta_i\;$$ so it reduces to $$\theta\;$$ in the equation above. You would normally have two indexes for multicelled section $$(i,j)\;$$.

$$\theta =\frac{1}{2 G \bar{A}}\sum_{j=1}^{3}{\frac{q_jl_j}{t_j}}$$

Note: $$\int$$ is and elongated S, standing for continuous summation.

$$\sum{}$$ is a discrete sum.

=Individual Contributions= Christopher PergolaEAS4200C.F08.WIKI.A 18:14, 8 October 2008 (UTC)

Braden Barnes Eas4200c.f08.wiki.b 19:54, 7 October 2008 (UTC)

Eric Viale Eas4200c.f08.wiki.c 22:15, 7 October 2008 (UTC)

Michael Rolle Eas4200c.f08.wiki.d 18:29, 8 October 2008 (UTC)

Jeff Huenink Eas4200c.f08.wiki.f 03:10, 8 October 2008 (UTC)

Chris Fontana Eas4200c.f08.WIKI.E 13:11, 8 October 2008 (UTC)

title: User:EAS4200C.F08.WIKI.A/Homework3


text:
User:EML4500.f08.JAMAMA/Homework3

September 22, 2008
$$ \mathbf{q}_{2 \times 1}^{(e)} = \mathbf{T}_{2 \times 4}^{(e)} * \mathbf{d}_{4 \times 1}^{(e)} \qquad \qquad \qquad \qquad \qquad \qquad (1) $$

Derivation of element F-D with respect to global coordinate sstem (p. 6-1)

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

(p. 4-5) $$ \ \ k^{(e)} \begin{bmatrix}

1 & -1 \\ -1 & 1

\end{bmatrix}

\begin{Bmatrix}

q_1^{(e)} \\ q_2^{(e)}

\end{Bmatrix}=

\begin{Bmatrix}

P_1^{(e)} \\ P_2^{(e)}

\end{Bmatrix} \qquad \qquad (2)$$

$$ q_1^{(e)}\ $$ = axial displacement of element $$ \mathbf{e}\ $$ at local node

$$ P_1^{(e)}\ $$ = axial force of element $$ \mathbf{e}\ $$ at local node

Global: Derive eq. $$ (1) $$ from eq. $$(2)$$ $$ \rightarrow $$ equation $$(2)$$ was already derived in Mtg.4

Want to find relation between $$ \begin{matrix} \mathbf{q}_{2 \times 1}^{(e)} & and & \mathbf{d}_{4 \times 1}^{(e)} \\ \mathbf{P}_{2 \times 1}^{(e)} & and & \mathbf{f}_{4 \times 1}^{(e)} \end{matrix} $$

These relations can be expressed in the form: $$ \mathbf{q}_{2 \times 1}^{(e)} = \mathbf{T}_{2 \times 4}^{(e)} * \mathbf{d}_{4 \times 1}^{(e)} $$

Consider the displacement vector of local node, denoted by $$ \mathbf{d}_{(1)}^{(e)} $$

$$ \Rightarrow $$ $$ \vec {\mathbf{d}}_{(1)}^{(e)} =  \mathbf{d}_1^{(e)}\vec {i} +  \mathbf{d}_2^{(e)} \vec {j}$$

$$\Rightarrow \ q_{1}^{(e)} $$ = axial displacement of node is the orthogonal projection of the dispalcement. Vector $$ \vec {\mathbf{d}}_{1}^{(e)}$$ of node on the axis of $$ \tilde{i} $$ of element $$ \mathbf{e} $$

$$ q_{1}^{(e)}= \vec {\mathbf{d}}_{(1)}^{(e)}* \vec{\tilde{i}} \Rightarrow (\mathbf{d}_1^{(e)} \vec {i} +  \mathbf{d}_2^{(e)} \vec {j})* \vec{\tilde{i}}= \mathbf{d}_1^{(e)}( \vec{i} * \vec{\tilde{i}} ) + \mathbf{d}_2^{(e)} (\vec {j}* \vec{\tilde{i}}) \Rightarrow $$

$$ \Rightarrow \begin{matrix} \vec{i} * \vec{\tilde{i}} & = & \cos(\theta^{(e)}) & = & l^{(e)} \\ \vec{j} * \vec{\tilde{i}} & = & \sin(\theta^{(e)}) & = & m^{(e)} \end{matrix} \Rightarrow \ $$ $$q_{1}^{(e)}= l^{(e)} \mathbf{d}_1^{(e)} + m^{(e)}\mathbf{d}_2^{(e)} = {\llcorner l^{(e)} m^{(e)} \lrcorner}_{1 \times 2}

{\begin{Bmatrix}

d_{1}^{(e)} \\ d_{2}^{(e)}

\end{Bmatrix}}_{2 \times 1}$$

Similarly for node, $$ q_{2}^{(e)}= {\llcorner l^{(e)} m^{(e)} \lrcorner}_{1 \times 2} {\begin{Bmatrix}

d_{3}^{(e)} \\ d_{4}^{(e)}

\end{Bmatrix}}_{2 \times 1}$$

$$ \begin{matrix}

\begin{Bmatrix}

q_{1}^{(e)} \\ q_{2}^{(e)} \end{Bmatrix} & = &

\begin{bmatrix}

l^{(e)} & m^{(e)} & 0 & 0 \\ 0 & 0 & l^{(e)} & m^{(e)} \end{bmatrix} &

\begin{Bmatrix}

d_{1}^{(e)} \\ d_{2}^{(e)} \\ d_{3}^{(e)} \\ d_{4}^{(e)} \end{Bmatrix} \\

& & & \\

\mathbf{q}_{2 \times 1}^{(e)} & = & \mathbf{T}_{2 \times 4}^{(e)} & \mathbf{d}_{4 \times 1}^{(e)} \end{matrix}

$$

September 26, 2008
Similarly (same argument)

$$\begin{Bmatrix} P{_{1}}^{(e)}\\ P{_{2}}^{(e)}\\ \end{Bmatrix} =  T^{(e)} \begin{Bmatrix} f{_{1}}^{(e)}\\ f{_{2}}^{(e)}\\ f{_{3}}^{(e)}\\ f{_{4}}^{(e)} \end{Bmatrix} $$  2X1 2X4 4X1

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

Recall element axial force displacement relationship

$$\hat{k}^{(e)}q^{(e)}=P^{(e)}$$ 2X2 2X12X1

$$\Rightarrow \hat{k}^{(e)} \underbrace{\left(T^{(e)}d^{(e)} \right)}_{\displaystyle{q^{(e)}}}= \underbrace{\left(T^{(e)}f^{(e)} \right)}_{\displaystyle{p^{(e)}}}$$

Goal: Want to have $$k^{(e)}d^{(e)}=f^{(e)}$$ so "move" $$T^{(e)}$$ from right to left hand side by pre-multiplying equation by $$T^{(e)-1}$$ (inverse of $$T^{(e)}$$ ). Unfortunately, $$T^{(e)}$$ is a rectangular matrix (2X4) so it cannot be inverted.

So the actual answer is: $$\begin{bmatrix} T{^{(e)T}}\bullet \hat{k}^{(e)}T^{(e)} \end{bmatrix} d^{(e)}=f^{(e)}$$ [4X22X22X4] 4X14X1 [ equal to a 4X4 ]

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

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

Justification for the above equation: Principle Virtual Work (PVW)

(See notes from the September1 17th for the first application of PVW. This is shown through Reduction of Global FD relationship.)

$$k d = F \Leftrightarrow  \bar{k} \bar{d} = \bar{F}$$ 6X6 6X1   6X6    2X2 2X1   2X1

Question: Why would you not solve as follows? $$d = k^{-1} F$$ ?

It is because the matrix k is a singularity, and the det k=0. Therefore $$k^{-1}$$ is not invertible. This is because

to find $$k^{-1}$$ you need to compute $$\frac{1}{\det k}$$

why? For an unconstrained structural system, there are three possible rigid body motions in 2-D. (2-translational and 1-rotational)

Dynamic eigenvalue problems: $$K v=\lambda Mv $$

K=Stiffness matrix

$$\lambda$$=Eigen value (related to the vibrational frequency)

M=Mass matrix

Zero eigenvalue corresponds to zero stored elastic energy, so rigid body modes. Modes, also called mode shapes, are used to describe the vibrations in an object. These correspond to the number of half wave vibrations in a body.

September 29, 2008
The previous simplification that was applied to the Stiffness Matrix was done by application of the boundary conditions as follows:

$$\left[\begin{array}{cccccc} K_{11}&K_{12}&K_{13}&K_{14}&K_{15}&K_{16}\\ K_{21}&K_{22}&K_{23}&K_{24}&K_{25}&K_{26}\\ K_{31}&K_{32}&K_{33}&K_{34}&K_{35}&K_{36}\\ K_{41}&K_{42}&K_{43}&K_{44}&K_{45}&K_{46}\\ K_{51}&K_{52}&K_{53}&K_{54}&K_{55}&K_{56}\\ K_{61}&K_{62}&K_{63}&K_{64}&K_{65}&K_{66}\end{array}\right] \left[\begin{array}{c}0\\0\\d_3\\d_4\\0\\0\end{array}\right]= \left[\begin{array}{cc} K_{13}&K_{14}\\K_{23}&K_{24}\\K_{33}&K_{34}\\K_{43}&K_{44}\\K_{53}&K_{54}\\K_{63}&K_{64} \end{array}\right] \left[\begin{array}{c}d_3\\d_4\end{array}\right]= \left[\begin{array}{c}F_1\\F_2\\F_3\\F_4\\F_5\\F_6\end{array}\right]$$

F3 and F4 being the known applied loads, the computation for rows 3 and 4 is no longer needed for it would only return F3 and F4. The only needed computations would be for rows 1, 2, 5, 6 to get F1, F2, F5, and F6

To close the loop between Statics and Finite Element Methods, Virtual Displacement can be used as follows:

Two bar Truss syst:

By statics, the reactions forces are known leading to the member forces P1(1), P2(1), component axial displacement degrees of freedom:

$$\begin{array}{l} q_1^{(1)}=\dfrac{P_1^{(1)}}{k^{(1)}}=\dfrac{P_2^{(1)}}{k^{(1)}}=AC\\ q_1^{(1)}=0,\, fixed\, node\, (1)\\ q_1^{(2)}=-\dfrac{P_2^{(2)}}{k^{(2)}}\\ q_2^{(2)}=0,\, fixed\, node\, (3)\\ \end{array}$$

Such a situation makes us wonder about how it would possible to back out the displacement degrees of freedom of node (2) from the results above?

Uy = R sin(α) ≈ Rα (if α small)

Ux = R(1 - cos(α)) ≈ 0 (first Order)

October 1, 2008
Closing the Loop (Continued)

Infinitesimal displacement (Relative to virtual displacement)

$$AC = \frac{P_2^1}{K^1} = \frac{5.1243}{\frac{3}{4}} = 6.8324$$

$$ AB = \frac{P_1^2}{K^2}$$

The process for finding the coordinates of $$( X_B, Y_B)$$ and $$( X_C , Y_C)$$ is as follows:

$$X_C$$ = 6.8324 cosθ$$ ^1$$

$$Y_C$$ = 6.8324 sinθ$$ ^1$$

By adding the values obtained for $$(X_C, Y_C)$$ and $$(X_B , Y_B)$$ to the values for the coordinate of A, you will obtain the actual values for points C & B.

Once B & C are known you have two unkowns remaining $$(X_D, Y_D)$$. To get these unkowns the equations for lines AB and BC are needed.

$$\overrightarrow {PQ}\ = (PQ)\hat i\ = (PQ) [cos\theta \hat i\ + \sin\theta \hat j] = (x - x_p) \hat i\ + (y - y_p) \hat j$$

$$X - X_p = (PQ) \cos \theta$$ $$Y - Y_p = (PQ) \sin \theta$$

$$\frac {Y - Y_p}{X - X_p} = \tan \theta \Rrightarrow Y - Y_p = (\tan \theta)(X - X_p) $$

Equation for a line perpendicular to the above line passing P:

$$ Y - Y_p = \tan(\theta + \frac{\pi}{2})(X - X_p)$$

$$ \overrightarrow {AD}\ = (X_d - X_a)\hat i\ + (Y_d - Y_a)\hat j$$

Due to A being at the origin A(0,0). $$\Rrightarrow X_a = 0 ,  Y_a = 0$$

By definition: $$\overrightarrow {AD} = d_3 \hat i\ + d_4 \hat j$$

$$\overrightarrow AD$$ is the displacement vector of A

$$d_3, d_4$$ are acquired from FEM $$\Rrightarrow d_3 = X_d , d_4 = Y_d$$

3 Bar Truss System

Convenient Local Node Numbering

$$\underline {P} = 30$$

$$\begin{matrix}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 \\ \theta^1 = 30^\circ   &   \theta^1 = -30^\circ   &   \theta^3 = 45^\circ \end{matrix}$$

October 3, 2008
The following schematic shows 2 equations, but 3 unknowns.

For the summation of the moment about $$A$$, the following free body diagram is used:

$$\begin{align} & \sum{F_{x}}=0 \\ & \sum{F_{y}=0} \\ & \sum{M_{A}=0} \\ \end{align}$$

The summation of the moment about the point A provides $$0=0$$, which makes it a trivial answer.

The same goes for the $$\sum{M_{B}}$$, as it also reveals a solution of $$0=0$$. Because the perpendicular distances are equivalent, and all of the forces act on a common line of action, the sum of the moments about B will also equal zero.

 * Note: This goes for both 2D and 3D schematics.

For just the bar with the Reaction Force labeled R1, the moment about the point B will become:

$$\sum{M_{B}=\overrightarrow{BA}\times \overrightarrow{F}}$$

For $${A}'$$ on the line of action labeled $$\vec{F}$$, $$\sum{M_{B}=\overrightarrow{B{A}'}\times \overrightarrow{F}}$$.

Where, $$\overrightarrow{B{A}'}=\overrightarrow{BA}+\overrightarrow{A{A}'}$$, and $$\overrightarrow{M_{B}}=(\overrightarrow{BA}+\overrightarrow{A{A}'})\times \vec{F}=\overrightarrow{BA}\times \vec{F}+\overrightarrow{A{A}'}\times \vec{F}$$

 * Note: $$\overrightarrow{A{A}'}\times \vec{F}=\vec{0}$$

$$\Rightarrow $$ Back to the 3-bar truss...

For the truss, Node A is in equilibrium: $$\sum\limits_{i=0}^{3}{\vec{F}_{i}=\vec{0}}$$

$$\sum\limits_{i}{\overrightarrow{M_{B,i}}=\sum\limits_{i}{\overrightarrow{B{A}'_{i}}}\times \overrightarrow{F_{i}}}$$

Where, $$A_{i}^{\prime }$$ is equal to any point on the line of action of $$\overrightarrow{F_{i}}$$.

Therefore, by combining the knowledge of $$A_{i}^{\prime }$$ with the knoledge of $$\sum\limits_{i}{\overrightarrow{M_{B,i}}=\sum\limits_{i}{\overrightarrow{B{A}'_{i}}}\times \overrightarrow{F_{i}}}$$ provides: $$\sum\limits_{i}{\overrightarrow{M_{B,i}}=\sum\limits_{i}{BA\times \overrightarrow{F_{i}}=\overrightarrow{BA}\times \sum\limits_{i}{\overrightarrow{F_{i}}}}}$$

 * Note that $$\sum\limits_{i}{\overrightarrow{F_{i}}=\vec{0}}$$.

Thus, the schematic above is statically indeterminant, because it has 2 equations and 3 unknowns.

Using the above diagram of the global matrix $$K$$, which combines the 3 element stiffness matrices: $$k_{ij}^{(1)}$$, $$k_{ij}^{(2)}$$, $$k_{ij}^{(3)}$$, the entries into the global matrix will be as follows:

$$\underline{K}=$$

Contributing Team Members
Arunas Janulevicius --Eml4500.f08.jamama.jan 20:55, 7 October 2008 (UTC)

William Mueller -- Eml4500.f08.jamama.mueller 00:02, 8 October 2008 (UTC)

Cedric Adam --Eml4500.f08.jamama.adam 05:46, 8 October 2008 (UTC)

Justin McIntire--Eml4500.f08.jamama.justin 16:04, 8 October 2008 (UTC)

Megan Alvarez --EML4500.f08.jamama.megan 18:34, 8 October 2008 (UTC)

Chris Alford --Eml4500.f08.jamama.chris 19:40, 8 October 2008 (UTC)

title: User:EML4500.f08.JAMAMA/Homework3


text:
User:Eml4500.f08.team.figuerrez/hw3

Question: How about ΣM B ?

$$ BA' = BA + AA' $$

$$ M_{B} = (BA + AA') x F $$

$$ M_{B} = BA x F + AA' x F $$

$$ (AA'x F) $$ is equal to zero.

Therefore the moment at point B is equal to $$ M_{B} = BA x F $$

=3-bar truss system=

Node A is in equilibrium:

$$\sum F_{i} = 0$$

$$\sum M_{B,i} = \sum BA'_{i} * F_{i}$$

$$A'_{i}$$ equals to any point on line of action of $$F_{i}$$ therefore  $$\sum M_{B,i} = \sum BA_{i} * F_{i}$$ $$\sum M_{B,i} = BA * \sum F_{i} = 0$$

=Stiffness Matrix=

The corresponding stiffness matrix of the 3-bar truss system is:

$$ \begin{bmatrix} K_{11}^{(1)} & K_{12}^{(1)} & K_{13}^{(1)} & K_{14}^{(1)} & - & - & - & -\\

K_{21}^{(1)} & K_{22}^{(1)} & K_{23}^{(1)} & K_{24}^{(1)} & - & - & - & -\\

K_{31}^{(1)} & K_{32}^{(1)} & (K_{33}^{(1)} + K_{11}^{(2)} + K_{11}^{(3)}) & (K_{34}^{(1)} + K_{12}^{(2)} + K_{12}^{(3)}) & K_{13}^{(2)} & K_{14}^{(2)} & K_{13}^{(3)} & K_{14}^{(3)}\\

K_{41}^{(1)} & K_{42}^{(1)} & (K_{43}^{(1)} + K_{21}^{(2)} + K_{21}^{(3)}) & (K_{44}^{(1)} + K_{22}^{(2)} + K_{22}^{(3)}) & K_{23}^{(2)} & K_{24}^{(2)} & K_{23}^{(3)} & K_{24}^{(3)}\\

- & - & K_{31}^{(2)} & K_{32}^{(2)} & K_{33}^{(2)} & K_{34}^{(2)} & - & -\\

- & - & K_{41}^{(2)} & K_{42}^{(2)} & K_{43}^{(2)} & K_{44}^{(2)} & - & -\\

- & - & K_{31}^{(3)} & K_{32}^{(3)} & - & - & - & -\\

- & - & K_{41}^{(3)} & K_{42}^{(3)} & - & - & - & -\\

\end{bmatrix} $$

title: User:Eml4500.f08.team.figuerrez/hw3


text:
User:EML4500.f08.delta 6.ramirez/pointcoordinates

The following are the unknown variables that will be determined:

(XD, YD) (XB, YB) (XC, YC)

We will begin by finding the coordinates of point D. Note: Point A has been designated as the origin.

Using the definition of line AD, (XD, YD) can be found.

$$\vec{AD} = (X_{D} - X_{A}) \vec{i} + (Y_{D} - Y_{A}) \vec{j}$$

Since point A is referenced as the origin, XA and YA are equal to zero.

By definition of the displacement vector of A,

$$\vec{AD} = d_{3}\vec{i} + d_{4}\vec{j}$$

Thus,

(XD, YD) = (d3, d4) = (4.352, 6.1271)

Next, using the relationship of two points on the same line gives,

$$y - y _{C} = tan(\theta) (x - x_{C})$$

where θ = 30°, (XC, YC) = (5.92, 3.42)

The relationship between point B and point D gives,

$$y - y _{B} = tan(\theta + \frac{\pi }{2}) (x - x_{B})$$

where θ = 135°,(XB, YB) = (-0.88, 0.88)

Results:

Point D: (4.35, 6.13) Point C: (5.92, 3.42) Point B: (-0.88, 0.88)


text:
User:Eml4500.f08.bottle.barnes/HW3

Using the global FD relation, where displacements $$\displaystyle d_1, d_2, d_5, d_6 $$ all equal zero, the first two columns and the last two columns for the $$\displaystyle K $$ matrix can be ignored, yielding the following relation:

 * $$ \begin{bmatrix} K_{13} & K_{14}

\\ K_{23} & K_{24} \\ K_{33} & K_{34} \\ K_{43} & K_{44} \\ K_{53} & K_{54} \\ K_{63} & K_{64} \end{bmatrix} \begin{Bmatrix} d_3 \\d_4 \\ \end{Bmatrix} = \begin{Bmatrix} F_1 \\F_2 \\ F_3 \\ F_4 \\ F_5 \\ F_6  \end{Bmatrix}$$

which can also be expressed as

 * $$ \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} \begin{Bmatrix} d_3 \\d_4 \\ \end{Bmatrix} = \begin{Bmatrix} F_1 \\F_2 \\ F_3 \\ F_4 \\ F_5 \\ F_6 \end{Bmatrix} $$

Inputting the values for each $$\displaystyle K $$ and $$\displaystyle d $$ yields

 * $$ \begin{bmatrix} -{9 \over 16} & -{3 \sqrt 3 \over 16}

\\ -{3 \sqrt 3 \over 16} & -{3 \over 16} \\ {9 \over 16} + {5 \over 2} & {3 \sqrt 3 \over 16} + -{5 \over 2} \\ {3 \sqrt 3 \over 16} + -{5 \over 2} & {3 \sqrt 3 \over 16} + {5 \over 2} \\ -{5 \over 2} & {5 \over 2} \\ {5 \over 2} & -{5 \over 2} \end{bmatrix} \begin{Bmatrix} 4.352 \\ 6.1271 \\ \end{Bmatrix} = \begin{Bmatrix} F_1 \\F_2 \\ 0 \\ 7 \\ F_5 \\ F_6 \end{Bmatrix} $$

Note that since $$\displaystyle F_3 $$ and $$\displaystyle F_4 $$ are known (because they are the applied loads, 0 and 7 respectively) rows three and four do not need to be computed.

Thus, the computed solution is

 * $$ F = \begin{Bmatrix} -4.4378 \\ -2.5622 \\ 0 \\ 7 \\ 4.4378 \\ -4.4378 \end{Bmatrix} $$

Notice that these results are the same as those determined from Method 1.

title: User:Eml4500.f08.bottle.barnes/HW3


text: File:HW3 Forces.JPG

title: File:HW3 Forces.JPG


text: File:Pqpicture.JPG

title: File:Pqpicture.JPG


text: File:HW3 B.JPG

title: File:HW3 B.JPG


text: File:Plot for truss.GIF

title: File:Plot for truss.GIF


text: File:Compare Beam2.png

title: File:Compare Beam2.png


title: File:Picture3.png


text:
User:Eml4500.f08.wiki1.aguilar/TeamHW3

Derivation of Elemental Free-body Diagram with respect to the Global Coordinate System

 * $$ \textbf{k}_{4x4}^{(e)}\textbf{d}_{4x1}^{(e)}=\textbf{f}_{4x1}^{(e)}

$$

$$ {k}^{(e)} \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix} \begin{Bmatrix} {q}_{1}^{(e)}\\ {q}_{2}^{(e)} \end{Bmatrix} = \begin{Bmatrix} {P}_{1}^{(e)}\\ {P}_{2}^{(e)} \end{Bmatrix} $$

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

We want to find the relationship between $$\textbf{q}^{(e)}$$ and $$\textbf{d}^{(e)}$$, and $$\textbf{P}^{(e)}$$ and $$\textbf{F}^{(e)}$$.

These relationships can be expressed in the form: $$ \textbf{q}_{2x1}^{(e)}\textbf{T}_{2x4}^{(e)}=\textbf{d}_{4x1}^{(e)} $$

Consider the displacement vector of the local node 1 denoted by $$\vec{d}_{1}^{(e)}$$.

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

$${q}_{1}^{(e)}$$ = axial displacement of node 1 is the prthagonal projection of the displacement vector $$\vec{d}_{1}^{(e)}$$ of node 1 on the $$\tilde{x}$$ axis of element 'e'.

$${q}_{1}^{(e)}=\vec{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}\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{i}\cdot\vec{\tilde{i}})=\sin\theta^{(e)}=m^{(e)}$$

$${q}_{1}^{(e)}=l^{(e)}{d}_{1}^{(e)}+m^{(e)}{d}_{2}^{(e)}$$

$${q}_{1}^{(e)}= \begin{bmatrix}l^{(e)} & m^{(e)}\end{bmatrix}_{1x2}\begin{Bmatrix}{d}_{1}^{(e)}\\{d}_{2}^{(e)}  \end{Bmatrix}_{2x1}$$

Here we can see that $${q}_{1}^{(e)}$$ is a 1x1 scalar.

We do this for node 2 as well. $${q}_{2}^{(e)}= \begin{bmatrix}l^{(e)} & m^{(e)}\end{bmatrix}_{1x2}\begin{Bmatrix}{d}_{3}^{(e)}\\{d}_{4}^{(e)}  \end{Bmatrix}_{2x1}$$

which leads us to:

$$\begin{Bmatrix} {q}_{1}^{(e)}\\{q}_{2}^{(e)}\end{Bmatrix}= \begin{bmatrix} l^{(e)}& m^{(e)}&0&0\\0&0& l^{(e)}& m^{(e)} \end{bmatrix} \begin{Bmatrix}{d}_{1}^{(e)}\\{d}_{2}^{(e)}\\{d}_{3}^{(e)}\\{d}_{4}^{(e)} \end{Bmatrix}$$

where this is the equation we set out to derive, $$ \textbf{q}_{2x1}^{(e)}\textbf{T}_{2x4}^{(e)}=\textbf{d}_{4x1}^{(e)} $$

Similarly, (same argument): $$ \begin{Bmatrix} P_1^{(e)} \\ P_2^{(e)} \end{Bmatrix} = I^{(e)} \begin{Bmatrix} f_1^{(e)} \\ f_2^{(e)} \\ f_3^{(e)} \\ f_4^{(e)}\end{Bmatrix} $$ b Where P is a 2x1 matrix, I is a 2x4 matrix, and f is a 4x1 matrix.

This relationship is the same as saying $$ \underline{P}^{(e)}= I^{(e)} \underline{f}^{(e)} $$ Recall the axial Fd relationship:

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

where k is a 2x2 matrix, q is a 2x1 matrix, and P is a 2x1 matrix.

$$ \hat k^{(e)} \underbrace{(\underline{T}^{(e)} \underline{d}^{(e)} )} = \underbrace{(\underline{T}^{(e)} \underline{f}^{(e)})} $$

$$ - \qquad \underline{q}^{(e)} \qquad \qquad \underline{P}^{(e)} $$

Goal: We want to have k (e) d (e)= f (e) so "move" T e from right side to the left side by pre multiplying equation by T (e)-1, the inverse of T (e). Unfortunately, T(e) is a rectangular matrix of the size 2x4 and cannot be inverted. Only square matricies can be inverted. To solve this issue, transpose T. Ans:

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

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

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

where k is on p6.1, k hat is on 12-2, and T is on 12-5. Justification for (1): PVW see 10-1 for 1st applications of PVW, reduction of global FD relationship.

$$ \underline {K}_{6x6} d_{6x1} = \underline{F}_{6x1} \rightarrow \underline{\bar K}_{2x2} \underline{\bar d}_{2x1} = \underline {\bar F}_{2x1} $$

Remember: Why not solve as follows?

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

Used mathcad to try to solve K-1 and could not, due to singularity of K ., i.e. the determinant of K =0 and thus K is not invertible. Recall that you need to computer $$\begin{matrix} \frac{1}{detK} \end{matrix}$$ to find K -1. Why? For an unconstrained structure system, there are 3 possible rigid body mostion in 2-D (2 translational and 1 rotational).

HW: Find the eigenvalues of K and make observations about the number of eigenvalues. Dyanmic eigenvalue problem K v = λ M  v  Where K is the stiffness matrix, M is the mass. Zero evaluation corresponds to zero stored elastic energy which means rigid body motion.

Using the Global Free-body Diagram Relationship
Using d3, and d4 in $$ \textbf{K}\textbf{d}=\textbf{F} $$ and using boundary conditions to eliminate rows in the global displacement matrix we get:

$$ \begin{bmatrix} K_{1,3} & K_{1,4} \\ K_{2,3} & K_{2,4} \\ K_{3,3} & K_{3,4} \\ K_{4,3} & K_{4,4} \\ K_{5,3} & K_{5,4} \\ K_{6,3} & K_{6,4} \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}$$

Note: We really only need to do the computations for rows 1,2,5, and 6 to get F1,F2,F5,and F6 because computation from rows 3 and 4 gives the applied load which is already known.

Closing the Loop between FEM and Statics: Virtual Displacement
Two-bar truss system:

Since our two-bar truss system is statically determinant because we can view element 1 and 2 as two force bodies. By statics, the reactions are known and therefore, so are the member forces $$P^{(1)}_1$$, and $$P^{(1)}_1$$.

We then compute axial displacement degrees of freedom. This is the amount of extension within the bars.

$${q}_{2}^{(1)} = \frac{{P}_{1}^{(1)}}{{k}^{(1)}}= \frac{{P}_{2}^{(1)}}{{k}^{(1)}}=AC$$

$${q}_{1}^{(1)} = 0$$    (Node 1 is fixed)

$${q}_{2}^{(1)} = \frac{{-P}_{2}^{(2)}}{{k}^{(2)}}=AB$$

$${q}_{2}^{(2)} = 0$$    (Node 2 is fixed)

How do we back out the displacement DOFs of node 2 from above results?

Note: The displacement of node 2 is in the direction of vector D.

Infinitesimal displacement (related to virtual displacement)
Global Free-Body Diagram showing Infinitesimal Displacements

$$ \overline{AC}= {|P_2^{(1)}| \over k^{(1)}}={5.1243 \over 3/4}=6.8324 $$

$$ \overline{AB}= {|P_1^{(1)}| \over k^{(2)}}={6.2762 \over 5}=1.2552 $$

Find (x,y) coordinates of B and C:

Point B:

$$ x= \overline{AB} cos{135^\circ}=1.2552*-{\sqrt 2 \over 2}=-0.8876 $$

$$ y= \overline{AB} sin{135^\circ}=1.2552*{\sqrt 2 \over 2}=0.8876 $$

Point C:

$$ x= \overline{AC} cos{30^\circ}=6.8324*{\sqrt 3 \over 2}=5.917 $$

$$ y= \overline{AC} sin{30^\circ}=6.8324*{1 \over 2}=3.4162 $$

We now have two unknowns, (XD,YD)

We need the equations for line $$\overline{AB}$$ and $$\overline{BC}$$

Method for the Determination of the Slope between Two Arbitrary Points
PQ Line Segment Diagram

$$ \overline{PQ}=(PQ) \tilde{i}=(PQ)[cos \theta \vec i +sin \theta \vec j]=(x-x_P)\vec i +(y-y_P)\vec j $$

$$ x-x_P=(PQ)cos \theta$$

$$ y-y_P=(PQ)sin \theta$$

$$ {y-y_P \over x-x_P}=tan \theta$$

$$ y-y_P=tan \theta(x-x_P)$$

Equation for line perpendicular to PQ passing through P:

$$ y-y_P=tan( \theta + {\pi \over 2})(x-x_P)$$

Determination of (XD,YD)
Line Perpendicular to Point B:

$$y-0.8876=tan{225^\circ}(x+0.8876)$$

$$y=x+1.7752$$

Line Perpendicular to Point C:

$$y-3.4162=tan{120^\circ}(x-5.917)$$

$$y=-1.732x + 13.665$$

Intersection at Point (XD,YD):

$$y_D+1.732(y_D-1.7752)=13.665$$

$$y_D=6.127$$

$$(x_D+1.7752)=-1.732x_D + 13.665$$

$$x_D=4.352$$

$$\vec {AD} = (x_D-x_A) \vec i + (y_D-y_A) \vec j$$

Simplification: $$x_A=0, y_A=0$$

$$\vec {AD} = x_D \vec i + y_D \vec j$$

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

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

$$ d_3=4.352$$

$$d_4=6.127$$

3-bar Truss System
3-bar Truss System

$$ E^{(1)}=2, A^{(1)}=3, L^{(1)}=5$$

$$E^{(2)}=4, A^{(2)}=1, L^{(2)}=5$$

$$E^{(3)}=3, A^{(3)}=2, L^{(3)}=10$$

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

$$P=30$$

Convenient Local Node Numbering:

3-bar Truss System Local Node Numbering

3-Bar Truss System

ΣFx = 0

ΣFy = 0

ΣMA = 0 (Trivial)

→ 2 equations, 3 unknowns → Statically Indeterminate

Question: How about MB? (3-D Explanation)

$$ \overline{M_B} = \overline{BA}\times\overline{F} = \overline{BA'} \times\overline{F} $$

(for A' on line of action of \overline{F}

$$ \overline{M_B} = (\overline{BA} + \overline{AA'}) \times\overline{F} = \overline{BA'} \times\overline{F} + \overline{AA'} \times\overline{F} $$

Back to 3-bar truss:

Node A is in equilibrium:

$$\sum_{i=0}^3 {\overline{F_i}} = \overline{0}$$

$$\sum_{i} {\overline{M_B,i}} = \sum_{i} {\overline{BA'_i}} \times \overline{F_i} $$

Ai' = any point on line of action of Fi

$$\sum_{i} {\overline{M_B,i}} = \sum_{i} {\overline{BA'_i}} \times \overline{F_i} = \overline{BA} \times \sum_{i=0} {\overline{F_i}} = \overline{0}$$

\begin{bmatrix} K_{11} & K_{12} & K_{13} & K_{14} & K_{15} & K_{16} & K_{17} & K_{18}\\ K_{21} & K_{22} & K_{23} & K_{24} & K_{25} & K_{26} & K_{27} & K_{28}\\ K_{31} & K_{32} & K_{33} & K_{34} & K_{35} & K_{36} & K_{37} & K_{38}\\ K_{41} & K_{42} & K_{43} & K_{44} & K_{45} & K_{46} & K_{47} & K_{48}\\ K_{51} & K_{52} & K_{53} & K_{54} & K_{55} & K_{56} & K_{57} & K_{58}\\ K_{61} & K_{62} & K_{63} & K_{64} & K_{65} & K_{66} & K_{67} & K_{68}\\ K_{71} & K_{72} & K_{73} & K_{74} & K_{75} & K_{76} & K_{77} & K_{78}\\ K_{81} & K_{82} & K_{83} & K_{84} & K_{85} & K_{86} & K_{87} & K_{88} \end{bmatrix}

$$

= \begin{bmatrix} k_{11}^{(1)} & k_{12}^{(1)} & k_{13}^{(1)} & k_{14}^{(1)} & 0 & 0 & 0 & 0\\ k_{21}^{(1)} & k_{22}^{(1)} & k_{23}^{(1)} & k_{24}^{(1)} & 0 & 0 & 0 & 0\\ k_{31}^{(1)} & k_{32}^{(1)} & (k_{33}^{(1)}+k_{11}^{(2)}+k_{11}^{(3)}) & (k_{34}^{(1)}+k_{12}^{(2)}+k_{12}^{(3)}) & k_{13}^{(2)} & k_{14}^{(2)} & k_{15}^{(3)} & k_{16}^{(3)}\\ k_{41}^{(1)} & k_{42}^{(1)} & (k_{43}^{(1)}+k_{21}^{(2)}+k_{21}^{(3)}) & (k_{44}^{(1)}+k_{22}^{(2)}+k_{22}^{(3)}) & k_{23}^{(2)} & k_{24}^{(2)} & k_{25}^{(3)} & k_{26}^{(3)}\\ 0 & 0 & k_{31}^{(2)} & k_{32}^{(2)} & k_{33}^{(2)} & k_{34}^{(2)} & 0 & 0\\ 0 & 0 & k_{41}^{(2)} & k_{42}^{(2)} & k_{43}^{(2)} & k_{44}^{(2)} & 0 & 0\\ 0 & 0 & k_{51}^{(3)} & k_{52}^{(3)} & 0 & 0 & k_{55}^{(3)} & k_{56}^{(3)}\\ 0 & 0 & k_{61}^{(3)} & k_{62}^{(3)} & 0 & 0 & k_{65}^{(3)} & k_{66}^{(3)} \end{bmatrix}

$$

MatLab Homework
Part of this assignment required the team to develop a MATLAB code that plots the deformed and un-deformed configurations of the two-bar truss system solved in class. The code is based on two models that were created by Dr. Vu-Quoc and X.G. Tan.

The blue dashed line represents the un-deformed truss, while the solid red line shows the truss after it has been deformed.

Authors
Eml4500.f08.wiki1.oatley 20:02, 8 October 2008 (UTC)

Eml4500.f08.wiki1.brannon 01:31, 8 October 2008 (UTC)

Eml4500.f08.wiki1.ambrosio 21:27, 7 October 2008 (UTC)

Eml4500.f08.wiki1.aguilar 20:45, 8 October 2008 (UTC)

Eml4500.f08.wiki1.schaet 19:42, 8 October 2008 (UTC)

Eml4500.f08.wiki1.handy 20:25, 8 October 2008 (UTC)

title: User:Eml4500.f08.wiki1.aguilar/TeamHW3


title: File:HW3b.JPG


title: File:Diagram16-1.1.png


text: File:Hw3 matlab.jpg

title: File:Hw3 matlab.jpg


title: File:3bartruss.JPG


text: File:Rotation of the cross-section of a bar2.jpg


text:
User:Eas4200c.f08.carbon.guidry

Beau Guidry

title: User:Eas4200c.f08.carbon.guidry


title: File:Diagram16-1.3.png


title: File:3bartrusslocalnodes.JPG


title: File:Diagram16-1.4.png


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


text:
User:Eas4200c.f08.carbon.orear/HW3/Contributors

Brian Taylor Eas4200c.f08.carbon.taylor 22:28, 7 October 2008 (UTC)

Jeff Wright Eas4200c.f08.carbon.w 04:31, 8 October 2008 (UTC)

Aaron Clausen Eas4200c.f08.carbon.clausen 06:53, 8 October 2008 (UTC)

Abby Booth Eas4200c.f08.carbon.booth 14:12, 8 October 2008 (UTC)

Jeff O'Rear Eas4200c.f08.carbon.orear 15:53, 8 October 2008 (UTC)

Beau Guidry Eas4200c.f08.carbon.guidry 15:56, 8 October 2008 (UTC)


text:
User:Eml4500.f08.delta 6.guzman/MATLAB Code (plot of undeformed two-bar truss structure)

MATLAB Plot of a Two-bar Truss Structure (undeformed and deformed)
The following code plots the analysis of an undeformed and deformed two-bar truss system using MATLAB MATLAB Code

Note: this MATLAB code was created using the principles obtained from:

 * The MATLAB script for the analysis of the Two Bar Truss System obtained from Dr. L. Vu-Quoc website http://clesm.mae.ufl.edu/~vql/courses/fead/2008.fall/codes/twoBarTrussEx.txt

 * The MATLAB script for the plot of a crab-leg structure created by X.G. Tan (2003) and obtained from Dr. L. Vu-Quoc website http://apollo.mae.ufl.edu/~vql/courses/fead/2006.fall/codes/matlab_plot_example2.m

 * The MATLAB script for the plot of an electric pylon (undeformed five-bar truss structure) created by Team MeshWorks (2003) and obtained from Dr. L. Vu-Quoc website http://apollo.mae.ufl.edu/~vql/courses/fead/2006.fall/code.web.pages/electric_pylon.m

Results


text: User:Cie

title: User:Cie


title: File:Meeting15loop.JPG


text:
User talk:Cie

 Hello Cie, and welcome to Wikiversity! If you need help, feel free to visit my talk page, or contact us and ask questions. After you leave a comment on a talk page, remember to sign and date; it helps everyone follow the threads of the discussion. The signature icon in the edit window makes it simple. To get started, you may

 * Take a guided tour and learn to edit.
 * Visit a (kind of) random random project.
 * Browse everything on Wikiversity, or visit a portal corresponding to your educational level: pre-school, primary, secondary, tertiary, non-formal education.
 * Find out about research activities on Wikiversity.

 * Read an introduction for teachers and find out how to write an educational resource for Wikiversity.
 * Give feedback about your initial observations
 * Discuss Wikiversity issues or ask questions at the colloquium.
 * Chat with other Wikiversitans on #wikiversity-en.

And don't forget to explore Wikiversity with the links to your left. Be bold to contribute and to experiment with the sandbox or your userpage, and see you around Wikiversity! Jtneill - Talk - c 22:43, 7 October 2008 (UTC)

 Hi! My name is WikiversityJack. I am a pumpkin and a Wikiversity mascot. My aunt told me that fun + learning = Wikiversity. Let me know if I can help you out.

title: User talk:Cie


text: File:Matlab graph.JPG

title: File:Matlab graph.JPG