VEM L5

Megjegyzések:
- lehet, hogy az elmozdulásra jó pontos eredményt kapunk, de a hajlítónyomaték ( amiből aztán feszültség számítható ) nagyon eltér a két deriválás miatt;
- a megoszló terheléssel terhelt elemen az Mh lineárisan közelít;
- a csomóponti elmozdulásokból a formafüggvények segítségével számítható az elemen belüli elmozdulásmező;

Adatok:

data={a->0.8,b->0.4,d2->0.02,d1->2d2,E1->50 (10^9),ν->0.3,F1->2500,M1->-500,p0->-5000,E2->4E1,I1->(((d1^4)π)/64),I2->(((d2^4)π)/64)};

Referencia-megoldás: Rugalmas szál differenciálegyenlete

Hajlítónyomatéki függvények:

Mh1=-FA*x+MA-(p0/2)(x^2);
Mh2=-FA*x+MA-(p0*a)(x-(a/2))-F1(x-a);

Rugalmas szál differenciálegyenletei:

rugszal1=D[v1[x],{x,2}]==-(Mh1/(I1*E1));
rugszal2=D[v2[x],{x,2}]==-(Mh2/(I2*E2));

Peremfeltételek és illesztési feltételek:

pf1=v1[0]==0;
pf2=v1'[0]==0;
pf3=v1[a]==v2[a];
pf4=v1'[a]==v2'[a];

Megoldás:

megoldas=FullSimplify[DSolve[{rugszal1,rugszal2,pf1,pf2,pf3,pf4},{v1,v2},x]][[1]]

vx1=FullSimplify[v1[x]/.megoldas]
vx2=FullSimplify[v2[x]/.megoldas]

Reakciók

reakciok=Solve[((v2[a+b]/.megoldas)==0)&&(FA+p0*a+F1+FB==0)&&((Mh2/.{x->(a+b)})==-M1),{FA,FB,MA}][[1]]

FullSimplify[FA/.reakciok]//.data
FullSimplify[MA/.reakciok]//.data
FullSimplify[FB/.reakciok]//.data

Ezek visszahelyettesítése

vx1=FullSimplify[vx1/.reakciok]
vx2=FullSimplify[vx2/.reakciok]

vx1//.data//Expand

vx2//.data//Expand

v1[x_]=vx1//.data;
v2[x_]=vx2//.data;

Show[{
Plot[v1[x],{x,0,a/.data},PlotStyle->Red],
Plot[v2[x],{x,a/.data,a+b/.data},PlotStyle->Green]
},PlotRange->All]

A szögelfordulás függvények ábrázolása:

