3.029 Spring 2022
Additional Reading Week 01: 01/31 - 02/04

Euler’s Equations of Rigid Bodies

Remember this quote from the (1^st) lecture?
“Undergraduate [physics] students come in expecting that the hardest thing they’ll have to learn will be relativity or quantum mechanics. 
Actually those are the most novel topics. The hardest thing that an undergraduate physics student must learn is the classical dynamics of spinning tops”

This (optional) reading will introduce the relevant concepts and equations in modeling rigid bodies such as spinning tops

The exposition here follows that in the excellent book “Structure and Interpretation of Classical Mechanics” by Prof. Sussman and Wisdom

Prof. Sussman and Wisdom teach an advanced class on these topics, 12.650J. If you find this reading to be of interest, you might enjoy taking that class in the future.

Rotation Kinematics

We start by examining the kinetic energy of a rigid body

Note: by rigid-body we mean an entity made up of large number of particles with mass , position , and velocities  with fixed positional constraints between them.

The kinetic energy is the simply given by

To proceed, it’s convenient to express the position of each particle as the relative position away from a fixed reference - here chosen to be the center of mass :

and express the kinetic energy as:

We can simplify this by noting that the average relative velocities of the rigid body vanish

and thus drop the second term to give

This can be identified as the sum of:

kinetic energy of motion of the total center of mass

kinetic energy of rotation about the center of mass

Finally, as we saw with the double pendulum, it’s convenient to express the velocities in terms of generalized coordinates.
To do this, we introduce the following formula:

where   is the angular velocity vector, and is related to the rotation matrix acting on the rigid body
(Note: a detailed derivation is outside the scope of this reading - but can be found here)

Moments of Inertia

We can now express the rotational kinetic energy as

Decomposing the angular velocity vector to its (Cartesian) coordinate system components

where   is the inertia tensor defined as

The moment of inertia tensor is built-in the Wolfram Language

So let’s explore some well-known objects

Solid Sphere (Ball)

MatrixForm[
MomentOfInertia[Ball[{0,0,0},r]]
]

MatrixForm[
MomentOfInertia[Ball[{0,0,0},r]]/Volume[Ball[{0,0,0},r]]
]

Hollow Sphere (Ball)

MatrixForm[
MomentOfInertia[Sphere[{0,0,0},r]]
]

MatrixForm[
MomentOfInertia[Sphere[{0,0,0},r]]/Area[Sphere[{0,0,0},r]]
]

Cylinder

MatrixForm[
Simplify[MomentOfInertia[Cylinder[{{0,0,0},{0,0,h}},r]],Assumptions->h>0]
]

MatrixForm[
Simplify[MomentOfInertia[Cylinder[{{0,0,0},{0,0,h}},r]]/Volume[Cylinder[{{0,0,0},{0,0,h}},r]],Assumptions->h>0]
]

Euler’s Equations

We now wish to explore the dynamics of rigid bodies using our knowledge of the rotational kinetic energy, and its decomposition into the moment of inertia tensor

To do this, we start by introducing Euler angles

Euler angles are an intuitive way to represent arbitrary rotations matrices M, as the composition of three simple rotations:

A counter-clockwise rotation of  around

rotationPsi=RotationTransform[ψ,{0,0,1}]

A counter-clockwise rotation of  around

rotationTheta=RotationTransform[θ,{0,1,0}]

A counter-clockwise rotation of  around

rotationPhi=RotationTransform[ϕ,{0,0,1}]

eulerRotation[ϕ_,θ_,ψ_]=Composition[rotationPhi,rotationTheta,rotationPsi]

Note: the Wolfram Language has built-in `EulerMatrix` functionality

AffineTransform[EulerMatrix[{ϕ,θ,ψ}]]

Euler’s equations of motion describe the rotation of rigid bodies using the principal axes of inertia

I.e. the eigenvalues of the moment of inertia,

Again, the derivation of Euler’s equations is outside the scope of this reading, but can be found here

eulersEqs[{I1_,I2_,I3_}]={
I1 ω1'[t]+(I3-I2)ω2[t] ω3[t]==0,
I2 ω2'[t]+(I1-I3)ω3[t] ω1[t]==0,
I3 ω3'[t]+(I2-I1)ω1[t] ω2[t]==0
};

where the angular velocity vector is specified in terms of Euler angles

eulerConstraints={
ω1[t]==ϕ'[t]Sin[θ[t]]Sin[ψ[t]]+θ'[t]Cos[ψ[t]],
ω2[t]==ϕ'[t]Sin[θ[t]]Cos[ψ[t]]-θ'[t]Sin[ψ[t]],
ω3[t]==ϕ'[t]Cos[θ[t]]+ψ'[t]
};

Intermediate Axis Theorem

Let’s use Euler’s equations to simulate an interesting phenomenon:

Try taking a (table-) tennis racket or book will also do and throw it in the air along one of its three main axes trying to make a full rotation

You should find that you can do this easily along the shortest and longest axes

But not along the intermediate axis, the racket keeps flipping!

First, let’s import a table-tennis paddle as a MeshRegion

and normalize such that its centroid is at {0,0,0} and its longest axis has length 1

paddle=RegionResize[Quiet[Import["https://cdn.thingiverse.com/assets/09/15/b4/5a/f1/PaddleFlat.stl"]],{1}];
Rasterize[
paddle=With[{centroid=RegionCentroid[paddle]},
TransformedRegion[paddle,RotationTransform[π/2,{0,1,0}]@*TranslationTransform[-centroid]]
],ImageSize->200]

