User:Eml4500.f08.bike.mcdonald/homework report 3

Element Force-Displacement Relation With Respect To The Global Coordinate System
The element force displacement (FD) relation can be determined with respect to the global coordinate system. The following equation shows the general force displacement relation for an element member.

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

The local free body diagram for the element member is pictured in the figure below.



The local axial force displacement diagram for the element member is pictured in the figure below.



The following force displacement equation can be constructed from the above axial force displacement diagram.

$$\mathbf{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}$$
 * $$ \mathbf{\hat{\; k}}^{(e)}_{2x2}\,\;\;\;\;\;\;\;\; \mathbf{q}^{(e)}_{2x1} \;\; = \;\;\;\; \mathbf{p}^{(e)}_{2x1}$$

Where:
 * qi(e) = axial displacement of element e at local node i.
 * pi(e) = axial force of element e at local node i.

The goal is to find the relations between $$\begin{cases} \mathbf{q}^{(e)}_{2x1} & \text{ and}\;\;\;\;\;   \mathbf{d}^{(e)}_{4x1}  \\ \mathbf{p}^{(e)}_{2x1} & \text{ and}\;\;\;\;\;  \mathbf{f}^{(e)}_{4x1} \end{cases}$$

One of the relations can be expressed in the form below.

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

The figure below depicts the displacement vector of local node 1 as well as a portion of the element member e.



In the figure above, the displacements $$d^{(e)}_1$$ and $$d^{(e)}_2$$ in the $$\vec{i}\,$$ and $$\,\vec{j}$$ directions are shown as components of $$\mathbf{\hat{\; d}}\,$$1(e). The equation of the $$\vec{i}\,$$ and $$\,\vec{j}\,$$ components of $$\mathbf{\hat{\; d}}\,$$1(e) can be seen in the following equation.

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

q1(e) = axial displacement of node 1 is the orthagonal projection of the displacement vector $$\bar{d}^{(e)}_1$$ of node 1 on the $$\tilde{x}$$ axis of element e. This is represented in the equation below.

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

The displacement vector of node 1 can be broken down into components. $$q^{(e)}_1= (d^{(e)}_1\vec{i}+d^{(e)}_2\vec{j}) \cdot \vec{\tilde{i}}$$

Which is equivalent to. $$q^{(e)}_1= d^{(e)}_1(\vec{i}\cdot \vec{\tilde{i}})+d^{(e)}_2(\vec{j} \cdot \vec{\tilde{i}})$$
 * $$ \cos \theta _e \;\;\;\;\; \sin \theta _e$$
 * $$ l^{(e)} \;\;\;\;\;\;\;\;\;\; m^{(e)}$$

q1(e) is a 1x1 scalar which can be calculated with the displacements of node 1 and the director cosines as shown in the following equations.

$$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{Bmatrix}d_1^{(e)}\\ d_2^{(e)}\end{Bmatrix}$$ $$1x1 \;\;\;\;\;\;\;\;\; 1x2 \;\;\;\;\;\;\;\;\; 2x1$$

The matrices that result from combining the displacement equations for nodes 1 and 2 are listed below. $$\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}^{(e)}=\mathbf{T}^{(e)}\mathbf{d}^{(e)}$$

Now, the equation below can be used to relate the axial forces to the internal forces and the transfer function. $$\textbf{P}_{2x1}^{(e)}=\textbf{T}^{(e)}\textbf{f}_{4x1}^{(e)}$$

Expanding the matrices out, it is seen that $$\textbf{P}^{(e)}$$ is comprised of the two axial forces while $$\textbf{f}^{(e)}$$ is comprised of the four internal forces.

$$\begin{Bmatrix}P_1^{(e)}\\P_2^{(e)}\end{Bmatrix}=\textbf{T}^{(e)}\begin{Bmatrix}f_1^{(e)}\\ f_2^{(e)} \\ f_3^{(e)} \\f_4^{(e)} \end{Bmatrix}$$

It can be observed that in order for the equation to be computable, the transfer matrix must be a 2 x 4 rectangular matrix.

Recalling the elemental axial force-displacement relation, the transfer matrix can be used to create an equation that replaces the axial force and displacement matrices with the elemental displacement and force matrices.

$$\hat{\textbf{k}}^{(e)}\textbf{q}^{(e)}=\textbf{p}^{(e)}\ = \hat{\textbf{k}}^{(e)}(\textbf{T}^{(e)}\textbf{d}^{(e)})\ = (\textbf{T}^{(e)}\textbf{f}^{(e)})$$