Show[{
Plot[v1'[x],{x,0,a/.data},PlotStyle->Red],
Plot[v2'[x],{x,a/.data,a+b/.data},PlotStyle->Green]
},PlotRange->All]

Hajlítónyomatéki függvények

Show[{
Plot[Mh1/.reakciok//.data,{x,0,a/.data},PlotStyle->Red,Filling->0,FillingStyle->Directive[{Opacity[0.2],Red}]],
Plot[Mh2/.reakciok//.data,{x,a/.data,a+b/.data},PlotStyle->Green,Filling->0,FillingStyle->Directive[{Opacity[0.2],Green}]]
},PlotRange->All]

VEM megoldás - n+m elemmel a két szakaszon

Elemek merevségi mátrixa

KE[L_,Iz_,E_]:=((Iz E)/(L^3))(({{12, 6*L, -12, 6*L}, {6*L, 4*(L^2), -6*L, 2*(L^2)}, {-12, -6*L, 12, -6*L}, {6*L, 2*(L^2), -6*L, 4*(L^2)}}))

Mátrix vagy Lista megtoldása balról, jobbról nullákkal

ArrayPad[KE[L,Iz,E1],{2,4}]//MatrixForm

Globális merevségi mx

makeK[n_,m_]:=Module[{L1,L2},
L1=a/n;
L2=b/m;
Plus@@(Table[ArrayPad[KE[L1,I1,E1],{2(i-1),2(n+m-i)}],{i,1,n}]~Join~Table[ArrayPad[KE[L2,I2,E2],{2(i-1),2(n+m-i)}],{i,n+1,n+m}])
]
makeK[1,1]//MatrixForm

Globális tehervektor (erre a feladatra)

FE[L_,p_]:={p L/2,p (L^2)/12,p L/2,-p (L^2)/12}

makeF[n_,m_]:=Module[{L1,L2,F},
L1=a/n;
L2=b/m;
F=Total@Table[ArrayPad[FE[L1,p0],{2(i-1),2(n+m-i)}],{i,1,n}];
F[[1]]+=FA;
F[[2]]+=MA;
F[[-2]]+=FB;
F[[-1]]+=M1;
F[[2(n+1)-1]]+=F1;
F
]
makeF[1,1]//MatrixForm

Egyenletrendszer kondenzációja, megoldása

Utolsó előtti sor elhagyása

(makeK[1,1])[[{1,2,3,4,6},{1,2,3,4,6}]]//MatrixForm

A megkötött szabadsági fokok: első, második és utolsó előtti 2(n+m+1)-1

solve[n_Integer,m_Integer]:=Module[{K,F,Kcond,Fcond,constrainedDOF,freeDOF,Ucond,U,Fsol,xi,vi,φi,Mi},
(* Teljes globális merevségi mx és terhelés vektor *)
K=makeK[n,m];
F=makeF[n,m];
(* Kötött szabadsági fokok *)
constrainedDOF={1,2,2(n+m+1)-1};
(* Szabad szabadsági fokok *)
freeDOF=Complement[Array[#&,{2(n+m+1)}],constrainedDOF];
Kcond=K[[freeDOF,freeDOF]];
Fcond=F[[freeDOF]];
(* KU=F megoldása *)
Ucond=LinearSolve[Kcond//.data,Fcond//.data];
(* Teljes csp-i elmozdulás vektor, 0-ákkal kitöltve *)
U=Array[0&,{2(n+m+1)}];
(* Szabad DOF helyekre a megoldás visszaírása *)
Table[U[[freeDOF[[i]]]]+=Ucond[[i]],{i,1,Length@freeDOF}];
(* Reakcióerők kiszámítása *)
Fsol=First@Solve@Thread[((K.U)[[constrainedDOF]]//.data)==F[[constrainedDOF]]//.data];
xi=Table[(i-1) a/n/.data,{i,1,n+1}]~Join~Table[a+i b/m/.data,{i,1,m}];
vi=U[[Table[i,{i,1,2(n+m+1),2}]]];
φi=U[[Table[i,{i,2,2(n+m+1),2}]]];
F=F/.Fsol//.data;
Mi=F[[Table[i,{i,2,2(n+m+1),2}]]];
{Transpose@{xi,vi},Transpose@{xi,φi},
Fsol}
]
teszt=solve[3,2]

A visszatérési lista első eleme az elmozdulásokat, a második az elfordulásokat tartalmazza,
a 3. pedig a reakcióerőkre vonatkozó megoldást

Összehasonlítás a referenciával:

Manipulate[
Module[{vemsol},
vemsol=solve[m,n];
Row[{
(* Lehajlás *)
Show[{
Plot[v1[x],{x,0,a/.data},PlotStyle->Red],
Plot[v2[x],{x,a/.data,a+b/.data},PlotStyle->Red],
ListPlot[vemsol[[1]],Joined->True],ListPlot[vemsol[[1]]]
},PlotRange->All,AxesLabel->{"x","v(x)"},ImageSize->350],
(* Szögelfordulás *)
Show[{
Plot[v1'[x],{x,0,a/.data},PlotStyle->Red],
Plot[v2'[x],{x,a/.data,a+b/.data},PlotStyle->Red],
ListPlot[vemsol[[2]],Joined->True],ListPlot[vemsol[[2]]]
},PlotRange->All,AxesLabel->{"x","v(x)"},ImageSize->350]
}]
],
{{m,2},1,10,1,Appearance->"Labeled"},
{{n,2},1,10,1,Appearance->"Labeled"}]

Formafüggvények: v(ξ)=(Subscript[N, 1])(ξ)(Subscript[u, 1])+(Subscript[N, 2])(ξ)(Subscript[u, 2])+(Subscript[N, 3])(ξ)(Subscript[u, 3])+(Subscript[N, 4])(ξ)(Subscript[u, 4])

NU1[ξ_,L_]:=(1/4)(((1-ξ))^2)(2+ξ);
NU2[ξ_,L_]:=(1/8)L(((1-ξ))^2)(1+ξ);
NU3[ξ_,L_]:=(1/4)(((1+ξ))^2)(2-ξ);
NU4[ξ_,L_]:=-(1/8)L(((1+ξ))^2)(1-ξ);

Hajlítónyomatéki függvény: (Subscript[Mh, VEM])(x)=-I E v''(x)

HF: Formafüggvények használatával teljes v(x) előállítása, majd abból (Subscript[Mh, VEM])(x) kiszámítása.