FEA code in Matlab for a
Truss Structure
By: Shiwei Zhou
Outline:
MATLAB programs for 2D truss
Verification of the MATLAB code by
Strand7
2D Truss transmitter tower
Strand7 example for 2D truss
transmitter tower
Input Geometrical Model
1. Node coordinate
2. Element connection
3. Force vector
4. Displacement vector
5. Boundary conditoins
Element Stiffness Matrix
Element Stiffness Matrix
Assemble Element Stiffness Matrix to
Globe Matrix
Boundary Conditions
Calculate Ax=b
Calculate Node Force
Internal force
Internal Force
Exercise
Node x y
1 4 0
2 4 3
3 0 0
4 0 3
Elem P1 P2
1 1 2
2 2 3
3 4 2
D=[0
0
D2x
D2y
0
0
0
0]
2 kN
Boundary
Conditions:
Node 1, 3 ,4 are
F=[F1x
fixed:
F1y
D=[0
0
-0.025 Q=[Q1x
0
D2x Q1y
F3x
D2y 2
F3y
0 0
F4x
0 Q3x
F4y]
0 Q3y
0]; Q4x
Q4y]
Flowchart for Exercise
clc;clear; close all;
AE=8e6;
% Determine the node coordinates % Determine the truss element
N=[]; E=[];
%The definitation of loads
F=zeros(2*NN,1);
D=zeros(2*NN,1);
D(2) = -0.025;
U(2) = -0.025;
syms U3 U4
U =[0 -0.025 U3 U4 0 0 0 0].';
%Assembling the global matrix
K=zeros(2*NN,2*NN);
for n=1:NE
L = sqrt((N(E(n,1),1)-N(E(n,2),1))^2 + (N(E(n,1),2)-N(E(n,2),2))^2);
lx = (N(E(n,2),1) - N(E(n,1),1))/L;
ly = (N(E(n,2),2) - N(E(n,1),2))/L;
%element matrix
Ke = AE*[ lx*lx lx*ly -lx*lx -lx*ly
lx*ly ly*ly -lx*ly -ly*ly
-lx*lx -lx*ly lx*lx lx*ly
-lx*ly -ly*ly lx*ly ly*ly]/L;
i = E(n,1); j = E(n,2);
K([2*i-1 2*i 2*j-1 2*j],[2*i-1 2*i 2*j-1 2*j]) = K([2*i-1 2*i 2*j-1 2*j],[2*i-1 2*i 2*j-1 2*j]) + Ke;
end
Flowchart for Exercise
T=K*U;
Q=solve(T(3),T(4));
D(3) = subs(Q.U3,'x',1);
D(4) = subs(Q.U4,'x',1);
F=K*D;
%Calculate internal force
for n=1:NE
L = sqrt((N(E(n,1),1)-N(E(n,2),1))^2 + (N(E(n,1),2)-
N(E(n,2),2))^2);
lx = (N(E(n,2),1) - N(E(n,1),1))/L;
ly = (N(E(n,2),2) - N(E(n,1),2))/L;
P1 = E(n,1); P2 = E(n,2);
q(n) = AE*[-lx -ly lx ly]*D([2*P1-1 2*P1 2*P2-1 2*P2])/L;
end
for i=1:2*NN
fprintf(['U(%i) = %5.3f\n'],i,D(i));
end
for i=1:2*NN
fprintf(['F(%i) = %5.3f\n'],i,F(i));
end
for i=1:NE
fprintf(['q(%i) = %5.3f\n'],i,q(i));
end
Go to Computer Lab
Develop a MATLAB code for the Exercise!
Verify Homework in Strand7
Create Beam Set BCs Apply load
Create Node
Set up Young’s Modulus Set up area of cross section
Linear Static Analysis
Log File
Model
Results
Output Strand7 Model in Txt Format
/ ______________________________________________________________________________
/ Strand7 MODEL EXCHANGE FILE
/ TIMESTAMP: 10:50:08 am, 18 August 2014
/ ______________________________________________________________________________
/ MODEL INFORMATION
FileFormat Strand7.2.4.4
ModelName "HOMEWORK"
Title ""
Project ""
Author ""
Reference ""
Comments ""
/ ______________________________________________________________________________
/ UNITS
LengthUnit m
MassUnit kg
EnergyUnit J
PressureUnit Pa
ForceUnit N
TemperatureUnit K
/ ______________________________________________________________________________
/ GROUP DEFINITIONS
Group 1 16711680 "\\Model"
/ ______________________________________________________________________________
/ FREEDOM CASE DEFINITIONS
FreedomCase 1 0 1 "Freedom Case 1"
/ ______________________________________________________________________________
/ LOAD CASE DEFINITIONS
LoadCase 1 0 "Load Case 1"
LCInclude 3
/ ______________________________________________________________________________
/ COORDINATE SYSTEM DEFINITIONS
CoordSys 1 "Global XYZ" GlobalXYZ
/ ______________________________________________________________________________
/ NODE COORDINATES
Node 1 0 4.00000000000000E+0 0.00000000000000E+0 0.00000000000000E+0
Node 2 0 4.00000000000000E+0 3.00000000000000E+0 0.00000000000000E+0
Node 3 0 0.00000000000000E+0 0.00000000000000E+0 0.00000000000000E+0
Node 4 0 0.00000000000000E+0 3.00000000000000E+0 0.00000000000000E+0
/ ______________________________________________________________________________
/ BEAM ELEMENTS
Beam 1 0 1 1 1 2
Beam 2 0 1 1 2 3
Beam 3 0 1 1 4 2
/ ______________________________________________________________________________
/ NODE RESTRAINTS (ROTATION AS RADIAN)
/ Freedom Case 1
NdFreedom 1 1 1 DX
NdFreedom 1 3 1 DX DY
NdFreedom 1 4 1 DX DY
/ ______________________________________________________________________________
/ NODE FORCES
/ Load Case 1
NdForce 1 2 -1.92000000000000E+4 -1.44000000000000E+4 0.00000000000000E+0
/ BEAM PROPERTIES
TrussProp 1 16737843 "Beam Property 1"
MaterialName "Unknown Material - Modified"
Modulus 8.00000000000000E+10
UsePoisson TRUE
InstantAlpha FALSE
Area 1.00000000000000E-4
MomentJ 1.40700000000000E-9
SectionType SolidRect
B 1.00000000000000E-2
D 1.00000000000000E-2
CT FALSE
IncludeTorsion FALSE
NonLinType Elasticplastic
Hardening Isotropic
________________
/ LINEAR STATIC SOLVER DATA
LoadFreedomSetLSA 1 ON
1
/ LINEAR BUCKLING SOLVER DATA
BuckNumModes 4
BuckShift 0.00000000000000E+0
/ LOAD INFLUENCE SOLVER DATA
LoadFreedomSetLIA 1 ON
1
/ GENERAL SOLVER DATA
SolverTempDependence None
SolverLoadCaseTempDependence 0
SolverActiveStage 0
SturmCheck FALSE
SolverFreedomCase 1
A Matlab Code for 2D Transmitter tower
Flow Chart
clc;clear;close all;
% read node and element information
N = textread('node.txt'); E = textread('elem.txt');
FN = [11 12 17 18 23 24]; F(2*FN-1) = -920*sin(20/180*pi); F(2*FN) = -920*cos(20/180*pi);
%Assembling the global matrix
K=zeros(2*NN,2*NN);
for n=1:NE
L = sqrt((N(E(n,1),1)-N(E(n,2),1))^2 + (N(E(n,1),2)-N(E(n,2),2))^2);
i = E(n,1); j = E(n,2);
K([2*i-1 2*i 2*j-1 2*j],[2*i-1 2*i 2*j-1 2*j]) = K([2*i-1 2*i 2*j-1 2*j],[2*i-1 2*i 2*j-1 2*j]) + Ke;
end
fixeddofs = [1 2 3 4]; % Fix end supports with zero deformation
alldofs = [1:2*NN];
freedofs = setdiff(alldofs,fixeddofs);
U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:);
U(fixeddofs,:)= 0;
for n=1:NE
P1 = E(n,1); P2 = E(n,2);
P(n) = A0*E0*[-lx -ly lx ly]*U([2*P1-1 2*P1 2*P2-1 2*P2])/L;
end
%Calculate the stress and strain
stress = P/A0;
strain = stress/E0;
Displacement
Strain Distribution
Result Comparison Between Matlab and
Strand 7
Strand7 Matlab
1.000000000000000 0 0
2.000000000000000 0 0
3.000000000000000 -0.056705187259539 -0.057696826067855
4.000000000000000 -0.020358528001983 0.012559782244823
5.000000000000000 -0.045999647171697 -0.172733192404131
6.000000000000000 -0.053109054594353 0.030202227373512
7.000000000000000 -0.263157875158521 -0.237774833979667
8.000000000000000 -0.251268324026876 0.058836063548212
9.000000000000000 -0.521945256702415 -0.290906572903267
10.000000000000000 -0.530747504449517 0.075960932937291
11.000000000000000 -0.693619429698391 -1.008756922796817
12.000000000000000 -0.681789204532035 0.521922821794403
13.000000000000000 -0.848748038328659 -0.332094346373710
14.000000000000000 -0.825251678313140 0.084272133505326
15.000000000000000 -1.165271291297121 -0.360053496810145
16.000000000000000 -1.183802187560931 0.087538142318972
17.000000000000000 -1.358028554905666 -1.543541703201050
18.000000000000000 -1.352392698081489 0.676291639035159
19.000000000000000 -1.540139025658684 -0.378790293549355
20.000000000000000 -1.514214220162403 0.086885916978769
21.000000000000000 -1.877909261927200 -0.389693981600771
22.000000000000000 -1.894365603525534 0.083788518739874
23.000000000000000 -2.069597330694618 -1.180928615989611
24.000000000000000 -2.067850950761376 0.612043780492507
25.000000000000000 -2.254750932371863 -0.394174398394259
26.000000000000000 -2.236519903204394 0.081969337109746
Debug and Helps for MATLA Code
clc;clear;close all;
% read node and element information
N = textread('node.txt');
E = textread('elem.txt');
%the numbers of nodes and elements
NN = length(N);
NE = length(E);
%Young's Modulus and cross section area
E0 = 2.1e6;
A0 = 6.91;
%plot the element
hold on;
for i=1:NE
plot([N(E(i,1),1),N(E(i,2),1)],[N(E(i,1),2),N(E(i,2),2)],'k','linewidth',2);
text((N(E(i,1),1)+N(E(i,2),1))/2,(N(E(i,1),2)+N(E(i,2),2))/2,num2str(i),'color','b','fontsize',12);
end
%plot the nodes
for i=1:NN
if i<3
plot(N(i,1),N(i,2),'r^','markersize',12,'markerfacecolor','r'); %plot fixed supports
else
plot(N(i,1),N(i,2),'ko','markersize',8,'markerfacecolor','k');
end
t t(N(i 1) N(i 2) 2 t (i) ' l ' ' ' 'f t i ' 12)
