# Kapitel 6: 
# Freie ungedämpfte Schwingungen einer schubstarren Rechteckplatte
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
> restart: with(LinearAlgebra): with(plots): unprotect(gamma); 
# Beispiel 6-1_a: 
# 
# Für die freien ungedämpften Transversalschwingungen einer an den beiden gegenüberliegenden Rändern x = 0 und x = a gelenkig gelagerten und sonst gelenkig bzw. freien schubstarren Rechteckplatte (s.h. obige Abbildung) ermittle man die Eigenkreisfrequenzen ohne Berücksichtigung der Drehträgheit nach der Methode von M. Lévy (1899). 
# Schubstarre Platte ohne Berücksichtigung der Drehträgheit. Die Ränder x = 0 und x = a sind gelenkig gelagert. 
# Die zeitfreien Ortsfunktionen Y und X sind:
> Y:=C1*cosh(kappa[j]*y)+C2*sinh(kappa[j]*y)+C3*cos(eta[j]*y)+C4*sin(eta[j]*y);
> X:=sin(alpha[j]*x);
# Für die zeitfreie Verschiebung folgt dann:
> W:=unapply(X*Y,x,y);
# Wir benötigen die folgenden partiellen Ableitungen der Funktion W(x,y):
> W2x:=diff(W(x,y),x$2);
> W2y:=diff(W(x,y),y$2);
> Wxy:=diff(W(x,y),x,y);
> DW :=simplify(W2x+W2y);
# Schnittlasten (M: Momente, Q: Querkräfte)
> Mxx:=unapply(-N*(W2x+nu*W2y),x,y);
> Myy:=unapply(-N*(W2y+nu*W2x),x,y);
> Mxy:=unapply(-(1-nu)*N*Wxy,x,y);
> Qx :=unapply(-N*diff(DW,x),x,y);
> Qy :=unapply(-N*diff(DW,y),x,y);
# Für die Randquerkräfte erhalten wir:
> Qqx:=unapply(-N*diff(W2x+(2-nu)*W2y,x),x,y);
> Qqy:=unapply(-N*diff(W2y+(2-nu)*W2x,y),x,y);
# 1. Randbedingungen für den Rand y = 0:
> gl1:= W(x,0); gl2:=simplify(Myy(x,0));
> solve([gl1,gl2],[C1,C3]); assign(%);
# 2. Randbedingungen für den Rand y = b:
> gl3:=eval(Myy(x,b)); gl4:=eval(Qqy(x,b));
# Berechnung der Systemdeterminante
> gl:=[gl3,gl4]:unb:=[C2,C4];
> A,c:= GenerateMatrix(gl,unb);
> DA:=simplify(Determinant(A));
> B:=simplify(GaussianElimination(A));
> DB:=simplify(Determinant(B));
# Die Konstanten
> C4:=-1; C2:=-C4*B[1,2]/B[1,1];
# Damit ist die zeitfreie Verschiebung
> W(x,y);
> #simplify(Myy(x,b));simplify(Myy(x,0)); simplify(Qqy(x,b));
# Wir führen noch folgende Abkürzungen ein:
> lambda:=z/b; alpha[j]:=gamma[j]/b;
> eta[j]:=sqrt(lambda^2-alpha[j]^2);kappa[j]:=sqrt(lambda^2+alpha[j]^2);
# Die Determinate DB (Frequenzgleichung) wird zur besseren nummerischen Bearbeitung noch manuell umgeformt. Wir erhalten nach einigen Rechenschritten:
> v:=(((1-nu)*gamma[j]^2+z^2)^2*sqrt(z^2-gamma[j]^2))/(((nu-1)*gamma[j]^2+z^2)^2*sqrt(z^2+gamma[j]^2));
> EWG:=-tanh(sqrt(z^2+gamma[j]^2))*cos(sqrt(z^2-gamma[j]^2))+sin(sqrt(z^2-gamma[j]^2))/v;
# Mit den Werten des Beispiels:
> a:=0.6; b:=0.4; beta:=b/a; h:=0.0014; E:=2.1E11; nu:=0.3; rho:=8400.0; mu:=rho*h;
> N:=E*h^3/(12*(1-nu^2));
> gamma[j]:=Pi*beta;
> EWG;
# Wir stellen die Frequenzgleichung grafisch dar:
> plot(EWG,z=0..15,y=-3..2,axes=boxed,gridlines);
# Der ersten fünf  Eigenwerte sind:
> zn:=[fsolve(EWG,z=2..15,maxsols=5)]; 
# Wir beschränken uns im Folgenden auf den kleinsten Eigenwert. Damit sind:
> z:=zn[1];lambda:=z/b; omega:=(N*lambda^4/mu)^(1/2);
> f:=omega/2/Pi;
# Für die Konstante C2 errechnen wir:
> C2;
# Damit haben wir die zeitfreie Auslenkung:
> W:=eval(X*Y);
> plot3d(W,x=0..a,y=0..b,orientation=[135, -65, 180],axes=framed);
> T:=sin(omega*t);
> w(x,y,t):= W*T;
> animate(plot3d, [w(x,y,t), x = 0 .. a, y = 0 .. b], t = 0 .. 0.15,orientation=[55, 75, 0],style = patchcontour,axes=framed,frames = 101);
> 
;
