User:Eml4500.f08.Ateam.barry

Finite Element Analysis

MatLab Tutorial Summary/Explanation


MatLab is developed by the MathWorks Company. Below is their link for further information outside the scope of the document below.

About MatLab: MatLab is a matrix-based enginering computation program that is user-interactive. By using matrices to enter and solve complex equations, MatLab can do in a fraction of the time it would take older programming languages to do. When dealing with finite element analysis, depending on the number of nodes, there can be numerous partial differential equations that need to be solved.

Starting Matlab: MatLab can be started 2 ways. The first is using the command window and typing run/matlab. The second is through the start menu >> programs >> matlab.

Matrices: MatLab works solely with matrices to handle all computations. These matrices can range in size from 1 by 1 all the way up to thousands of elements. Entering a matrix into matlab is easy and is done with the following notation: A = [1 2 3; 4 5 6; 7 8 9] This creates a 3 by 3 matrix with the values shown and a new row started after each semi-colon. Additionally, matrices do not have to be entered with just real numbers. Complex numbers (i) can also be entered into matrices. For example: A = [1+3i; 2-4i; 3]

When it comes to large matrices, it may be very time consuming to enter in all the matrix values by hand. In such a case, the values can be read in from an ASCII file. As long as the file has a rectangular matrix in it, simply type: >> load numbers.ext This will load the matrix in the file numbers into the workspace.

