Skręcanie swobodne dowolnego przekroju

Clear["Global`*"];Off[Inverse::luc];tu=TimeUsed[];

Dane wejściowe: Region oraz oczko siatki
Należy opatrzeć komentarzem (*...*) wszystkie prócz jednej definicji

(*DWUTEOWNIK ZWARTY*)
oczko=0.0003;h=1;b=0.8;tw=0.3;tf=0.3;
ℛ=BoundaryMeshRegion[{{-b/2,-h/2},{b/2,-h/2},{b/2,-h/2+tf},{tw/2,-h/2+tf},{tw/2,h/2-tf},{b/2,h/2-tf},{b/2,h/2},{-b/2,h/2},{-b/2,h/2-tf},{-tw/2,h/2-tf},{-tw/2,-h/2+tf},{-b/2,-h/2+tf}},Line[{1,2,3,4,5,6,7,8,9,10,11,12,1}]];
(*ELIPSA*)
(*a=2;b=1;oczko=0.001;
ℛ=BoundaryDiscretizeRegion[ImplicitRegion[((x^2)/(a^2))+((y^2)/(b^2))≤1,{x,y}]];*)
(*PROSTOKĄT*)
(*a0=3;b0=1;oczko=0.001;
ℛ=BoundaryMeshRegion[{{0,0},{a0,0},{a0,b0},{0,b0}},Line[{1,2,3,4,1}]];*)
(*TRÓJKĄT*)
(*b=1;a=((b(Sqrt[3]))/2);oczko=0.0001;
ℛ=BoundaryMeshRegion[{{0,0},{b,0},{b/2,a}},Line[{1,2,3,1}]];*)

siatka=TriangulateMesh[ℛ,MaxCellMeasure->oczko]
le=Length[MeshCells[siatka,2]](*ilość elementów*)

Równania dla pojedynczego ES:

Definicja trójkątnego ES: funkcje kształtu L1,L2,L3; Pole oraz a,b,c wczytane później dla operowania na prostszych wzorach

(*Δ=Abs[(1/2)Det[(({{x1, x2, x3}, {y1, y2, y3}, {1, 1, 1}}))]];*)
(*(({{a1, b1, c1}, {a2, b2, c2}, {a3, b3, c3}}))=(({{x2 y3-x3 y2, y2-y3, x3-x2}, {x3 y1-x1 y3, y3-y1, x1-x3}, {x1 y2-x2 y1, y1-y2, x2-x1}}));*)
(({{L1}, {L2}, {L3}}))=(({{a1+b1 x +c1 y}, {a2+b2 x+c2 y}, {a3+b3 x+c3 y}}))/(2Δ);

Definicja poszukiwanej funkcji naprężeń Prandtl’a i funkcjonału opisującego skręcanie swobodne. 

Równanie swbodonego skręcania: Δϕ=-2Gθ

Poniżej funkcjonał dla zagadnienia Poissona z innymi oznaczeniami:


