> restart:
# HF, Übung 3: Hohlzylinder, rotationssymmetrischer Zustand
# 
# Spannungen aus Erwärmung und Rotation
# 
# Achtung: Wegen der Ungewissheit der Nummerierung der 
# Integrationskonstanten kann diese Sitzung nur Befehl für
# Befehl abgearbeitet werden
# 
> 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: 
> id:=diag(1,1,1);
# 
# Geometrisch gemessene Verzerrungen:
# 
> e_geom:=matrix(3,3):
# 
# Verzerrungsanteil, der durch Spannungen hervorgerufen ist:
# 
> e_spg:=evalm(e_geom-alpha*temp*id);
# 
# Isotropes Hookesches Gesetz:
# 
> t:=evalm(2*G*(e_spg+nu/(1-2*nu)*trace(e_spg)*id));
# 
# Rotationssymmetrischer ebener Verschiebungsansatz:
# 
> u:=vector([ur(r),0,0]);
> alias(zyl=([r,phi,z],coords=cylindrical)):
> grad_u := vecgrad(u,zyl);
# 
# Die (linearisierten) geometrischen Verzerrungen sind der symmetrische
# Teil des Verschiebungsgradienten:
# 
> e_geom:=evalm((grad_u + transpose(grad_u))/2);
# 
# Rotationssymmetrischer Temperaturansatz:
# 
> temp:=theta(r);
> la:=laplacian(temp,zyl);
# 
# Lösung der stationären Wärmeleitungsgleichung mit den Randbedingungen:
# Außentemperatur = 0, Innentemperatur = ti
# 
> loestemp:=dsolve({la=0,theta(ra)=0,theta(ri)=ti},theta(r));
> theta:=unapply(rhs(loestemp),r);
# 
# Aufstellung und Lösung der Gleichgewichtsbedingungen:
# 
# Volumenkräfte sind Zentrifugalkräfte
# 
> divt:=tensdiv(t,zyl):
> divt1:=expand(divt[1]);
> divt[2];
> divt[3];
> loesu:=dsolve(divt1+rho*omega^2*r=0,ur(r));
> # 
> ur:=unapply(rhs(loesu),r);
# 
# Bestimmung der zwei Integrationskonstanten aus den Bedingungen, dass 
# Innen- und Außenoberfläche spannungsfrei sein sollen:
# 
> loes:=solve({subs(r=ri,t[1,1])=0,subs(r=ra,t[1,1]=0)},
> {_C1,_C2});
> assign(loes);
# 
# Beispiel: Stahl
# 
# Einheiten:  N, mm, s, K
# 
> G:=80000.; rho:=7.85E-9; nu:=0.3; alpha:=0.000012;
> ra:=100.;
# 
# Nur Lastfall Temperatur, innen 100 K
# 
> omega:=0.; ti:=100.;
# 
# Großer Innenradius:
# 
> ri:=80.;
> ptemp1:=plot(temp,r=ri..ra):
> ptemp1; # Temperatur
> pu1:=plot(ur(r),r=ri..ra,0..0.08):
> 
> pu1; # Radialverschiebung
> pt1:=plot([t[1,1],t[2,2]],r=ri..ra,color=[red,blue]): 
> 
> pt1; # Radial- und Umfangsspannung
# Kleiner Innenradius:
# 
> ri:=5.;
> ptemp2:=plot(temp,r=ri..ra):
> ptemp2; # Temperatur
> pu2:=plot(ur(r),r=ri..ra,0..0.04):
> pu2; # Radialverschiebung
> pt2:=plot([t[1,1],t[2,2]],r=ri..ra
> ,color=[red,blue]): 
> pt2; # Radial- und Umfangsspannung
# Vergleich der beiden Fälle:
# 
> with(plots):
> display({ptemp1,ptemp2}); # Temperatur
> display({pu1,pu2}); # Radialverschiebung
> display({pt1,pt2}); # Radial- und Umfangsspannung
# Nur Lastfall Zentrifugalkraft, 6000 Umdrehungen/min
# 
> ti:=0.; omega:= evalf(200*Pi);
# 
# Großer Innenradius:
# 
> ri:=80.;
> pu3:=plot(ur(r),r=ri..ra,0..0.015):
> 
> pu3; # Radialverschiebung
> pt3:=plot([t[1,1],t[2,2]],r=ri..ra,color=[red,blue]): 
> 
> pt3; # Radial- und Umfangsspannung
# Kleiner Innenradius:
> ri:=5.;
> pu4:=plot(ur(r),r=ri..ra,0..0.0022):
> 
> pu4; # Radialverschiebung
> pt4:=plot([t[1,1],t[2,2]],r=ri..ra,color=[red,blue]): 
> pt4; # Radial- und Umfangsspannung
# 
# Vollzylinder (Also Innenradius = 0)
# 
> _C1:='_C1': _C2:='_C2':
> ur(r);
# 
# Statt Randbedingung bei r=0: Forderung, dass die Verschiebungen
# beschränkt bleiben.
# Der Ausdruck  r*u(r) soll also bei r=0  den Wert  0  besitzen.
# 
> loes2:=solve({subs(r=0,expand(r*ur(r)))=0},{_C2});
> assign(loes2):
# 
# Spannungsfreiheit an der Außenfläche:
# 
> loes1:=solve({subs(r=ra,t[1,1])=0},{_C1});
> assign(loes1);
> 
> pu5:=plot(ur(r),r=0..ra,0..0.0022):
> 
> pu5; # Radialverschiebung
> pt5:=plot([t[1,1],t[2,2]],r=0..ra,color=[red,blue]):
> 
> pt5; # Radial- und Umfangsspannung 
# 
# 
> display({pu3,pu4,pu5});  # Radialverschiebung
> display({pt3,pt4,pt5}); # Radial- und Umfangsspannung
> 
>