# 
# Kapitel 2: 
# Schwingungsverhalten eines Stabes aus Kelvin-Material
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(plots): unprotect(gamma):
# Beispiel 2-5:
# 
# 
# Der Stab in der obigen Abbildung wird zum Zeitpunkt t = 0 an der Stelle x = a um das Maß u0 ausgelenkt und dann sich selbst überlassen. Gesucht werden sämtliche Zustandsgrößen. Animieren Sie den Bewegungsvorgang. 
> X:=sin(kappa*l*xi); T:=An*exp(zeta1*t)+Bn*exp(zeta2*t);
> un:=X*T;
# Berechnung der Konstanten I1
> I1:=u0/alpha*l*int(xi*X,xi=0..alpha) + u0*l*int(X,xi=alpha..1);
> gl:=cos(kappa*l)=0;
> I1:=simplify(I1,{gl});
> I2:=0;
> An:= 2/l*(I2-I1*zeta2)/(zeta1-zeta2);
> Bn:=-2/l*(I2-I1*zeta1)/(zeta1-zeta2);
# Berechnung der charakteristischen Exponenten
> zeta1:=-delta + sqrt(delta^2-omega^2); zeta2:=-delta - sqrt(delta^2-omega^2);
> cL:=sqrt(E/rho); omega:=kappa*cL;
> kappa:=(2*n-1)*Pi/2/l; 
# Mit den Werten des Beispiels
> rho:=3500.; E:=2.5*10^(9); nu:= 2.8*10^(5); l:=5.; A:=7.57*10.^(-4); alpha:= 0.5; u0:=0.005;
> gamma:=nu/E; delta:=0.5*gamma*omega^2;
> u:=0: nE:=50:
> for n to nE do
>   u:=u + un;
> end do:
# Normierung der Zustandsgrößen
> v:=diff(u,t):
> N:=E*A/l*(diff(u,xi) + gamma*diff(u,xi,t)):
# Animation der Zustandsgrößen
> tE:=0.05:
> 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 = ["¾ = x/", "u(¾,t)"]],t=0..tE):
> p2:=plots[animate]( plot,[v,para,title = "\nGeschwindigkeit\n",labels = ["¾ = x/", "v(¾,t)"]],t=0..tE):
> p3:=plots[animate]( plot,[N,para,title = "\nNormalkraft\n",    labels = ["¾ = x/", "N(¾,t)"]],t=0..tE):
> plots[display](Matrix(1,3,[p1,p2,p3]));
# Von Interesse sind noch die Verschiebung des rechten Endes
> plot(subs(xi=1,u),t=0..tE,thickness=2,gridlines=true,title = "\nVerschiebung des rechten Endes\n",labels = ["t", "u(l,t)"],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed);
# und der zeitliche Verlauf der Einspannkraft:
> plot(subs(xi=0,N),t=0..tE,thickness=2,gridlines=true,title = "\nEinspannkraft\n",labels = ["t", "N(x=0,t)"],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed);
> 
;
