User:Eml4500.f08.qwiki.bishop/HW3

Comment:
The pictures that were previously used in this homework submission have been changed. The previous images (HW_3.GIF and HW_3.2.GIF) were based on images that were not originally authored by Team Qwiki, however I did not intentionally plagiarize Team Echo's work.

As Lee Harris, fellow Team Qwiki member, also stated in Homework 4, while working on Homework 2, both Team Qwiki and Team Echo created two images, with the same name: Element1.jpg and Element2.jpg. Due to the nature of MediaWiki, if two different teams create an image with the same name, only the latest version is displayed when called upon in an article. Therefore, even though Team Qwiki uploaded two of our own images (Element1.jpg and Element2.jpg) during homework 2, Team Echo's images were displayed because their upload was the most recent. In the future, we will make sure to upload files that have a unique name.

For Homework 3 , I created HW_3.GIF and HW_3.2.GIF on October 5 because it was the two bar truss used in Homework 2, and to keep the images showing the same two bar truss with identical angles and lengths, I based the mentioned images on them. I was under the impression that Element1.jpg was the work of Team Qwiki. The problem has been resolved and Team Qwiki files are now displayed here.

In the future Team Qwiki and myself will be sure to:

• Create unique names when uploading files
• Check the file history of all images before I right click and save to my computer

A comparison between the current version of Homework 3 and our prior submission can be found here.
Andrea Booher EML4500.f08.qwiki.booher 00:17, 13 November 2008 (UTC)

Lecture 12 - Monday, September 22, 2008

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

Recall from Page 6-1, $k^{(e)}d^{(e)}=f^{(e)}$ (Equation 1)

Therefore, the element force displacement matrix can be written as follows:

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

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

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

The overall goal is to derive equation 1 from equation 2 (already derived in Meeting 4). We want to find the relationship between:

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

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

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

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

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

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

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

HW From September 22

HW Problem: Find $q_2^{(e)}$ for node [2] using the same procedure above that was used to find $q_2^{(e)}$. $q_{2}^{(e)}=d_{[3]}^{(e)}*\vec{\tilde{i}}$ $q_{1}^{(e)}=(d_{3}^{(e)}\vec{i}+d_{4}^{(e)}\vec{j})\vec{\tilde{i}}$ $q_{1}^{(e)}=d_{3}^{(e)}(\vec{i}*\vec{\tilde{i}})+d_{4}^{(e)}(\vec{j}*\vec{\tilde{i}})$ $q_{1}^{(e)}=l^{(e)}d_{3}^{(e)}+m^{(e)}d_{4}^{(e)}$ Answer: $q_{1}^{(e)}=\begin{bmatrix} l^{(e)} & m^{(e)}\\ \end{bmatrix}\begin{pmatrix} d_{3}^{(e)}\\ d_{4}^{(e)}\\ \end{pmatrix}$

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

$\begin{pmatrix} q_{1}^{(e)}\\ q_{2}^{(e)}\\ \end{pmatrix}=\begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)}\\ \end{bmatrix}\begin{pmatrix} d_{1}^{(e)}\\ d_{2}^{(e)}\\ d_{3}^{(e)}\\ d_{4}^{(e)}\\ \end{pmatrix}$

Lecture 13 - Wednesday, September 24, 2008

Test day, no lecture notes given. Exam 1 Solutions

Lecture 14 - Friday, September 26, 2008

Similar to our objective in Meeting 12, today we would like to transform from a 4x2 matrix to a 2x2. The transformation relationship

$\begin{bmatrix} P_1^{(e)}\\ P_2^{(e)}\\ \end{bmatrix} =T_{(2x4)}^{(e)} \begin{bmatrix} f_1^{(e)}\\ f_2^{(e)}\\ f_3^{(e)}\\ f_4^{(e)}\\ \end{bmatrix}$ derived in Meeting 12 can also be expressed in a more general form:

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

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

$\mathbf{T}^{(e)}_{(2x4)}$ is the matrix used to transform force components to an axial force

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

The overall goal is to find the relationship between:

• $\mathbf{q}^{(e)}_{(2x2)}$ and $\mathbf{d}^{(e)}_{(4x2)}$, (axial displacements of each element node and their respective displacement component)
• $\mathbf{P}^{(e)}_{(2x2)}$ and $\mathbf{f}^{(e)}_{(4x2)}$, (axial forces at each element node and their respective force components)

