# Kapitel  3_i:  
# Das mathematische und das physische Pendel
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(plottools): with(plots):
# Beispiel 3-18:
# 
# Das mathematische Pendel (s.h. obige Skizze) ist ein idealisiertes Pendel, bei dem in der Modellvorstellung eine konzentrierte Masse m (idealerweise eine Punktmasse) an einem masselosen starren Stab befestigt ist. Auf die freigeschnittene Masse, die in der (x1,x2)-Ebene eine Kreisbewegung mit dem Radius  durchführt, wirken die Gewichtskraft G und die Stabkraft S. Die nichtlineare Differenzialgleichung 2. Ordnung lautet
> GL:=diff(phi(t),t,t)+omega^2*sin(phi(t));
# Zur Abkürzung wurde 
> omega:=sqrt(g/l);
# gesetzt. Der Differenzialgleichung entnehmen wir die Winkelbeschleunigung
> phippN:=-omega^2*sin(phi(t));
# Es sind folgende Anfangsbedingungen vorgegeben:
> AB:=phi(0)=phi0, D(phi)(0)=phip0;
# Das mathematische Pendel startet aus der vertikalen Ruhelage (Æ = 0, s.h. obige Skizze) mit der Anfangswinkelgeschwindigkeit phip0.
> phi0:=0.; phip0:=3.;
> g:=9.81; l:=3.; m:=1.;
# Wir lösen die Aufgabe nummerisch:
> N:=200; #Anzahl der Zeitschritte
;
> TE:=5.; #Schwingungszeit
;
> val:=Array(1..N+1): for i to N+1 do val[i]:=(i-1)*TE/N end do:
> LN:=dsolve({GL,AB},numeric,output=val); #Lösung der nichtlinearen Differenzialgleichung
;
# Die Stabkraft S (als Zugkraft positiv) errechnet sich zu (s.h. Buch TM3)
> SN:=m*l*(diff(phi(t),t)^2+omega^2*cos(phi(t)));
> para:= titlefont = ["ARIAL", 15],labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],gridlines = true,axes=boxed:
> q1:=plots[odeplot](LN,[t,phi(t)],0..TE, para,title = "Drehwinkel\n",labels=["Zeit t", "phi"],color=blue):
> q2:=plots[odeplot](LN,[t,diff(phi(t),t)],0..TE,para,title ="Winkelgeschwindigkeit\n",labels = ["Zeit t","phi'"],color=blue):
> q3:=plots[odeplot](LN,[t,phippN],0..TE,para,title = "Winkelbeschleunigung\n",labels = ["Zeit t", "phi''"],color=blue):
> q4:=plots[odeplot](LN,[t,SN],0..TE,para,title = "Stabkraft S\n",labels = ["Zeit t", "S"],color=blue):
# Im Fall kleiner Auslenkungen darf linearisiert werden. Dann gilt wegen sin Æ H Æ  die linearisierte Gleichung
> GLL:=diff(phi(t),t,t)+omega^2*phi(t);
# Der obigen Differenzialgleichungen entnehmen wir die Winkelbeschleunigung
> phippL:=-omega^2*phi(t);
# Wir lösen auch diese Aufgabe nummerisch
> LL :=dsolve({GLL,AB},numeric,output=val); #Lösung der linearen Differenzialgleichung
;
# Die Stabkraft S (als Zugkraft positiv) errechnet sich zu (s.h. Buch TM3)
> SL:=m*l*(diff(phi(t),t)^2+omega^2*cos(phi(t)));
> q5:=plots[odeplot](LL,[t,phi(t)],0..TE, para,labels=["Zeit t", "phi"],color=red):
> q6:=plots[odeplot](LL,[t,diff(phi(t),t)],0..TE,para,labels = ["Zeit t","phi'"],color=red):
> q7:=plots[odeplot](LL,[t,phippN],0..TE,para,labels = ["Zeit t", "phi''"],color=red):
> q8:=plots[odeplot](LL,[t,SL],0..TE,para,labels = ["Zeit t", "S"],color=red):
# Wir geben die Lösungen grafisch aus: 
> #farbe:=nichtlinear (blau), linear (rot):
> p1:=plots[display]([q1,q5],title = "Drehwinkel\nnichtlinear (blau), linear (rot)\n"):
> p2:=plots[display]([q2,q6],title = "Winkelgeschwindigkeit\nnichtlinear (blau), linear (rot)\n"): 
> p3:=plots[display]([q3,q7],title = "Winkelbeschleunigung\nnichtlinear (blau), linear (rot)\n"): 
> p4:=plots[display]([q4,q8],title = "Stabkraft\nfarbe\n"):
> A:= Array(1 .. 4): A[1]:= p1: A[2]:= p2: A[3]:= p3: A[4]:= p4:
> plots[display](A);
# Die obigen Grafiken zeigen, dass bis t = 0.5 die Kurven recht gut übereinstimmen. Bei großen Drehwinkeln weichen die Lösungen jedoch erheblich voneinander ab. Bei großen Drehwinkeln ist deshalb immer die Lösung der nichtlinearen Differenzialgleichung zu verwenden. 
# Wir animieren den Bewegungsvorgang des mathematischen Pendels:
> Scheibe:=plottools[disk]([0, -l],0.05*l,color = blue):
> Linie  :=plottools[line]([0,0],[0,-l],color = blue, thickness = 3):
> T      :=plots[display]([Scheibe,Linie]):
> NT:=LinearAlgebra[RowDimension](LN[2,1]): R:=NULL:
> for j to NT do    #Schleife über jeden Zeitschritt
>   S := plottools[rotate](T,LN[2,1][j,2],[0,0]):
>   P := plots[display](S): 
>   R := R,P:
> end do:
> AN:=plots[display](R,insequence=true,scaling=constrained,gridlines=true):
> AN;
# In den Umkehrpunkten kommt das Pendel mit phip = 0 zur Ruhe, und der größte Ausschlagswinkel im nichtlinearen Fall wird erreicht mit
> phimax:=2*arcsin(phip0/2/omega); evalf(convert(phimax,degrees));
# Beispiel 3-19:
# 
> restart: with(plottools): with(plots):
# Ein homogener dünner Stab ist im Punkt A reibungsfrei drehbar gelagert. Ein solches Pendel wird physisches Pendel genannt. Gesucht werden die Schnittlasten N, Q und M sowie die Lagerreaktionslasten bei A. Die Bestimmung der Schnittlasten in einem schwingenden Stab gehört in die Problemklasse der Kinetostatik. 
> GL:=diff(phi(t),t,t)+omegar^2*sin(phi(t));
# Der Differenzialgleichung entnehmen wir die Winkelbeschleunigung
> phipp:=-omegar^2*sin(phi(t));
> l:=3;      #Stablänge
;
> lr:=2/3*l; #Reduzierte Stablänge 
;
> omegar:=evalf(sqrt(9.81/lr));
# Das Stabpendel startet aus der vertikalen Ruhelage mit der Anfangswinkelgeschwindigkeit phip0.
> phi0 :=0.;  #Anfangsauslenkung
;
> phip0:=4.;  #Anfangswinkelgeschwindigkeit
;
> kappa:=2*omegar/phip0;
# Wegen º2 > 1 kommt das Pendel bei einem maximalen Ausschlagwinkel ± zum Stehen und kehrt dann seine Bewegungsrichtung um. Es liegt ein hin- und herschwingendes Pendel vor.
> AB:=phi(0)=phi0,D(phi)(0)=phip0;
> N:=200; #Anzahl der Zeitschritte
;
> TE:=5.; #Schwingungszeit
;
> val:=Array(1..N+1): for i to N+1 do val[i]:=(i-1)*TE/N end do:
> L:=dsolve({GL,AB},numeric,output=val);
> para:= titlefont = ["ARIAL", 15],labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],gridlines = true,axes=boxed:
> q1:=plots[odeplot](L,[t,phi(t)],0..TE, para,title = "Drehwinkel\n",labels=["Zeit t", "phi"],color=blue):
> q2:=plots[odeplot](L,[t,diff(phi(t),t)],0..TE,para,title ="Winkelgeschwindigkeit\n",labels = ["Zeit t","phi'"],color=blue):
> q3:=plots[odeplot](L,[t,phipp],0..TE,para,title = "Winkelbeschleunigung\n",labels = ["Zeit t", "phi''"],color=blue):
> A:= Array(1 .. 3): A[1]:= q1: A[2]:= q2: A[3]:= q3:
> plots[display](A);
# Wir bestimmen noch den maximalen Ausschlagwinkel, dort muss phip verschwinden.
> LM:=dsolve({GL,AB},numeric,output=Array([1.02,1.03,1.04,1.05]));
# Der maximale Ausschlagwinkel ± nach t =1.03 s hat nummerisch näherungsweise den Wert 
> alpha:=2.2535;
> alphagrad:=evalf(convert(alpha,degrees));
# Animation des Bewegungsvorganges 
> Stab:=rectangle([-0.1, 0.],[0.1, -l], color = "Nautical GrayViolet"):
> NT:=LinearAlgebra[RowDimension](L[2,1]): R:=NULL:
> for j to NT do    #Schleife über jeden Zeitschritt
>   P :=plottools[rotate](Stab,L[2,1][j,2],[0,0]):
>   Q :=plots[display](P): R:= R,Q:
> end do:
> AN:=plots[display](R,insequence=true,scaling=constrained,axes=boxed,gridlines=true):
> AN;
# Berechnung der Schnittlasten
# Systemwerte
> l:=3; lr:=2/3*l; omegar:=evalf(sqrt(9.81/lr));
# Wir wählen folgende Anfangsbedingungen (AB)
> phi0:=Pi/2;   #Anfangsauslenkung, der Stab wird aus der Horizontalen ohne Anfangsgeschwindigkeit losgelassen
;
> phip0:=0.;    #Anfangswinkelgeschwindigkeit
;
# Für diese Anfangsbedingungen erfolgt die Schnittlastenberechnung (s.h. Buch TM 3): Normalkraft und Querkraft werden auf das Stabgewicht G und das Biegemoment auf G bezogen
# 1. Die Normalkraft
> NS:=unapply(1/4*(1-xi)*(3*(1+xi)*(phip0^2/omegar^2 +2*(cos(phi)-cos(phi0)))+4*cos(phi)),xi,phi);
# 2. Die Querkraft
> QS:=unapply(1/4*(1-xi)*(1-3*xi)*sin(phi),xi,phi);
# 3. Das Biegemoment
> MS:=unapply(1/4*xi*(1-xi)^2*sin(phi),xi,phi);
# Das Biegemoment hat seinen Extremwert, wo die Querkraft verschwindet. Das ist an der Stelle xi = 1/3
> max_MS:=MS(1/3,phi);
# Wir geben im Folgenden die normierten Schnittlasten grafisch als 3D-Plots aus. 
> p1:=plot3d(NS(xi,phi), xi = 0 .. 1, phi = -phi0 .. phi0,titlefont = ["ROMAN", 15], title="\nNormierte Normalkraft\n",orientation=[145,-75,175]):
> p2:=plot3d(QS(xi,phi), xi = 0 .. 1, phi = -phi0 .. phi0, color = "SkyBlue",titlefont = ["ROMAN", 15], title="\nNormierte Querkraft\n",orientation=[145,-75,175]):
> p3:=plot3d(MS(xi,phi), xi = 0 .. 1, phi = -phi0 .. phi0, colorscheme = ["Blue", "LimeGreen"],titlefont = ["ROMAN", 15], title="\nNormiertes Biegemoment\n",orientation=[125,-75,180]):
> plots[display](Matrix(1,3,[p1,p2,p3]));
# Die Größtwerte der Schnittlasten können den obigen Grafiken entnommen werden.
> NS(0,0);
> QS(0,Pi/2); QS(0,-Pi/2);
# Das betragsmäßig größte Biegemoment tritt im Feld dort auf, wo die Querkraft verschwindet. Das ist bei ¾ = 1/3
> MS(1/3,Pi/2); MS(1/3,-Pi/2);
# Wir geben die Schnittlasten für einen festen Winkel Æ in Abhängigkeit von der Stabkoordinate ¾ aus:
> para:= xi=0..1,axes=boxed,gridlines=true,titlefont = ["ROMAN", 15]:
> q1:=plot(NS(xi,   0),para, title="\nNormierte Normalkraft N/G für phi = 0\n"):
> q2:=plot(QS(xi,Pi/2),para, title="\nNormierte Querkraft Q/G für phi = Pi/2\n"):
> q3:=plot(MS(xi,Pi/2),para, title="\nNormiertes Biegemoment M/(Gl) für phi = Pi/2\n"):
> plots[display](Array([q1, q2, q3]));
> 
;
#  
