# Kapitel  3_h: 
# Reines Rollen einer inhomogenen Scheibe, das Rollpendel
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): unprotect(gamma):
# Beispiel 3-17:
# 
# Eine inhomogene Scheibe der Masse m bewegt sich durch reines Rollen unter dem Einfluß der Gewichtskraft G längs einer schiefen Ebene mit dem Neigungswinkel ± abwärts. Der Schwerpunkt der Scheibe mit dem Radius a liegt um die Strecke cS außermittig. Zum Zeitpunkt t = 0 befindet sich die Scheibe mit Æ = 0 in der gestrichelt dargestellten Ruhelage. Gesucht wird das Bewegungsgesetz, wenn sich die Scheibe nach dem Loslassen unter dem Einfluß der äußeren Gewichtskraft weiterbewegt.
> PDEtools[declare](phi(t), prime = t);
# Wir formen die im Buch TM 3 angegebene Differenzialgleichung 3.47 für Õ(t) noch etwas um und erhalten:
> GL:=(Thetaq + 1 + gamma^2 - 2*gamma*sin(phi(t)))*diff(phi(t),t,t) - gamma*diff(phi(t),t)^2*cos(phi(t)) - gq*(gamma*cos(alpha+phi(t)) + sin(alpha));
# Für die grafische Darstellung der Winkelbeschleunigung lösen wir die Differenzialgleichung nach phi'' auf und erhalten:
> phipp:=solve(GL,diff(phi(t),t,t));
# Zur Abkürzung haben wir folgende Größen eingeführt:
> Thetaq:=Theta/(m*a^2);gamma:=cS/a; gq:=g/a;
# Die Systemwerte des Beispiels sind (d: Scheibendicke)
> alpha:=evalf(Pi/9); a:=0.5; cS:=a/3;  rho:=1.; d:=a/10 ; m:=evalf(rho*a^2*Pi*d); Theta:= evalf(m/2*a^2); g:=9.81;
# Mit diesen Werten lautet die Differenzialgleichung
> GL; phipp;
# Die Scheibe startet aus der Ruhelage mit den folgenden Anfangsbedingungen
> AB:=phi(0) = 0, D(phi)(0) = 0;
> N:=200; #Anzahl der Zeitschritte
;
> TE:=3;  #Simulationszeit
;
> 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);
# Für die Darstellung der Winkelbeschleunigung lösen wir die Differenzialgleichung nach phi'' auf 
> phipp:=(gamma*(diff(phi(t), t))^2*cos(phi(t))+gq*(gamma*cos(alpha+phi(t))+sin(alpha)))/(Thetaq+1+gamma^2-2*gamma*sin(phi(t)));
> 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"]):
> q2:=plots[odeplot](L,[t,diff(phi(t),t)],0..TE,para,title ="Winkelgeschwindigkeit\n",labels = ["Zeit t","phi'"]):
> q3:=plots[odeplot](L,[t,phipp],0..TE,para,title = "Winkelbeschleunigung\n",labels = ["Zeit t", "phi''"]):
> A := Array(1 .. 3): A[1]:= q1: A[2]:= q2: A[3]:= q3:
> plots[display](A);
# Animation des Bewegungsvorganges für die inhomogene Scheibe auf schiefer Ebene
> Scheibe:=plottools[disk]([0,a],a,color="Nautical GrayViolet"):
> NT:=LinearAlgebra[RowDimension](L[2,1]): R:= NULL:
> for j to NT do    #Schleife über jeden Zeitschritt
>   P:=plottools[translate](Scheibe,a*L[2,1][j,2],0);
>   Q:=plots[display](P): R:= R,Q:
> end do:
> AN:=plots[display](R,insequence = true,scaling = constrained):
> AN;
# Das Rollpendel
> restart: with(LinearAlgebra): unprotect(gamma): with(plottools):
# Die im Buch TM 3 angegebene Differenzialgleichung wird noch etwas umgeformt.
> PDEtools[declare](psi(t), prime = t);
> GL:=(Thetaq + 1 + gamma^2 -2*gamma*cos(psi(t)))*diff(psi(t),t,t) + gamma*(diff(psi(t),t)^2+gq)*sin(psi(t));
# Zur Abkürzung haben wir folgende Größen eingeführt:
> Thetaq:=Theta/(m*a^2); gamma:=cS/a; gq:=g/a;
# Für die grafische Darstellung der Winkelbeschleunigung lösen wir die Differenzialgleichung nach psi'' auf und erhalten:
> psipp:=solve(GL,diff(psi(t),t,t));
# Systemwerte:
> a:=0.5; cS:=a/3;  rho:=1.; d:=a/10 ; m:=evalf(rho*a^2*Pi*d); Theta:= evalf(m/2*a^2); g:=9.81;
> GL;
# Das Rollpendel startet laut Aufgabenstellung aus der Ruhelage mit der Anfangsauslenkung È(t = 0) = À/9 = 20°.
> AB:=psi(0)=Pi/9, D(psi)(0)=0;
> N:=100;   #Anzahl der Zeitschritte
;
> TE:=4.8; #Simulationszeit
;
> 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, psi(t)],0..TE, para,title = "Drehwinkel\n",labels=["Zeit t", "psi"]):
> q2:=plots[odeplot](L,[t,diff(psi(t),t)],0..TE,para,title ="Winkelgeschwindigkeit\n",labels = ["Zeit t","psi'"]):
> q3:=plots[odeplot](L,[t,psipp],0..TE,para,title = "Winkelbeschleunigung\n",labels = ["Zeit t", "psi''"]):
> A := Array(1 .. 3): A[1]:= q1: A[2]:= q2: A[3]:= q3:
> plots[display](A);
# Animation des Bewegungsvorganges für das Rollpendel
> Scheibe:= plottools[disk]([0,a],a,color=red,transparency=0.5):
> NT:=LinearAlgebra[RowDimension](L[2,1]): R:=NULL:
> for j to NT do    #Schleife über jeden Zeitschritt
>   P := plottools[translate](Scheibe,a*L[2,1][j,2],0):
>   M := plottools[point]([a*L[2,1][j,2],a],symbol=circle,symbolsize =15,color=black):
>   S := plottools[point]([cS*sin(L[2,1][j,2]),a-cS*cos(L[2,1][j,2])],symbol=circle,symbolsize =15,color=white):
>   Q := plots[display]([P,M,S]): 
>   R := R,Q:
> end do:
> AN:=plots[display](R,insequence=true,scaling=constrained,gridlines=true):
> AN;
# Wir geben noch die Lage des Scheibenschwerpunktes aus:
> q4:=plots[odeplot](L,[t,  cS*sin(psi(t))],0..TE,color=blue): #x-Komponente
;
> q5:=plots[odeplot](L,[t,a-cS*cos(psi(t))],0..TE,color= red): #y-Komponente
;
> plots[display]([q4,q5],axes=boxed,gridlines=true,title = "\nKomponenten des Schwerpunktes\nx: blau, y: rot\n",titlefont = ["ARIAL", 15], labels = ["Zeit t", "x,y-Komponenten"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8]);
# Hinweis: Für kleine Drehwinkel È(t)  und kleine erste Ableitungen dÈ(t)/dt  geht die Gleichung des Rollpendels über in eine linearisierte Schwingungsdifferenzialgleichung mit der Eigenkreisfrequenz É0 und der Schwingungsdauer T
> omega0:=sqrt(m*g*cS/(Theta+m*(a^2+cS^2-2*a*cS)));
> T:=evalf(2*Pi/omega0);2*T;
> 
;
