Fundamentals of FEM Laboratory No. 5

Creation of the element stiffness matrices in the local coordinate system

Ke1=A1 E1/L1 {{1,-1},{-1,1}};
Ke2=A2 E2/L2 {{1,-1},{-1,1}};
Ke3=A3 E3/L3 {{1,-1},{-1,1}};

Definition of the transformation matrices and their transposed matrices

T1={{c1,0},{s1,0},{0,c1},{0,s1}};
T1t=Transpose[T1];
T2={{c2,0},{s2,0},{0,c2},{0,s2}};
T2t=Transpose[T2];
T3={{c3,0},{s3,0},{0,c3},{0,s3}};
T3t=Transpose[T3];

Checking the orthogonality of the transformation matrices

T1t.T1
T2t.T2
T3t.T3

Transformation of stiffness matrices into the global system

Ke1g=T1.Ke1.T1t;
Ke2g=T2.Ke2.T2t;
Ke3g=T3.Ke3.T3t;

Definition of matrices containing zeros

K1=Table[0,{6},{6}];
K2=K1;K3=K1;

Expansion of the stiffness matrices in accordance with the DOFs

K1[[{1,2,3,4},{1,2,3,4}]]=K1[[{1,2,3,4},{1,2,3,4}]]+Ke1g;
K2[[{3,4,5,6},{3,4,5,6}]]=K2[[{3,4,5,6},{3,4,5,6}]]+Ke2g;
K3[[{1,2,5,6},{1,2,5,6}]]=K3[[{1,2,5,6},{1,2,5,6}]]+Ke3g;

Creation of the structural stiffness matrix

Ks=K1+K2+K3;

Definition of boundary conditions and the structural nodal displacement vector

v1=0;u2=0;v2=0;
U={u1,v1,u2,v2,u3,v3};

Definition of the structural load vector

F={0,Ay,Bx,By,F1,0};

Definition of the geometrical and material properties

a=0.6; b=0.4; A1=100 (10^(-6)); E1=210 (10^9); A2=120 (10^(-6)); E2=180 (10^9); A3=90 (10^(-6)); E3=195 (10^9); F1=100 (10^3);

Nodal coordinates and beam element lengths

x1=0; y1=0; x2=a; y2=0;x3=a;y3=b;
L1=(Sqrt[(((x2-x1))^2)+(((y2-y1))^2)]);
L2=(Sqrt[(((x3-x2))^2)+(((y3-y2))^2)]);
L3=(Sqrt[(((x3-x1))^2)+(((y3-y1))^2)]);

Display of the structure by polygonplot

L={{x1,y1},{x2,y2},{x3,y3}};
Poly=Graphics[{White,Opacity[0], EdgeForm[Directive[Thick,Blue]],Polygon[L]}]

Definition of the sine and cosine of the angles

s1=(y2-y1)/L1; s2=(y3-y2)/L2; s3=(y3-y1)/L3;
c1=(x2-x1)/L1; c2=(x3-x2)/L2; c3=(x3-x1)/L3;

The structural equilibrium equation

actdof={1,5,6};
fixdof={2,3,4};
Kc=Ks[[actdof,actdof]];
Fc=F[[actdof]];
U[[actdof]]=LinearSolve[Kc,Fc];
forces=Ks.U;
displacement=U;
 u1=displacement[[1]]
u3=displacement[[5]]
v3=displacement[[6]]
Ay=forces[[2]]
Bx=forces[[3]]
By=forces[[4]]

Definition of a scale factor

sf=5;

Definition of a polygon

Ldef=L+sf Partition[displacement,2];

Display the deformed shape

Defpoly=Graphics[{{White,Opacity[0], EdgeForm[Thick],Polygon[Ldef]},{Text[Style["Deformed shape of the TRUSS2D structure",15,Bold],{0.3,0.4}]}}]
Show[Poly,Defpoly]

Animation of the deformed shape

Animate[Graphics[{{White,Opacity[0], EdgeForm[Directive[Thick,Blue]],Polygon[L]},{White,Opacity[0],EdgeForm[Directive[Thick,Black]],Polygon[L+sf t Partition[displacement,2]]},{Text[Style["Deformed shape of the TRUSS2D structure",15,Bold],{0.3,0.45}]}},PlotRange->{{-0.05,0.7},{-0.05,0.5}}],{t,0,1},AnimationRunning->False]

Calculation of the beam diagrams, element displacement vector, force vector, stiffness matrix

ue1={u1,v1,u2,v2};
Fe1=Ke1g.ue1;
Fe1L=T1t.Fe1;
ue2={u2,v2,u3,v3};
Fe2=Ke2g.ue2;
Fe2L=T2t.Fe2;
ue3={u1,v1,u3,v3};
Fe3=Ke3g.ue3;
Fe3L=T3t.Fe3;

Reduction factor for the plot of the beam diagram

rf=1 (10^(-6));

Creation of the points of the beam diagram

Le1={{x1,y1},{x2,y2},{x2,y2+rf Fe1[[3]]},{x1,y1-rf Fe1[[1]]}};
Le2={{x2,y2},{x3,y3},{x3-rf Fe2[[4]],y3},{x2+rf Fe2[[2]],y2}};
Le3={{x1,y1},{x3,y3},{x3-rf Fe3[[4]],y3+rf Fe3[[3]]},{x1+rf Fe3[[2]],y1-rf Fe3[[1]]}};

Display of the beam diagram

Graphics[{{Gray,EdgeForm[Thick],Polygon[Le1]},{Gray,EdgeForm[Thick],Polygon[Le2]},{Gray,EdgeForm[Thick],Polygon[Le3]},{Text[Style["Beam diagram - normal force",15,Bold],{0.3,0.55}]}},Axes->True]