We can find its Moment of Inertia tensor

and diagonalize it to obtain the principal moments and axes

centroid=RegionCentroid[paddle];
MatrixForm[ℐ=MomentOfInertia[paddle,centroid]]
{principalMoments,principalAxes}=Eigensystem[ℐ];
principalMoments=Diagonal[ℐ]

This allows us to visualize its inertia ellipsoid and 3 principal axes

inertiaellipsoid=Ellipsoid[centroid,Inverse[ℐ]/5000];
axes=Graphics3D[{
{Black,inertiaellipsoid},
MapThread[{#2,Arrowheads[{-0.05,0.05}],Arrow[Tube[{centroid-#/5,centroid+#/5}]]}&,{principalAxes,{Red,RGBColor[0,1/3,0],Blue}}]}];
Rasterize[Show[paddle,axes,BaseStyle->Opacity[0.25],ImageSize->200]]

Let’s try the experiment numerically

First, we’ll flip along the stable red axis

initialConditions["stable"]={ϕ[0]==0.001,θ[0]==0.001,ψ[0]==0.001,ω3[0]==0.,ω1[0]==10,ω2[0]==0.025};
{ϕSolution["stable"],θSolution["stable"],ψSolution["stable"]}=NDSolveValue[{eulerConstraints,eulersEqs[principalMoments],initialConditions["stable"]},{ϕ,θ,ψ},{t,0,1}];

Rasterize[
Multicolumn[
Table[
Show[
TransformedRegion[paddle,eulerRotation[
ϕSolution["stable"][t],
θSolution["stable"][t],
ψSolution["stable"][t]]
],Boxed->True,ImageSize->{Automatic,150},PlotRange->{{-3/4,3/4},{-3/4,3/4},{-3/4,3/4}}],{t,Subdivide[0,1,17]}],6,Appearance->"Horizontal"]
]

As expected, the paddle undergoes a clean rotation

Let’s try the same along the unstable blue axis

initialConditions["unstable"]={ϕ[0]==0.001,θ[0]==0.001,ψ[0]==0.001,ω3[0]==10,ω1[0]==0.,ω2[0]==0.025};
{ϕSolution["unstable"],θSolution["unstable"],ψSolution["unstable"]}=NDSolveValue[{eulerConstraints,eulersEqs[principalMoments],initialConditions["unstable"]},{ϕ,θ,ψ},{t,0,2π}];

Rasterize[
Multicolumn[
Table[
Show[
TransformedRegion[paddle,eulerRotation[
ϕSolution["unstable"][t],
θSolution["unstable"][t],
ψSolution["unstable"][t]]
],Boxed->True,ImageSize->{Automatic,150},PlotRange->{{-3/4,3/4},{-3/4,3/4},{-3/4,3/4}}],{t,Subdivide[0,2π,17]}],6,Appearance->"Horizontal"]
]

The paddle does in-fact flip before it makes a full rotation!

Separatrix

Let’s try and get a more intuitive picture of what’s happening!

First, let’s consider our degrees of freedom:

We have the three Euler angles

and their three derivatives

Now, let’s count our constants of motion (i.e. constraints)

First, we know the total energy must be conserved

For rigid bodies the Cartesian components of angular momentum,  are also conserved

We therefore have 6 degrees of freedom and 4 constraints

As such, we expect our trajectories to trace 2D-curves in phase-space!

Let’s move to the body’s principal frame of reference, and express the total angular momentum as

Similarly, we can express the energy in terms of the principal angular momentum components:

If we plot our trajectories in the principal angular momentum phase-space, our trajectories must respect both constraints

The first constraint, is that of a sphere

The second constraint, is that of an ellipsoid

As such, the only allowed trajectories are the curves of their intersection

sphere[{la_,lb_,lc_}]=(la^2)+(lb^2)+(lc^2);
ellipsoid[{a_,b_,c_}][{la_,lb_,lc_}]:=((((la/a)))^2)+((((lb/b)))^2)+((((lc/c)))^2)-2;

Show[
Table[
ContourPlot3D[{sphere[{la,lb,lc}]==a,ellipsoid[{1,0.5,0.75}][{la,lb,lc}]==0},{la,-2,2},{lb,-2,2},{lc,-2,2},
MeshFunctions->{Function[{la,lb,lc,f},(sphere[{la,lb,lc}]-a)-ellipsoid[{1,0.5,0.75}][{la,lb,lc}]]},MeshStyle->{{Opacity[.5],Dashed,Black}},PlotPoints->20,Mesh->{{0}},ContourStyle->{None,None},Boxed->False,Axes->False],{a,.001,4,.1}],
Graphics3D[{Opacity[1],Arrowheads[{-0.04,0.04}],{Specularity[Gray,10]},MapThread[{RGBColor[#2],Arrow[Tube[{centroid-1.5 #,centroid+1.5#},.02]]}&,{{{1,0,0},{0,1,0},{0,0,1}},IdentityMatrix[3]}]}],PlotRange->All
]

The above figure is called a separatrix and very insightful

It says that for small variations about two of the principal axes it's stable (i.e. closed trajectories exist in its vicinity)

About the blue axis however, variations are unstable and cause the body to flip!