In order to eliminate all axial forces, the transfer matrix must be moved from the right hand side of the equation to the left hand side. In order to do this, the inverse of the transfer matrix must be determined.

Unfortunately, the transfer matrix was previously determined to be a $$2 x 4$$ rectangular matrix and thus an inverse cannot be taken. To fix this problem, the transpose will be used instead of the inverse to find the elemental forces.

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

The principle of virtual work can also be used to reduce this matrix down to one that would be easier to compute. The principle is based on the reduction of the global force-displacement relation and can be reduced to what is shown below by assuming that $$d_1 = d_2 = d_5 = d_6 = 0$$.

$$\textbf{k}_{6x6}\textbf{d}_{6x1}=\textbf{F}_{6x1}$$ --> $$\bar{\textbf{k}}_{2x2}\bar{\textbf{d}}_{2x1}=\bar{\textbf{F}}_{2x1}$$

The matrix cannot be solved by simply moving the stiffness matrix over because the stiffness matrix is singular. Singular matrices are defined as having a determinant which equals zero. Since the inverse of a matrix requires one over the determinant, this creates an invertible situation.

In an unconstrained structural system, there are three possible rigid body movements: two translations and one rotation.

Using dynamics and vibrations, the eigenvalues found can be used in the equation shown below where $$\textbf{k}$$ is the stiffness matrix, $$\lambda$$ is the eigenvalue, and $$\textbf{M}$$ is the mass matrix. This equation relates the stiffness and the mass matrices.

The four zero eigenvalues correspond to the system having no stored elastic energy. This becomes a rigid body mode. Modes can also be referred to as the eigenvectors.

Alternative Method For Determining Reactions of Two Member Truss System
When solving for the reactions of the two member truss system, the only method described so far was to use the Element FD Relation in order to relate the Element Displacement Matrices with the Element Force Matrices. Another method that can be used to solve for the reactions of the system is to use the Global FD Relation.

The use of the Global FD for this purpose was assigned as a homework assignment, and is illustrated below.

Closing the Loop Between the Finite Element Method and Statics
The key to closing the loop between the Finite Element Method and Statics is to apply the Principle of Virtual Work. This principle has been used in Homework Report 2 when attempting to reduce the Global Stiffness Matrix.



The image below shows a Free Body Diagram of a two-force member along with the relationship between the axial forces in the member.



Due to the axial forces that exist in the two force members in the figure below, a certain amount of axial deformation, $$q$$ will occur. The figure below shows the two member truss shape before the given load, $$P$$ is applied to the system causing axial forces and displacements.



Because the reactions in both two-force members in the two member truss system are now known, the methods of statics allow for the axial member forces $$P_{1}^{(1)}$$ and $$P_{2}^{(2)}$$ to be determined. Knowing $$P_{1}^{(1)}$$ and $$P_{2}^{(2)}$$, the axial displacement DOFs (the amount of extension of each member) can be determined as shown below.

$$

q_{1}^{(1)}=0,\text{ fixed node 1}

$$

$$

q_{2}^{(1)}=\dfrac{P_{1}^{(1)}}{k^{(1)}}=\dfrac{P_{2}^{(1)}}{k^{(1)}}

$$

$$

q_{1}^{(2)}=\dfrac{-P_{2}^{(2)}}{k^{(2)}}

$$

$$

q_{2}^{(2)}=0,\text{ fixed node 2}

$$

Knowing the axial displacements in each member, the overall displacement DOFs of Node 2 can be determined using a graphical process shown in the figure below.



In the above figure, the red dotted line is the deformed shape of the two member truss system, caused by application of the load, $$P$$. The geometric length  corresponds to the value of axial displacement in Element 1, $$q_{2}^{(1)}$$, whereas the geometric length   corresponds to the value of axial displacement in Element 2, $$ q_{1}^{(2)}$$.

NOTE: The reason that segments  and   are perpendicular to   and , respectively, is due to the Principle of Virtual Work and the use of extremely small angles of deflection in the principle. The below figure shows the derivation of this principle by means of the small angle theorem.



The Small Angle Theorem states that for small angles, $$\alpha$$, $$\sin {\alpha}=\alpha$$ and $$\cos {\alpha}=1$$. Applying this theorem to the relations for $$u_{x}$$ and $$u_{y}$$ shown in the above figure, the following values are determined.

$$

