> restart:
# HF, Übung 6: Zylindrischer Flüssigkeitsbehälter
# als Beispiel einer rotationssymmetrisch belasteten Zylinderschale
# 
> with(linalg):
> vecgrad:=proc()
>  local i,v,q,co,h,vq,hq,hinv,dia,vdia,ddia;
>  description ` Vektorgradient:  v dyadisch nabla`;
>  v:=args[1];
>  q:=args[2];
>  h:=[1,1,1];
>  if nargs>2 then 
>    co:=args[3];
>    if type(co,`=`)then
>      if rhs(co)=cylindrical then
>        h:=[1,q[1],1] fi;
>      if rhs(co)=spherical then
>        h:=[1,q[1],q[1]*sin(q[2])] fi
>    else h:=co
>    fi
>  fi;
>  vdia:=array(sparse,1..3,1..3):
>  for i to 3 do vdia[i,i]:=v[i] od;
>  vq:=jacobian(v,q);
>  hq:=jacobian(h,q);
>  hinv:=array(sparse,1..3,1..3):
>  for i to 3 do hinv[i,i]:=1/h[i] od;
>  dia:=evalm(hq&*hinv&*v);
>  ddia:=array(sparse,1..3,1..3):
>  for i to 3 do ddia[i,i]:=dia[i] od;
>  evalm((vq-hinv&*transpose(hq)&*vdia
>  +ddia)&*hinv);
>  end:
>  tensdiv:=proc()
>  local i,l,t,q,co,h,hinv,hq,vol,tz,tzeile,di;
>  description ` Tensordivergenz:  t punkt nabla`;
>  t:=args[1];
>  q:=args[2];
>  h:=[1,1,1];
>  if nargs>2 then 
>    co:=args[3];
>    if type(co,`=`)then
>      if rhs(co)=cylindrical then
>        h:=[1,q[1],1] fi;
>      if rhs(co)=spherical then
>        h:=[1,q[1],q[1]*sin(q[2])] fi
>    else h:=co
>    fi
>  fi;
>  hinv:=array(sparse,1..3,1..3):
>  for i to 3 do hinv[i,i]:=1/h[i] od;
>  vol:=h[1]*h[2]*h[3];
>  tz:=evalm(t&*hinv*vol);
>  hq:=jacobian(h,q);
>  for i to 3 do
>      di[i]:=diverge(row(tz,i),q)/vol +1/h[i]*
>  add((t[l,i]*hq[i,l]-t[l,l]*hq[l,i])/h[l],
>  l=1..3)
>  od;
>  vector([di[1],di[2],di[3]])
>  end: 
> alias(zyl=([r,phi,z],coords=cylindrical)):
# 
# Krümmungstensor:
# 
> n:=vector([-1,0,0]);
> grad_n:=vecgrad(n,zyl);
> k:=subs(r=a,evalm(-grad_n));
# 
# Aufstellung der Gleichgewichtsbedingungen:
# 
> s:=matrix([[0,0,0],[0,sphi(z),0],[0,0,sz(z)]]);
> div_s_pre:=tensdiv(s,zyl);
> div_s:=subs(r=a,eval(div_s_pre));
> qp:=vector([0,0,qz(z)]);
> div_qp:=diverge(qp,zyl);
# 
# Kraftgleichgewichtsbedingungen normal zur Schale und in Richtung der
# Mantellinie:
# 
> gl1:=div_s[1]-div_qp+gam*(l-z)=0;
> gl2:=div_s[3]=0;
# 
# 
> m:=matrix([[0,0,0],[0,mphi(z),0],[0,0,mz(z)]]);
> div_mpre:=tensdiv(m,zyl);
> div_m:=subs(r=a,eval(div_mpre));
# 
# Momentengleichgewichtsbedingung um eine Achse tangential zum
# Breitenkreis:
# 
> gl3:=div_m[3]-qz(z)=0;
# 
# Stoffgesetz: Isotrop-elastisch
# 
> sphi:=z->2*G*h/(1-nu)*(ephi(z)+nu*ez(z)):
> sz:=z->2*G*h/(1-nu)*(ez(z)+nu*ephi(z)):
> sphi(z);sz(z);
> mphi:=z->2*G*h^3/12/(1-nu)*(bphi(z)+nu*bz(z)):
> mz:=z->2*G*h^3/12/(1-nu)*(bz(z)+nu*bphi(z)):
> mphi(z); mz(z);
# 
# Kirchhoffsche Hypothese:
# 
> omegaz:=-diff(w(z),z);
# Zusammenhänge zwischen Verzerrungen und Weggrößen
# 
> up:=vector([0,0,uz(z)]);
> grad_u:=vecgrad(up,zyl);
> ephi:=z->-w(z)*k[2,2]:
> ez:=z->grad_u[3,3]:
> ephi(z); ez(z);
> omegap:=vector([0,0,omegaz]);
> grad_om:=vecgrad(omegap,zyl);
> bphi:=z->grad_om[2,2]:
> bz:=z->grad_om[3,3]:
> bphi(z); bz(z);
# 
# Die Gleichgewichtsbedingungen, ausgedrückt als drei gekoppelte
# gewöhnliche Differentialgleichungen in den Weggrößen und der
# Querkraft:
# 
> fak:=1/2/G/h;
> g1:=expand(gl1*fak*2);
> g2:=expand(gl2*fak*(1-nu));
> g3:=expand(gl3*fak*2);
# 
# Der Membranspannungszustand als partikuläre Lösung
# dieser drei Gleichungen
# 
# Wir versuchen folgenden Ansatz:
# 
> uz_Membran(z):=-nu/2/a*last*(l-z)^2+D1*z+D2;
> uz(z):=uz_Membran(z):
> w_Membran(z):=nu*a*D1+last*(l-z);
> w(z):=w_Membran(z):
> qz(z):=0;
# 
# Im Folgenden sieht man, daß die drei Gleichungen durch diesen Ansatz
# in der Tat erfüllt werden,
# wenn für   last  ein geeigneter Wert gesetzt wird,
# und daß es nur Membrankräfte gibt, aber keine Biegemomente:
# 
> g2;g3;
> lastwert:=last=-gam*a^2/2/G/h/(1+nu);
> assign(lastwert);
> simplify(eval(g1));
> simplify(eval(s[2,2])); # Das gibt die bekannte Kesselformel
> simplify(eval(s[3,3]));
> simplify(m[2,2]);
> simplify(m[3,3]);
> last:='last';
# 
# Randbedingungen für Schnittgrößen am oberen Behälterrand ( z=l):
# m[3,3]=0, qz=0 sind bereits erfüllt
# s[3,3]=0 läßt sich erfüllen durch die Wahl:
# 
> D1:=0:
# 
# Randbedingungen für Weggrößen an der Einspannung  (z=0):
# 
# 
# uz=0 ließe sich erfüllen durch die Wahl:
# 
> D2=-1/4*nu*a*gam*l^2/G/h/(1+nu):
> eval(w(z));
> omegaz;
# w=0, omegaz=0 an der Einspannung sind mit der Membranlösung nicht zu
# erfüllen.
# 
# 
# Es muß eine Biegelösung überlagert werden. Diese muß dem homogenen
# Differentialgleichungssystem genügen. Da das System linear ist und
# konstante Koeffizienten besitzt, führt ein Exponentialansatz zum Ziel:
> qz(z):=Q*exp(-lambda*z/L);
> w(z):=W*a*exp(-lambda*z/L);
> uz(z):=U*a*exp(-lambda*z/L);
> gln:=[subs(gam=0,g1), g2, g3];
# 
# Wir dividieren diese Gleichungen durch den Exponentialausdruck. Das
# Ergebnis ist dasselbe, als wenn z=0 gesetzt wird.
# 
> glnneu:=eval(subs(z=0,gln));
# 
# Die Matrix dieses homogenen Gleichungssystems wird wie folgt
# generiert:
# 
> A:=genmatrix(glnneu,[Q,W,U]);
# Nichttriviale Lösungen der homogenen Gleichungen erfordern das
# Verschwinden der Determinante:
> pol:=det(A);
> Lwahl:=L=sqrt(a*h/sqrt(3*(1-nu^2)));
> assign(Lwahl);
> pols:=simplify(pol,symbolic);
> wurzeln:=solve(pols,lambda);
# 
# Bedeutung der komplexen Exponenten für die Lösungen: 
# 
> L:='L':
> f:=exp(-(1+I)*z/L);
> evalc(f);
# Nur die Wurzeln mit positivem Realteil liefern Lösungen, die mit
# wachsendem z rasch abklingen und  den oberen Behälterrand, wo die
# Randbedingungen bereits erfüllt sind, nicht mehr beeinflussen. Deren
# reelle Bestandteile haben die Gestalt - mit s=z/L -
# 
> f1:=exp(-s)*cos(s); f2:=exp(-s)*sin(s);
# 
# 
> plot({f1,f2},s=0..4);
# Man sieht: Die Randstörungen sind abgeklungen bei etwa  s=4, also
# z=4*L (<a, wenn h<a/10)
# 
> qz(z):=exp(-z/L)*(Q1*cos(z/L)+Q2*sin(z/L));
> w(z):=exp(-z/L)*(W1*cos(z/L)+W2*sin(z/L));
> 
> 
> uz(z):=exp(-z/L)*(U1*cos(z/L)+U2*sin(z/L));
# 
# Einsetzen in die homogenen Differentialgleichungen:
# 
> gln:=expand([subs(gam=0,g1)*exp(z/L), g2*exp(z/L), g3*exp(z/L)]);
# 
# Erfüllung der Randbedingungen bei z=0:
# 
> randb:=eval(subs(z=0,w(z)+w_Membran(z)))=0,
> eval(subs(z=0,diff(w(z)+w_Membran(z),z)))=0,
> eval(subs(z=0,uz(z)+uz_Membran(z)))=0;
> 
> lo:=solve({randb},{W1,W2,U1});
> assign(lo);
# 
# 
# Jetzt sind noch die Differentialgleichungen zu erfüllen.
# Sie haben die Gestalt
# 
# g1c*cos(z/L)+g1s*sin(z/L) = 0
# g2c*cos(z/L)+g2s*sin(z/L) = 0
# g3c*cos(z/L)+g3s*sin(z/L) = 0
# 
# Sie sind für alle z erfüllt, wenn gilt 
# 
# g1c = 0,  g1s = 0, g2c  = 0, g2s = 0,  g3c = 0,  g3s = 0
# 
# 
> g3c:=eval(subs(z=0,gln[3]));
> g3s:=eval(subs(z=L*Pi/2,gln[3]));
> loesunga:=solve({g3c,g3s},{Q1,Q2});
> assign(loesunga);
> g2c:=eval(subs(z=0,gln[2]));
> g2s:=eval(subs(z=L*Pi/2,gln[2]));
> loesungb:=solve({g2c,g2s},{U2,D2});
> assign(loesungb);
> g1c:=eval(subs(z=0,gln[1]));
> g1s:=eval(subs(z=L*Pi/2,gln[1]));
> g1csimp:=simplify(g1c);
> g1ssimp:=simplify(g1s);
> assign(Lwahl);
> simplify(g1csimp);
> simplify(g1ssimp);
# 
# Ergebnisse
# 
# Verlauf der Verformungen
# 
> L:='L';
> we:=eval(w(z)+w_Membran(z));
> wp:=collect(we,last);
> dwp:=diff(wp,z);
> plot(subs(last=1,L=1,l=50,wp),z=0..10);
> plot(subs(last=1,L=1,l=50,wp),z=0..1);
> plot(subs(last=1,L=1,l=50,dwp),z=0..10);
# 
# Verlauf der Schnittgrößen
# 
> me:=expand(eval(m[3,3]/G/h^3/last));
> plot(subs(L=1,l=50,nu=0.3,me),z=0..10);
> qe:=expand(eval(qz(z)/G/h^3/last));
> plot(subs(L=1,l=50,nu=0.3,qe),z=0..10);
>