Berechnung von mehrschichtigen Verbundrohrstrukturen bei rotationssymmetrischen Belastungen
Lars Ulke-Winter
BTU Cottbus-Senftenberg
Fachgebiet Polymerbasierter Leichtbau
lars.ulke-winter@b-tu.de

Clear["Global`*"]

Voraussetzungen:
	1. Keine Schubspannungen auf den Mantelflächen
	2. Ansätze für die Verschiebungsfunktionen mit konstanten Längsdehnungen (Subscript[K, c]) und konstanten Verdrillungen (Subscript[K, d])
	   der gesamten Zylinderstruktur

Lösen des rotationssymmetrischen Problems

Eingabe

Clear[f,m,pi,pa,ω,E1, E2,nu21, nu23,G12,ρ,oris,n,ti,ts,rInnen,laenge];

Belastungen

Statisch

f=3000; m=40000;pi=5; pa=0;

Dynamisch durch Rotation um die Längsachse (z.B. in rad/s)

ω=0;

Materialeigenschaften im System mit 1Faserlängsrichtung, 2 Faserquerichtung in der intralaminaren Ebene, 3 Faserquerrichtung in der interlaminaren Ebene

Ingenieurskonstanten Beispiel aus VDI 2014 Blatt-3, CFK S.27
Bei transversaler Isotropie gilt:

({{({{(Subscript[E, 2])=(Subscript[E, 3])}, {(Subscript[G, 12])=(Subscript[G, 13])}, {(Subscript[ν, 21])=(Subscript[ν, 31])}, {(Subscript[G, 23])=((Subscript[E, 2])/(2(1+(Subscript[ν, 23]))))}})}})

Die Indizierung, insbesondere der Querkontraktionszahlen, richtet sich nach der VDI 2014 Blatt-3 (Subscript[ν, ij]) charakterisiert die Querdehnung in i - Richtung bei Beanspruchung in j -Richtung (bei vielen FEM-Systemen bspw. ABAQUS, wird ist die Reihenfolge  umgekehrt angegeben)

E1=128600; E2=9500;
nu21=0.27 ; nu23=0.471;
G12=5400;
ρ=1.8 (10^(-9));

Orientierungen der Einzelschichten (im Uhrzeigersinn positiv)

Achtung: In einigen FEM-Programmen werden die Faserorientierungen gegen den Uhrzeigersinn positiv angegeben (Vorzeichenwechsel beachten!)

oris=-1*{-45°,-30°,15°,60°}; n=Length[oris];

Schichtdicken

wenn konstant:

ti=0.208975; ts=ConstantArray[ti,n]

wenn beliebig:

(*ts={0.8359}*)

Geometrie

rInnen=20; laenge=250;

Hilfsfunktionen zum aufstellen und der Transformation der gemischtvarianten Elastizitätstensoren 4-Stufe (Subscript[(C^ij), kl])

Clear[createStiffTensorKartesisch,transformStiffTensor];

Stellt den 4-stufigen Steifigkeitstensor im kartesischen System bezüglich seiner Hauptrichtung auf (1-Faserichtung, 2-Faserquerrichtung, 3-Dickenrichtung), wobei von einem transversal isotropen Materialgesetz ausgegangen wird.

createStiffTensorKartesisch[e1_,e2_,nu21_,nu23_,g12_]:=Module[{s11=1/e1,s12=-(nu21/e1),s13=-(nu21/e1),s22=1/e2,s23=-(nu23/e2),s33=1/e2,s44=(2(1+nu23))/e2,s55=1/g12,s66=1/g12,
                     sij,cij,stiff},
												        sij=({ 
                                                                                                           {s11, s12, s13, 0, 0, 0},
                                                                                                           {s12, s22, s23, 0, 0, 0},
                                                                                                                       {s13, s23, s33, 0, 0, 0},
                                                                                                                       {0, 0, 0, s44, 0, 0},
                                                                                                                       {0, 0, 0, 0, s55, 0},
                                                                                                                       {0, 0, 0, 0, 0, s66}
                                                                                                                       });cij=Inverse[sij];
                                    (* 4-stufigen Steifigkeitstensor im ueblichen Kartesischen System aufstellen *)
				stiff=SymmetrizedArray[{{1,1,1,1}->cij[[1,1]],{1,1,2,2}->cij[[1,2]],{1,1,3,3}->cij[[1,3]],
										  {2,2,2,2}->cij[[2,2]],{2,2,3,3}->cij[[2,3]],
					                                                 {3,3,3,3}->cij[[3,3]],
					                                                 {2,3,2,3}->cij[[4,4]],
					                                                 {1,3,1,3}->cij[[5,5]],
					                                                 {1,2,1,2}->cij[[6,6]],{_,_,_,_}->0},{3,3,3,3},{{{2,1,3,4},1},{{1,2,4,3},1},{{3,4,1,2},1}}]]