u_{y}=R\sin {\alpha }\cong R\alpha

$$

$$

u_{x}=R(1-\cos \alpha )\cong 0

$$



The method of infinitesimal displacement, which is related to virtual displacement, is used to analyze the two bar truss system.

Some sample calculations for the lengths of AC and AB were performed as follows:

$$AC = \frac{|P_{2}^{(1)}|}{k^{(1)}} = \frac{5.1243}{.75} = 6.8324$$

$$ AB = \frac{|P_{1}^{(2)}|}{k^{(2)}} = \frac{6.276}{5} = 1.255$$

The (x,y) coordinates $$(x_{B}, y_{B})$$ and $$(x_{C}, y_{C})$$ must then be found. This can be done using simple trigonometry and the values of the magnitude of the line and the values of theta obtained in HW 2. This method is further explained using the figure on the left, along with the intermediate calculation steps shown below.

2 unknowns $$(x_{D}, y_{D})$$

The equations for lines AB and BC are then needed.

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

$$= (x - x_{p})\vec{i} + (y - y_{p})\vec{j}$$

$$\Rightarrow x - x_{p} = (PQ)cos\Theta$$

and $$y - y_{p} = (PQ)sin\Theta$$

$$\Rightarrow \frac{y - y_{p}}{x - x_{p}} = tan\Theta$$

$$y - y_{p} = (tan\Theta )(x - x_{p})$$

Using this equation along with $$\Theta^{1} = 30^{o}$$, it can be found that $$(x_{C}, y_{C})$$ = $$(5.92, 3.416)$$

Also, the equation for a line perpendicular to the above line, passing P is as follows:

$$y - y_{p} = (tan(\Theta + \frac{\pi}{2}))(x - x_{p})$$

This is useful in finding the value of $$(x_{B}, y_{B})$$, and using the value of $$\Theta^{2} = 135^{o}$$ it can be found that $$(x_{B}, y_{B})$$ = $$(-.877, .877)$$

To find the value of $$(x_{D}, y_{D})$$ one can use the following equation:

$$ \vec{AD} = (x_{D} - x_{A})\vec{i} + (y_{D} - y_{A})\vec{j}$$

Where the value of $$(x_{A}, y_{A})$$ are zero due to choosing the axis' origin to be point A.

By definition:

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

This means that, using the values obtained in HW 2, the line AD can be expressed as:

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

Which is the displacement vector of A from FEM. This ends closing the loop between statics and the Finite Element Method. The next topic will cover the analysis of a three-bar truss system.

Analysis of the Three-Bar Truss System
The three bar truss system shown above will now be analyzed using the values below as parameters.

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

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

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

$$P = 30$$ $$\Theta ^{(1)} = 30^{o}$$ $$\Theta ^{(2)} = 30^{o}$$ $$\Theta ^{(3)} = 45^{o}$$

Care should be taken when numbering the local nodes in order to make the assembly of the $$K$$ matrix convenient.

The 3-bar truss system is shown with a known force P and three unknown axial reaction forces $$R_{1}, $$$$R_{2}, $$$$ R_{3}$$.

In this system, the statics method yields the following equations. $$

\sum{F_{X}=0} $$ $$ \sum{F_{Y}=0}$$ $$ \sum{M_{A}=0}$$

The moment equation about point A does not solve for any of the unknown force values because all the forces' lines of action go through point A. This means that only the first two equations, the sum of the forces in the X and Y directions can be used. In order to solve this system statically, a third equation must be found. A moment equation around another point could have been the answer, but as shown in the following 3-D explanation, the moment equation is proven to be unusable because it is also equal to zero.

