# Kapitel 5: 
# Freie gedämpfte Schwingung einer allseits eingespannten Rechteckmembran
# mit Anfangsauslenkung ¦(x,y)
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(plots): with(Student[MultivariateCalculus]): unprotect(gamma):
# Beispiel 5-2
# 
# Für eine allseits eingespannte Membran sind die gedämpften Schwingungen w(x,y,t) zu ermitteln, wenn sich die Membran in einem medium der Dichte ñ bewegt. Zum Zeitpunkt t = 0 besitzt die Membran die Anfangsauslenkung ¦(x,y). Der Bewegungsvorgang ist zu animieren. 
# Die Eigenkreisfrequenzen  Émn  einer allseits eingespannten Rechteckmembran sind:
> cM:=sqrt(n0/mu); omega:=unapply(cM*Pi/a*sqrt(m^2+(alpha*n)^2),m,n);
# Mit den Werten des Beispiels sind:
> mu:=0.05; rho:= 0.05; a:=0.4; b:=0.3 ; n0:=2.; u:=a/3; v:=b/4; c:=a/20; d:= b/30;  
> alpha:=a/b; mue:=Pi*u/a; gamma:=Pi*c/a; nu:=Pi*v/b; deltav:=Pi*d/b;
# 
>  Phi:= 0.: me:= 2; ne:= me; w0:=1.;
# Festlegung der Anfangsauslenkung
> for m to me do
>   for n to ne do
>     Phi:= Phi + sin(m*mue)*sin(m*gamma)*sin(n*nu)*sin(n*deltav)*sin(m*Pi*x/a)*sin(n*Pi*y/b):
>   end do:
> end do:
> Phi:=evalf(w0*Phi);
# Grafische Ausgabe der Anfangsgeschwindigkeit
> p1:=plot3d(Phi , x=0..a, y=0..b,axes=normal,orientation=[-50, 65, 0], style = surfacecontour, scaling = unconstrained):
> p2:=contourplot(Phi, x=0..a,y=0..b,axes=boxed,contours = 10, view = [0..a, 0..b], gridlines,scaling = constrained):
> plots[display](Matrix(1,2,[p1,p2]));
# Bestimmung der Integrationskonstanten Amn und Bmn
> delta:=rho/2/mu; m:='m'; n:='n';
> zeta1:=unapply(-delta + sqrt(delta^2-omega(m,n)),m,n);
> zeta2:=unapply(-delta - sqrt(delta^2-omega(m,n)),m,n);
> wB:=0.;
> for m from 1 to me do
>  for n from 1 to ne do
>     omegamn:=cM*Pi/a*sqrt(m^2 + (alpha*n)^2);
>     Wmn    :=sin(m*Pi*x/a)*sin(n*Pi*y/b);
>     I1mn   :=MultiInt(Phi*Wmn, y = 0..b, x = 0..a);
>     dzeta  :=zeta1(m,n)-zeta2(m,n);
>     Amn    :=-I1mn*zeta2(m,n)/dzeta;
>     Bmn    := I1mn*zeta1(m,n)/dzeta;
>      wB    := wB + Wmn*(Amn*exp(zeta1(m,n)*t) + Bmn*exp(zeta2(m,n)*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:=1: #Simulationszeit 1 Sekunde
;
> p3:=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):
> p4:=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,[p3,p4]));
# 
# Berechnung von Contourplots
> maxi:=3; t:= 0; Dt:=te/maxi; 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("\nContourplot 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);
> 
;