Let us consider the force transformation relationship:

$\mathbf{P}^{(e)}_{(2x2)}=\mathbf{T}^{(e)}_{(2x4)}\mathbf{f}^{(e)}_{(4x2)}$

Recall that the element axial FD relationship is given by:

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

By introducing a transformation matrix into the FD relationship, we obtain the following relationship:

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

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

What is an invertible matrix?

A given Matrix A is said to be invertible, or non-singular, if and only if there exists a matrix B such that $\mathbf{AB} = \mathbf{BA} = \mathbf{I}_n \$ where $\mathbf{I_n}$ is the nxn identity matrix. If this relationship between A and B holds true, then B is determined from A and is said to be the inverse of A. The inverse of A is denoted as A-1, and $\mathbf{AA^{-1}}=\mathbf{I}$. NOTE: Non-square matrices do not have an inverse.

Since we cannot solve for the inverse of the transformation ($\mathbf{T}^{(e)}$ is not an nxn matrix), we must instead make use of the transpose:

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

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

Homework: Verify that $\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$

Let us begin with the given expression: $\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$, where $\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}$ $\mathbf{T}^{(e)T}= \begin{bmatrix} l^{(e)} & 0\\ m^{(e)} & 0 \\ 0 & l^{(e)}\\ 0 & m^{(e)}\\ \end{bmatrix}$ $\mathbf{\hat{k}}^{(e)}_{(2x2)}=k^{(e)}\begin{bmatrix} 1 & -1 \\ -1 & 1\\ \end{bmatrix}$ $\mathbf{T}^{(e)}= \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)}\\ \end{bmatrix}$ Since we are dealing with matrices, we can first carry out the operation for the first two matrices $\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}$ and multply the resultant matrix by $T^{(e)}$ $\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$ $\mathbf{k}^{(e)}=(\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)})\mathbf{T}^{(e)}$ $\mathbf{k}^{(e)}=\Biggl( \begin{bmatrix} l^{(e)} & 0\\ m^{(e)} & 0 \\ 0 & l^{(e)}\\ 0 & m^{(e)}\\ \end{bmatrix} k^{(e)}\begin{bmatrix} 1 & -1 \\ -1 & 1\\ \end{bmatrix} \Biggr)\mathbf{T}^{(e)}$ $\mathbf{k}^{(e)}=k^{(e)}\Biggl( \begin{bmatrix} l^{(e)} & -l^{(e)}\\ m^{(e)} & -m^{(e)} \\ -l^{(e)} & l^{(e)}\\ -m^{(e)} & m^{(e)}\\ \end{bmatrix} \Biggr) \mathbf{T}^{(e)}$ $\mathbf{k}^{(e)}=k^{(e)} \begin{bmatrix} l^{(e)} & -l^{(e)}\\ m^{(e)} & -m^{(e)} \\ -l^{(e)} & l^{(e)}\\ -m^{(e)} & m^{(e)}\\ \end{bmatrix} \begin{bmatrix} l^{(e)} & m^{(e)} & 0 & 0\\ 0 & 0 & l^{(e)} & m^{(e)}\\ \end{bmatrix}$ $\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}$ Therefore, it holds true that $\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$

The relationship $\mathbf{k}^{(e)}=\mathbf{T}^{(e)T}\mathbf{\hat{k}}^{(e)}\mathbf{T}^{(e)}$ can be justified by the principle of virtual work, which was first introduced in Meeting 10. Additionally, applying the principle of virtual work to the relationship $K_{(6x6)}d_{(6x1)}=F_{(6x1)}$ reduces it to $K_{(2x2)}d_{(2x1)}=F_{(2x1)}$.

One may wonder why we cannot simply solve $d=K^{(-1)}F$
Answer: The equation cannot be solved in the current form due to the singularity of $K$, i.e, $detK=0$ and thus $K$ is not invertible. (See above for definition of interbile matrix) Recall that in order to find $K^{(-1)}$ we must first compute $1/{(detK)}$. If a matrix is singular then $1/{(detK)}$ is undefined and cannot be solved.
For our case of an unconstrained structural system, there are three possible rigid body motions in 2-dimensions (2 translations and 1 rotation) for each node.

HOMEWORK: Compute the eigenvalues of $K_{(6x6)}$

