# Kapitel 5: 
# Freie Schwingungen einer indischen Trommel (Tabla)
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
> restart: with(LinearAlgebra): with(plots):
> addcoords(ZK,[z,r,phi],[r*cos(phi),r*sin(phi),z]);
# Beispiel 5-6:
# 
# Die indische Trommel Tabla, die auch Dayan genannt wird, besteht in ihrer Bespannung aus drei übereinandergeschichteten Ziegenfellen, wobei untere und obere Haut im Zentrum kreisringförmig ausgeschnitten sind. Dieser Ausschnitt wird mehrfach mit Stimmpaste bestrichen, poliert und anschließend kontrolliert gebrochen. Als Modell zur Berechnung der Eigenfre-quenzen der freien Schwingungen gehen wir von einer Membran mit bereichsweise konstanten Massenverteilungen aus. Im mittleren Bereich der Bespannung ist ¼ = ¼1, im Außenbereich dagegen ¼ = ¼2 . Die Membran ist mit der konstanten Kraft n0 vorgespannt. In Ermangelung realer Materialdaten für die Dichten werden diese fiktiv angenommen. 
# Für die Massenverhältnisse ¼2/¼1 = (0.5, 1.0, 2.0) sind die Eigenfrequenzen zu berechnen.
# Geg.: a = 0,15 m, b = 0,075 m, ¼1 = 0,2 kg/m2, n0 = 1200 N/m. 
# Lösung: Die in der Umgebung von r = b unstetig verlaufende Massenbelegung ¼ hat zur Folge, dass die Lösung bereichsweise beschafft werden muss. An der Stelle r = b sind Übergangsbedingungen zu formulieren, die beide Lösungen zusammenführen. 
> R :=unapply(A*BesselJ(n,lambda*r)+B*BesselY(n,lambda*r),A,B,n,lambda,r);
# Wir benötigen zur Erfüllung einer Übergangsbedingung an der Stelle r = a die  erste Ableitung von R nach r
> DR:=D[5](R);
# Lösung im Bereich I (0 d r < b)
> R1 :=unapply( R(A,0,n,lambda1,r),A,n,lambda1,r);
> DR1:=unapply(DR(A,0,n,lambda1,r),A,n,lambda1,r);
# Lösung im Bereich II (b < r < a)
> R2 :=unapply( R(B,C,n,lambda2,r),B,C,n,lambda2,r);
> DR2:=unapply(DR(B,C,n,lambda2,r),B,C,n,lambda2,r);
# Rand- und Übergangsbedingungen: 
# 1.) Bei r = a ist die Membran eingespannt: 
> gl1:=R2(B,C,n,kappa*zn/a,a)=0;
# Übergangsbedingungen bei r = b
# 2.) Gleiche Auslenkung der Membran an der Übergangsstelle
> gl2:= R1(A,n,zn/a,beta*a) = R2(B,C,n,kappa*zn/a,beta*a);
# 3.)  Gleiche Vertikalkomponenten V der schrägen Normalkraft n0
> gl3:=(DR1(A,n,zn/a,beta*a) = DR2(B,C,n,kappa*zn/a,beta*a));
# Aufbau der Koeffizientenmatrix AM
> gl:=[gl1,gl2,gl3]: unb:=[A,B,C];
> AM,c:=GenerateMatrix(gl,unb);
# Determinante der Koeffizientenmatrix
> DAM:=simplify(Determinant(AM));
> lambda1:=omega*sqrt(mu1/n0);
> lambda2:=omega*sqrt(mu2/n0);
# Mit den Werten des Beispiels:
> a:=0.15; b:=a/2;  mu1:= 0.2; mu2:= 0.5*mu1; kappa:=sqrt(mu2/mu1); beta:=b/a; n0:=1200.;
> lambda1; lambda2;
# Rotationssymmetrische Lösungen (n = 0: Ordnung der Besselfunktionen,  m: Nummer der Nullstelle)
> plot([seq(DAM,n=0..2)],zn=0..20,y= -5..10,title = "\nEigenwertgleichungen (mu2/mu1 = 1/2)\n", legend = ["n=0", "n=1","n=2"], titlefont = ["ARIAL", 15], labels = ["zn", ""], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], linestyle = [solid, longdash, dashdot], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],axes = boxed, gridlines=true);
# Wir berechnen für die Ordnungszahlen n = (0,1,2) jeweils genau 5 Nullstellen zn. In der ersten Spalte stehen die Nullstellen für die Ordnungszahl n = 0.
> ZN:=Matrix(5,3):
> for n from 0 to 2 do
> loe:=fsolve(DAM,zn=1..25,maxsols = 5):
>     for i to 5 do
>       ZN[i,n+1]:=loe[i];
>     end do:
> end do:
> ZN;
> OMEGAN:=sqrt(n0/mu1)/a*ZN; #Eigenkreisfrequenzen
;
> FREQN:=OMEGAN/2./Pi;       #Eigenfrequenzen
;
> 
;
