Laboratory 3

Initialization

Geometrical and material properties (in SI units)

p=750;
L=1.2;
d=0.046;
Ex=200 (10^9);
s=2 (10^5);
Iz=(d^4) Pi/64;

Base shape function

ω=(x^2);

Approximation number

Nappr=4;

Calculation

Vector of basis polynomials

Table[{(x^(i-1))},{i,1,5}]

B=ω Table[{(x^(i-1))},{i,1,Nappr+1}]

MatrixForm[B]

Transpose of B

Bt=Transpose[B]

MatrixForm[Bt]

(2^nd) derivative of B

ddB=D[B,{x,2}];
ddBt=Transpose[ddB];

MatrixForm[ddB]

ddB.ddBt//MatrixForm

Construction of the stiffness matrix and the vector of external forces

S=Integrate[Iz Ex ddB.ddBt,{x,0,L}]+s (B.Bt)/.x->L;
S//MatrixForm

Q=Integrate[-p B,{x,0,L}];
Q//MatrixForm

Calculation of the vector of unknown coefficients

“or use A=LinearSolve[S,Q]; instead”;

A=Inverse[S].Q;
A//MatrixForm

Transpose[B].A

Deflection and angle of rotation functions

v[x_]=Transpose[B].A;
ϕ[x_]=D[v[x],x];

Plot[{v[x],ϕ[x]},{x,0,L},PlotLegends->"Expressions"]

Calculate the deflection and the rotation at cross section B

vb=v[L][[1]]
ϕb=ϕ[L][[1]]

Calculation for a series of Nappr

Maximum value of the approximation number

Nmax=6;

Vb=ConstantArray[0,Nmax];
Φb=ConstantArray[0,Nmax];
V=ConstantArray[0,Nmax];
Φ=ConstantArray[0,Nmax];
Do[
Nappr=j;
B=ω Table[{(x^(i-1))},{i,1,Nappr+1}];
Bt=Transpose[B];
ddB=D[B,{x,2}];
ddBt=Transpose[ddB];
S=Integrate[Iz Ex ddB.ddBt,{x,0,L}]+s (B.Bt)/.x->L;
Q=Integrate[-p B,{x,0,L}];
A=Inverse[S].Q;
"or use A=LinearSolve[S,Q]; instead"; 
v[x_]=Transpose[B].A;
ϕ[x_]=D[v[x],x];
"storing the functions in a vector";
V[[j]]=v[x][[1]];
Φ[[j]]=ϕ[x][[1]];
"storing in vectors the deflection and rotation values at section B ";
Vb[[j]]=v[L][[1]];
Φb[[j]]=ϕ[L][[1]],{j,Nmax}
]

Φb

Vb

V//Column

Φ//Column

Plotting the results

Convergence of deflection and rotation values at section B

ListPlot[ArrayFlatten[Vb,1],PlotRange->{{0,Nmax},{0,-0.0013}},PlotStyle->Red,PlotMarkers->{Automatic,Medium},AxesOrigin->{0,0},AxesLabel->{Style["(Subscript[N, appr])",14],Style["v(L) [m]",14]},AxesStyle->10]

ListPlot[ArrayFlatten[Φb,1],PlotRange->{{0,Nmax},{0,-0.001}},PlotStyle->Red,PlotMarkers->{Automatic,Medium},AxesOrigin->{0,0},AxesLabel->{Style["(Subscript[N, appr])",14],Style["φ(L) [rad]",14]},AxesStyle->10]

Convergence of deflection and angle of rotation functions

PlotV=Table[Plot[V[[i,1]],{x,0,L},PlotRange->Full,AxesOrigin->{0,0},AxesLabel->{Style["x [m]",14],Style["v(x) [m]",14]},AxesStyle->10],{i,Nmax}];
PlotΦ=Table[Plot[Φ[[i,1]],{x,0,L},PlotRange->Full,AxesOrigin->{0,0},AxesLabel->{Style["x [m]",14],Style["φ(x) [rad]",14]},AxesStyle->10],{i,Nmax}];

Manipulate[(i^2),{i,1,5,1}]

Manipulate[PlotV[[i]],{{i,1,"(Subscript[N, appr])"},1,Nmax,1}]

Manipulate[PlotΦ[[i]],{i,1,Nmax,1}]

Calculation of bending moment and shear force

Bending moment functions

Mb=Iz Ex D[V,{x,2}]

Shear force functions

Vs=-D[Mb,x];
Vs//Column

PlotMb=Table[Plot[Mb[[i,1]],{x,0,L},PlotRange->Full,AxesOrigin->{0,0},AxesLabel->{Style["x [m]",14],Style["Mb(x) [Nm]",14]},AxesStyle->10],{i,Nmax}];
PlotVs=Table[Plot[Vs[[i,1]],{x,0,L},PlotRange->Full,AxesOrigin->{0,0},AxesLabel->{Style["x [m]",14],Style["V(x) [N]",14]},AxesStyle->10],{i,Nmax}];

Manipulate[PlotMb[[i]],{{i,1,"(Subscript[N, appr])"},1,Nmax,1}]
Manipulate[PlotVs[[i]],{{i,1,"(Subscript[N, appr])"},1,Nmax,1}]