(żródło: online coursebook for Advanced Variational Methods In Mechanics, Department of Aerospace Engineering Sciences
University of Colorado at Boulder, http://www.colorado.edu/engineering/CAS/courses.d/AVMM.d/)

Jeśli przyjąć k=1, s=-2Gθ, u=ϕ to mamy równanie i funkcjonał:
Δϕ=-2Gθ
Π=(1/2)∫[((((∂_x)ϕ))^2)+((((∂_y)ϕ))^2)]dA-2Gθ∫ϕdA 
	// korzystamy, że ∫((Subscript[L, 1])^a)((Subscript[L, 2])^b)((Subscript[L, 3])^c)dA=((a!b!c!)/((a+b+c+2)!))2Δ , więc ((2Δ)/(3!))= (Δ/3)  (Zienkiewicz, Vol.1, wzór 8.38)
	∫ϕdA=∫(L1 ϕ1+L2 ϕ2+L3 ϕ3)dA=(Δ/3)(ϕ1+ϕ2+ϕ3)
	//zuważamy, że (∂_x)ϕ oraz (∂_y)ϕ są stałe, więc (1/2)∫[((((∂_x)ϕ))^2)+((((∂_y)ϕ))^2)]dA=(1/2)Δ(((((∂_x)ϕ))^2)+((((∂_y)ϕ))^2))
Ostatecznie więc: 

Π=(1/2)Δ(((((∂_x)ϕ))^2)+((((∂_y)ϕ))^2))-(2/3)G θ Δ (ϕ1+ϕ2+ϕ3)

Φ=(({{ϕ1}, {ϕ2}, {ϕ3}}));ϕ=L1 ϕ1+L2 ϕ2+L3 ϕ3;
Π=(1/2)Δ(((((∂_x)ϕ))^2)+((((∂_y)ϕ))^2))-(2/3)G θ Δ (ϕ1+ϕ2+ϕ3)

Przekształcenie funkcjonału do  postaci macierzowej

Można łatwo wyprowadzić K i f tak, że Π=(1/2)Transpose[Φ].K.Φ-Transpose[f].Φ, gdziek K,f i dowód poniżej:

K=(1/(4Δ))((({{b1 b1, b1 b2, b1 b3}, {b2 b1, b2 b2, b2 b3}, {b3 b1, b3 b2, b3 b3}}))+(({{c1 c1, c1 c2, c1 c3}, {c2 c1, c2 c2, c2 c3}, {c3 c1, c3 c2, c3 c3}})));
f=(2/3) G Δ θ (({{1}, {1}, {1}}));
Simplify[Π-((1/2)Transpose[Φ].K.Φ-Transpose[f].Φ)]

Uzupełnienie definicji ES

Δ=Abs[(1/2)Det[(({{x1, x2, x3}, {y1, y2, y3}, {1, 1, 1}}))]];(({{a1, b1, c1}, {a2, b2, c2}, {a3, b3, c3}}))=(({{x2 y3-x3 y2, y2-y3, x3-x2}, {x3 y1-x1 y3, y3-y1, x1-x3}, {x1 y2-x2 y1, y1-y2, x2-x1}}));

Dla pojedynczego elementu musi być zatem spełnione równanie: K.Φ=f (co wynika ze zróżniczkowania funkcjonału po Φ)

K//MatrixForm

Równania dla całego obszaru:

Podział na elementy skończone prostokąta
Punkty[[1..lw]]		-zawieraja współrzędne wierzchołka
V[[1..le]]		-zawiera numery wierzchołków danego ES

brzegpts=MeshCoordinates[BoundaryMesh[siatka]];
brzeg={};
Punkty=MeshCoordinates[siatka];
lw=Length[Punkty];
V=Table[Flatten[Cases[MeshCells[siatka,2][[i]],pts_->pts]],{i,le}];
Do[	If[ MemberQ[brzegpts,Punkty[[i]]],AppendTo[brzeg,i]];,{i,lw}];

Sprawdzenie czy rzeczywiście wszystkie punkty na brzegu zostały uwzględnione do przypisania warunków brzegowych. Raczej potrzebne, ponieważ porównanie punktów jest przez współrzędne, na co być może wpłynąć jakiś błąd numeryczny

Length[brzegpts]==Length[brzeg]

Wygenerowanie macierzy iw,która zawiera informację o tym ile elementów schodzi się w danym wierzchołku-potrzebne przy ustaleniu wartości naprężcenia/odkształcenia w wierzchołku:

iw=Table[0,{i,lw}];
Do[k=V[[i,j]];iw[[k]]=iw[[k]]+1,{i,le},{j,3}];

Rysunek z numeracją węzłów

(*rysunek1={};
For[i=1,i≤le,	rysunek1=Append[rysunek1,Line[{Punkty[[V[[i,1]]]],Punkty[[V[[i,2]]]],Punkty[[V[[i,3]]]],Punkty[[V[[i,1]]]]}]];i++];
For [i=1,i≤ Length[Punkty],	rysunek1=Append[rysunek1,Text[i,Punkty[[i]]]];i++];
Graphics[rysunek1]*)

K i f dla każdego elementu:

Ktab=Table[K/.{G->1,θ->1,
x1-> Punkty[[V[[i,1]],1]],y1-> Punkty[[V[[i,1]],2]],
x2-> Punkty[[V[[i,2]],1]],y2-> Punkty[[V[[i,2]],2]],
x3-> Punkty[[V[[i,3]],1]],y3-> Punkty[[V[[i,3]],2]]
},{i,1,le}];
ftab=Table[f/.{G->1,θ->1,
x1-> Punkty[[V[[i,1]],1]],y1-> Punkty[[V[[i,1]],2]],
x2-> Punkty[[V[[i,2]],1]],y2-> Punkty[[V[[i,2]],2]],
x3-> Punkty[[V[[i,3]],1]],y3-> Punkty[[V[[i,3]],2]]
},{i,1,le}];

Agregacja do macierzy globalnych

Kglob=Table[0,{i,lw},{j,lw}];
fglob=Table[0,{i,lw}];
Do[
fglob[[   V[[el,wezel]]   ]]+=ftab[[el,wezel]];
,{wezel,1,3},{el,1,le}]; 
Do[
Kglob[[   V[[el,w1]],V[[el,w2]]   ]]+=Ktab[[el,w1,w2]];
,{w1,1,3},{w2,1,3},{el,1,le}];

Warunki brzegowe i wyznaczenie naprężeń w elemetach a potem węzłach

DL=N[1.0 (10^100)];
Do[Kglob[[brzeg[[i]],brzeg[[i]]]]=DL,{i,Length[brzeg]}];
wynik=Inverse[Kglob].fglob//Flatten;
Txze=Tyze=Table[0,{i,le}];(*naprężenia w elementach*)
Txzw=Tyzw=Table[0,{i,lw}];(*naprężenia w węzłach*)
wynike=Table[0,{i,le}];(*średnia wartość funkcji ϕ dla elementu*)
Do[Txze[[i]]=(1/(2Δ))(wynik[[V[[i,1]]]]c1+wynik[[V[[i,2]]]]c2+wynik[[V[[i,3]]]]c3)/.
{x1-> Punkty[[V[[i,1]],1]],y1-> Punkty[[V[[i,1]],2]],
x2-> Punkty[[V[[i,2]],1]],y2-> Punkty[[V[[i,2]],2]],
x3-> Punkty[[V[[i,3]],1]],y3-> Punkty[[V[[i,3]],2]]}
,{i,le}];
Do[Tyze[[i]]=-(1/(2Δ))(wynik[[V[[i,1]]]]b1+wynik[[V[[i,2]]]]b2+wynik[[V[[i,3]]]]b3)/.
{x1-> Punkty[[V[[i,1]],1]],y1-> Punkty[[V[[i,1]],2]],
x2-> Punkty[[V[[i,2]],1]],y2-> Punkty[[V[[i,2]],2]],
x3-> Punkty[[V[[i,3]],1]],y3-> Punkty[[V[[i,3]],2]]}
,{i,le}];
Do[wynike[[i]]=((wynik[[V[[i,1]]]]+wynik[[V[[i,2]]]]+wynik[[V[[i,3]]]])/3),{i,le}];
Do[
	Txzw[[V[[e,w]]]]+=((Txze[[e]])/(iw[[V[[e,w]]]]));	
	Tyzw[[V[[e,w]]]]+=((Tyze[[e]])/(iw[[V[[e,w]]]]));
,{w,3},{e,le}]

Przedstawienie wyników (wykres wypadkowych naprężeń stycznych):

Wykres przedstawia wartości naprężeń przemnożone przez ((Subscript[I, s])/((Subscript[M, s])j)), gdzie (Subscript[I, s])to policzony poniżej moment bezwładności, (Subscript[M, s])to moment skręcający, zaś j to jednostka w jakiej zadano wymiary. Aby otrzymać naprężenie, należy podzielić wartość z wykresu przez ten ułamek.

dataT=Table[{Punkty[[i,1]],Punkty[[i,2]],(Sqrt[((Txzw[[i]])^2)+((Tyzw[[i]])^2)])},{i,lw}];
MinX=Min[Table[Punkty[[i,1]],{i,lw}]];MaxX=Max[Table[Punkty[[i,1]],{i,lw}]];
MinY=Min[Table[Punkty[[i,2]],{i,lw}]];MaxY=Max[Table[Punkty[[i,2]],{i,lw}]];
Show[ListContourPlot[dataT,AspectRatio->Automatic,PlotLegends->Automatic,PlotLabel->"τ((Subscript[I, s])/((Subscript[M, s]) j))"],RegionPlot[!({x,y}∈ℛ),{x,MinX,MaxX},{y,MinY,MaxY},PlotStyle->White]]

Moment bezwładności na skręcanie (w [jednostka w jakiej zadano ((wymiary])^4)). Aby uzyskać sztywność na skręcanie należy pomnożyć przez G

2*Sum[Δ*wynike[[i]]/.
{x1-> Punkty[[V[[i,1]],1]],y1-> Punkty[[V[[i,1]],2]],x2-> Punkty[[V[[i,2]],1]],y2-> Punkty[[V[[i,2]],2]],x3-> Punkty[[V[[i,3]],1]],y3-> Punkty[[V[[i,3]],2]]}
,{i,le}]

Czas w sekundach potrzebny na wykonanie obliczeń

TimeUsed[]-tu