The 6x6 global stiffness matris is given as: $\mathbf{K}_{(6x6)}=\begin{bmatrix} 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\\ \end{bmatrix}$ The best way to solve for the eigenvalues of a 6x6 matrix is by using MATLAB. For the global stiffness, the input and output MATLAB code are provided below: >>clear all; >>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]; >> eig(K) ans = -0.0000 -0.0000 0.0000 0.0000 1.4706 10.0294 It can be seen that the two of the six eigenvalues are non-zero, with values of 1.47 and 10.02. It appears that the number of non-zero eigenvalues indicates that there are only two degrees of freedom in the system that can absorb energy. They are attributed to the 2D motion (x- and y-plane) motion of global node 2. If we solved for three non-zero eigenvalues, one could infer that global node 2 also stores rotational energy. A system with zero eigenvalues will not allow any energy storage in the specifically analyzed mode. Zero eigenvalues corresponds to no stored elastic energy, which implies rigid body motion. In general, the number of vibrational mode shapes of a given system is directly related to its eigenvalues and eigenvectors. Note that dynamic equilibrium problems are written at $K \nu = \lambda M \nu$

References

Lecture 15 - Monday, September 29, 2008

After attacking the basics of Finite Element Analysis, and reviewing Statics, it is time now to close the gap between the two subjects.

Closing the Loop

To do so, one must look to find the virtual displacement. Closing The Loop

The problem then arises of how to check the displacements. Using the below free body diagrams of the 2 force bodies (elements 1 and 2).

By Statics, the reactions are known and therefore the forces P₁⁽¹⁾, P₂⁽²⁾. Compute the axial displacement degrees of freedom, q (the amount of extension of the bars).

Computing the Displacements

Displacement Equations
$q_{2}^{(1)}=\frac{P_{1}^{(1)}}{k^{(1)}}=\frac{P_{2}^{(1)}}{k^{(1)}}$ $q_{1}^{(2)}=\frac{-P_{2}^{(2)}}{k^{(2)}}$ Knowing that q1(1)=0 because it is a fixed node (1). Also q2(2)=0 because it is a fixed node (3).

Question: How to back out from above results the displacement dof's of node (2)?

Lecture 16 - Wednesday, October 1, 2008

Closing the Loop: continued

• The infinitesimal displacement of point global node 2 (point A to point D) is related to the concept of virtual displacement from statics.
• The figure below shows the displacement of global node 1:
  • The displacement in the y direction, y_D, is equivalent to the displacement, d₄
  • The displacement in the x direction, x_D, is equivalent to the displacement, d₃
  • The vector summation of these two displacement provide the displacement vector $\vec{AD}$

• To determine the axial forces, P, we must first find the axial displacements, q
  • $q^{(1)}_{2}=l*d^{(1)}_{3}+m*d^{(1)}_{4}=cos(30)*4.352+sin(30)*6.1871=6.8625$
  • $q^{(2)}_{1}=l*d^{(2)}_{1}+m*d^{(2)}_{2}=cos(-45)*4.352+sin(-45)*6.1871=-1.2976$
• Then the axial forces
  • $P^{(1)}_{2}=k^{(1)}*q^{(1)}_{2}=(0.75)*(6.8625)=5.145$
  • $P^{(2)}_{1}=k^{(2)}*q^{(2)}_{1}=(5)*(-1.2976)=-6.488$

• The length of line AC can be found by:
  • $AC=\frac{\left|P^{(1)}_2 \right|}{k^{(1)}}=\frac{5.145}{.75}=6.8625$
• The length of line AB can be found by:
  • $AB=\frac{\left|P^{(2)}_1 \right|}{k^{(2)}}=\frac{-6.488}{5}=-1.2976$
    • Note 1: The negative value the length of AB is due to the assigned local axes of element 2.
    • Note 2: The lengths of AB and AC match the displacements for Elements 2 and 1 respectively.
• So the x and y coordinates for points B and C are (if point A is the origin):
  • $\left(x_{B},y_{B} \right)=(-.92,.92)$
  • $\left(x_{C},y_{C} \right)=(5.94,3.43)$

• So, we still have the 2 unknowns (x_D, y_D)
  • The figure below shows the xy axes and also the $\tilde{x}\tilde{y}$ axes
  • Let's analyze the vector PQ

• • The equation derived from this figure is: $y-y_{P}=tan\theta (x-x_{P})$
  • And the equation for the line perpendicular to this line is: $y-y_{P}=tan(\theta +\frac{\Pi }{2}) (x-x_{P})$
  • So, by applying some trigonometry and algebra $(x_D, y_D)=(4.352, 6.187)$
  • And the vector AD by definition is:
    • $\vec{AD}=d_{3}\vec{i}+d_{4}\vec{j}$