Analyzing the system shown on the left, the following equations are formulated. $$ \sum{M_{B}}=\vec{BA}\times \vec{F}=\vec{BA'}\times \vec{F}$$

$$A'$$ is assumed to be on the line of action of $$\vec{F}$$. The vector $$BA'$$ can be rewritten as shown. $$ \vec{BA'}=\vec{BA}\times \vec{AA'}$$ Therefore, the moment around point B can be rewritten as shown. $$ \sum{M_{B}}=(\vec{BA}+ \vec{AA'})\times \vec{F}= \vec{AA'}\times \vec{F}+\vec{BA}\times \vec{F} $$

$$\vec{AA'}\times \vec{F}$$ becomes zero because $$\vec{AA'} $$ and $$\vec{F}$$ are in the same direction.

Now the moment equation around B is as shown: $$\sum{M_{B}}=\vec{BA}\times \vec{F} $$

Now, back to analyzing a specific 3-bar truss system, the diagram is shown below. Node A in the 3-bar truss system is in equilibrium. $$\sum_{i=0}^{3}{F_{i}=\vec{0}}$$ The moment equation around point B, which was derived above, is rewritten as shown. $$\sum_{i}^{}{M_{B,i}=\sum_{i}^{}{\vec{BA'_{i}}}}\times \vec{F_{i}}$$ $$A'_{i}$$= any point on the line of action of $$\vec{F_{i}}$$. The moment equation around point B can be rewritten as shown. $$\sum_{i}^{}{M_{B,i}=\sum_{i}^{}{\vec{BA}}}\times \vec{F_{i}}$$ $$\sum_{i}^{}{M_{B,i}=\vec{BA}}\times \sum{\vec{F_{i}}}$$ As stated above, the sum of the forces is equal to zero due to point A being in equilibrium. The moment around B is proved to be zero.

$$\sum_{i}^{}{M_{B,i}=\vec{0}}$$

With this conclusion, it was determined that the system is statically indeterminate. The Finite Element Analysis Method must be used to solve the 3-bar truss system.

To begin the Finite Element Method, the stiffness matrix is created. This matrix is shown below.

Plot of Undeformed and Deformed Two-Bar Truss System Using MATLAB
The plot below shows the deformed and undeformed states of the two-bar truss system analyzed during class. This plot was created through the use of a MATLAB script developed by Andrew J. McDonald, using the MATLAB script, developed by X.G. Tan and hosted on Professor Loc Vu Quoc's University of Florida website.



The blue dotted lines represent the undeformed state of the two-member truss system, and the solid red lines represent the deformed state of the two-member truss system.

The code of the MATLAB script developed to produce this plot is given below.

% Matlab Script %********************************************************************* % filename: two_bar_truss_def_and_und.m  % % PURPOSE: %   Plot the two-bar truss system for both deformed %   and undeformed cases % % AUTHOR: Andrew J. McDonald % % Modified on: ______________________________ % Created on : Sat, 04 Oct 2008, 11:12:54 EST % % DEPENDENCIES: %  call: % % REMARKS: % %*********************************************************************

clear; close;

% model with 2-D beam elements with 2 dofs per node dof = 2;        %  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

% two-bar truss system data %    L_1 = 4; L_2 = 2; theta_1 = 30*(pi/180); theta_2 = -45*(pi/180);

position(:, 1) = [ 0; 0]; position(:, 2) = [ L_1*cos(theta_1); L_1*sin(theta_1)]; position(:, 3) = [ L_1*cos(theta_2)+L_2*cos(theta_2); L_1*sin(theta_1)+L_2*sin(theta_2)];

% set up the nodal coordinate arrays x, y    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', '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,'Global 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 ,'Global Element 2',...     'HorizontalAlignment', 'center')

hold on

d_3 = 4.352; d_4 = 6.1271;

position_d(:, 1) = [ 0; 0]; position_d(:, 2) = [ L_1*cos(theta_1)+ d_3; L_1*sin(theta_1)+ d_4 ]; position_d(:, 3) = [ L_1*cos(theta_2)+L_2*cos(theta_2); L_1*sin(theta_1)+L_2*sin(theta_2)];

% 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,'-r') % 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', 'center') text(x(node_connect_d(2,1)),y(node_connect_d(2,1)),'Global Node 2',...      'HorizontalAlignment', 'center') 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,'Global Element 1',...      'HorizontalAlignment', 'center') 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 ,'Global Element 2',...     'HorizontalAlignment', 'center')

% put the title on the figure title('Deformed and Undeformed States of Two-Bar Truss System') % label x axis xlabel('x') % label y axis ylabel('y')

Contributing Team Members
Andrew McDonald - Eml4500.f08.bike.mcdonald 01:31, 5 October 2008 (UTC) Garrett Pataky - (Eml4500.f08.bike.pataky 00:55, 7 October 2008 (UTC)) Sam Bernal - Eml4500.f08.bike.bernal 15:11, 7 October 2008 (UTC) Bobby Sweeting - Eml4500.f08.bike.sweeting 22:51, 6 October 2008 (UTC) Shawn Gravois - Eml4500.f08.bike.gravois 00:26, 30 September 2008 (UTC) Eric Viale - Eml4500.f08.bike.viale 02:50, 7 October 2008 (UTC)