A finite element program for rectilinear beams
Author
Maria Girardi
Title
A finite element program for rectilinear beams
Description
The notebook contains a simple fe program to model rectilinear beams. Beam elements with cubic shape functions are implemented in the code, which can be easily modified to account for different distributed and concentrated loads. The code is for educational purposes only.
Category
Educational Materials
Keywords
beam element, cubic shape functions, finite elements, linear elastic solutions
URL
http://www.notebookarchive.org/2021-03-0hx8a52/
DOI
https://notebookarchive.org/2021-03-0hx8a52
Date Added
2021-03-01
Date Last Modified
2021-03-01
File Size
170.61 kilobytes
Supplements
Rights
Redistribution rights reserved
Download
Open in Wolfram Cloud
In[]:=
"Finite element program: beam elements with Hermite polynomials shape functions. Linear elastic solution and modal analysis. By Maria Girardi (maria.girardi@isti.cnr.it).For educational purposes only";
In[]:=
Off[General::spell1]
In[]:=
"Nodal coordinates (abscissa along the beam's axis)";
In[]:=
c:={0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.,3.2,3.4,3.6,3.8,4.,4.2,4.4,4.6,4.8,5.,5.2,5.4,5.6,5.8,6.}
In[]:=
"Number n of nodes";
In[]:=
n=Length[c]
Out[]=
31
In[]:=
Ltot=Extract[c,n]-Extract[c,1]
Out[]=
6.
In[]:=
"Length of the elements";
In[]:=
l[i_]:=Table[Extract[c,i+1]-Extract[c,i]]
In[]:=
"Number m of elements";
In[]:=
m:=n-1
In[]:=
m
Out[]=
30
In[]:=
"Width of the elements";
In[]:=
b:={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1}
In[]:=
Length[b]
Out[]=
30
In[]:=
"Height of the elements";
In[]:=
h:={0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4}
In[]:=
Length[h]
Out[]=
30
In[]:=
"Young's moduli";
In[]:=
Y:={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3}
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
9
10
In[]:=
Length[Y]
Out[]=
30
In[]:=
"mass density";
In[]:=
ρ:={1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800}
In[]:=
Length[ρ]
Out[]=
30
In[]:=
J[i_]:=Extract[b,i]12
3
(Extract[h,i])
In[]:=
J[1]
Out[]=
0.00533333
In[]:=
"distributed transverse load";
In[]:=
p[x_]=45000;
In[]:=
Plot[p[x],{x,0,Ltot}]
Out[]=
In[]:=
"Shape functions";
In[]:=
ψ
10
1
4
3
ξ
2y-l[i]
l[i]
ψ
20
1
4
3
ξ
2y-l[i]
l[i]
ψ
11
1
4
2
ξ
3
ξ
2y-l[i]
l[i]
ψ
21
1
4
2
ξ
3
ξ
2y-l[i]
l[i]
In[]:=
Plot[[9],{y,0,l[9]}]
ψ
20
Out[]=
In[]:=
"local load vector";
In[]:=
Clear[B]
In[]:=
B[i_]:=p[x]([i]/.y(x-Extract[c,i]))x,p[x]([i]/.y(x-Extract[c,i]))x,p[x]([i]/.y(x-Extract[c,i]))x,p[x]([i]/.y(x-Extract[c,i]))x;
Extract[c,i+1]
∫
Extract[c,i]
ψ
10
Extract[c,i+1]
∫
Extract[c,i]
ψ
11
Extract[c,i+1]
∫
Extract[c,i]
ψ
20
Extract[c,i+1]
∫
Extract[c,i]
ψ
21
In[]:=
B[16]
Out[]=
{4500.,1500.,4500.,-1500.}
In[]:=
B[18]
Out[]=
{4500.,1500.,4500.,-1500.}
In[]:=
p[0]
Out[]=
45000
In[]:=
F
In[]:=
MatrixForm[[11]]
F
Out[]//MatrixForm=
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
4500. |
1500. |
4500. |
-1500. |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
In[]:=
"Global load vector: assembly";
In[]:=
M[0]:=[0];
F
In[]:=
Table[M[s]=M[s-1]+[s],{s,1,m,1}];f=M[m];
F
In[]:=
MatrixForm[f]
Out[]//MatrixForm=
4500. |
1500. |
9000. |
-5.91172× -11 10 |
9000. |
-2.91038× -11 10 |
9000. |
8.14907× -10 10 |
9000. |
4.65661× -10 10 |
9000. |
9.31323× -10 10 |
9000. |
2.79397× -9 10 |
9000. |
-5.58794× -9 10 |
9000. |
-6.51926× -9 10 |
9000. |
1.86265× -9 10 |
9000. |
-2.6077× -8 10 |
9000. |
-3.35276× -8 10 |
9000. |
8.19564× -8 10 |
9000. |
9.68575× -8 10 |
9000. |
8.19564× -8 10 |
9000. |
-1.49012× -8 10 |
9000. |
1.11759× -7 10 |
9000. |
-1.49012× -7 10 |
9000. |
5.96046× -8 10 |
9000. |
-3.8743× -7 10 |
9000. |
1.49012× -8 10 |
9000. |
-5.36442× -7 10 |
9000. |
8.9407× -8 10 |
9000. |
-5.06639× -7 10 |
9000. |
-1.43051× -6 10 |
9000. |
-2.98023× -7 10 |
9000. |
-2.38419× -7 10 |
9000. |
1.78814× -7 10 |
9000. |
-1.13249× -6 10 |
9000. |
-8.9407× -7 10 |
4500. |
-1500. |
In[]:=
Dimensions[f]
Out[]=
{62}
In[]:=
"Global load vector: boundary conditions (simply supported beam)";
In[]:=
fv:=ffv=ReplacePart[fv,0,1];fv=ReplacePart[fv,0,2n-1];
In[]:=
MatrixForm[fv];
In[]:=
(*Massmatrix*)
In[]:=
V[i_]:=Extract[ρ,i]Extract[b,i]Extract[h,i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y,[i][i]y;
l[i]
∫
0
2
ψ
10
l[i]
∫
0
ψ
11
ψ
10
l[i]
∫
0
ψ
20
ψ
10
l[i]
∫
0
ψ
21
ψ
10
l[i]
∫
0
ψ
11
ψ
10
l[i]
∫
0
2
ψ
11
l[i]
∫
0
ψ
11
ψ
20
l[i]
∫
0
ψ
21
ψ
11
l[i]
∫
0
ψ
20
ψ
10
l[i]
∫
0
ψ
11
ψ
20
l[i]
∫
0
ψ
20
ψ
20
l[i]
∫
0
ψ
21
ψ
20
l[i]
∫
0
ψ
10
ψ
21
l[i]
∫
0
ψ
11
ψ
21
l[i]
∫
0
ψ
20
ψ
21
l[i]
∫
0
ψ
21
ψ
21
In[]:=
MatrixForm[V[1]]
Out[]//MatrixForm=
53.4857 | 15.0857 | 18.5143 | -8.91429 |
15.0857 | 5.48571 | 8.91429 | -4.11429 |
18.5143 | 8.91429 | 53.4857 | -15.0857 |
-8.91429 | -4.11429 | -15.0857 | 5.48571 |
In[]:=
"Global mass matrix: assembly";
In[]:=
V
In[]:=
S[0]:=[0];
V
In[]:=
Table[S[s]=S[s-1]+[s],{s,1,m,1}];mass=S[m];
V
In[]:=
MatrixForm[mass];
In[]:=
Dimensions[mass]
Out[]=
{62,62}
In[]:=
"Boundary conditions on the mass matrix: simply supported beam";
In[]:=
mv:=massDo[mv=ReplacePart[mv,If[j==1,1,0],{1,j}],{j,1,2n,1}];Do[mv=ReplacePart[mv,If[l==1,1,0],{l,1}],{l,1,2n,1}];Do[mv=ReplacePart[mv,If[j2n-1,1,0],{2n-1,j}],{j,1,2n,1}];Do[mv=ReplacePart[mv,If[l2n-1,1,0],{l,2n-1}],{l,1,2n,1}];
In[]:=
MatrixForm[mv];
In[]:=
Dimensions[mv]
Out[]=
{62,62}
In[]:=
"Element stiffness matrix";
In[]:=
ATel[i_]:=Extract[Y,i]J[i]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y,D[[i],{y,2}]D[[i],{y,2}]y;
l[i]
∫
0
2
D[[i],{y,2}]
ψ
10
l[i]
∫
0
ψ
11
ψ
10
l[i]
∫
0
ψ
20
ψ
10
l[i]
∫
0
ψ
21
ψ
10
l[i]
∫
0
ψ
11
ψ
10
l[i]
∫
0
2
D[[i],{y,2}]
ψ
11
l[i]
∫
0
ψ
11
ψ
20
l[i]
∫
0
ψ
21
ψ
11
l[i]
∫
0
ψ
10
ψ
20
l[i]
∫
0
ψ
11
ψ
20
l[i]
∫
0
ψ
20
ψ
20
l[i]
∫
0
ψ
21
ψ
20
l[i]
∫
0
ψ
10
ψ
21
l[i]
∫
0
ψ
11
ψ
21
l[i]
∫
0
ψ
20
ψ
21
l[i]
∫
0
ψ
21
ψ
21
In[]:=
"Global stiffness matrix: assembly"
Out[]=
Global stiffness matrix: assembly
In[]:=
Kel
Kel
Kel
In[]:=
"Global stiffness matrix: boundary conditions (simply supported beam)";
In[]:=
kel=Kel[m];kvel:=kelDo[kvel=ReplacePart[kvel,If[j==1,1,0],{1,j}],{j,1,2n,1}];Do[kvel=ReplacePart[kvel,If[l==1,1,0],{l,1}],{l,1,2n,1}];Do[kvel=ReplacePart[kvel,If[j2n-1,1,0],{2n-1,j}],{j,1,2n,1}];Do[kvel=ReplacePart[kvel,If[l2n-1,1,0],{l,2n-1}],{l,1,2n,1}];
In[]:=
"System solution";
In[]:=
U=Flatten[Table[,,{j,1,n,1}]];
u
j
du
j
In[]:=
Ured=Drop[Drop[U,1],{2n-2}]
Out[]=
{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,}
du
1
u
2
du
2
u
3
du
3
u
4
du
4
u
5
du
5
u
6
du
6
u
7
du
7
u
8
du
8
u
9
du
9
u
10
du
10
u
11
du
11
u
12
du
12
u
13
du
13
u
14
du
14
u
15
du
15
u
16
du
16
u
17
du
17
u
18
du
18
u
19
du
19
u
20
du
20
u
21
du
21
u
22
du
22
u
23
du
23
u
24
du
24
u
25
du
25
u
26
du
26
u
27
du
27
u
28
du
28
u
29
du
29
u
30
du
30
du
31
In[]:=
fvred=Drop[Drop[fv,1],{2n-2}]
Out[]=
{1500.,9000.,-5.91172×,9000.,-2.91038×,9000.,8.14907×,9000.,4.65661×,9000.,9.31323×,9000.,2.79397×,9000.,-5.58794×,9000.,-6.51926×,9000.,1.86265×,9000.,-2.6077×,9000.,-3.35276×,9000.,8.19564×,9000.,9.68575×,9000.,8.19564×,9000.,-1.49012×,9000.,1.11759×,9000.,-1.49012×,9000.,5.96046×,9000.,-3.8743×,9000.,1.49012×,9000.,-5.36442×,9000.,8.9407×,9000.,-5.06639×,9000.,-1.43051×,9000.,-2.98023×,9000.,-2.38419×,9000.,1.78814×,9000.,-1.13249×,9000.,-8.9407×,-1500.}
-11
10
-11
10
-10
10
-10
10
-10
10
-9
10
-9
10
-9
10
-9
10
-8
10
-8
10
-8
10
-8
10
-8
10
-8
10
-7
10
-7
10
-8
10
-7
10
-8
10
-7
10
-8
10
-7
10
-6
10
-7
10
-7
10
-7
10
-6
10
-7
10
In[]:=
kred=Drop[Drop[kvel,{1},{1}],{2n-2},{2n-2}];
In[]:=
Dimensions[kred]
Out[]=
{60,60}
In[]:=
MatrixForm[kred];
In[]:=
da=Flatten[N[Ured/.Solve[kred.Uredfvred,Ured]]]
Out[]=
{0.00253125,0.00505144,0.00251475,0.010038,0.00246675,0.0148989,0.0023895,0.019578,0.00228525,0.0240234,0.00215625,0.028188,0.00200475,0.0320289,0.001833,0.035508,0.00164325,0.0385914,0.00143775,0.04125,0.00121875,0.0434589,0.0009885,0.045198,0.00074925,0.0464514,0.00050325,0.047208,0.00025275,0.0474609,-1.6364×,0.047208,-0.00025275,0.0464514,-0.00050325,0.045198,-0.00074925,0.0434589,-0.0009885,0.04125,-0.00121875,0.0385914,-0.00143775,0.035508,-0.00164325,0.0320289,-0.001833,0.028188,-0.00200475,0.0240234,-0.00215625,0.019578,-0.00228525,0.0148989,-0.0023895,0.010038,-0.00246675,0.00505144,-0.00251475,-0.00253125}
-15
10
In[]:=
da=Insert[Insert[da,0,1],0,{2n-1}]
Out[]=
{0,0.00253125,0.00505144,0.00251475,0.010038,0.00246675,0.0148989,0.0023895,0.019578,0.00228525,0.0240234,0.00215625,0.028188,0.00200475,0.0320289,0.001833,0.035508,0.00164325,0.0385914,0.00143775,0.04125,0.00121875,0.0434589,0.0009885,0.045198,0.00074925,0.0464514,0.00050325,0.047208,0.00025275,0.0474609,-1.6364×,0.047208,-0.00025275,0.0464514,-0.00050325,0.045198,-0.00074925,0.0434589,-0.0009885,0.04125,-0.00121875,0.0385914,-0.00143775,0.035508,-0.00164325,0.0320289,-0.001833,0.028188,-0.00200475,0.0240234,-0.00215625,0.019578,-0.00228525,0.0148989,-0.0023895,0.010038,-0.00246675,0.00505144,-0.00251475,0,-0.00253125}
-15
10
In[]:=
N[5p[0]/(J[1]Extract[Y,1]384)]
4
Ltot
Out[]=
0.0474609
In[]:=
Extract[da,31]
Out[]=
0.0474609
In[]:=
Np[0](24J[1]Extract[Y,1])
3
Ltot
Out[]=
0.0253125
In[]:=
2Extract[da,2]/l[1]
Out[]=
0.0253125
In[]:=
deformedshape=Tablel,Extract[da,j],{j,1,2n,2}
j-1
2
j-1
2
Out[]=
{{0,0},{0.2,0.00505144},{0.4,0.010038},{0.6,0.0148989},{0.8,0.019578},{1.,0.0240234},{1.2,0.028188},{1.4,0.0320289},{1.6,0.035508},{1.8,0.0385914},{2.,0.04125},{2.2,0.0434589},{2.4,0.045198},{2.6,0.0464514},{2.8,0.047208},{3.,0.0474609},{3.2,0.047208},{3.4,0.0464514},{3.6,0.045198},{3.8,0.0434589},{4.,0.04125},{4.2,0.0385914},{4.4,0.035508},{4.6,0.0320289},{4.8,0.028188},{5.,0.0240234},{5.2,0.019578},{5.4,0.0148989},{5.6,0.010038},{5.8,0.00505144},{6.,0}}
In[]:=
ListPlot[deformedshape,PlotJoinedTrue]
Out[]=
In[]:=
"Elastic curvature";
In[]:=
chiel[i_,x_]:=-If[Extract[c,i]≤x≤Extract[c,i+1],((D[[i],{y,2}]Extract[da,(1+2(i-1))]+D[[i],{y,2}]Extract[da,(3+2(i-1))]+D[[i],{y,2}]Extract[da,(2+2(i-1))]+D[[i],{y,2}]Extract[da,(4+2(i-1))])/.y(x-Extract[c,i])),0]
ψ
10
ψ
20
ψ
11
ψ
21
In[]:=
Plot[{chiel[1,x],chiel[2,x],chiel[3,x],chiel[4,x],chiel[5,x],chiel[6,x],chiel[7,x],chiel[8,x],chiel[9,x],chiel[10,x],chiel[11,x],chiel[12,x],chiel[13,x],chiel[14,x],chiel[15,x],chiel[16,x],chiel[17,x],chiel[18,x],chiel[19,x],chiel[20,x]},{x,0,Ltot}]
Out[]=
In[]:=
chimax=
p[0]
2
Ltot
8J[1]Extract[Y,1]
Out[]=
0.0126562
In[]:=
chiel[15,3]
Out[]=
0.0126656
In[]:=
"---------------------Modal Analysis-------------------";
In[]:=
mred=Drop[Drop[mv,{1},{1}],{2n-2},{2n-2}];
In[]:=
MatrixForm[%];
In[]:=
Dimensions[mred]
Out[]=
{60,60}
In[]:=
"Natural frequencies of the beam";
In[]:=
N
Eigenvalues[{kred,mred}]
Out[]=
{187083.,186334.,184132.,180600.,175921.,170317.,164017.,157240.,150180.,143000.,135827.,128758.,121863.,115191.,108770.,102617.,96738.6,91135.4,85802.,80730.5,75910.7,71331.9,66982.8,62852.4,58930.5,55208.3,51679.9,48345.1,45222.9,42422.6,40824.8,35591.2,33214.9,30801.1,28462.2,26221.,24082.7,22047.4,20113.4,18278.6,16540.7,14897.4,13346.6,11886.5,10515.3,9231.75,8034.47,6922.45,5894.85,4950.94,4090.18,3312.13,2616.47,2002.96,1471.43,1021.77,653.913,367.821,163.475,40.8687}
In[]:=
freq=Sort[%/(2π)]
Out[]=
{6.50446,26.0179,58.5405,104.074,162.62,234.186,318.781,416.425,527.143,650.972,787.967,938.194,1101.74,1278.73,1469.28,1673.57,1891.79,2124.18,2371.,2632.54,2909.14,3201.15,3508.95,3832.88,4173.2,4529.9,4902.15,5286.31,5664.52,6497.47,6751.77,7197.45,7694.36,8225.1,8786.68,9379.08,10003.3,10660.6,11352.8,12081.6,12848.7,13655.8,14504.6,15396.4,16331.9,17311.2,18333.2,19395.2,20492.5,21617.5,22759.1,23902.,25025.5,26104.1,27106.8,27998.8,28743.3,29305.5,29656.,29775.2}
In[]:=
Dimensions[%]
Out[]=
{60}
In[]:=
fund=Extract[freq,1]
Out[]=
6.50446
In[]:=
"Fundamental frequency [Hz]";
In[]:=
fund_elastica=N
2
π
2
Ltot
Extract[Y,1]J[1]
Extract[ρ,1]Extract[b,1]Extract[h,1]
2π
Out[]=
6.50446
Cite this as: Maria Girardi, "A finite element program for rectilinear beams" from the Notebook Archive (2021), https://notebookarchive.org/2021-03-0hx8a52
Download
