Gelöste Aufgaben/GYRO

Aufgabenstellung

Kreisel und Ihre Bewegungsgleichungen sind immer wieder eine Herausforderung für uns Ingenieure. Hier nähern wir uns hier dem Thema mit dem Aufstellen und Lösen der Bewegungsgleichungen für große Winkel - und damit nichtlinearen Differentialgleichungen. Die sonst üblichen Linearisierungen führen wir nicht durch - und schauen, ob wir die Gleichungen noch gelöst bekommen.

Kreisel mit gelenkiger Fußpunktlagerung.

Gesucht ist die Präzessions-Bahn eines Kreisels mit Euler-Drehwinklen. Unser Kreisel ist ein Kegel der Höhe H und Radius R. Die Bewegungsgleichungen sollen mit dem Prinzip der virtuellen Verrückungen aufgestellt werden. Wir suchen nach der Trajektorie des Mittelpunktes des Kreisels für große Kippwinkel.

Lösung mit Maxima und Matlab^®

Wir arbeiten mit Maxima und Matlab.

Maxima brauchen wir zum Aufstellen der Bewegungsgleichungen, die zunächst mal recht komplex aussehen. Dabei nehmen wir die Bewegungsgleichungen wie sie aus den Gleichgewichtsbedingungen kommen und lösen sie als Anfangswertproblem.

Geomertrie des Kreisels

Header

Wir lösen hier das Anfangswertproblem mit den nichtlinearen Bewegungsgleichungen zu den Koordinaten

$\underset{\_}{\phi }=\left(\begin{array}{c}{\phi }_1(t)\\ {\phi }_2(t)\\ {\phi }_3(t)\end{array}\right)$,

wobei φ₁, φ₂ Kippwinkel bezüglich der x₁- und x₂-Achsen sind, φ₃ ist der Rotationswinkel um die Kreisel-Symmetrieachse.

Euler-Kippwinkel φ₁, φ₂

Die Kippwinkel φ₁, φ₂ sind im nebenstehenden Bild eingetragen, die aus einer Transformation hervorgehenden Koordinatensysteme erhalten jeweils einen >'< zur Kennzeichnung. Dabei steht am Anfang das blue Koordinatensystem, das durch Kippung um φ₂ in das grüne und dann durch Kippung um φ₁ in das rote überführt wird. Nicht dargestellt ist die folgende Drehung um den Rotationwinkel φ₃.

Die Gleichgewichtsbedingungen kommen aus dem Prinzip der virtuellen Verrückungen. Für diese Gleichungen brauchen wir die virtuelle Arbeit der D'Alembert'schen Trägheitskräfte und die virtuelle Arbeit der Gewichtskraft des Kreisels.

📑 toggle code

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 21.05.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2022-03-27                            */
/* ref: gyroscopic precession                          */
/* description: solve using principle of virt. work    */
/*******************************************************/

Declarations

Wir brauchen zunächst die Transformationsmatrizen der Euler-Rotation, also die Drehmatrizen, die unser inertiales Koordinatensystem in das Kreiselfeste Koordinatensystem überführt. Wir kippen zunächst um die x₂-Achse, dann um die resultierende x'₁-Achse und schließlich um die Rotationsachse x"₃. Die Koordinatentransformation ist dann

${\underset{\_}{\underset{\_}{D}}}_{0\to K}={\underset{\_}{\underset{\_}{D}}}_3({\phi }_3(t))\cdot {\underset{\_}{\underset{\_}{D}}}_1({\phi }_1(t))\cdot {\underset{\_}{\underset{\_}{D}}}_2({\phi }_2(t))$.

Und wir führen noch Abkürzungen für die Systemparameter

$I_y=\frac{3(R^2+4H^2)m}{20},I_z=\frac{3R^{2}m}{10},m=\frac{\pi H{R}^2\rho }{3}$

und die Zustandsgrößen

$\frac{d}{dt}{\phi }_i(t)={\omega }_i(t)$

