# 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_b: (Rechteckplatte allseits gelenkig gelagert)
# 
# Für die freien ungedämpften Transversalschwingungen einer allseits gelenkig gelagerten schubstarren Rechteckplatte ermittle man die Eigenkreisfrequenzen ohne Berücksichtigung der Drehträgheit nach der Methode von M. Lévy (1899). 
# Die zeitfreien Ortsfunktionen Y und X sind 
> Y:=C1*cosh(kappa[k]*y)+C2*sinh(kappa[k]*y)+C3*cos(eta[k]*y)+C4*sin(eta[k]*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);
# 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(W(x,b));gl4:=eval(Myy(x,b));
# Berechnung der Systemdeterminante (Frequenzgleichung)
> gl:=[gl3,gl4]: unb:=[C2,C4];
> A,c:= GenerateMatrix(gl,unb);
> DA:=simplify(Determinant(A));
> B:=simplify(GaussianElimination(A));
# Wegen B[2,2] = 0 setzen wir 
> C2:=0; C4:=1; 
> alpha[j]:=j*Pi/a; lambda:= mu/N*omega^2; eta[k]:=sqrt(lambda^2-alpha[j]^2);
> loe:=solve(eta[k]*b-k*Pi,omega);
> omega:=unapply(simplify(loe[1]),j,k);
> eta[k]:=simplify(subs(omega=omega(j,k),eta[k]),sqrt, symbolic);
# 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));
> W:=unapply(eval(X*Y),j,k);
# Wir untersuchen die kleinste Eigenkreisfrequenz:
> j:=1; k:=1;
# Das liefert die Grundkreisfrequenz
> omega(j,k);
# Wir stellen die zugehörige zeitfreie Eigenfunktion grafisch dar
> plot3d(W(j,k),x=0..a,y=0..b,orientation=[135, -65, 180],axes=framed);
# und animieren den Bewegungsvorgang 
> T:=sin(omega(j,k)*t);
> w:= W(j,k)*T;
> animate(plot3d, [w,x = 0 .. a,y = 0 .. b],t = 0 .. 1.,orientation=[55, 75, 0],style = patchcontour,axes=framed,frames = 51);
> 
;