Transformiert die Komponenten eines gemischtvarianten 4-stufigen Tensors (stiff) (Subscript[(C^ij), kl]) auf eine neue Basis die sich nach transMatrix  transformiert

transformStiffTensor[stiff_,transMatrix_]:=Module[{stiffTrans=TensorContract[stiff⊗Inverse[transMatrix]⊗Inverse[transMatrix]⊗transMatrix⊗transMatrix,{{1,6},{2,8},{3,9},{4,11}}]},
stiffTrans]

Aufbauen der Elastizitätstensoren (Subscript[(C^ij), kl]) der Einzelschichten

Kartesischer Steifigkeitstensor im Hauptachsensystem mit der “üblichen” Richtungsnummierung s.o.

Clear[stiffKart,transformMatrixChangeAxes,transformMatrixRotUm1,α,stiffZyl]

stiffKart=createStiffTensorKartesisch[E1,E2,nu21,nu23,G12]
stiffKart//MatrixForm

Transformieren der Steifigkeitstensoren auf: 
Achsen tauschen 1->3, 3->1; Drehen um die neue 1-Achse um (Subscript[α, i]), in Zylinderkoordinaten überführen

Transformationsmatrizen

Achsen tauschen 1->3, 3->1

transformMatrixChangeAxes=(({{0, 0, 1}, {0, 1, 0}, {1, 0, 0}}));

Rotation um α um 1-Achse

transformMatrixRotUm1=RotationMatrix[α,{1,0,0}] ;
transformMatrixRotUm1//MatrixForm

Transformieren der Steifigkeitstensoren auf das Rohr-Bezugssystem mit Faserlängsrichtung z-Achse und Dickenrichtung r-Achse (Bauteilsystem)

stiffZyl=Table[Fold[transformStiffTensor,stiffKart,{transformMatrixChangeAxes,(transformMatrixRotUm1/.{α->oris[[i]]})}],{i,n}];

Zwischenrechnungen

Clear[u1,u2,var,us,ris,kc,kd]

Allgemeiner Verschiebungsansatz für die Einzelschichten

Für (Subscript[C, 11])≠ (Subscript[C, 22]) bei allen Schichten mit Orientierung ≠0°

u1[r_,k1_,k2_]:=k1 (r^(Sqrt[(c22/c11)]))+k2 (r^(-(Sqrt[(c22/c11)])))+kc r((c23-c13)/(c11-c22))+kd (r^2)((c24-2c14)/(4c11-c22))-(ω^2)(r^3)(ρ/(9 c11-c22))

Für (Subscript[C, 11])= (Subscript[C, 22]) bei allen Schichten mit Orientierung =0° (Axiale Verstärkung oder Isotropie)

u2[r_,k1_,k2_]:=k1 r+k2 (1/r)-(ω^2)(r^3)(ρ/(8 c11))

Zu bestimmende unbekannte Variablen und Integrationskonstanten

Tiefgestellter Index kennzeichnet zugehörige Schicht (gezählt von Rohrinnenseite)

