# Kapitel 5: 
# Die Rechteckmembran mit Anfangsgeschwindigkeit ¨(x,y)
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(plots): with(Student[MultivariateCalculus]): unprotect(gamma):
# Beispiel 5-1
# 
# 
# Eine rechteckige Membran wird am Punkt P(x = u,y = v) in der Ruhelage zum Zeitpunkt  t = 0 von einem starren Körper getroffen und so zu transversalen Schwingungen angeregt. Aus Messungen ergibt sich die zweidimensionale Verteilung der Anfangsgeschwindigkeit ¨(x,y) . Es sind der Ausschlag w(x,y,t) und die Geschwindigkeit diff(w(x, y, t), t);zu berechnen. 
# Mit den Werten des Beispiels sind 
> n0:= 2.; mue:= 0.005; cM:=sqrt(n0/mue);
> a:=0.3; b:=0.4; u:=a/3; v:=b/4; c:=a/20; d:= b/30;  
> alpha:=a/b; mu:=Pi*u/a; gamma:=Pi*c/a; nu:=Pi*v/b; delta:=Pi*d/b;
>  PSI:= 0.: me:= 2; ne:= me; v0:= 500.;
# Festlegung der Anfangsgeschwindigkeit (s.h. die Darstellung in der Aufgabe)
> for m to me do
>   for n to ne do
>     amn:= sin(m*mu)*sin(m*gamma)*sin(n*nu)*sin(n*delta); 
>     kmn:= m*n*(m^2+(alpha*n)^2);
>     PSI:= PSI + amn/kmn*sin(m*Pi*x/a)*sin(n*Pi*y/b):
>   end do:
> end do:
> PSI:=evalf(v0*PSI);
# Grafische Ausgabe der Anfangsgeschwindigkeit
> plot3d(PSI , x=0..a, y=0..b,axes=normal,orientation=[130, -65, 178], style = hidden, scaling = unconstrained,color = black);
> contourplot(PSI, x=0..a,y=0..b,axes=boxed,contours = 10, view = [0..a, 0..b], gridlines);
# Bestimmung der Integrationskonstanten Bmn
> wB:=0.; me:= 2; ne:= me;
> for m from 1 to me do
>  for n from 1 to ne do
>     omegamn:=cM*Pi/a*sqrt(m^2 + (alpha*n)^2);
>     Bmn:= 1./omegamn*MultiInt(PSI*sin(m*Pi*x/a)*sin(n*Pi*y/b), y = 0..b, x = 0..a);
>      wB:= wB + Bmn*sin(m*Pi*x/a)*sin(n*Pi*y/b)*sin(omegamn*t);
>   end do;
> end do;
> wB:=unapply(4.*wB/a/b,x,y,t);  #Auslenkung
;
> vB:=diff(wB(x, y, t),t);       #Geschwindigkeit
;
# Animation der Bewegungsgrößen
> te:=0.1:
> p1:=animate(plot3d, [wB(x,y,t), x = 0 .. a, y = 0 .. b], t = 0 .. te,orientation=[130, -65, 178], style = patchcontour,axes=framed,frames = 101):
> p2:=animate(plot3d, [vB,        x = 0 .. a, y = 0 .. b], t = 0 .. te,orientation=[130, -65, 178], style = patchcontour,axes=framed,frames = 101):
> plots[display](Matrix(1,2,[p1,p2]));
# Berechnung von Contourplots
> maxi:=3; t:= 0; Dt:=0.04; HL:=Matrix(1,maxi);
> for i to maxi do
>   t:= t + Dt;
>   HL(1,i):=contourplot(wB(x,y,t), x=0..a,y=0..b,title =typeset("\nVerschiebungsfeld zum Zeitpunkt t = ", t," s\n"), titlefont = [], axes=boxed,contours = 10, view = [0..a, 0..b], gridlines=true,titlefont = ["ARIAL", 15],labelfont = ["Helvetica", 12]):
> end do:
> display(HL);
# Rotierende Knotenlinien für eine quadratische Membran (s.h. Kap. 5.2)
> a:=1; b:=1; alpha:=a/b; m:='m'; n:='n'; t:='t'; cM;
> omega:=unapply(cM*Pi/a*sqrt(m^2+(alpha*n)^2),m,n);
> W:=unapply(sin(m*Pi*x/a)*sin(n*Pi*y/b),m,n);
> w1:=W(m,n)*cos(omega(m,n)*t); w2:= W(n,m)*sin(omega(m,n)*t);
> w:=unapply(w1+w2,m,n);
# Wir prüfen, ob diese Lösung die Bewegungsdifferenzialgleichung erfüllt
> simplify(diff(w(m,n),t,t) - cM^2*(diff(w(m,n),x,x) + diff(w(m,n),y,y)));
> m:=2;n:=3;
> p3:=animate(plot3d, [w(m,n), x = 0 .. a, y = 0 .. b], t = 0 .. te, contours=[0], orientation=[130, -65, 178], style = patchcontour,axes=framed,frames = 101):
> p4:=animate(plot3d, [w(m,n), x = 0 .. a, y = 0 .. b], t = 0 .. te, contours=[0], orientation=[ 90, 0, 0], style = patchcontour,axes=framed,frames = 101):
> plots[display](Matrix(1,2,[p3,p4]));
> 
;