ein.

📑 toggle code

/* Euler-Drehmatrizen */
D[1](φ) := matrix([1,0,0],[0,cos(φ),-sin(φ)],[0,sin(φ),cos(φ)]);
D[2](φ) := matrix([cos(φ),0,sin(φ)],[0,1,0],[-sin(φ),0,cos(φ)]);
D[3](φ) := matrix([cos(φ),-sin(φ),0],[sin(φ),cos(φ),0],[0,0,1]);

coord: [[ φ[1](t), φ[2](t), φ[3](t)],
        [δφ[1]   ,δφ[2]   ,δφ[3]   ]];
q: funmake('matrix,makelist([coord[1][i]],i,1,3));

abbr: makelist(diff(φ[i](t),t)=ω[i](t),i,1,3);

params: [
  I[y] =      m       * 3/20*( R^2+4*H^2),
  I[z] =      m       * 3/10*  R^2,
   m   = ρ*%pi*R^2*H/3           ];
params: append(subst(params[3], [params[1], params[2]]),[params[3]]);
params: [params, makelist(solve(params[i], [H^3,R^4,R^2][i])[1],i,1,3)];

Equations of Motion

Die Gleichgewichtsbedingungen für den starren Kreisel (die virtuelle Formänderungenergie $\delta \mathrm{\Pi }$ ist Null) kommen aus

$\begin{array}{lll}\delta W& =& \sum _i({\overrightarrow{F}}_i\cdot \delta r_i)-\rho \underset{V}{\int }{\ddot{\overrightarrow{r}}}_P\cdot \delta {\overrightarrow{r}}_{P}dV\\ & \stackrel{!}{=}& 0\end{array}$,

wobei wir nur eine Kraft F - die Gewichtskraft $mg(-{\overrightarrow{e}}_3)$ - haben und V das Kreisel-Volumen ist.

Mit dem Ortsvektor zu einem Punkt des Kreisels

${\overrightarrow{r}}_P=\left(r0z\right)\cdot {\underset{\_}{\underset{\_}{D}}}_3(\alpha )\cdot {\underset{\_}{\underset{\_}{D}}}_{0\to K}({\phi }_1,{\phi }_2,{\phi }_3)\cdot \left(\begin{array}{c}{\overrightarrow{e}}_1\\ {\overrightarrow{e}}_2\\ {\overrightarrow{e}}_3\end{array}\right)$

können wir nun die virtuelle Gesamtarbeit des Systems ausrechnen. Die Parameter mit R und H ersetzen wir durch die oben definierten Abkürzungen und erhalten die gekoppelten Bewegungsgleichungen

$\left(\begin{array}{ccc}{I}_y& 0& 0\\ 0& (1-{\cos \left({\phi }_1\right)}^2)I_z+{\cos \left({\phi }_1\right)}^2{I}_y& -\sin \left({\phi }_1\right)I_z\\ 0& -\sin \left({\phi }_1\right)I_z& I_z\end{array}\right)\cdot \ddot{\underset{\_}{\phi }}+\left(\begin{array}{c}({{\omega }_2}^2\cos \left({\phi }_1\right)\sin \left({\phi }_1\right)-{\omega }_2{\omega }_3\cos \left({\phi }_1\right))I_z-{{\omega }_2}^2\cos \left({\phi }_1\right)\sin \left({\phi }_1\right)I_y\\ ({\omega }_1{\omega }_3\cos \left({\phi }_1\right)-2{\omega }_1{\omega }_2\cos \left({\phi }_1\right)\sin \left({\phi }_1\right))I_z+2{\omega }_1{\omega }_2\cos \left({\phi }_1\right)\sin \left({\phi }_1\right)I_y\\ {\omega }_1{\omega }_2\cos \left({\phi }_1\right)I_z\end{array}\right)=\frac{3}{4}Hmg\left(\begin{array}{c}\sin \left({\phi }_1\right)\cos \left({\phi }_2\right)\\ \cos \left({\phi }_1\right)\sin \left({\phi }_2\right)\\ 0\end{array}\right)$

Dieses System von gewöhnlichen Differentialgleichung müssen wir nach den Winkelbeschleunigungen lösen - und das geht nur, wenn die Massenmatrix regulär ist. Wir prüfen das anhand der Systemdeterminante

$\det \underset{\_}{\underset{\_}{M}}=\cos ({\varphi }_1{)}^2{I}_y^2{I}_z$

Für den Kippwinkel muss also cos(φ₁)≠0 oder 0 < φ₁ < 90° gelten. Das können wir für unsere Problemstellung gut einhalten.

📑 toggle code

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* local vector to some material point with coordinates φ[i](t), i=1,2,3 */
/* trafo, order of rotations: y, x', z'' */
r[p] : matrix([r,0,z]).D[3](α).D[3](φ[3](t)).D[1](φ[1](t)).D[2](φ[2](t));

/* variation of local vector */
δr[p] : sum(subst([ε=0],diff(subst([coord[1][i]=coord[1][i]+ε*coord[2][i]],r[p]),ε)),i,1,3);

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* virtual work of d'Alembert forces */ 
δW[a] :-ρ*diff(r[p],t,2).δr[p]$
/* integrate over body               */
δW[a] : integrate(integrate(integrate(r*δW[a], α,0,2*%pi), r,0,R*(z/H)), z,0,H)$
δW[a] : expand(expand(expand(expand(expand(expand(expand(δW[a])))))))$
δW[a] : expand(subst(params[2],δW[a]));

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* virtual work of external (gravitational) force */ 
r[m] : subst([r=0, z=3*H/4], r[p]);
δr[m] : sum(subst([ε=0],diff(subst([coord[1][i]=coord[1][i]+ε*coord[2][i]],r[m]),ε)),i,1,3);
δW[g] : m*matrix([0,0,-g]).δr[m];
δW[g] : expand(δW[g]);

Solving

Mit der Spaltenmatrix der Zustandsgrößen

$\underset{\_}{q}(t)=\left(\begin{array}{c}{\phi }_1(t)\\ {\phi }_2(t)\\ {\phi }_3(t)\\ {\omega }_1(t)\\ {\omega }_2(t)\\ {\omega }_3(t)\\ \end{array}\right)$

lösen wir die gewöhnlichen Differentialgleichungen

$\underset{\_}{\dot{q}}(t)=\underset{\_}{f}(\underset{\_}{q}(t))$

Dazu formulieren wir die Bewegungsgleichungen zunächst in Maxima und übertragen sie dann in ein Matlab^©-Script. Die Bewegungsgleichungen machen wir mit der Zeit

$T=\sqrt{\frac{3Hmg}{4I_y}}$

dimensionslos (die Winkel selbst sind ja bereits dimensionslos). Die dimensionslose Zeit $\tau =t/T$.

📑 toggle code

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* equations of motion */
eom: makelist(coeff(-δW[a],coord[2][row]),row,1,3);

/* mass matrix */
M: makelist(makelist(coeff(eom[row],diff(coord[1][col],t,2)),col,1,3),row,1,3);
M: trigsimp(funmake('matrix,M));
/* M^(-1) */
Mi : trigsimp(invert(M));

/* the rest ... */
gyr: trigsimp(expand(eom - M.diff(q,t,2)));
P:  makelist(coeff(-δW[g],coord[2][row]),row,1,3);

Matlab^©-files

Folder Structure

Die Datei-Struktur zeit das Skript GYRO.m. Classes und Functions sind in den jeweiligen Ordnern. Die Excel-Datei hält alle System-Parameter.

Den komplette Quellcode zu diesem Programm können Sie über dieses ZIP-File rechts herunterladen.

archive

download compressed archive →

Postprocessing

Die Ergebnisse der numerischen Integration tragen wir für die Winkel φ₁(t), φ₂(t) sowie für die Winkelgeschwindigkeiten ω₁(t), ω₂(t) auf. Es folgen noch die Trajektorie des Mittellinien-Endpunkts des Kreisels und die Phasendiagramme.

Winkel	Winkelgeschwindigkeit
Trajektorie	Phasendiagramme

📑 toggle code

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* if you wanted to linearize w.r.t. φ[1], φ[2] and ω[1], ω[2]            */
gyrs: trigsimp(subst(abbr,gyr));
smalls: [ω[1](t)=ε*ω[1](t),ω[2](t)=ε*ω[2](t)];
gyrs: subst([ε^2=0,ε=1],expand(subst(smalls,gyrs)));

G: makelist(makelist(coeff(expand(gyrs[row][1]/ω[3](t)),ω[col](t)),col,1,3),row,1,3);
G: trigsimp(funmake('matrix,G));

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* prepare output for Matlab                                              */
toNumerics: append(makelist(sin(φ[i](t))=s[i],i,1,3),makelist(cos(φ[i](t))=c[i],i,1,3), [I[z]=i*I[y]], solve(3/4*H*m*g=j*I[y],H));
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

MiG: ratsimp(subst(abbr,subst(toNumerics,trigsimp(Mi.gyr))));
MiP: ratsimp(subst(abbr,subst(toNumerics,trigsimp(Mi.P))));

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
rhs1: -MiG;
/*
[-c(1)*i*omega(2)*omega(3)-(c(1)*s(1)-c(1)*s(1)*i)*omega(2)^2;
 (i*omega(1)*omega(3)+(2*s(1)-s(1)*i)*omega(1)*omega(2))/c(1);
 (s(1)*i*omega(1)*omega(3)+(-s(1)^2*i+s(1)^2+1)*omega(1)*omega(2))/c(1)];
*/
rhs2: -MiP;
/*
s(1)*c(2)*obj.j
(s(2)*obj.j)/c(1)
(s(1)*s(2)*obj.j)/c(1))
*/

/* Center-Line*/

r[c]: subst(toNumerics,subst([r=0, z=1],r[p]));
/*
[c[1]*s[2];
-s[1];
c[1]*c[2]];

*/

Nachtrag

Meist linearisiert man die Bewegungsgleichungen, indem man von kleinen Kippwinkeln ausgeht (φ₁<<1, φ₂<<1). Für den Kreisel gilt ω₁<<ω₃, ω₂<<ω₃, so dass man außerdem ω₃ = constant setzten darf. Damit erhalten wir die "übliche" vereinfachte Bewegungsgleichung

$\left(\begin{array}{ccc}{I}_y& 0& 0\\ 0& I_y& 0\\ 0& 0& I_z\end{array}\right)\cdot \ddot{\underset{\_}{\phi }}+{\omega }_3\cdot \left(\begin{array}{ccc}0& -I_z& 0\\ I_z& 0& 0\\ 0& 0& 0\end{array}\right)\cdot \dot{\underset{\_}{\phi }}=\frac{3}{4}Hmg\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right)\cdot \underset{\_}{\phi }$

📑 toggle code

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* simplify */
MS: subst([cos(φ[1](t))=1, sin(φ[1](t))=0], M);
gyrS: subst([ε^2=0, ε=1],expand(subst([cos(φ[1](t))=1, sin(φ[1](t))=φ[1](t), ω[1](t)=ε*ω[1](t), ω[2](t)=ε*ω[2](t)],subst(abbr,gyr))));
PS: subst([cos(φ[1](t))=1, cos(φ[2](t))=1, sin(φ[1](t))=φ[1](t), sin(φ[2](t))=φ[2](t)], P);

Dieses lineare System von gewöhnlichen Bewegungsgleichungen kann man dann hervorragend als Eigenwertproblem analysieren.

---

Links

• ...

Literature

• ...