EN130: Mechanics of Structures
FEM Truss Analysis Program Listing.
Non-Graphics Version
> # Finte Element Truss
Analysis Program. Text-only version
#
Written for EN130 at Brown Univ.
# Adapted from AF Bower's 2-D
FEM code by JA Blume
# restart:
#
#********** GIVE THE NAME OF THE INPUTFILE IN THE LINE
BELOW ***********
#
Use frontslash/ in place of backslash\
#
Enclose the name in single backquotes`
# End the line with a
colon:
#
infilename:=`C:/Documents
and Settings/laptop306/Desktop/tutorial.inp.txt`:
#
#*************************** THEN HIT RETURN
***************************
#
#***********************************************************************
# load linear
algebra maple package
#
with(linalg):
#
# define an output file name
#
lifn:=length(infilename):
outfilename:=cat(substring(infilename,1..lifn-8),`.out.txt`):
#
# ================= Read data from the input file
==================
#
infile:=fopen(infilename,READ):
#
# Title:
#
runtitle := readline(infile):
#
# No. nodes
and nodal coordinates
#
nnode := fscanf(infile,`%d`)[1]:
coord := array(1..nnode,1..2):
for i from 1
to nnode do
coord[i,1] := fscanf(infile,`%f`)[1]:
coord[i,2] := fscanf(infile,`%f`)[1]:
od:
#
# No.
elements connectivity, and material data.
#
connect[i,1] and connect[i,2] hold the two joint numbers for the
# ith member. connect[i,3] holds EA for that member
#
nelem := fscanf(infile,`%d`)[1]:
connect :=
array(1..nelem,1..3):
for i from
1 to nelem do
for j
from 1 to 2 do
connect[i,j] := fscanf(infile,`%d`)[1]:
od:
connect[i,3] := fscanf(infile,`%e`)[1]:
od:
#
# No. degrees
of freedom with prescribed displacements; prescribed
#
displacements
#
nfix := fscanf(infile,`%d`)[1]:
fixnodes := array(1..nfix,1..3):
for i from
1 to nfix do
#
# fixnodes[i,1] holds the
number of the ith fixed node
# fixnodes[i,2] holds the
direction (1 or 2)
# fixnodes[i,3] holds the valus of the prescribed displacement
#
fixnodes[i,1] := fscanf(infile,`%d`)[1]:
fixnodes[i,2] := fscanf(infile,`%d`)[1]:
fixnodes[i,3] := fscanf(infile,`%e`)[1]:
od:
#
# No. loaded
nodes, with the loads
#
nload := fscanf(infile,`%d`)[1]:
loads :=
array(1..nload,1..3):
#
# Node number
of ith loaded node is in loads[i,1]
# F1 and F2
for the ith loaded node are in loads[i,2] and
loads[i,3]
#
for i
from 1 to nload do
loads[i,1] := fscanf(infile,`%d`)[1]:
loads[i,2] := fscanf(infile,`%e`)[1]:
loads[i,3] := fscanf(infile,`%e`)[1]:
od:
fclose(infile):
# Print input data to the screen:
print(`number of nodes and elements:`,nnode,
nelem);
nodecoords:=array(1..nnode,1..3,
[seq([nn,coord[nn,1],coord[nn,2]],nn=1..nnode)]):
print(`coords: [node #, x,y]:`,nodecoords);
elcon:=array(1..nelem,1..4, [seq([nn,connect[nn,1],connect[nn,2],connect[nn,3]],nn=1..nelem)]):
print(`Connectivity: [el.#, node 1, node 2, EA]:`,elcon);
print(`[fixed node #, dof,
value]:`,fixnodes);
print(`[loaded node #, Px, Py]:`,loads);
#
> #
#================ ELEMENT STIFFNESS MATRIX
===================
#
elstif := proc (xa,ya,xb,yb,EA)
local elstif, nx, ny,
length;
#
# Define the
element length and direction vector
#
length := sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb));
#
nx := (xb-xa)/length;
ny := (yb-ya)/length;
#
# Define the
element stiffness
#
elstif := evalm((EA/length)*array(1..4,1..4,
[[nx*nx,nx*ny,-nx*nx,-nx*ny],
[ny*nx,ny*ny,-ny*nx,-ny*ny],
[-nx*nx,-nx*ny,nx*nx,nx*ny],
[-ny*nx,-ny*ny,ny*nx,ny*ny]]));
end:
#
#===========Assemble the global stiffness
matrix==============
#
Stif :=
array(1..2*nnode,1..2*nnode):
for i from 1 to 2*nnode do
for j from
1 to 2*nnode do
Stif[i,j] := 0.0:
od:
od:
#
# Loop over
all the elements
#
for lmn from 1 to nelem do
#
# Set up the stiffness
for the current element
#
a :=
connect[lmn,1]:
b :=
connect[lmn,2]:
EA :=
connect[lmn,3]:
k :=
array(1..4,1..4):
k:=elstif(coord[a,1],coord[a,2],coord[b,1],coord[b,2],EA):
#
# Add the
current element stiffness to the global stiffness
#
for i from 1
to 2 do
for ii
from 1 to 2 do
for j
from 1 to 2 do
for jj from 1 to 2 do
rw := 2*(connect[lmn,i]-1)+ii:
cl := 2*(connect[lmn,j]-1)+jj:
Stif[rw,cl] := Stif[rw,cl] +
k[2*(i-1)+ii,2*(j-1)+jj]:
od:
od:
od:
od:
od:
#
# ==================== Assemble global residual vector
============
#
# Define the
residual
#
resid := array(1..2*nnode):
for i from 1 to 2*nnode do
resid[i] := 0.:
od:
#
# Loop loaded noads
#
for i from 1
to nload do
a :=
loads[i,1]:
f1 :=
loads[i,2]:
f2 :=
loads[i,3]:
resid[2*a-1] := resid[2*a-1] +
f1:
resid[2*a] := resid[2*a] + f2:
od:
#
# Modify the global
stiffness and residual to include constraints
#
for i from 1
to nfix do
rw := 2*fixnodes[i,1]+fixnodes[i,2]-2:
for j
from 1 to 2*nnode do
resid[j]
:= resid[j] - Stif[j,rw]*fixnodes[i,3]:
Stif[rw,j] := 0:
Stif[j,rw] := 0:
od;
Stif[rw,rw] := 1.0:
resid[rw] := fixnodes[i,3]:
od:
#
# ================== Solve the FEM equations
===================
#
u := linsolve(Stif,resid):
#
# Define a
procedure to calculate the element strains
#
elstrain := proc (xa,ya,xb,yb,uax,uay,ubx,uby)
local elstrain, nx, ny,
length, uel;
#
# length and
direction vector for the element
#
length := sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb));
#
nx := (xb-xa)/length;
ny := (yb-ya)/length;
#
# Element
displacement vector
#
uel := array(1..4,[uax,uay,ubx,uby]);
#
# Element
strains
#
elstrain := ((uax-ubx)*nx+(uay-uby)*ny)/length;
end:
#
# compute
element forces and strains:
#
for lmn from 1 to nelem do
a :=
connect[lmn,1]:
b := connect[lmn,2]:
EA :=
connect[lmn,3]:
xa := coord[a,1]:
ya := coord[a,2]:
xb := coord[b,1]:
yb := coord[b,2]:
uxa := u[2*a-1]:
uya := u[2*a]:
uxb := u[2*b-1]:
uyb := u[2*b]:
Strain[lmn] := elstrain(xa,ya,xb,yb,uxa,uya,uxb,uyb):
Force[lmn]:=EA*Strain[lmn];
od:
#
# Print nodal
displacements, element strains and forces to a file
#
outfile:=fopen(outfilename,WRITE):
fprintf(outfile, `%s\n`,runtitle):
fprintf(outfile,`%s %3d %3d\n` ,`Number of Nodes and Elements:`,nnode, nelem):
fprintf(outfile,`%s\n`,`Node#
x1 x2`):
for i from 1
to nnode do
fprintf(outfile,`%5d %9.4f
%9.4f\n`,
i,coord[i,1],coord[i,2]):
od:
fprintf(outfile,`%s\n`,`El#
Node 1 Node 2 EA`):
for i from 1
to nelem do
fprintf(outfile,`%4d %10d
%10d %9.4f\n`,
i,connect[i,1],connect[i,2],connect[i,3]):
od:
fprintf(outfile,`%s\n`,`Fixed Node#
DOF Value`):
for i from 1
to nfix do
fprintf(outfile,`%9d %12d %10.4f\n`,
fixnodes[i,1],fixnodes[i,2],fixnodes[i,3]):
od:
fprintf(outfile,`%s\n`,`Loaded Node# P1
P2`):
for i from 1
to nload do
fprintf(outfile,`%9d %13.4f %9.4f\n`,
loads[i,1],loads[i,2],loads[i,3]):
od:
fprintf(outfile,`%s\n`,` `):
fprintf(outfile,`%s\n`,`*********************************`):
fprintf(outfile,`%s\n`,` `):
#
fprintf(outfile,`%s\n`,`Nodal Displacements:`):
fprintf(outfile,`%s\n`,` Node
u1 u2`):
for i from 1
to nnode do
fprintf(outfile,`%3d %8.4f %8.4f\n`,i,u[2*i-1],u[2*i]):
od:
fprintf(outfile,`%s\n`,`Strains and Forces:`):
fprintf(outfile,`%s\n`,`Element
Strain Force`):
for i from 1
to nelem do
fprintf(outfile,`%3d %9.4f %9.4f\n`,i,Strain[i],Force[i]):
od:
fclose(outfile):
print(`Finished!`);