> restart:
# HF, Übung 4: Bohrung im unendlichen Medium, Kerbfaktor
# 
# 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);
> e:=matrix(3,3): # Matrix des Verzerrungstensors
# 
# Spannung, auf  2G  bezogen,nach dem isotropen Hooke'schen Gesetz:
# 
> t:=evalm(e+nu/(1-2*nu)*trace(e)*id);
# 
# Ebener Verschiebungszustand in Zylinderkoordinaten:
# 
> u:=vector([ur(r,phi),uphi(r,phi),0]);
> alias(zyl=([r,phi,z],coords=cylindrical)):
> grad_u := vecgrad(u,zyl);
# 
# Verzerrungen, ausgedrückt durch die Verschiebungsableitungen:
# 
> e:=evalm((grad_u + transpose(grad_u))/2);
# 
# Gleichgewicht:
# 
> glgw:=tensdiv(t,zyl):
# 
# Verschiebungsansatz unter Beachtung der Symmetrie des Problems:
# Beide Verschiebungskomponenten haben die Periode  Pi, die
# Radialverschiebung ist bezüglich
# der Achsen x und y symmetrisch, die Tangentialverschiebung
# antimetrisch
> ur:=(r,phi)-> a0(r)+ a2(r)*cos(2*phi);
> uphi:=(r,phi)->b2(r)*sin(2*phi);
# 
# Auswirkung auf die drei Kraft-Gleichgewichtsbedingungen:
# 
> glgw[3];
> collect(glgw[1],cos);
# 
# Diese Bedingung hat die Gestalt   c12(r)*cos(2*phi) + c11(r) = 0  mit
# 
> c12:=coeff(glgw[1],cos(2*phi)):
> c11:=eval(subs(phi=Pi/4,glgw[1])):
> collect(glgw[2],sin);
# 
# Diese Bedingung hat die Gestalt   c22(r)*sin(2*phi) = 0  mit
# 
> c22:=coeff(glgw[2],sin(2*phi)):
# 
# Lösung der drei gekoppelten Differentialgleichungen   c11(r) = 0, 
# c12(r) = 0,  c22(r) = 0:
# 
> loesneu:=dsolve({c11=0,c12=0,c22=0},{a0(r),a2(r),b2(r)});
> assign(loesneu);
> a0:=unapply(a0(r),r);
# 
# Wegen der Beschränktheit der Verzerrungen im Unendlichen dürfen die
# Verschiebungen nur mit der
# ersten Potenz von r beim Übergang ins Unendliche anwachsen. Daher
# müssen gewisse Konstanten 
# in a2 und b2 gleich 0 sein. Ohne vorher ihre Nummer zu kennen,
# erreicht man das so:
# 
> li:=limit(a2(r)/r^3,r=infinity);
> a2:=unapply(simplify(a2(r)-li*r^3),r);
> li:=limit(b2(r)/r^3,r=infinity);
> b2:=unapply(simplify(b2(r)-li*r^3),r);
# 
# Spannungszustand im Unendlichen, also praktisch in hinreichender
# Entfernung von der Bohrung:
# 
> t11inf:=limit(t[1,1],r=infinity);
> t22inf:=limit(t[2,2],r=infinity);
> sigmax0inf:=eval(subs(phi=0,t11inf));
> sigmay0inf:=eval(subs(phi=0,t22inf));
> sigmay90inf:=eval(subs(phi=Pi/2,t11inf));
> sigmax90inf:=eval(subs(phi=Pi/2,t22inf));
# Kontrolle:
> simplify(sigmax0inf-sigmax90inf);
> simplify(sigmay0inf-sigmay90inf);
# Vorsicht: Nummerierung der Konstanten:  
> loes:=solve({sigmax0inf=sx,sigmay0inf=sy},{_C3,_C5});
> assign(loes);
> t[1,1];
> t[1,2];
> t[2,2];
# 
# Bestimmung der verbliebenen drei Integrationskonstanten aus
# Spannungsfreiheit am Lochrand:
# 
# Vorsicht: Nummerierung der Konstanten:  
> loes:=solve({eval(subs(r=a,phi=0,t[1,1]))=0,eval(subs(r=a,phi=Pi/2,t[1
> ,1]))=0,eval(subs(r=a,phi=Pi/4,t[1,2]))=0},{_C1,_C4,_C6});
> assign(loes);
# 
# Ergebnis:Umfangsspannung am Lochrand:
# 
> sigmarand:=simplify(eval(subs(r=a,t[2,2])));
# 
# Sonderfall: allseitiger Zug (im Unendlichen):
# 
> plot(subs(sx=1,sy=1,sigmarand),phi=0..Pi/2);
# Abgelesen: Kerbfaktor  2
# 
# Sonderfall: einachsiger Druck in x-Richtung  (im Unendlichen):
# 
> plot(subs(sx=-1,sy=0,sigmarand),phi=0..Pi/2);
# Abgelesen: Kerbfaktor   3
# 
# Außerdem Spaltzugspannung bei phi=0 !!
# 
# 
# Im Folgenden die Spannungsverteilungen und Verschiebungen für den Fall
# des einachsigen Druckes:
# 
# 
> sy:=0; sx:=-1; a:=1; 
> ta:=simplify(eval(subs(phi=0,t[1,1])));
> tc:=simplify(eval(subs(phi=Pi/2,t[1,1])));
> # 
> plot([ta,tc],r=1..8,color=[red,blue]); # Radialspannungen
> td:=simplify(eval(subs(phi=0,t[2,2])));
> te:=simplify(eval(subs(phi=Pi/2,t[2,2])));
> # 
> plot([td,te],r=1..8,color=[blue,red]); # Tangentialspannungen
> tb:=simplify(eval(subs(phi=Pi/4,t[1,2])));
> # 
> plot(tb,r=1..8,color=black); # Schubspannungen
> # 
> uuur:=combine(ur(r,phi),trig);
> # 
> uur:=collect(uuur,cos);
> # 
> uuphi:=simplify(uphi(r,phi));# 
# 
# Man erkennt: Die Verschiebungen hängen von der Querkontraktionszahl 
# nu ab, die Spannungen nicht
# 
# Auswertung der Verschiebungen für  nu=0.3:
# 
> nu:=0.3;
> uur;
> plot(subs(phi=0,uur),r=1..5,-4..0);
> plot(subs(phi=Pi/2,uur),r=1..5,0..2);
> uuphi;
> plot(subs(phi=Pi/4,uuphi),r=1..5,0..3);
> 
>