var=Join[{kc,kd},Array[{(Subscript[k1, #]),(Subscript[k2, #])}&,{n}]]//Flatten

Verschiebungsansätze der Einzelschichten parametrisiert mit ihren zugehörigen Integrationskonstanten

Je nachdem ob eine rein axiale Verstärkung (0°) vorliegt oder nicht, ergibt sich eine unterschiedliche Differentialgleichungen mit den jeweiligen Lösungen u1(r,(Subscript[k, 1i]),(Subscript[k, 2i])) bzw. u2(r,(Subscript[k, 1i]),(Subscript[k, 2i]))

wenn oris[[i]]≠0 dann: Ansatz u1 sonst: Ansatz u2

us=Table[If[oris[[i]]≠0°,u1[r,(Subscript[k1, i]),(Subscript[k2, i])],u2[r,(Subscript[k1, i]),(Subscript[k2, i])]],{i,n}];
us//TableForm

(*us=Array[u1[r,(Subscript[k1, #]),(Subscript[k2, #])]&,{n}]*)

Grenzen der Einzelschichten

ris=FoldList[Plus,rInnen,ts]

Gleichungssystem aufstellen

Clear[gl1,gl2,gl3,gl4,glUeb1,glUeb2]

Gleichgewichtsbedingungen

Resultierende Längskraft

gl1=Plus@@Table[(2π (∫_(ris[[i]])^(ris[[i+1]]))(D[us[[i]],r]stiffZyl[[i]][[3,3,1,1]]+((us[[i]])/r)stiffZyl[[i]][[3,3,2,2]]+kc stiffZyl[[i]][[3,3,3,3]]+kd r stiffZyl[[i]][[3,3,2,3]])r\[DifferentialD]r/.{c11->stiffZyl[[i]][[1,1,1,1]],c22->stiffZyl[[i]][[2,2,2,2]],c23->stiffZyl[[i]][[3,3,2,2]],c13->stiffZyl[[i]][[3,3,1,1]],c14->stiffZyl[[i]][[2,3,1,1]],c24->stiffZyl[[i]][[2,3,2,2]]}),{i,n}]==f

Resultierendes Torsionsmoment

gl2=Plus@@Table[
(2π (∫_(ris[[i]])^(ris[[i+1]]))(D[us[[i]],r]stiffZyl[[i]][[2,3,1,1]]+((us[[i]])/r)stiffZyl[[i]][[2,3,2,2]]+kc stiffZyl[[i]][[2,3,3,3]]+kd r stiffZyl[[i]][[2,3,2,3]])(r^2)\[DifferentialD]r)/.{c11->stiffZyl[[i]][[1,1,1,1]],c22->stiffZyl[[i]][[2,2,2,2]],c23->stiffZyl[[i]][[3,3,2,2]],c13->stiffZyl[[i]][[3,3,1,1]],c14->stiffZyl[[i]][[2,3,1,1]],c24->stiffZyl[[i]][[2,3,2,2]]},{i,n}]==m

Randbedingungen an der Rohrinnen- und -außenseite

gl3=(D[us[[1]],r]stiffZyl[[1]][[1,1,1,1]]+((us[[1]])/r) stiffZyl[[1]][[1,1,2,2]]+kc stiffZyl[[1]][[1,1,3,3]]+kd r stiffZyl[[1]][[1,1,2,3]]/.{r->ris[[1]],c11->stiffZyl[[1]][[1,1,1,1]],c22->stiffZyl[[1]][[2,2,2,2]],c23->stiffZyl[[1]][[3,3,2,2]],c13->stiffZyl[[1]][[3,3,1,1]],c14->stiffZyl[[1]][[2,3,1,1]],c24->stiffZyl[[1]][[2,3,2,2]]})==-pi
gl4=(D[us[[n]],r]stiffZyl[[n]][[1,1,1,1]]+((us[[n]])/r) stiffZyl[[n]][[1,1,2,2]]+kc stiffZyl[[n]][[1,1,3,3]]+kd r stiffZyl[[n]][[1,1,2,3]]/.{r->ris[[-1]],c11->stiffZyl[[n]][[1,1,1,1]],c22->stiffZyl[[n]][[2,2,2,2]],c23->stiffZyl[[n]][[3,3,2,2]],c13->stiffZyl[[n]][[3,3,1,1]],c14->stiffZyl[[n]][[2,3,1,1]],c24->stiffZyl[[n]][[2,3,2,2]]})==-pa

Übergangsbedingungen

Gleiche Radialspannungen zwischen 2 Schichten

glUeb1=Table[(D[us[[i]],r]stiffZyl[[i]][[1,1,1,1]]+((us[[i]])/r) stiffZyl[[i]][[1,1,2,2]]+kc stiffZyl[[i]][[1,1,3,3]]+kd r stiffZyl[[i]][[1,1,2,3]]/.{r->ris[[i+1]],c11->stiffZyl[[i]][[1,1,1,1]],c22->stiffZyl[[i]][[2,2,2,2]],c23->stiffZyl[[i]][[3,3,2,2]],c13->stiffZyl[[i]][[3,3,1,1]],c14->stiffZyl[[i]][[2,3,1,1]],c24->stiffZyl[[i]][[2,3,2,2]]})==(D[us[[i+1]],r]stiffZyl[[i+1]][[1,1,1,1]]+((us[[i+1]])/r) stiffZyl[[i+1]][[1,1,2,2]]+kc stiffZyl[[i+1]][[1,1,3,3]]+kd r stiffZyl[[i+1]][[1,1,2,3]]/.{r->ris[[i+1]],c11->stiffZyl[[i+1]][[1,1,1,1]],c22->stiffZyl[[i+1]][[2,2,2,2]],c23->stiffZyl[[i+1]][[3,3,2,2]],c13->stiffZyl[[i+1]][[3,3,1,1]],c14->stiffZyl[[i+1]][[2,3,1,1]],c24->stiffZyl[[i+1]][[2,3,2,2]]}),{i,n-1}]

Gleiche Radialverschiebungen zwischen 2 Schichten

glUeb2=Table[(us[[i]]/.{r->ris[[i+1]],c11->stiffZyl[[i]][[1,1,1,1]],c22->stiffZyl[[i]][[2,2,2,2]],c23->stiffZyl[[i]][[3,3,2,2]],c13->stiffZyl[[i]][[3,3,1,1]],c14->stiffZyl[[i]][[2,3,1,1]],c24->stiffZyl[[i]][[2,3,2,2]]})==(us[[i+1]]/.{r->ris[[i+1]],c11->stiffZyl[[i+1]][[1,1,1,1]],c22->stiffZyl[[i+1]][[2,2,2,2]],c23->stiffZyl[[i+1]][[3,3,2,2]],c13->stiffZyl[[i+1]][[3,3,1,1]],c14->stiffZyl[[i+1]][[2,3,1,1]],c24->stiffZyl[[i+1]][[2,3,2,2]]}),{i,n-1}]

Lösen des Gleichungssystems

Clear[lsg,kc,kd,upw,u,v,w,σrpw,σϕpw,σzpw,τϕzpw,σr,σϕ,σz,τϕz]

Lösen

lsg=Solve[{gl1,gl2,gl3,gl4,glUeb1,glUeb2}//Flatten,var]

Lösung ausgeben

Gesamtdehnung in Axialrichtung und Gesamtverdrillung

kc=(kc/.lsg)[[1]]
kd=(kd/.lsg)[[1]]

Verschiebungsfunktionen

Für die jeweiligen Einzelschichten in radialer Richtung

uins=us/.lsg//Flatten;
uins//TableForm

einsetzen der jeweiligen Steifigkeiten

upw=Table[uins[[i]]/.{c11->stiffZyl[[i]][[1,1,1,1]],c22->stiffZyl[[i]][[2,2,2,2]],c23->stiffZyl[[i]][[3,3,2,2]],c13->stiffZyl[[i]][[3,3,1,1]],c14->stiffZyl[[i]][[2,3,1,1]],c24->stiffZyl[[i]][[2,3,2,2]]},{i,n}]

Damit ergibt sich in radialer Richtung eine abschnittsweise definierte Verschiebungsfunktion

u=Piecewise[Table[{upw[[i]],ris[[i]]<=r<ris[[i+1]]},{i,n}]~Join~{{upw[[-1]],r==ris[[-1]]}}]

Für den Gesamtverbund in Umfangs - und Rohrlängsrichtung

v=kd laenge r
w=kc z

Spannungsfunktionen

Da die Spannung von den radialen Verschiebungen abhängen, sind diese ebenfalls abschnittsweise definiert

σrpw=Table[D[upw[[i]],r]stiffZyl[[i]][[1,1,1,1]]+((upw[[i]])/r)stiffZyl[[i]][[1,1,2,2]]+kc stiffZyl[[i]][[1,1,3,3]]+kd r stiffZyl[[i]][[1,1,2,3]],{i,n}];
σϕpw=Table[D[upw[[i]],r]stiffZyl[[i]][[2,2,1,1]]+((upw[[i]])/r) stiffZyl[[i]][[2,2,2,2]]+kc stiffZyl[[i]][[2,2,3,3]]+kd r stiffZyl[[i]][[2,2,2,3]],{i,n}];
σzpw=Table[D[upw[[i]],r]stiffZyl[[i]][[3,3,1,1]]+((upw[[i]])/r) stiffZyl[[i]][[3,3,2,2]]+kc stiffZyl[[i]][[3,3,3,3]]+kd r stiffZyl[[i]][[3,3,2,3]],{i,n}];
τϕzpw=Table[D[upw[[i]],r]stiffZyl[[i]][[2,3,1,1]]+((upw[[i]])/r) stiffZyl[[i]][[2,3,2,2]]+kc stiffZyl[[i]][[2,3,3,3]]+kd r stiffZyl[[i]][[2,3,2,3]],{i,n}];

σr=Piecewise[Table[{σrpw[[i]],ris[[i]]<=r<ris[[i+1]]},{i,Length[oris]}]~Join~{{σrpw[[-1]],r==ris[[-1]]}}];
σϕ=Piecewise[Table[{σϕpw[[i]],ris[[i]]<=r<ris[[i+1]]},{i,Length[oris]}]~Join~{{σϕpw[[-1]],r==ris[[-1]]}}];
σz=Piecewise[Table[{σzpw[[i]],ris[[i]]<=r<ris[[i+1]]},{i,Length[oris]}]~Join~{{σzpw[[-1]],r==ris[[-1]]}}];
τϕz=Piecewise[Table[{τϕzpw[[i]],ris[[i]]<=r<ris[[i+1]]},{i,Length[oris]}]~Join~{{τϕzpw[[-1]],r==ris[[-1]]}}];

Grafische Ausgabe (Bauteilsystem, Rohr)

{"f="<>ToString[f]<>"; ","m="<>ToString[m]<>"; ","pi="<>ToString[pi]<>"; ","pa="<>ToString[pa]<>"; ","ω="<>ToString[ω]}//Row
"ti="<>ToString[ti]
Row[oris,"/",Frame->True]

GraphicsRow[{Show[Plot[Evaluate[u],{r,ris[[1]],ris[[-1]]},PlotStyle->Orange,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[r],HoldForm[u]},PlotLabel->HoldForm[Radialverschiebung u über Dicke],LabelStyle->{13,GrayLevel[0]}],
Show[Plot[Evaluate[v],{r,ris[[1]],ris[[-1]]},PlotStyle->Brown,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[r],HoldForm[ϕ]},PlotLabel->HoldForm[Verdrehung v über Dicke (am Rohrende)],LabelStyle->{13,GrayLevel[0]}],
Show[Plot[Evaluate[w],{z,0,laenge},PlotStyle->Blue,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[z],HoldForm[w]},PlotLabel->HoldForm[Verschiebung w über Länge ],LabelStyle->{13,GrayLevel[0]}]},Frame->All]

GraphicsGrid[{
{Show[Plot[Evaluate[σr],{r,ris[[1]],ris[[-1]]},PlotStyle->Orange,Filling->Automatic,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[r],HoldForm[Subscript[σ, r]]},PlotLabel->HoldForm[Radialspannung Subscript[σ, r] über Dicke],LabelStyle->{13,GrayLevel[0]}],Show[Plot[Evaluate[σϕ],{r,ris[[1]],ris[[-1]]},PlotStyle->Brown,Filling->Automatic,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[r],HoldForm[Subscript[σ, ϕ]]},PlotLabel->HoldForm[Umfangsspannung Subscript[σ, ϕ] über Dicke],LabelStyle->{13,GrayLevel[0]}]},

{Show[Plot[Evaluate[σz],{r,ris[[1]],ris[[-1]]},PlotStyle->Blue,Filling->Automatic,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[r],HoldForm[Subscript[σ, z]]},PlotLabel->HoldForm[Längsspannung Subscript[σ, z] über Dicke],LabelStyle->{13,GrayLevel[0]}],Show[Plot[Evaluate[τϕz],{r,ris[[1]],ris[[-1]]},PlotStyle->Red,Filling->Automatic,GridLines->Automatic],ImageSize->Large,AxesLabel->{HoldForm[r],HoldForm[Subscript[τ, ϕz]]},PlotLabel->HoldForm[Schubspannung Subscript[τ, ϕz] über Dicke],LabelStyle->{13,GrayLevel[0]}]}},Frame->All]

Tabellarische Ausgabe Schichtunter-, -mitte, -oberseite (Bauteilsystem, Rohr)

valuesSectionP=Table[{ris[[i]]//Chop,ris[[i]]+((ris[[i+1]]-ris[[i]])/2)//Chop,ris[[i+1]]//Chop},{i,n}]//Flatten;
valuesSigmaR=Table[(σr/.r->valuesSectionP[[i]])//Chop,{i,Length[valuesSectionP]}]//Partition[#,3]&;
valuesSigmaϕ=Table[(σϕ/.r->valuesSectionP[[i]])//Chop,{i,Length[valuesSectionP]}]//Partition[#,3]&;
valuesSigmaZ=Table[(σz/.r->valuesSectionP[[i]])//Chop,{i,Length[valuesSectionP]}]//Partition[#,3]&;
valuesTauϕZ=Table[(τϕz/.r->valuesSectionP[[i]])//Chop,{i,Length[valuesSectionP]}]//Partition[#,3]&;
Table[TableForm[MapThread[List,{Partition[valuesSectionP,3][[i]],valuesSigmaR[[i]],valuesSigmaϕ[[i]],valuesSigmaZ[[i]],valuesTauϕZ[[i]]}],TableAlignments->Center,TableHeadings->{{"Schicht "<>ToString[i],oris[[i]]},{"r","(Subscript[σ, r])","(Subscript[σ, ϕ])","(Subscript[σ, z])","(Subscript[τ, ϕz])"}}],{i,n}]//Column