User:Eml4500.f08.delta 6.guzman/MATLAB Two Bar Truss Example Modified

From Wikiversity
Jump to navigation Jump to search

Two Bar Truss Example (MATLAB)[edit]

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

Take this truss example with the following:
Data

Length Cross-sectional Area Young's Modulus Inclination Angle
Element 1 L(1) = 4

A(1) = 1

E(1) = 3

θ(1) = pi/3

Element 2 L(2) = 2

A(2) = 2

E(2) = 5

θ(2) = pi/4


MATLAB Code

% 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

Note: this MATLAB code was obtained from Dr. L. Vu-Quoc website http://clesm.mae.ufl.edu/~vql/courses/fead/2008.fall/codes/twoBarTrussEx.txt

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

function k = PlaneTrussElement(e, A, coord)
% PlaneTrussElement(e, A, coord)
% Generates stiffness matrix for a plane truss element
% e = modulus of elasticity
% A = area of cross-section
% coord = coordinates at the element ends

x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
k = e*A/L*[ls^2, ls*ms,-ls^2,-ls*ms;
ls*ms, ms^2,-ls*ms,-ms^2;
-ls^2,-ls*ms,ls^2,ls*ms;
-ls*ms,-ms^2,ls*ms,ms^2];

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

function [d, rf] = NodalSoln(K, R, debc, ebcVals)
% [nd, rf] = NodalSoln(K, R, debc, ebcVals)
% Computes nodal solution and reactions
% K = global coefficient matrix
% R = global right hand side vector
% debc = list of degrees of freedom with specified values
% ebcVals = specified values
dof = length(R);
df = setdiff(1:dof, debc);
Kf = K(df, df);
Rf = R(df) - K(df, debc)*ebcVals;
dfVals = Kf\Rf;
d = zeros(dof,1);
d(debc) = ebcVals;
d(df) = dfVals;
rf = K(debc,:)*d - R(debc);

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

function [eps, sigma, force] = PlaneTrussResults(e, A, coord, disps)
% [eps, sigma, force] = PlaneTrussResults(e, A, coord, disps)
% Compute plane truss element results
% e = modulus of elasticity
% A = Area of cross-section
% coord = coordinates at the element ends
% disps = displacements at element ends
% The output quantities are eps = axial strain
% sigma = axial stress and force = axial force.

x1=coord(1,1); y1=coord(1,2);
x2=coord(2,1); y2=coord(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
ls=(x2-x1)/L; ms=(y2-y1)/L;
T=[ls,ms,0,0; 0,0,ls,ms];
d = T*disps;
eps= (d(2)-d(1))/L;
sigma = e.*eps;
force = sigma.*A;

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

Results

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