> restart:
# HF, Übung 8: 
# Plastizierung einer Hohlkugel unter Innendruck
# 
# Das Beispiel stammt von Reuss (ZAMM 1930)
# 
# 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:
>  
> alias(sph=([r,theta,phi],coords=spherical)):
# 
# Der Spannungszustand ist rotationssymmetrisch zu jeder Geraden durch
# den Kugelmittelpunkt
# 
> t:=diag(sigmar,sigmat,sigmat);
> id:=diag(1,1,1);
> tdev:=evalm(t-trace(t)/3*id);
> sigmar:=sr(r); sigmat:=st(r);
> tev:=map(eval,t);
> div_t:=tensdiv(tev,sph);
> glgw:=simplify(div_t[1]);
# Das ist nur eine Differentialgleichung für zwei unbekannte
# Spannungskomponenten.
# 
# Im elastischen Bereich müssen daher auch die Verformungen betrachtet
# werden.
# (Das Problem ist also statisch unbestimmt.)
# 
# Im plastischen Bereich steht die Fließbedingung als zweite Gleichung
# zur Verfügung.
# 
# 
# Zunächst der plastische Bereich
# 
# Es soll die isotrope Fließbedingung nach von Mises benutzt werden
# 
# Einschub: Betrag einer Matrix
> s:=matrix(3,3):
> trace(s&*transpose(s));
> norm(s,frobenius);
# Ende Einschub
# 
# Jetzt Anwendung auf den Spannungsdeviator
> tdevbetragquadrat:=trace(tdev&*transpose(tdev));
> tdbq:=factor(tdevbetragquadrat);
> fliessbed:=tdbq-2/3*sigmav^2;
> solve(fliessbed=0,st(r));
> st:=unapply(vorz*sigmav+sr(r),r);
# Der Wert von  vorz  ist   +1  oder  -1.
> glgw;
# 
# Randbedingung innen
# 
> randi:=sr(ri)=-p_i;
> loesung:=dsolve({glgw=0,randi},sr(r));
> sr:=unapply(rhs(loesung),r);
> sr_plast:=sr(r); st_plast:=st(r);
> st_plast-sr_plast; # Kontrolle
# 
# Bei voll durchplastizierter Kugel gilt diese Lösung von r=ri bis r=ra
# 
# Am Außenrand ist dann
# 
> randa:=sr(ra)=-p_a;
> vorzeichen:=solve(randa,vorz);
# als der Außendruck p_a ist. Die maximal aufnehmbare Druckdifferenz 
# (bei der ungehindertes Fließen beginnt) ist offenbar:
> abs(p_i-p_a) = 2*sigmav*ln(ra/ri):
# 
# Qualitativer Verlauf der Spannungen im voll plastizierten Zustand:
# 
> plot({2*ln(r),2*ln(r)+1},r=0.6..0.8,-1.2..0.8);
# Verformung im elastischen und im plastischen Bereich
# 
# Der Kugeltensoranteil der Stoffgleichung ist rein elastisch
# und kann in beiden Gebieten gleichermaßen angesetzt werden
# 
> u:=vector([ur(r),0,0]);
> grad_u:=vecgrad(u,sph);
> tre:=trace(grad_u);
> kugel:=trt-3*K*tre;
# 
# Zunächst im plastischen Gebiet:
# 
> trt:=trace(t);
> loesung:=dsolve(kugel=0,ur(r));
> ur_plast:=subs(_C1=U1,rhs(loesung));
> sr:='sr': st:='st':
# 
# Nun im elastischen Gebiet
# Hier muss zusätzlich der Deviatoranteil der elastischen Stoffgleichung
# herangezogen werden.
# 
> tdevev:=map(eval,tdev);
> edev:=evalm(grad_u-tre/3*id);
> deviatorgl:=tdevev[1,1]-2*G*edev[1,1];
> trt:=trace(t);
> kugel;
> glgw;
# Es werden  sr  und   st  durch die Verschiebung    ur  ausgedrückt und
# dann in die Gleichgewichtsbedingung eingesetzt, die so zu einer
# Differentialgleichung 
# 2. Ordnung  in   ur  wird.
> loes:=solve({deviatorgl=0,kugel=0},{sr(r),st(r)});
> assign(loes);
> g:=expand(r^2*glgw/G);
> loesung:=dsolve({g=0},{ur(r)});
> ur_elast:=subs(_C1=U2,_C2=U3,rhs(op(loesung)));
> sr_elast:=simplify(subs(ur(r)=ur_elast,sr(r)));
> st_elast:=simplify(subs(ur(r)=ur_elast,st(r)));
# Reservierung für später:
> sr_elast_res:=subs(U2=V2,U3=V3,sr_elast);
> st_elast_res:=subs(U2=V2,U3=V3,st_elast);
# Ende der Reservierung
# 
> differenz:=simplify(st_elast-sr_elast);
> plot(1/r^3,r=1..2,0..1);
# 
# Im elastischen Bereich ist der Betrag dieser Spannungsdifferenz
# kleiner als 
# die Vergleichsspannung  sigmav. Der größte Wert ergibt sich für den
# kleinsten Nenner,
# also, wenn die elastische Lösung in der gesamten Kugel gilt, für 
# r=ri. 
# 
# Das Fließen beginnt also immer an der Innenfläche der Hohlkugel.
# 
# 
# Nun soll Teilplastizierung betrachtet werden.
# Das Innengebiet von ri bis rg ist schon plastisch, das Außengebiet von
# rg bis ra noch elastisch.
# 
# Vorsicht mit der Benennung der Konstanten
# 
> loes:=solve(subs(r=rg,differenz)=vorz*sigmav,{U2});
> assign(loes);
# 
# Übergangsbedingungen an der Grenzfläche:
# 
> ueber1:=subs(r=rg,ur_elast-ur_plast);
> ueber2:=subs(r=rg,sr_elast-sr_plast);
> loes:=solve({ueber1=0,ueber2=0},{U1,U3});
> assign(loes);
# 
# 
> sr_plast;
> st_plast;
> ur_plast;
> sr_elast;
> st_elast;
> ur_elast;
# 
# Der Druck auf der Außenfläche ist:
# 
> p_a:=-subs(r=ra,sr_elast);
> deltap:=factor(p_i-p_a);
# 
# Druckdifferenz bei voller Plastizierung (schon oben ermittelt):
# 
> deltap_Ende:=simplify(eval(subs(rg=ra,deltap)));
# 
# Druckdifferenz bei Beginn der Plastizierung:
# 
> deltap_Anfang:=simplify(eval(subs(rg=ri,deltap)));
# Verhältnis:
> q:=simplify(deltap_Ende/deltap_Anfang);
> qq:=simplify(subs(ra=psi*ri,q),symbolic);
# 
# Verhältnis der Druckdifferenzen (Tragreserve)
# in Abhängigkeit vom geometrischen Parameter  psi =  ra/ri
# 
> plot(qq,psi=1..4,0..5);
# 
# Zahlenbeispiel:
# 
> K:=8/3*G; vorz:=1; sigmav:=0.001*G; p_i:=0.8*sigmav;ri:=1; rg:=1.25;
> 
> ur_ges:=ur_plast*Heaviside(rg-r)+ur_elast*Heaviside(r-rg):
> plot(ur_ges,r=1..1.6,0..0.0004);
> sr_ges:=(sr_plast*Heaviside(rg-r)+sr_elast*Heaviside(r-rg))
> /sigmav:
> st_ges:=(st_plast*Heaviside(rg-r)+st_elast*Heaviside(r-rg))
> /sigmav:
> plot([sr_ges,st_ges],r=1..1.6,color=[red,blue]);
# Man entnimmt dem Bild, dass im vorliegenden Falle der Außenrand,
# wenn er spannungsfrei sein soll, etwa bei r=1.6 liegen muss.
# 
> plot(st_ges-sr_ges,r=1..1.6,0..1.1);
# 
# Entlastung
# Die Entlastung geschieht durch Überlagerung eines in der ganzen Kugel
# elastischen Zustandes.
# Einen solchen Spannungszustand, der noch  2 freie Konstanten V2 und V3
# enthält, hatten wir unter
# dem Namen sr_elast_res und st_elast_res reserviert.
# 
# Auch nach Entlastung muss der Außenrand bei r=1.6 spannungsfrei sein,
# außerdem auch der Innenrand r=1.
# 
> rba:=subs(r=1.6,sr_elast_res);
> rbi:=subs(r=ri,sr_elast_res+sr_plast);
> loes:=solve({rba=0,rbi=0},{V2,V3});
> assign(loes);
> sr_ent:=sr_ges+sr_elast_res/sigmav:
> st_ent:=st_ges+st_elast_res/sigmav:
> plot([sr_ent,st_ent],r=1..1.6,color=[red,blue]);
> plot(st_ent-sr_ent,r=1..1.6,-1.1..1.1);
> # Das ist der nach Entlastung verbleibende Eigenspannungszustand.
# 
# Bei Wiederbelastung tritt keine Plastizierung mehr auf.
# Durch einmalige Überlastung kann also der elastische
# Belastungsspielraum vergrößert werden.
# Von dieser Möglichkeit macht man bei Behältern Gebrauch
# (Autofrettage).. 
> 
>