Matrix Operations: MatLab uses the same basic mathematical operators plus a few additional ones. Below are the operators: * + addition * - subtraction * * muntiplication * / right division * \ left division * ^ power * ` conjugate The left division operator is special in that it allows for division of matrices that are not square. If you use the right division operator, MatLab will give an error saying that the matrices are not square.

Statements, Expressions, and Variables:

Matlab is an expression language meaning the user types an expression, usually composed of variables, functions, and operators, and Matlab interprets it and performs the expression. After the expression is evaluated, the output is put into a matrix. If no matrix name is specified by the user, matlab puts the result in a matrix called ans. Normally after each expression is evaluated, Matlab outputs the result of the expression. To supress(hold) the output, put a semi-colon at the end of the expression. Variables can have virtually any name you choose, as long as it is not a reserved variable or function name. Matlab is case sensitive but has no restriction on variable name length. Some Helpful Commands: * whos - displays a list of variables currently in the workspace Also, there are come commands that can be used to help the user build matrices Matrix Building Commands: * eye - identity matrix * zeros - matrix of zeros * ones - matrix of ones Each of the commands needs to be followed by the following notation: (# of rows, # on columns) Example: zeros(3,3) yeilds a square 3 by 3 matrix of zeros.

Loops:

Loops are used to repeat a series of operations until a certain condition is met. Below are the two types of loops that are used frequently in Matlab.

For:

For loops exectue continually once a condition is met up until the condition is no longer met.

function problem_1_37

% This function will plot 3 conditional functions

dx = 0.1; x = [-2:dx:6]; for k = 1:length(x) if x(k) >= 5 y(k) = 10*(x(k)-5) + 1; elseif x(k) >= -1 y(k) = 2 + cos(pi*x(k)); else y(k) = exp(x(k)+1); end end

plot(x,y), xlabel('Time x (seconds)'), ylabel('Height y (kilometers)')

While:

While loops are similar to for loops except for while loops are guaranteed to execute at least once. It will exectue the statement, then check it against the while statement. If it is true, it will run back through the loop. If it is false, it will exit the loop.

function e = macheps

% This function determines the smallest number that can be added to 1 in % Matlab and still get a number greater than 1.

e = 1; while (1+e >1) e = (e/2); end

e = e*2;

'''In addition to loops, one can have an operation performed only IF a condition is met. That type of operator is an if statement. Below is an example of an IF statement.'''

function y = problem_1_36(x) % This function will check conditional statements to evaluate a function given different input values. if x >= 5 y = 10*(x-5)+1; elseif x >= -1 y = 2+cos(pi*x); else y = exp(x+1); end

NOTE: All code was taken from my previous homework assignments in Numerical Methods

Relations:

The relational operators are used to compare the value of a variable to another variable or to a constant value. Below are some common relational operators: * == equal * ~= not equal * < less than * > greater than * >= greater than or equal to    * <= less than or equal to     * & and * | or    * ~ not

Relational operators are always used in for, while,and if statements to see when to execute, or stop executing expressions within the loop.

Scalar Functions:

Certain matlab functions operate on scalar numbers and when applied to a matrix, they act on 1 matrix element at a time. Some scalar Matlab functions are listed below: sin, cos, tan, asin, acos, atan exp, log, rem, abs, sqrt, sin, round, floor, ceil

Vector Functions:

Vector functions act on a vector matrix column by column to produce a row matrix as the solution. These functions can only be used when the matrix is larger than 1 by 1. Below are some examples of vector functions: max, min, sort, sum, prod, median, mean, std, any, all

Command Line Editing:

The command line in Matlab is denoted by 2 arrows (>>). The command line is the line that the user uses to type in commands. To keep the user from having to re-enter previously used commands, Matlab allows the user to scroll through the old commands used during the session by hitting the up-arrow and down-arrow keys.

Also, it is important to remember that instead of entering a series of commands by hand every time, create an m-file to execute. This will save loads of time and needless typing.

M-Files

M-files are special Matlab files that have a .m extension. These files contain a series of statements that can be exectued at one time, thus saving the user time from having to exectue expressions 1 at a time. There are two types of m-files, script files and function files.

Script files contain a series of normal Matlab scripts. They are traditionally used to enter large amounts of data into matrices. Any variable or matrix is global, meaning the parameter can be used in any other portion of the workspace.

Function files are a type of m-file that allow users to essentially create new Matlab functions. The functions can be specific to a particuar problem that the user is working on.

Text Strings

You can assign text strings to variables by using single quotes around the text to be entered into the variable. Additionally, text can be displayed to the screen using the disp function. The error function also outputs a message to the screen but also aborts the code in an m-file when it is executed. Below are examples of the 2 functions: disp('Hello World') error('Warning!!!')

Homework # 2 Lecture Notes
Assembly of Force Displacement Matrix The global force displacement matrix is assembled from the 2 element displacement matrices.

The Global Stiffness Matrix with symbolic representation:
 * $$\begin{bmatrix}

K^1_{11}&K^1_{12}&K^1_{13}&K^1_{14}&0&0\\ K^1_{21}&K^1_{22}&K^1_{23}&K^1_{24}&0&0\\ K^1_{31}&K^1_{32}&K^1_{33}+K^2_{11}&K^1_{34}+K^2_{12}&K^2_{13}&K^2_{14}\\ K^1_{41}&K^1_{42}&K^1_{43}+K^2_{21}&K^1_{44}+K^2_{22}&K^2_{23}&K^2_{24}\\ 0&0&K^2_{31}&K^2_{32}&K^2_{33}&K^2_{34}\\ 0&0&K^2_{41}&K^2_{42}&K^2_{43}&K^2_{44} \\ \end{bmatrix} $$

Note: The overlap in the 2 matrices comes from the fact that global node 2 is shared between both elements.

Calculating the Matrix Values:

K11 = k111 = 0.5625 K12 = k112 = 0.3247 K33 = k133 + k211 = 3.0625 K34 = k134 + k212 = -2.1752 K43 = k143 + k221 = -2.1752 K44 = k144 + k222 = 2.6875

Elimination of the known degrees of freedom reduces the global force displacement relationship. d1 = d2 = d5 = d6 = 0.

Two Bar Truss Matlab Example

Here is the Matlab code run for this example:

% 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

Here are the results for the code:

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

Homework # 3 Lecture Notes 9/29/08
Use the global force-displacement relationship

HW: Use the force displacement relationship: (6x2)x(2x1)=(6x1) Note: Only need to do computation for rows 1, 2, 5, 6 to get reactions F1. F2, F5, F6.

Take the 6x2 matrix:


 * $$\begin{bmatrix}

K^1_{13}&K^1_{14}\\ K^1_{23}&K^1_{24}\\ K^1_{33}+K^2_{11}&K^1_{34}+K^2_{12}\\ K^1_{43}+K^2_{21}&K^1_{44}+K^2_{22}\\ K^2_{31}&K^2_{32}\\ K^2_{41}&K^2_{42}\\ \end{bmatrix} $$

multiplied by the 2x1 matrix:


 * $$\begin{bmatrix}

d_3&d_4\\ \end{bmatrix} $$

and that yeilds the 6x1 reaction force matrix:


 * $$\begin{bmatrix}

F_1\\ F_2\\ F_3\\ F_4\\ F_5\\ F_6\\ \end{bmatrix} $$

Closing the loop between the FEM and Statics
-To close the loop between the finite element method and the method of statics, we use the principle of virtual displacement to solve for the reactions and thus the member forces in the 2 force body.

Two Bar Truss System

To solve the two bar truss system and thus close the loop here is the progression of steps:

FEM -> Compute Displacement -> Compute Reactions Statics -> Compute Reactions -> Compute Displacement When solved, both of these solution steps should yeild the same values for the reactions and displacement. This is how we can check our work.

We are given the 2-Force bodies for each element(1,2). By statics, the reactions are known and thus the member forces are known. Member Forces: P1 and P2 Next, we compute the axial displacement degrees of freedom. (This gives us the amount of extension in the bars) q2 = P1/K1 = BC  q1 = 0  (fixed node 1) q1 = -P2/K2 = AB  q2 = 0   (fixed node 3)

Question: How do we back out from the above results, the displacement degrees of freedom of node 2?

The dotted lines represent the deformed object and the solid lines represent the original object. The distances AB and CD are the displacements of node 2. Given the forces P1 and P2 we can back out the displacement of node 2.

Below is an explanation of the displacement in the figure above. This is why we can approximate the displacement in 1 direction as being zero and being the length of the bar multiplied by the displacement angle.

Homework Report # 4 Lecture Notes 10/15/08
Justification of assembly of element stiffness martix into global stiffness matrix.

k(e) = element stiffness matrix k = global stiffness matrix

Consider the example of the 2 bar truss: Recall the element force displacement (FD) relationship: k(e)*d = f(e)

Use the euler cut principle (Method 2) to solve for the forces. The euler cut principle states that node 2 should be in equilibrium. We take this and solve for the forces in each of the 2 bars.

Free Body Diagrams of elements 1 and 2.



The element degrees of freedom are given in the matrix: d(e). This is a 4x1 matrix.

For node 2, identify the global degrees of freedom to element degrees of freedom for both of the bars.



Now we can use the equilibrium concept from the euler cut method to solve for the forces.

Sum of forces in x direction: F = 0 = -f(1)3 - f(2)1 Sum of forces in y direction: F = 0 = P -f (1)4 - f(2)2

Next we use the element force displacement relationship to obtain the displacement of node 2: k(e)*d = f(e)

Lecture Notes 11/12/08
FEM Via Principle of Virtual Work (Continued from 11/10/08)

The finite element method uses a weighting matrix W to solve for the solution. The first set of solutions is for a bar element i. The PVW solution followes a similar pattern with the displacement matrix being replaced by the weighting matrix W.

$$ U(x_{i+1}) = d_{i+1} $$

Apply the same interpretation for the W(x) matrix: For example,

$$ W(x) = N_iW_i + N_{i+1}W_{i+1} $$

Element Stiffness Matrix for element i:

$$ \int_{x_i}^{x_{i+1}}{} [N_iW_i + N_{i+1}W_{i+1}](EA)[N_id_i + N_{i+1}d_{i+1}]dx $$

The Definition of $$ N_i and N_{i+1} $$ can be seen below:

$$ N_i = \frac{dN_i(x)}{dx} $$       $$ N_{i+1} = \frac{dN_{i+1}(x)}{dx} $$

Note:

$$ u(x) = \N_i(x) \N_{i+1}(x)| \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} $$

The N(x) matrix is a 1x2 matrix and the d matrix is 2x1.

Taking the derivative yeilds the following:

$$ \frac{du(x)}{dx} = \N_i'(x) \N'_{i+1}(x)| \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} $$

The N'(x) matrix is a 1x2 matrix and is called the B(x) matrix.

Similarly for the W matrix from the Principle of Virtual Work (PVW):

$$ W(x) = N(x) \begin{Bmatrix} W_i\\ W_{i+1} \end{Bmatrix} $$

$$ \frac{dW(x)}{dx} = B \begin{Bmatrix} W_i\\ W_{i+1} \end{Bmatrix} $$

Recall the Element DOFs:



$$ \begin{Bmatrix} d_i\\ d_{i+1} \end{Bmatrix} = \begin{Bmatrix} d_1^i\\ d_2^i \end{Bmatrix} d^{(i)} $$

$$ \begin{Bmatrix} W_i\\ W_{i+1} \end{Bmatrix} = \begin{Bmatrix} W_1^i\\ W_2^i \end{Bmatrix} W^{(i)} $$

From these relations, we can find the integral that contains both W and d.

$$ \int_{x_i}^{x_{i+1}}{} (BW^i)(EA)(Bd^i)dx $$

The B matrix relates d to W. (Principle of Virtual Work)

Our Goal now is to derive $$ W^i(k_id_i) $$ from the previous integral.

Re-write the integral with the EA in the front:

$$ \int_{x_i}^{x_{i+1}}{} (EA)(BW^i)(Bd^i)dx $$

Take the Transpose of the BW factor to eliminate the B matrix:

$$ \int_{x_i}^{x_{i+1}}{} (EA)(BW^i)^T(Bd^i)dx $$

Distribute the Transpose inside of the Parentheses:

$$ \int_{x_i}^{x_{i+1}}{} (EA)(B^TW^{iT})(Bd^i)dx $$

The W matrix is equal to its transpose so it can be removed:

$$ \int_{x_i}^{x_{i+1}}{} (EA)(B^TW^{i})(Bd^i)dx $$

What we are left with is ( With the W matrix pulled outside the integrand because it does not vary with x):

$$ W^{i} \int_{x_i}^{x_{i+1}}{}(B^T)(EA)(B)dx d^i$$

The B matrix is a function of x and is equal too:

$$ B(x) = \|\frac{-1}{L^i} \|\frac{-1}{L^i}| $$

Note: $$ L^{i} = x_{i+1} - x_i $$