# Kapitel 2: 
# Freie Stabschwingungen bei veränderlicher Querschnittsfläche
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots):
# Beispiel 2-4:
# 
# Der links eingespannte und rechts freie Stab in der obigen Abbildung mit linear veränderlicher Querschnittsfläche wird am rechten Rand zum Zeitpunkt t = 0 um das Maß u0 ausgelenkt und sodann -ohne Anfangsgeschwindigkeit- sich selbst überlassen. Für den dann einsetzenden Schwingungsvorgang sind sämtliche Zustandsgrößen für die Flächenverhältnisse ± = Ar/A =1/ 2 und ± = 2  zu berechnen.
# Wir lösen die ortsabhängige Differenzialgleichung
> ode2:=diff(Z(z),z,z)+1/z*diff(Z(z),z)+kappa^2*Z(z)=0;
> dsolve(ode2); assign(%);
> Z:=unapply(Z(z),z);
# Die Besselfunktionen erster und zweiter Gattung 0-ter Ordnung besitzen folgende Ableitungen:
> Z1str:=simplify(diff(Z(z),z));
> Z2str:=collect(diff(%,z),[_C1,_C2]);
# Wir lassen uns den Verlauf der Besselfunktionen und deren Ableitungen über die Variable  Â = z º von Maple grafisch darstellen:
> J0:= BesselJ(0,zeta): Y0:=BesselY(0,zeta): J1:=BesselJ(1, zeta): Y1:=BesselY(1, zeta): 
> plot([J0,Y0,J1,Y1],zeta=0..10,y=-1..1,title = "\nBesselfunktionen erster und zweiter Gattung\n",titlefont = ["ARIAL", 15],labels = ["Â", "Besselfunktionen"], labeldirections = ["horizontal", "vertical"], linestyle = [solid, dot,dash,spacedash], axesfont = ["HELVETICA", "ROMAN", 8],labelfont = ["HELVETICA", 12],legend = ["J0(Â)", "Y0(Â)","J1(Â)", "Y1(Â)"],legendstyle = [font = ["HELVETICA", 9], location = bottom],color=[black,black,black,black],axes=boxed,gridlines=true);
# Wir beschaffen uns zunächst aus den Randbedingungen die Eigenwerte des Systems. Am linken eingeklemmten Rand bei x = 0 und damit z = 1 bzw. Â = º muss die Verschiebung verschwinden, was Z = C J0(º) + D Y0(º) = 0 erfordert. Das rechte freie Stabende x =  und damit z = ± = Ar / A ist normalkraftfrei. Das erreichen wir wie folgt:
> gl1:=C*BesselJ(0,kappa)+D*BesselY(0,kappa);
> gl2:=C*BesselJ(1,kappa*alpha)+D*BesselY(1,kappa*alpha);
# In Matrizenschreibweise erhalten wir
> A,b:=GenerateMatrix([gl1,gl2],[C,D]);
> EW:=unapply(Determinant(A),kappa,alpha);
> plot([EW(kappa,1/2),EW(kappa,2)],kappa=0..30,y=-0.2..0.2,title = "\nEigenfunktionen\n",titlefont = ["ARIAL", 15],labels = ["º", ""], labeldirections = ["horizontal", "vertical"], linestyle = [solid,dash], axesfont = ["HELVETICA", "ROMAN", 8],labelfont = ["HELVETICA", 12],legend = ["± = 1/2", "± = 2"],legendstyle = [font = ["HELVETICA", 9], location = bottom],color=[black,black],axes=boxed,gridlines=true);
# Zur Berechnung der Nullstellen schreiben wir folgende Prozeduren, wobei  F die Eigenwertgleichung auswertet.
> F:= proc(kappa)
> local alpha;
> description "Eigenfunktion des Stabes mit linear veränderlicher Querschnittsfläche";
> alpha:=1/2;
> evalf(BesselJ(0, kappa)*BesselY(1, alpha*kappa)-BesselY(0, kappa)*BesselJ(1, alpha*kappa));
> end proc:
# In der folgenden Prozedur wird die Funktion F im Intervall [a,b] mit der Schrittweite dx ausgewertet. Ist F(x)*F(x+dx)  < 0, dann hat F in diesem Intervall einen Nulldurchgang. In diesem Fall wird die Nullstelle mit  fsolve  nummerisch konkretisiert. 
> Proc_Calc_03:=proc(a,b,dx)
> local nmax, x, n;
> global Kappa;
> description "Berechnung der Eigenwerte º";
> nmax:=round((b-a)/dx);
> x:= a; Kappa:= [];
> for n to nmax do
>   if F(x)*F(x+dx) < 0 then
>     Kappa:= [op(Kappa),fsolve(F(z),z=x..x+dx)];
>   end if;
> x:= x + dx;
> end do;
> print(`Eigenwerte º =`,Kappa);
> end proc:
> 
;
# Wir rufen die Prozedur zur Berechnung der Eigenwerte º auf:
> Proc_Calc_03(1,150,1):
# Mit den bekannten Eigenwerten werden in der folgenden Prozedur die Eigenkreisfrequenzen É und die Eigenfrequenzen f = É/(2À) berechnet.
> Proc_Calc_05:=proc(rho,E,Al,Ar,l,Kappa)
> local cL,eta;
> global omega,f; 
> description "Berechnung der Eigenkreisfrequenzen É und Eigenfrequenzen f";
> cL:=sqrt(E/rho);
> eta:=Ar/Al-1;
> omega:=cL*sqrt((eta/l)^2)*Kappa;
> f:=omega/(2.*Pi);
> print(`Eigenkreisfrequenzen É = `,omega);
> print(`Eigenfrequenzen f = `,f);
> end proc:
> 
;
> rho:=7850.; E:=2.1*10^11; Al:=0.001; alpha:= 1/2; Ar:= alpha*Al; l:=5.0; u0:=0.001;
> Proc_Calc_05(rho,E,Al,Ar,l,Kappa);
> Z:= proc (kappa, alpha) options operator, arrow;
> BesselY(0, kappa)*BesselJ(0, kappa*(1+(alpha-1)*xi)) - BesselJ(0, kappa)*BesselY(0, kappa*(1+(alpha-1)*xi)) end proc;
> 
;
# Wir geben die Eigenfunktionen Z für die ersten drei Eigenwerte aus
> plot([seq(Z(Kappa[i],alpha),i=1..3)],xi=0..1,title = "\nEigenfunktionen\n",titlefont = ["ARIAL", 15],labels = ["¾ = x/", ""], labeldirections = ["horizontal", "vertical"], linestyle = [solid,dash,dashdot], axesfont = ["HELVETICA", "ROMAN", 8],labelfont = ["HELVETICA", 12],legend = ["º1", "º2", "º3"],legendstyle = [font = ["HELVETICA", 9], location = bottom],color=[black,black,black],axes=boxed,gridlines=true);
> Konst:=proc(alpha,Phi)
> local n,IZ,IN;
> global An;
> description "Berechnung der Konstanten An";
> An:=[];
> for n to nops(Kappa) do
>   IZ:= int((1+(alpha-1)*xi)*Phi*Z(Kappa[n],alpha),xi=0..1,numeric = true);
>   IN:= int((1+(alpha-1)*xi)*Z(Kappa[n],alpha)^2,xi=0..1,numeric = true);
>   An:= [op(An),IZ/IN];
> end do;
> end proc:
> 
;
> Phi:=ln(1+(alpha-1)*xi); #Auslenkung u(x) zum Startzeitpunkt t = 0.
;
> Konst(alpha,Phi):
> An:=evalf(An*u0/ln(alpha));
> u:=0:
> for n to nops(Kappa) do
>    u := u  + evalf(An[n]*Z(Kappa[n],alpha))*cos(omega[n]*t);
> end do:
> v:=diff(u,t):
> N:=E*Al/l*(1+(alpha-1)*xi)*diff(u,xi):
> para:=xi=0..1,t=0..tE,frames=150,thickness=2,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed:
> tE:=4/sqrt(E/rho/l^2);
> #Animation der Zustandsgrößen;
> para:=xi=0..1,t=0..tE,frames=150,thickness=2,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed:
> p1:=plots[animate](u,para,title = "\nAuslenkung in [m]\n",       labels = ["¾", "u(¾,t)"]):
> p2:=plots[animate](v,para,title = "\nGeschwindigkeit in [m/s]\n",labels = ["¾", "v(¾,t)"]):
> p3:=plots[animate](N,para,title = "\nNormalkraftin [N]\n",       labels = ["¾", "N(¾,t)"]):
> plots[display](Matrix(1,3,[p1,p2,p3]));
> 
;