3-Bar Truss System

Properties of Each Element

Element	E	A	L	Theta (degrees)
1	2	3	5	30
2	4	1	5	-30
3	3	2	10	-135

• The force, P, has a value of 50
• Make sure the local node numbering is selected so that it is convenient to assemble the stiffness matrix, K
• The image below shows each Element and its local nodes as well as the global nodes

Lecture 17 - Friday, October 3, 2008

If we sum the forces on the X and Y axis and take the moment about A we only get two useful equations with three unknowns because the moment about A gives 0=0 which is useless.

SO to find another equation and solve the system, we ask ourselves, what about the moment about point B?

3-D Explanation

$\sum M_b = \overrightarrow {BA}x \overrightarrow{F} = \overrightarrow {BA'} x \overrightarrow{F}$

This is true for all A' on F's line of action

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

$\sum M_b = (\overrightarrow {BA} + \overrightarrow{AA'})x \overrightarrow{F} = \overrightarrow {BA} x \overrightarrow{F} + \overrightarrow{AA'} x \overrightarrow{F}$

$\overrightarrow{AA'}x \overrightarrow{F} = \overrightarrow{0}$ because they are colinear

Back to the 3-Bar Truss

Node A is in equilibrium

So...

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

$\sum_{i} \overrightarrow{M_{b,i}} = \sum_{i} \overrightarrow{BA'_i}x\overrightarrow{F_i}$

where A' is a point on the line of action of $\overrightarrow{F_i}$

So...

$\sum_{i} \overrightarrow{M_{b,i}} = \overrightarrow{0}$

We find that the moment about B is also zero and we are left with still only two equations and three unknowns.

The problem is thus statically 'indeterminate' and we must use the Finite Element Method.

So here is our resulting stiffness matrix

$\mathbf{K}= \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_{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)}\\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_{31}^{(3)} & k_{32}^{(3)} & 0 & 0 & k_{33}^{(3)} & k_{34}^{(3)}\\0 & 0 & k_{41}^{(3)} & k_{42}^{(3)} & 0 & 0 & k_{43}^{(3)} & k_{44}^{(3)} \end{bmatrix}$

and

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

HW 3 Matlab Assignment

The coding was developed using matlab_plot_example2.m. This code was created by X.G. Tan, and can be found on Dr. Loc Vu-Quoc's website.

 % Matlab Script
 %*********************************************************************
 % filename: matlab_plot_example2.m
 %
 % PURPOSE:
 %    Example of plotting with Matlab
 %    Input nodes and connectivity information
 %    Plot the crab-leg structure
 %
 % AUTHOR: X.G. Tan
 %
 % Modified on: Tue, 08 Oct 2003, 16:08:58 EDT
 % Created on : Tue, 09 Feb 1999, 09:40:47 EST
 %
 % DEPENDENCIES:
 %   call:
 %
 % REMARKS:
 %
 %*********************************************************************
 clear;
 close;
 % model with 3-D beam elements with 6 dofs per node
   dof = 6;         %  dof per node
 % input the nodal coordinates
 %
   n_node = 3;             % total number of nodes
   n_elem = 2;             % total number of elements
   total_dof = 6 * n_node;  % total dofs of system
   Le1 = 4;   % beam length
   Le2 = 2;   % beam length
   thetae1 = pi/6;
   thetae2 = (-pi/4);
   position(:, 1) = [ 0; 0];
   position(:, 2) = [ Le1*cos(thetae1); Le1*sin(thetae1)];
   position(:, 3) = [ Le1*cos(thetae2) + Le2*cos(thetae2); Le1*sin(thetae1) + Le2*sin(thetae2)];
   % set up the nodal coordinate arrays x, y, x
   for i = 1 : n_node
      x(i) = position(1,i);
      y(i) = position(2,i);
   end
 % set up the element connectivity array node_connect, defined as
 % follows:
 % node_connect(local node number, element number) = global node number
   node_connect(1, 1) = 1;   % element 1
   node_connect(2, 1) = 2;
   node_connect(1, 2) = 2;   % element 2
   node_connect(2, 2) = 3;  
   % plot the whole truss system
   % loop over the elements
   for i =  1 : n_elem  
      % for element i, do the following:
      % node_1 = global node number corresponding to the local node 1 
      node_1 = node_connect(1,i);
      % node_2 = global node number corresponding to the local node 1 
      node_2 = node_connect(2,i);
      % xx : 1x2 array containing x coordinates of node_1 and node_2
      xx = [x(node_1),x(node_2)];
      % yy : 1x2 array containing y coordinates of node_1 and node_2
      yy = [y(node_1),y(node_2)];
      % plot the element i in 2-D using command plot
      % use axis command to avoid having the plot window too tight
      axis([0 5 0 5])
      % solid line
      plot(xx,yy,'-')
      % hold the current figure for the next element
      hold on
   end
    text(x(node_connect(1,1)),y(node_connect(1,1)),'Global Node 1', 'HorizontalAlignment', 'left')
    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', 'right')
    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', 'left')
 
    hold on
 
    d3 = 4.352;
    d4 = 6.1271;
 
    position_d(:, 1) = [ 0; 0];
    position_d(:, 2) = [ Le1*cos(thetae1)+ d3; Le1*sin(thetae1)+ d4 ];
    position_d(:, 3) = [ Le1*cos(thetae2)+Le2*cos(thetae2); Le1*sin(thetae1)+Le2*sin(thetae2)];
 
    % set up the nodal coordinate arrays x, y
    for i = 1 : n_node
       x(i) = position_d(1,i);
       y(i) = position_d(2,i);
    end
 
 % set up the element connectivity array node_connect, defined as
 % follows:
 % node_connect(local node number, element number) = global node number
 
 node_connect_d(1, 1) = 1;   % element 1
 node_connect_d(2, 1) = 2;
 
 node_connect_d(1, 2) = 2;   % element 2
 node_connect_d(2, 2) = 3;
 
    % plot the whole truss system
    % loop over the elements
    for i =  1 : n_elem  
       % for element i, do the following:
       % node_1 = global node number corresponding to the local node 1 
       node_1 = node_connect_d(1,i);
       % node_2 = global node number corresponding to the local node 1 
       node_2 = node_connect_d(2,i);
       % xx : 1x2 array containing x coordinates of node_1 and node_2
       xx = [x(node_1),x(node_2)];
       % yy : 1x2 array containing y coordinates of node_1 and node_2
       yy = [y(node_1),y(node_2)];
       % plot the element i in 2-D using command plot
       % use axis command to avoid having the plot window too tight
       axis([0 8 0 9])
       % solid line
       plot(xx,yy,'--g')
       % hold the current figure for the next element
       hold on
    end
 
    text(x(node_connect_d(1,1)),y(node_connect_d(1,1)),'Global Node 1', 'HorizontalAlignment', 'left')
    text(x(node_connect_d(2,1)),y(node_connect_d(2,1)),'Global Node 2', 'HorizontalAlignment', 'right')
    text(x(node_connect_d(2,2)),y(node_connect_d(2,2)),'Global Node 3', 'HorizontalAlignment', 'center')
    text(x(node_connect_d(2,1))/2,y(node_connect_d(2,1))/2,'Element 1', 'HorizontalAlignment', 'right')
    text( x(node_connect_d(2,1)) + (x(node_connect_d(2,2))-x(node_connect_d(2,1)))/2 , y(node_connect_d(2,1))...
        + (y(node_connect_d(2,2))-y(node_connect_d(2,1)))/2 ,'Element 2', 'HorizontalAlignment', 'left')
 
   % put the title on the figure
   title('Two-bar Truss System: Undeformed and Deformed')
   % label x axis
   xlabel('x')
   % label y axis
   ylabel('y')

The blue solid line represents the undeformed truss system, while the green dashed line represents the deformed truss system.

Contributing Team Members

The following students contributed to this report

• Michael Berry Eml4500.f08.qwiki.berry 10:41, 8 October 2008 (UTC)
• Andrea Booher Eml4500.f08.qwiki.booher 11:32, 8 October 2008 (UTC)
• Chris Bishop Eml4500.f08.qwiki.bishop 13:02, 8 October 2008 (UTC)
• David Nobles Eml4500.f08.qwiki.nobles 15:29, 8 October 2008 (UTC)
• Eric Gatch Eml4500.f08.qwiki.gatch 12:15, 8 October 2008 (UTC)
• Lee HarrisEml4500.f08.qwiki.harris 17:16, 8 October 2008 (UTC)