# 
# Kapitel 2: 
# Der Stab konstanten Querschnitts unter einer harmonischen Randlast F(t)
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
# 
# 
> restart: with(LinearAlgebra): with(plots):
# 
# 
# Der Stab in in der obigen Abbildung wird am rechten Rand durch eine äußere harmonische Kraft zu Schwingungen angeregt.  
# Beispiel 2-7:
> Proc_Calc_03:=proc(rho,E,A,l,F0,Omega,nE,tE)
> #------------------------------------------------------------------------
> #Eingabe: 
> #      rho: Dichte in [kg/m^3]
> #        E: Elastizitätsmodul in [N/m^2]
> #        A: Elastizitätsmodul in [N/m^2]
> #        l: Länge der Saite in [m]
> #       F0: Erregerkraft am Stabende in [N]
> #       ©: Erregerkraftfrequenz in [s-1]
> #       nE: Anzahl der zu berücksichtigenden Reihenglieder
> #       tE: Simulationszeit in [s]
> #Ausgabe:
> #        u: Auslenkung der Stabelemente
> #        v: Geschwindigkeit der Stabelemente
> #        N: Normalkraft
> #------------------------------------------------------------------------
> local cL,kqL,kLn,omn,arg1,arg2,I1qn,etan,SU,eps,n,Bn,kst,up0,para,p1,p2,p3;
> global u,v,N;
> description "Der Stab konstanten Querschnitts unter Randlast";
> cL:=sqrt(E/rho); kqL:=Omega/cL; kLn:=(2*n-1)*Pi/2/l; omn:=kLn*cL; etan:=Omega/omn;
> arg1:=(kqL-kLn)*l; arg2:=(kqL+kLn)*l;
> I1qn:=sin(arg1)/arg1 - sin(arg2)/arg2;
> SU:=0.; eps:=1.*10^(-6);
> for n to nE do
>   if abs(eval(etan)-1.) < eps then   #Resonanzfall
>     Bn:=1;
>   else;
>     Bn:=etan*I1qn;
>   end if;
>   SU:= SU + map(eval,Bn*sin(kLn*x)*sin(omn*t));
> end do; 
> kst:=E*A/l; up0:=F0/(kst*kqL*l*cos(kqL*l));
> u :=subs(x=xi*l,up0*(sin(kqL*x)*sin(Omega*t)-SU));
> v :=diff(u,t);
> N :=E*A/l*diff(u,xi);
> #Ausgabe der Ergebnisse
> para:=xi=0..1,frames=250,thickness=2,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed:
> p1:=plots[animate](plot,[u ,para,title = "\nAuslenkung\n",     labels = ["¾","u(¾,t)"]],t=0..tE):
> p2:=plots[animate](plot,[v ,para,title = "\nGeschwindigkeit\n",labels = ["¾","v(¾,t)"]],t=0..tE):
> p3:=plots[animate](plot,[N ,para,title = "\nNormalkraft\n",    labels = ["¾","N(¾,t)"]],t=0..tE):
> plots[display](Matrix(1,3,[p1,p2,p3]));
> end proc:
# Mit den Werten des Beispiels
> rho:=2100.; A:=0.01; E:=3.*10^(10); l:=8.; F0:=100000.; Omega:=600.; nE:=20;  tE:=0.05; 
# rufen wir die Prozedur auf:
> Proc_Calc_03(rho,E,A,l,F0,Omega,nE,tE);
# Von Interesse sind noch die Normalkräfte an den beiden Stabenden und die Auslenkung am rechten Stabende :
> p4:=plot([subs(xi=0,N),subs(xi=1,N)],t=0..tE,thickness=2,gridlines=true,title = "\nNormalkräfte an den Stabenden\n",labels = ["t", "N(t)"],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 8],axesfont = ["HELVETICA", "ROMAN", 10],legend = ["N(x=0,t)", "N(x=l,t)"],legendstyle = [font = ["HELVETICA", 9], location = bottom],color=[black,blue],axes=boxed,gridlines=true):
> p5:=plot(subs(xi=1,u),t=0..tE,thickness=2,gridlines=true,title = "\nAuslenkung des rechten Stabendes\n",labels = ["t", "u(1,t)"],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 8],axesfont = ["HELVETICA", "ROMAN", 10],axes=boxed,gridlines=true):
> plots[display](Matrix(1,2,[p4,p5]));
> 
;
