# 
# Kapitel 2: 
# Der Stab mit Endmasse
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots): unprotect(gamma):
> 
;
# 
# Der in der obigen Abbildung skizzierte Dehnstab trägt am rechten Ende eine konzentrierte Masse M. Zum Zeitpunkt t = 0 wird diese um das Maß u0 ausgelenkt und dann ohne Anfangsgeschwindigkeit sich selbst überlassen. 
# Es sind die Zustandsgrößen für das Massenverhältnis µ = M/mst zu berechnen. Stellen Sie dazu eine Prozedur zur Verfügung.
> Proc_Calc_03:=proc(rho,E,A,l,u0,eps,N,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]
> #     u0: Anfangsauslenkung des rechten Stabendes
> #    eps: Massenverhältnis
> #      N: Anzahl der zu berechnenden Eigenwerte
> #     tE: Dauer der Animationszeit
> #Ausgabe:
> #     EW: Eigenwerte » 
> #     OM: Eigenkreisfrequenzen É
> #     An: Konstanten
> #      u: Auslenkung 
> #      v: Geschwindigkeit
> #      N: Normalkraft
> #--------------------------------------------------------------------
> local lambda,g1,g2,ewg,n,i,j,us,para,p1,p2,p3;
> global EW,OM,FR,C,An,b,u,v,NK;
> description "Berechnung der Zustandsgrößen eines Stabes mit Endmasse";
> EW:=Vector(N);OM:=Vector(N); FR:=Vector(N); C:=Matrix(N,N); b:=Vector(N);
> g1:=cot(lambda);
> g2:=eps*lambda;
> ewg:=g1-g2;
> for n to N do
>   EW[n]:=fsolve(ewg,lambda=(n-1)*Pi...n*Pi);
>   OM[n]:=sqrt(E/rho)/l*EW[n];
>   FR[n]:=OM[n]/2/Pi;
> end do:
> for i to N do
>   b[i]:=i/N;
>   for j to N do
>     C[i,j]:=sin(EW[j]*i/N);
>   end do:
> end do:
> An:=LinearSolve(C,b);
> us:=0; 
> for n to N do
>   us:=us + An[n]*sin(EW[n]*xi)*cos(OM[n]*t):
> end do:
> u :=u0*us;
> v :=diff(u,t);
> NK:=E*A/l*diff(u,xi):
> #Ausgabe der Ergebnisse
> plot([g1,g2,g1-g2],lambda=0..4*Pi,y=-2..4,title = "\nAuswertung der Eigenfunktion\n",titlefont = ["ARIAL", 15], labels = ["»", ""], labeldirections = ["horizontal", "vertical"],linestyle = [dash,dash,solid], axesfont = ["HELVETICA", "ROMAN", 8],labelfont = ["HELVETICA", 12],legend = ["cot »", "µ»","cot »-µ»"],legendstyle = [font = ["HELVETICA", 9], location = bottom],color=[blue,red,black],axes=boxed,gridlines=true);
> print(`Eigenwerte » = `,Transpose(EW));
> print(`Eigenkreisfrequenzen É = `,Transpose(OM));
> print(`Eigenfrequenzen f = `,Transpose(FR));
> print(`Konstanten An = `,Transpose(An));
> #Animation der Bewegungsgrößen
> 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 in [m]\n",       labels = ["¾ = x/", "u(¾,t)"]],t=0..tE):
> p2:=plots[animate]( plot,[v ,para,title = "\nGeschwindigkeit in [m/s]\n",labels = ["¾ = x/", "v(¾,t)"]],t=0..tE):
> p3:=plots[animate]( plot,[NK,para,title = "\nNormalkraft in [N]\n",      labels = ["¾ = x/", "N(¾,t)"]],t=0..tE):
> plots[display](Matrix(1,3,[p1,p2,p3]));
> end proc:
# Beispiel 3-6:
>  rho:=3500.; A:=7.57*10^(-4); E:=2.5*10^(8); l:=5.; u0:=0.005; epsilon:=0.2; N:=10; tE:=0.2;
                             3500.
                        0.0007570000000
                                      8

                        2.500000000 10 
                               5.
                             0.005
                              0.2
                               10
                              0.2
;
> Proc_Calc_03(rho,E,A,l,u0,epsilon,N,tE);
 Eigenwerte &lambda; = ,Vector[row](%id = 18446745942416053718)
Eigenkreisfrequenzen &omega; = ,Vector[row](%id = 18446745942416\



  714262)
  Eigenfrequenzen f = ,Vector[row](%id = 18446745942416714534)
    Konstanten An = ,Vector[row](%id = 18446745942416714398)

;
> 
;
# Von Interesse sind noch die Normalkräfte an den beiden Stabenden:
> plot([subs(xi=0,NK),subs(xi=1,NK)],t=0..tE,thickness=2,gridlines=true,title = "\nNormalkräfte an den Stabenden in [N]\n",labels = ["t", "N(t)"],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],legend = ["N(0,t)", "N(l,t)"],legendstyle = [font = ["HELVETICA", 9], location = bottom],color=[black,blue],axes=boxed,gridlines=true);

> 
;
