# Kapitel 2: 
# Erzwungene Stabschwingung
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart:with(linalg):with(plots):
# Beispiel 2-12:
# 
# Der in der obigen Abbildung skizzierte Stab, der am linken Ende frei gelagert und rechts an einen (k-m-c)-Schwinger gekoppelt ist, wird durch die pulsierende Längskraft  F(t)  belastet. Gesucht werden die Auslenkung und Geschwindigkeit im eingeschwungenen Zustand. Animieren Sie sämtliche Zustandsgrößen. 
> 
;
> Proc_Calc_03:=proc(rho, E, A, l, kF, m , c, F0, omega, phiF)
> #---------------------------------------------------------------------------------
> # Eingabe Stab:
> #      rho: Dichte in [kg/m^3]
> #        E: Elastizitätsmodul in [N/m^2]
> #        A: Flächeninhalt des Stabquerschnitts in [m^2]
> #        l: Länge des Stabes in [m]
> # Eingabe (kF-m-c)-Schwinger:
> #       kF: Federsteifigkeit in [N/m]
> #        m: Masse in [kg]
> #        c: Dämpfungskonstante in [kg/s]
> # Eingabe Erregerkraft F(t) = F0*cos(omega*t+phiF)
> #       F0: Erregerkraftamplitude
> #    omega: Erregerkreisfrequenz
> #     phiF: Nullphasenverschiebungswinkel
> #  Ausgabe: Sämtliche Systemwerte und Zustandsgrößen
> #----------------------------------------------------------------------------------
> local cL,kappaL,lambda,ZL,Zl,Z0,zeta,vdx,Ndx,Nd0,udx,para,T,p1,p2,p3,par,p4,p5,p6;
> global u,v,N;
> description "";
> #Berechnung der Systemwerte
> cL:=sqrt(E/rho); 
> kappaL:=omega/cL;
> lambda:=kappaL*l;          
> ZL:=A*sqrt(rho*E);
> Zl:=I*(m*omega-kF/omega) + c;         
> zeta:=Zl/ZL;
> Z0:=-ZL*(zeta*cos(lambda) + I*sin(lambda))/(cos(lambda)+I*zeta*sin(lambda));
> vdx:=Nd0*(cos(lambda*xi)/Z0 + I/ZL*sin(lambda*xi));
> Ndx:=Nd0*(cos(lambda*xi)+I*ZL/Z0*sin(lambda*xi));
> Nd0:=F0*exp(I*phiF);
> udx:=-I/omega*vdx;
> u:=unapply(evalc(Re(udx*exp(I*omega*t))),xi,t);
> v:=unapply(evalc(Re(vdx*exp(I*omega*t))),xi,t);
> N:=unapply(evalc(Re(Ndx*exp(I*omega*t))),xi,t);
> # Ausgabe der Systemwerte und Zustandsgrößen
> print(`cL = `,cL);
> print(`ZL = `,ZL);
> print(`Zl = `,Zl);
> print(`Z0 = `,Z0);
> print(`u  = `,expand(u(xi,t)));
> print(`v  = `,expand(v(xi,t)));
> print(`N  = `,expand(N(xi,t)));
> #Animation der Zustandsgrößen
> T:=2*Pi/omega;
> para:=xi=0..1,frames=250,thickness=2,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "horizontal"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed:
> p1:=plots[animate](plot,[u(xi,t) ,para,title = "\nAuslenkung\n"             ,labels = ["¾ = x/","u(¾,t)"]],t=0..2*T):
> p2:=plots[animate](plot,[v(xi,t) ,para,title = "\nTeilchengeschwindigkeit\n",labels = ["¾ = x/","v(¾,t)"]],t=0..2*T):
> p3:=plots[animate](plot,[N(xi,t) ,para,title = "\nNormalkraft\n"            ,labels = ["¾ = x/","N(¾,t)"]],t=0..2*T):
> par:=t=0..2*T,thickness=2,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "horizontal"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,color=[black,blue],legend = ["¾ = 0", "¾ = 1"],legendstyle = [font = ["HELVETICA", 9], location = bottom]:
> p4:=plot([u(0,t),u(1,t)],par,title = "\nAuslenkungen der Stabenden\n"     ,labels = ["t","u"]):
> p5:=plot([v(0,t),v(1,t)],par,title = "\nGeschwindigkeiten der Stabenden\n",labels = ["t","v"]):
> p6:=plot([N(0,t),N(1,t)],par,title = "\nNormalkräfte an den Stabenden\n"  ,labels = ["t","N"]):
> plots[display](Matrix(2,3,[p1,p2,p3,p4,p5,p6]));
> end proc:
> 
;
# Mit den Werten des Beispiels:
> rho:=7850.   : E:=2.1*10^(11): A:=7.57*10^(-4): l:=10.0:   #Eingabewerte Stab
;
> kF :=7.9485E6: m:= 59.4245:    c:=30735.50514:             #Eingabewerte (kF-m-c)-Schwinger
;
> F0:=100000.: omega:=500: phiF:=-Pi/2:                      #Eingabewerte für die Erregerkraft
;
> Proc_Calc_03(rho, E, A, l, kF, m , c, F0, omega, phiF);
> 
;
