> restart:
# HF, Übung 5: Rechteckplatte, isotrop-elastisch
# 
# Schüttkraft       p= pmax*sin(Pi*x/a)*sin(Pi*y/b)
# 
# 
> 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: 
# 
# Verformungsansätze
# 
> w:=W*sin(Pi*x/a)*sin(Pi*y/b);
> omegax:=Psiy*cos(Pi*x/a)*sin(Pi*y/b);
> omegay:=-Psix*sin(Pi*x/a)*cos(Pi*y/b);
> omega:=vector([omegax,omegay,0]);
> grad_o:=vecgrad(omega,[x,y,z]);
# 
# Biegeverzerrung:
# 
> bm:=evalm((grad_o+transpose(grad_o))/2);
> idp:=diag(1,1,0); # Matrix des planaren Einheitstensors
# 
# Momententensor
# 
> m:=evalm(G*h^3/6*(bm+nu/(1-nu)*trace(bm)*idp));
> simplify(m[1,1]);
> simplify(m[2,2]);
> simplify(m[1,2]);
# 
# Planarer Querkraftoperator:
# 
> q:=tensdiv(m,[x,y,z]);
> div_q:=simplify(diverge(q,[x,y,z]));
# 
# Maximalwert der Schüttkraft:
# 
> pmax:=factor(eval(subs(x=a/2,y=b/2,-div_q)));
# 
# Schubwinkel:
# 
> gam:=evalm(grad(w,[x,y,z])+omega);
# 
# Verknüpfung der Drehungen mit den Durchbiegungen:
# 
> glg:=evalm(q-G*h*gam);
> map(factor,glg);
> c1:=eval(subs(x=0,y=b/2,simplify(glg[1])));
> c2:=eval(subs(x=a/2,y=0,simplify(glg[2])));
> loes:=solve({c1=0,c2=0},{Psix,Psiy});
> assign(loes);
> simplify(pmax);
# 
# Ergebnis für die Durchbiegung:
# 
> loes:=solve({pmax=pgeg},{W});
> assign(loes);
# 
# Ergebnisse für die Biege- und Torsionsmomente:
# 
# Diese hängen nicht von der Plattendicke, wohl aber von der
# Querkontraktionszahl ab.
# 
> simplify(m[1,1]);
> simplify(m[2,2]);
> simplify(m[1,2]);
# 
# Definition einer (dimensionslosen) bezogenen Durchbiegung in
# Plattenmitte:
# 
> W_dim:=W*G*h^3/pgeg/a^4;
# Durchbiegung der Quadratplatte:
# 
# Abhängigkeit von Querkontraktionszahl und Plattendicke
# 
> b:=a; h:=zeta*a;
> plot([seq(W_dim,nu=[0,1/3,1/2])],zeta=0..0.3,
> 0..0.02,color=[red,blue,green]);
> b:='b'; h:='h';
# 
# Zum Vergleich der Größenordnung der Biegemomente:
#  Streifen parallel zur x-Achse (Balkenlösung)
# 
> dgl:=diff(Mb(x),x,x)=-pgeg*sin(Pi*x/a)*sin(Pi*y/b);
> loes:=dsolve({dgl,Mb(0)=0,Mb(a)=0},{Mb(x)});
> M:=rhs(loes);
# Verhältnis der Plattenmomente zu den Streifenmomenten:
# 
# Abhängigkeit vom Seitenverhältnis der Platte
# 
> verh:=simplify(m[1,1]/M);
> b:=lambda*a;
> ver:=simplify(verh);
> plot([seq(ver,nu=[0,1/3,1/2])],lambda=0..5,
> color=[red,blue,green]);
# 
# Sonderfall: Quadratplatte
# 
> ver_quad:=subs(lambda=1,ver);
> 
>