Course:PHYS350/Rigid Bodies

Definition

• In a rigid body, no part of the body moves relative to another part of the body regardless of the force applied.
• A rigid body is equivlent to a system of particles restricted to have fixed distances apart.

There are two methods for solving for the motion of rigid bodies.

1. Write out the kinetic energy $T$ and use the Lagrangian to derive the equations of motion
2. Use Euler's equations to get the motion of the centre of mass and the rotation.

The Motion of a Rigid Body

The restrictions of rigid body motion require that distance between any two particles within the body remain constant. For example

$l_{ij}^{2}=\left({\vec {r}}_{i}-{\vec {r}}_{j}\right)^{2}=\left(x_{i}-x_{j}\right)^{2}+\left(y_{i}-y_{j}\right)^{2}+\left(z_{i}-z_{j}\right)^{2}$

will be constant during the motion. Let's take the time derivative of this quantity,

${\frac {dl_{ij}^{2}}{dt}}=2\left({\vec {r}}_{i}-{\vec {r}}_{j}\right)\cdot \left({\vec {v}}_{i}-{\vec {v}}_{j}\right)$

For this to vanish either ${\vec {v}}_{i}-{\vec {v}}_{j}$ must vanish or the dot product must vanish. The first corresponds to a translation of the object. For the dot product to vanish, the velocity difference must be perpendicular to the position difference. Specifically,

${\vec {v}}_{i}-{\vec {v}}_{j}={\vec {\omega }}_{ij}\times \left({\vec {r}}_{i}-{\vec {r}}_{j}\right).$

At first it might appear that you can have a different ${\vec {\omega }}_{ij}$ for each pair of particles. It turns out that this is not the case. Let's calculate

${\vec {v}}_{i}-{\vec {v}}_{k}=\left({\vec {v}}_{i}-{\vec {v}}_{j}\right)-\left({\vec {v}}_{j}-{\vec {v}}_{k}\right)={\vec {\omega }}_{ij}\times \left({\vec {r}}_{i}-{\vec {r}}_{j}\right)-{\vec {\omega }}_{jk}\times \left({\vec {r}}_{j}-{\vec {r}}_{k}\right).$

On the other hand we could have calculated it as follows,

${\vec {v}}_{i}-{\vec {v}}_{k}={\vec {\omega }}_{ik}\times \left({\vec {r}}_{i}-{\vec {r}}_{k}\right)={\vec {\omega }}_{ik}\times \left({\vec {r}}_{i}-{\vec {r}}_{j}\right)-{\vec {\omega }}_{ik}\times \left({\vec {r}}_{j}-{\vec {r}}_{k}\right).$

For these expressions to be equal for arbitrary values of ${\vec {r}}_{i}$, ${\vec {r}}_{j}$ and ${\vec {r}}_{k}$, the three values of ${\vec {\omega }}$ must be equal to each other so we will drop the subscripts and call this quantity the angular velocity of the rotation.

To summarize, a rigid body has two types of motion:

1. Translation of the entire body, and
2. Uniform rotation about a particular axis (the direction of ${\vec {\omega }}$ and an angular speed given by the magnitude of ${\vec {\omega }}$.

Angular Velocity

We can use some point within or outside the body as the origin and denote its position at any time by ${\vec {r}}_{0}$. The velocity of a particular part of the body is

${\vec {v}}={\vec {v}}_{0}+{\vec {\omega }}\times \left({\vec {r}}-{\vec {r}}_{0}\right)$

where ${\vec {r}}$ is the location of the part of the body and ${\vec {v}}_{0}$ is the velocity of the origin.

Let's move the origin to ${\vec {r}}_{1}$. The velocity of the particle does not change so we have

${\vec {v}}={\vec {v}}_{0}+{\vec {\omega }}\times \left({\vec {r}}-{\vec {r}}_{0}\right)={\vec {v}}_{1}+{\vec {\omega }}\times \left({\vec {r}}-{\vec {r}}_{1}\right)$

so

${\vec {v}}_{1}={\vec {v}}_{0}+{\vec {\omega }}\times \left({\vec {r}}_{1}-{\vec {r}}_{0}\right).$

The vector ${\vec {\omega }}$ is constant with respect to changes in the origin. A useful example of this is a wheel that rolls without slipping.

Moment of Inertia

Rigid bodies are such an important part of classical mechanics that we have developed special techniques to calculate their kinetic energies and angular momenta. Let's calculate the kinetic energy of a mass in rigid body motion. Its velocity is

${\vec {v}}_{i}={\vec {v}}_{0}+{\vec {\omega }}\times {\vec {r}}_{i}$

where we have written the coordinates relative to the origin (${\vec {r}}_{0}$).

Kinetic Energy

Its kinetic energy is

$T_{i}={\frac {1}{2}}m_{i}\left({\vec {v}}_{i}\cdot {\vec {v}}_{i}\right)={\frac {1}{2}}m_{i}\left({\vec {v}}_{0}+{\vec {\omega }}\times {\vec {r}}_{i}\right)\cdot \left({\vec {v}}_{0}+{\vec {\omega }}\times {\vec {r}}_{i}\right).$

and multiplying out the terms we get

$T_{i}={\frac {1}{2}}m_{i}\left[v_{0}^{2}+2{\vec {v}}_{0}\cdot \left({\vec {\omega }}\times {\vec {r}}_{i}\right)+\left({\vec {\omega }}\times {\vec {r}}_{i}\right)^{2}\right].$

We can simplify the result using the following fact

$\left({\vec {A}}\times {\vec {B}}\right)^{2}=A^{2}B^{2}\cos ^{2}\theta =A^{2}B^{2}\left(1-\sin ^{2}\theta \right)=A^{2}B^{2}\left[1-{\frac {\left({\vec {A}}\cdot {\vec {B}}\right)^{2}}{A^{2}B^{2}}}\right]=A^{2}B^{2}-\left({\vec {A}}\cdot {\vec {B}}\right)^{2}$

to get

$T_{i}={\frac {1}{2}}m_{i}\left[v_{0}^{2}+2{\vec {v}}_{0}\cdot \left({\vec {\omega }}\times {\vec {r}}_{i}\right)+{\vec {\omega }}^{2}{\vec {r}}_{i}^{2}-\left({\vec {\omega }}\cdot {\vec {r}}_{i}\right)^{2}\right].$

Angular Momentum

Let's look at the angular momentum of the particle, we have

${\vec {L}}_{i}=m_{i}{\vec {r}}_{i}\times {\vec {v}}_{i}=m_{i}{\vec {r}}_{i}\times {\vec {v}}_{0}+{\vec {r}}_{i}\times \left({\vec {\omega }}{\vec {r}}_{i}\right).$

Here we can use the vector triple product to simplify the result

${\vec {A}}\times \left({\vec {B}}\times {\vec {C}}\right)={\vec {B}}\left({\vec {A}}\cdot {\vec {C}}\right)-{\vec {C}}\left({\vec {A}}\cdot {\vec {B}}\right)$

which for us yields

${\vec {L}}_{i}=m_{i}{\vec {r}}_{i}\times {\vec {v}}_{i}=m_{i}{\vec {r}}_{i}\times {\vec {v}}_{0}+m_{i}\left[{\vec {\omega }}r^{2}-{\vec {r}}\left({\vec {\omega }}\cdot {\vec {r}}\right)\right].$

If ${\vec {v}}_{0}=0$ then

$T_{i}={\frac {1}{2}}{\vec {\omega }}\cdot {\vec {L}}_{i}.$

In this case, let's write out ${\vec {L}}_{i}$ component by component, we have

${\vec {L}}_{i}=m_{i}\left[{\begin{matrix}\omega _{x}r_{i}^{2}-x_{i}\left(\omega _{x}x_{i}+\omega _{y}y_{i}+\omega _{z}z_{i}\right)\\\omega _{y}r_{i}^{2}-y_{i}\left(\omega _{x}x_{i}+\omega _{y}y_{i}+\omega _{z}z_{i}\right)\\\omega _{z}r_{i}^{2}-z_{i}\left(\omega _{x}x_{i}+\omega _{y}y_{i}+\omega _{z}z_{i}\right)\end{matrix}}\right]=m_{i}\left[{\begin{matrix}y_{i}^{2}+z_{i}^{2}&x_{i}y_{i}&x_{i}z_{i}\\x_{i}y_{i}&x_{i}^{2}+z_{i}^{2}&y_{i}z_{i}\\x_{i}z_{i}&y_{i}z_{i}&x_{i}^{2}+y_{i}^{2}\end{matrix}}\right]\left[{\begin{matrix}\omega _{x}\\\omega _{y}\\\omega _{z}\end{matrix}}\right]\equiv {\vec {\vec {I}}}_{i}\cdot {\vec {\omega }}$

defining the moment of inertia matrix for that particular particle.

Putting It All Together

Because the quantities ${\vec {v}}_{0}$ and ${\vec {\omega }}$ are the same for every particle, we can sum up the kinetic energy and angular momentum for all of the particles to get

$T=\sum _{i}T_{i}={\frac {1}{2}}Mv_{0}^{2}+{\frac {1}{2}}{\vec {\omega }}\cdot {\vec {\vec {I}}}\cdot {\vec {\omega }}+M{\vec {v}}_{0}\cdot \left({\vec {\omega }}\times {\vec {r}}_{CM}\right)$

where $M$ is the total mass, ${\vec {r}}_{CM}$ is the location of the centre of mass and ${\vec {\vec {I}}}=\sum _{i}{\vec {\vec {I}}}_{i}$ is the moment of inertia of the body about the origin. The total angular momentum is

${\vec {L}}=M{\vec {r}}_{CM}\times {\vec {v}}_{0}+{\vec {\vec {I}}}\cdot {\vec {\omega }}$

Using the result for the scalar triple product we can show that

${\vec {L}}={\frac {\partial T}{\partial \omega }}$

so the angular momentum is the generalized momentum conjugate to the angular velocity; if the potential does not depend on angle, the angular momentum is a first integral of the motion.

Manupulating the Moment of Inertia

Principal Axes

The moment of inertia is a positive definite matrix. We know this because when we multiply it by any non-zero vector ${\vec {\omega }}$, we get the kinetic energy of the body that is necessarily positive. Such a matrix can always be diagonalized. From a physical point of view, that means that the eigenvector of the matrix are known as the principal axes of the body. If the angular velocity points along one of the principal axes, then the angular momentum is parallel to the angular velocity. Furthermore, if one chooses the principal axes of the body to be one's coordinate axes, the diagonal element of the moment of inertia matrix vanish, so for example,

${\vec {L}}=\left[{\begin{matrix}I_{xx}\omega _{x}\\I_{yy}\omega _{y}\\I_{zz}\omega _{z}\end{matrix}}\right]$

and the kinetic energy is

$T={\frac {1}{2}}Mv_{0}^{2}+M{\vec {v}}_{0}\cdot \left({\vec {\omega }}\times {\vec {r}}_{CM}\right)+{\frac {1}{2}}I_{xx}\omega _{x}^{2}+{\frac {1}{2}}I_{yy}\omega _{y}^{2}+{\frac {1}{2}}I_{zz}\omega _{z}^{2}$

Change of Origin

Sometimes it is easiest to calculate the moment of inertia about a particular point but you are interested in the moment of inertia about another point. There is a straightforward prescription to achieve this. First, let's use the centre of mass of the system as the origin, so we have

$T={\frac {1}{2}}Mv_{CM}^{2}+{\frac {1}{2}}{\vec {\omega }}\cdot {\vec {\vec {I}}}_{CM}\cdot {\vec {\omega }}$

Let's write this same equation with respect to a new origin ${\vec {r}}_{0}$ so we have

$T={\frac {1}{2}}M\left({\vec {v}}_{0}+{\vec {\omega }}\times {\vec {r}}_{CM}\right)^{2}+{\frac {1}{2}}{\vec {\omega }}\cdot {\vec {\vec {I}}}_{CM}\cdot {\vec {\omega }}$

      $={\frac {1}{2}}Mv_{0}^{2}+M{\vec {v}}_{0}\cdot \left({\vec {\omega }}\times {\vec {r}}_{CM}\right)+{\frac {1}{2}}M\left({\vec {\omega }}\times {\vec {r}}_{CM}\right)^{2}+{\frac {1}{2}}{\vec {\omega }}\cdot {\vec {\vec {I}}}_{CM}\cdot {\vec {\omega }}.$

If we compare this result to the earlier result we find

${\vec {\omega }}\cdot {\vec {\vec {I}}}{\vec {\omega }}={\vec {\omega }}\cdot {\vec {\vec {I}}}_{CM}\cdot {\vec {\omega }}+M\left({\vec {\omega }}\times {\vec {r}}_{CM}\right)^{2}.$

On a component by component basis we have

${\vec {\vec {I}}}={\vec {\vec {I}}}_{CM}+M\left[{\begin{matrix}y_{CM}^{2}+z_{CM}^{2}&x_{CM}y_{CM}&x_{CM}z_{CM}\\x_{CM}y_{CM}&x_{CM}^{2}+z_{CM}^{2}&My_{CM}z_{CM}\\x_{CM}z_{CM}&y_{CM}z_{CM}&x_{CM}^{2}+y_{CM}^{2}\end{matrix}}\right]$

To change from one origin to another you first use the equation above to get the moment of inertia relative to the centre of mass and then use it a second time to get the moment of inertia about the new origin.

Calculating a Moment of Inertia

Let's calculate the moment of inertia of an ellipsoid. Let's assume that the axes of the ellipsoid are line up along the coordinate axes. Along the $x-$axis the ellipsoid spans from $-a$ to $a$. Along the $y$ and $z-$axes, the bounds are $b$ and $c$ respectively.

Because the coordinate axes are the principal axes of the ellipsoid, the off-diagonal components of the matrix will vanish. We will calculate the sum of $z^{2}$ over the ellipsoid first. We have

$\sum _{i}m_{i}z_{i}^{2}=\int _{\rm {Ellipsoid}}\rho dVx^{2}=\int _{\rm {Ellipsoid}}\rho dxdydzz^{2}.$

The region of integration is rather complicate in the Cartesian coordinates but this gives us a place to start. We would like to simplify things a bit by defining

$x=au,y=bv,z=cw.$

In these new coordinates the ellipsoid becomes a sphere of unit radius. We also have to include the Jacobian of the coordinate transformation

${\frac {\partial \left(x,y,z\right)}{\partial \left(u,v,w\right)}}=abc.$

Now our integral looks like

$\sum _{i}m_{i}z_{i}^{2}=\int _{\rm {UnitSphere}}\rho abc(cw)^{2}dudvdw.$

Let's calculate the density of the ellipsoid in terms of its $M$ and volume. The volume of the ellipsoid is

$V=\int _{\rm {Ellipsoid}}dxdydz=\int _{\rm {UnitSphere}}abcdudvdw={\frac {4}{3}}\pi abc$

so we can rewrite the density in the moment of inertia integral

$\sum _{i}m_{i}z_{i}^{2}=\int _{\rm {UnitSphere}}{\frac {3M}{4\pi abc}}abc(au)^{2}dudvdw={\frac {3Mc^{2}}{4\pi }}\int _{\rm {UnitSphere}}w^{2}dudvdw.$

The easiest way to integrate over a unit sphere is to use spherical coordinates. Let's do it:

$\sum _{i}m_{i}z_{i}^{2}={\frac {3Mc^{2}}{4\pi }}\int _{0}^{1}dr\int _{0}^{\pi }d\theta \int _{0}^{2\pi }d\phi r^{2}\sin \theta r^{2}\cos ^{2}\theta ={\frac {3Mc^{2}}{4\pi }}{\frac {2}{5}}\pi \int _{-1}^{1}d\left(\cos \theta \right)\cos ^{2}\theta ={\frac {3Mc^{2}}{4\pi }}{\frac {4}{15}}\pi ={\frac {1}{5}}Mc^{2}.$

By symmetry we can calculate the results for the other two coordinates. We have

$\sum _{i}m_{i}x_{i}^{2}={\frac {1}{5}}Ma^{2},\sum _{i}m_{i}y_{i}^{2}={\frac {1}{5}}Mb^{2}$

and the moment of inertia matrix about the centre of mass have the following non-zero components:

$I_{xx}={\frac {M}{5}}\left(b^{2}+c^{2}\right),I_{yy}={\frac {M}{5}}\left(a^{2}+c^{2}\right),I_{zz}={\frac {M}{5}}\left(a^{2}+b^{2}\right).$

This is an example of reducing a general shape to a symmetric one through a transformation of variables. One important result of this technique is Routh's rule.

Ellipsoid of Inertia

We saw earlier that if the velocity of the origin vanished, then we could write $T={\frac {1}{2}}{\vec {L}}\cdot {\vec {\omega }}$. Now if there are no forces on the body, the kinetic energy of the body is conserved. If there are no forces on the body, there are no torques, so the angular momentum is conserved as well. Does this mean that the angular velocity is constant?

The answer is of course no because the equation only says that the component of the angular velocity along the direction of the angular momentum is constant. The other components can change. For example, if the angular velocity lies along one of the principal axes (eigenvectors of the moment of inertia matrix), then the angular momentum points in the same direction and the angular velocity must stay constant in direction and magnitude.

A usual picture to understand the force-free motion of a rigid body is the ellipsoid of inertia. Let's use the principal axes of the body as a coordinate frame to write the kinetic energy

$T={\frac {1}{2}}\left(I_{xx}\omega _{x}^{2}+I_{yy}\omega _{y}^{2}+I_{zz}\omega _{z}^{2}\right).$

The ellipsoid of inertia is the locus of values of ${\vec {\omega }}$ that give a kinetic energy of the body.

Ellipsoid of Inertia of an Ellipsoid

Let's calculate the shape of the ellipsoid of inertia of an ellipsoid. We have

$T={\frac {1}{2}}{\frac {M}{5}}\left[\left(b^{2}+c^{2}\right)\omega _{x}^{2}+\left(a^{2}+c^{2}\right)\omega _{y}^{2}+\left(a^{2}+b^{2}\right)\omega _{y}^{2}\right],$

so along the $x$, $y$ and $z$ directions, the axes of the ellipsoid of inertia are

${\sqrt {\frac {10T}{M(b^{2}+c^{2})}}},{\sqrt {\frac {10T}{M(a^{2}+c^{2})}}},{\sqrt {\frac {10T}{M(a^{2}+b^{2})}}}.$

The size of the ellipsoid of inertia increases as $T$ and decreases with the mass and size of the actual ellipsoidal body. The shape of the ellipsoid of inertia reflects the shape of the physical ellipsoid. The largest axis of the ellipsoid is also the largest axis of the corresponding ellipsoid of inertia. If two axes of the ellipsoid are the same size, the corresponding axes of the ellipsoid of inertia will be equal as well.

Motion of the Ellipsoid of Inertia

The Ellipsoid of Inertia Rolling on the Invariable Plane

To conserve energy the angular velocity must remain on the ellipsoid of inertia and the orientation of ellipsoid of inertia determines the orientation of the body. We can use the conservation of angular momentum to orient the ellipsoid. If we remember that ${\vec {L}}=\partial T/\partial {\vec {\omega }}$ we see that the normal to the ellipsoid of inertia is the angular momentum; therefore, as the ellipsoid of inertia rotates about the angular velocity, it must remain tangent to a plane perpendicular to the angular momentum, the invariable plane. Furthermore, the angular velocity is restricted to move along the ellipsoid of inertia along curves of constant $\left|{\vec {L}}\right|$. This curve is called the polhode, and it traces the path of the angular velocity through the body. The angular velocity also traces a path in the invariable plane called the herpolhode. The herpolhode is the path of the angular velocity through space. We can also write that $2T={\vec {\omega }}\cdot {\vec {L}}$ so the centre of the ellipsoid of inertia must remain a constant distance above the invariable plane. The height of the centre of the ellipsoid is $2T/\left|{\vec {L}}\right|$, so when the ellipsoid is down low the angular momentum is high relative to the kinetic energy.

Keeping the ellipsoid of inertia tangent to the invariable plane preserves the direction of the angular momentum. The magnitude of the angular momentum is also conserved. This defines a second ellipsoid through the equation

$\left|{\vec {L}}\right|^{2}=I_{xx}^{2}\omega _{x}^{2}+I_{yy}^{2}\omega _{y}^{2}+I_{zz}^{2}\omega _{z}^{2}$

so the polhodes are the intersections of the ellipsoid of inertia with the angular momentum ellipsoid. The angular momentum ellipsoid deviates more from a sphere than the ellipsoid of inertia; this gives an idea of where the polhodes lie for small and large angular momenta.

If the body has two principal axes with equal moments of inertia, the polhodes are circles centered on the axis with the unique moment of inertia, and the herpolhodes are circles in the invariable plane. If the unique moment of inertia is larger than the others, the body is called oblate. Otherwise it is prolate -- in analogy to the spheroids with the similar properties. In this case, the free rotation of the body consists of a constant precession of the angular velocity in a circle about the angular momentum in the space frame (the space cone) and about the unique principal axis in the body frame (the body cone). In the oblate case the body cone rolls within the space cone, and in the prolate case, the body cone rolls on the outside of the space cone.

The polhodes near the major axis

The polhodes near the intermediate axis

The polhodes near the minor axis

If no two moments of inertia are equal to each other, the polhodes are much more complicated. The polhodes near the axes with the largest and smallest moment of inertia are closed, but those near the intermediate axis are not; consequently, rotation about the intermediate axis is not stable as we shall see when we analyze the free rotation using Euler's equations.

Euler's Equations

So far we have tried to find a graphical description of the motion of an object without any torques. Although "the polhode rolls on the herpolhode without slipping" might paint a nice Victorian picture, it doesn't exploit the mathematical language of physics to which we are now accustomed. We know that the change in the angular momentum is equal to the torques on the body

${\frac {d{\vec {L}}}{dt}}={\vec {N}}.$

Let's write the angular momentum in terms of the principal axes of the body. Here we will be looking at how the principal axes move so we will use the unit vectors ${\hat {1}},{\hat {2}}$ and ${\hat {3}}$ to denote the current direction of the principal axes of the body. We can write

${\frac {d{\vec {L}}}{dt}}={\frac {\partial L_{1}}{\partial t}}{\hat {1}}+{\frac {\partial L_{2}}{\partial t}}{\hat {2}}+{\frac {\partial L_{3}}{\partial t}}{\hat {3}}+L_{1}{\frac {\partial {\hat {1}}}{\partial t}}+L_{2}{\frac {\partial {\hat {2}}}{\partial t}}+L_{3}{\frac {\partial {\hat {3}}}{\partial t}}$

        $={\frac {\partial L_{1}}{\partial t}}{\hat {1}}+{\frac {\partial L_{2}}{\partial t}}{\hat {2}}+{\frac {\partial L_{3}}{\partial t}}{\hat {3}}+L_{1}{\vec {\omega }}\times {\hat {1}}+L_{2}{\vec {\omega }}\times {\hat {1}}+L_{3}{\vec {\omega }}\times {\hat {3}}$

Let's define the angular momentum relative to the body axes as

${\vec {L}}_{b}=\left[{\begin{matrix}L_{1}\\L_{2}\\L_{3}\end{matrix}}\right]=\left[{\begin{matrix}\omega _{1}I_{1}\\\omega _{2}I_{2}\\\omega _{3}I_{3}\end{matrix}}\right]$

to get

${\frac {d{\vec {L}}}{dt}}={\frac {d{\vec {L}}_{b}}{dt}}+{\vec {\omega }}\times {\vec {L}}_{b}={\vec {N}}$

where the torque and angular velocity are written relative to the body axes. One can imagine using Euler's equations would be quite cumbersome if one had to include the torques because one would always have to transform the torques from the inertia frame to the body frame. Let's write out these equations component by component

$I_{1}{\dot {\omega }}_{1}-\omega _{2}\omega _{3}\left(I_{2}-I_{3}\right)=N_{1}$

$I_{2}{\dot {\omega }}_{2}-\omega _{3}\omega _{1}\left(I_{3}-I_{1}\right)=N_{2}$

$I_{3}{\dot {\omega }}_{3}-\omega _{1}\omega _{2}\left(I_{1}-I_{2}\right)=N_{3}$

These equations are most powerful when there are no torques, so the right-hand sides are zero (but check out Euler Equations with Torques). In this case we can see immediately that if any two moments of inertia differ than the angular velocity must lie along one of the principal axes to remain constant.

If two of the moments of inertia are equal to each other (let's take $I_{1}=I_{2}$), we can solve the torque-free Euler's equations exactly. We have

$I_{1}{\dot {\omega }}_{1}=\left(I_{1}-I_{3}\right)\omega _{3}\omega _{2}$

$I_{2}{\dot {\omega }}_{2}=-\left(I_{1}-I_{3}\right)\omega _{3}\omega _{1}$

$I_{3}{\dot {\omega }}_{3}=0$

We can solve this with

${\dot {\omega }}_{1}=-\Omega \omega _{2},{\dot {\omega }}_{2}=\Omega \omega _{1}$

with

$\Omega ={\frac {I_{3}-I_{1}}{I_{1}}}\omega _{3},$

yielding the solution

$\omega _{1}=A\cos \Omega t,\,\omega _{2}=A\sin \Omega t.$

This result coincides with the graphical picture that the polhodes of a symmetric body are circles centered on the axis with the unique moment of inertia.

Chandler Wobble

We can calculate the expected rate for the Earth to precess due to its oblatness. Let's estimate the ratio of the moments of inertia in terms of the radii of the Earth at the equator and at the poles, we have

${\frac {I_{3}-I_{1}}{I_{1}}}={\frac {2R_{\rm {eq}}^{2}-\left(R_{\rm {eq}}^{2}+R_{\rm {pole}}^{2}\right)}{R_{\rm {eq}}^{2}+R_{\rm {pole}}^{2}}}$

      $={\frac {\left(R_{\rm {eq}}-R_{\rm {pole}}\right)\left(R_{\rm {eq}}+R_{\rm {pole}}\right)}{R_{\rm {eq}}^{2}+R_{\rm {pole}}^{2}}}\approx {\frac {R_{\rm {eq}}-R_{\rm {pole}}}{R_{\rm {eq}}}}.$

where in the last step we have assume that the difference between the polar and equatorial radius is small. We have $R_{\rm {eq}}\approx 6378{\rm {km}},R_{\rm {pole}}\approx 6356{\rm {km}}$, so the ratio is about $1/290$. What is $\omega _{3}$ for the Earth? It is $2\pi /1{\rm {day}}$, so the period of the free precession of the Earth should be about 290 days. The period is somewhat longer about 433 days [1]. You can check out a plot of the wobble at the Paris Observatory. Why the Earth wobbles is a bit of a puzzle since because the Earth is not a solid body, without a driving force it would have damped long ago. Since the period of the wobble is similar to a year, seasonal changes are thought to be to blame.

I first learned about the wobble while spending a summer at the Pulkovo Astronomical Observatory near St. Petersburg. I saw a plot of the location of the Earth's pole on the surface of the Earth -- it moves in a irregular circular shape with a radius of 3 to 15 meters. The wobble can be measured by accurate measurements of the locations of stars relative to landmarks on the Earth. I asked my hosts why did they need to keep track the wobble so accurately. "To aim our missiles at your cities -- of course" came the reply.

Triaxial Body

Unfortunately, we can't find such a straightforward to the solution of the motion of a triaxial body. In the symmetric case, one of the components of the angular velocity was constant with time. This is not the case for a triaxial body. We have the following Euler's equations

$I_{1}{\dot {\omega }}_{1}=\left(I_{2}-I_{3}\right)\omega _{3}\omega _{2}$

$I_{2}{\dot {\omega }}_{2}=\left(I_{3}-I_{1}\right)\omega _{3}\omega _{1}$

$I_{3}{\dot {\omega }}_{3}=\left(I_{1}-I_{2}\right)\omega _{1}\omega _{2}.$

Let's take $I_{3}>I_{2}>I_{1}$ and define the following three positive ratios of moments of inertia

$A={\frac {I_{3}-I_{2}}{I_{1}}},B={\frac {I_{3}-I_{1}}{I_{2}}},C={\frac {I_{2}-I_{1}}{I_{3}}}$

giving the following equations

${\dot {\omega }}_{1}=-A\omega _{3}\omega _{2}$

${\dot {\omega }}_{2}=+B\omega _{3}\omega _{1}$

${\dot {\omega }}_{3}=-C\omega _{1}\omega _{2}$

Major Axis

Let's assume that the rotation is nearly about the major axis, so $\omega _{1}\gg \ \omega _{2},\omega _{3}$. We will only include the latter components of the angular velocity to first order, so we have

${\dot {\omega }}_{1}\approx 0,{\dot {\omega }}_{2}=B\omega _{3}\omega _{1},{\dot {\omega }}_{3}=-C\omega _{2}\omega _{1}.$

If we define $\Omega ={\sqrt {BC}}\omega _{1}$ we have the following solution for $\omega _{2}$ and $\omega _{3}$,

$\omega _{3}={\sqrt {\frac {B}{C}}}K\cos \Omega t,\omega _{2}=K\sin \Omega t,$

so the angular velocity travels in an ellipse centered on the major axis.

Minor Axis

Let's assume that the rotation is nearly about the major axis, so $\omega _{3}\gg \ \omega _{1},\omega _{2}$. We will only include the latter components of the angular velocity to first order, so we have

${\dot {\omega }}_{3}\approx 0,{\dot {\omega }}_{1}=-A\omega _{3}\omega _{2},{\dot {\omega }}_{2}=-B\omega _{3}\omega _{1}.$

If we define $\Omega ={\sqrt {AB}}\omega _{3}$ we have the following solution for $\omega _{1}$ and $\omega _{2}$,

$\omega _{1}={\sqrt {\frac {A}{B}}}K\cos \Omega t,\omega _{2}=K\sin \Omega t,$

so the angular velocity travels in an ellipse centered on the minor axis.

Intermediate Axis

Let's assume that the rotation is nearly about the intermediate axis, so $\omega _{2}\gg \ \omega _{1},\omega _{3}$. We will only include the latter components of the angular velocity to first order, so we have

${\dot {\omega }}_{2}\approx 0,{\dot {\omega }}_{1}=-A\omega _{3}\omega _{2},{\dot {\omega }}_{3}=-C\omega _{1}\omega _{2}.$

If we define $\alpha ={\sqrt {AC}}\omega _{2}$ we have the following solution for $\omega _{1}$ and $\omega _{2}$,

$\omega _{1}=K_{1}e^{\alpha t},\omega _{3}=K_{3}e^{\alpha t},\,$

so the angular velocity diverges exponetially away from the intermediate axis.

Lagrangian Treatment of Rigid-Body Motion

Although it is possible to include the effects of torques on the motion of a rigid body within the framework of Euler's equations, it is rather cumbersome. Here we will develop the Lagragnian treatment of a top -- a symmetric rigid body with a torque from the Earth's gravity.

Euler Angles

In a Lagrangian treatment it is necessary to find a set of independent coordinates to describe the position of the system at any time. For rigid bodies a convenient set of such coordinates are called the Euler angles. There are several different conventions for the Euler angles. The convention used here makes the analysis of the top simpler.

To get from one set of coordinates to the other, we have to do a series of three rotations. We imagine that the primed coordinates are the body coordinates and the unprimed coordinates are the inertia coordinates; therefore, it is straightforward to write the kinetic energy in the primed coordinates and the potential energy in the unprimed coordinates. Here are the three steps in words and in terms of a rotation matrix

1. Rotate by an angle $\phi$ about the initial $z-$axis. The final position of the $x-$axis is called the line of nodes.

$D=\left[{\begin{matrix}\cos \phi &\sin \phi &0\\-\sin \phi &\cos \phi &0\\0&0&1\end{matrix}}\right]$

1. Rotate about the line of nodes by an angle $\theta$. This changes the direction of the $z-$axis and pitches the $y-$axis out of the initial $x-y-$plane.

$C=\left[{\begin{matrix}1&0&0\\0&\cos \theta &\sin \theta \\0&-\sin \theta &\cos \theta \end{matrix}}\right]$

1. Finally rotate about the new $z-$axis by an angle $\psi$. This takes the $x-$axis out of the initial $x-y-$plane.

$B=\left[{\begin{matrix}\cos \psi &\sin \psi &0\\-\sin \psi &\cos \psi &0\\0&0&1\end{matrix}}\right]$

Euler proved that any rotation can be decomposed into these three rotations. One can see this mathematically by constructing the rotation matrix for these three successive transformations ($A=BCD$) and verifying that it is the most general rotation matrix.

The Kinetic Energy

To calculate the kinetic energy of the body, we have to find the angular velocity in terms of the coordinates, $\phi ,\theta$ and $\psi$ and their time derivatives. In any coordinate system, the angular velocity is directed along the axis of rotation. We can construct the total angular velocity vector using the following information:

1. Changes in $\phi$ produce rotations about the $z-$axis (inertial frame).
2. Changes in $\psi$ produce rotations about the $z'-$axis (body frame).
3. Changes in $\theta$ produce rotations about the line of nodes.

We can write the angular velocity relative to the inertial frame to get

${\vec {\omega }}=\left[{\begin{matrix}0\\0\\1\end{matrix}}\right]{\dot {\phi }}+\left[{\begin{matrix}\cos \phi \\\sin \phi \\0\end{matrix}}\right]{\dot {\theta }}+\left[{\begin{matrix}\sin \theta \sin \phi \\-\sin \theta \cos \phi \\\cos \theta \end{matrix}}\right]{\dot {\psi }}$

or relative to the body frame

${\vec {\omega }}=\left[{\begin{matrix}0\\0\\1\end{matrix}}\right]{\dot {\psi }}+\left[{\begin{matrix}\cos \psi \\-\sin \psi \\0\end{matrix}}\right]{\dot {\theta }}+\left[{\begin{matrix}\sin \theta \sin \psi \\\sin \theta \cos \psi \\\cos \theta \end{matrix}}\right]{\dot {\phi }}.$

The first form is useful to visualize what the body is doing, while the second form is good for calculating the kinetic energy because we have used the primed coordinate system to characterize the location of the principal axes of the body.

For simplicity let's specialize to symmetric bodies with $I_{x'}=I_{y'}$, so we have

$T={\frac {1}{2}}I_{x'}\left(\omega _{x'}^{2}+\omega _{y'}^{2}\right)+{\frac {1}{2}}I_{z'}\omega _{z'}^{2}.$

Because the primed frame is fixed to the body, the moments of inertia in the kinetic energy are constant with time. Let's first calculate

$\omega _{x'}^{2}+\omega _{y'}^{2}=\left(\cos \psi {\dot {\theta }}+\sin \psi \sin \theta {\dot {\phi }}\right)^{2}+\left(-\sin \psi {\dot {\theta }}+\cos \psi \sin \theta {\dot {\phi }}\right)^{2}$

    $=\cos ^{2}\psi {\dot {\theta }}^{2}+\sin ^{2}\psi \sin ^{2}\theta {\dot {\phi }}^{2}+2\sin \psi \cos \psi \sin \theta {\dot {\theta }}{\dot {\phi }}+\sin ^{2}\psi {\dot {\theta }}^{2}+\cos ^{2}\psi \sin ^{2}\theta {\dot {\phi }}^{2}-2\sin \psi \cos \psi \sin \theta {\dot {\theta }}{\dot {\phi }}$

    $={\dot {\theta }}^{2}+\sin ^{2}\theta {\dot {\phi }}^{2}.$

This simplication only obtains for a symmetric body. The $z'-$component is

$\omega _{z'}^{2}=\left({\dot {\psi }}+{\dot {\phi }}\cos \theta \right)^{2}.$

Combining these results yields the kinetic energy

$T={\frac {I_{x'}}{2}}\left({\dot {\theta }}^{2}+{\dot {\phi }}^{2}\sin ^{2}\theta \right)+{\frac {I_{z'}}{2}}\left({\dot {\psi }}+{\dot {\phi }}\cos \theta \right)^{2}$

The Potential Energy

The potential energy is pretty simple using the Euler angles. If $l$ is the distance between the pivot point and the center of mass of the top, we have

$V=Mgl\cos \theta$

The Lagrangian

The Lagrangian is the difference between the kinetic and the potential energy, so we have

$L=T-V={\frac {I_{x'}}{2}}\left({\dot {\theta }}^{2}+{\dot {\phi }}^{2}\sin ^{2}\theta \right)^{2}+{\frac {I_{z'}}{2}}\left({\dot {\psi }}+{\dot {\phi }}\cos \theta \right)^{2}-Mgl\cos \theta .$

We notice immediately that $\psi$ and $\phi$ don't appear in the Lagrangian. The system looks the same as in turns around the $z-$ and $z'-$axes. This means that we have two conserved quantities, reducing the three dimensional problem to a single dimension (we will choose this to be $\theta$). Because the Lagrangian does not depend on time, the Hamiltonian is conserved as well so we can use the techniques outlined in Linear Motion to understand the top.

Integrals of the Motion

To calculate the conserved momenta we have to calculate the partial derivatives of the Lagrangian with respect to the generalized velocities. We have

$p_{\psi }={\frac {\partial L}{\partial {\dot {\psi }}}}=I_{z'}\left({\dot {\psi }}+{\dot {\phi }}\cos \theta \right)=I_{z'}\omega _{z'}\equiv I_{x'}a$

where the last step defines the conserved quantity $a$ that has units of angular velocity. Similarly for the coordinate $\phi$, we have

$p_{\phi }={\frac {\partial L}{\partial {\dot {\phi }}}}=\left(I_{x'}\sin ^{2}\theta +I_{z'}\cos ^{2}\theta \right){\dot {\phi }}+I_{z'}\cos \theta {\dot {\psi }}\equiv I_{x'}b.$

Finally we have the Hamiltonian that is also conserved,

$H=T+V={\frac {I_{x'}}{2}}\left({\dot {\theta }}^{2}+{\dot {\psi }}^{2}\sin ^{2}\theta \right)+{\frac {I_{z'}}{2}}\left({\dot {\psi }}+{\dot {\phi }}\cos \theta \right)^{2}+Mgl\cos \theta .$

As we mentioned earlier, these three integrals are sufficient to solve for the motion. We can use the two conserved momenta to solve for the generalized velocities, ${\dot {\psi }}$ and ${\dot {\phi }}$, in terms of the conserved quantities $a$ and $b$ and the value of the third coordinate $\theta$. We have

${\dot {\phi }}={\frac {b-a\cos \theta }{\sin ^{2}\theta }},$

${\dot {\psi }}={\frac {I_{x'}a}{I_{z'}}}-\cos \theta {\frac {b-a\cos \theta }{\sin ^{2}\theta }}.$

Conservation of Energy

Since both the Hamiltonian and the angular velocity along the $z'-$axis are conserved, let's define the conserved quantity

$E'=E-{\frac {I_{z'}\omega _{z'}^{2}}{2}}={\frac {I_{x'}{\dot {\theta }}^{2}}{2}}+{\frac {I_{x'}}{2}}{\frac {\left(b-a\cos \theta \right)^{2}}{\sin ^{2}\theta }}+Mgl\cos \theta$

which we can rewrite as

$E'={\frac {I_{x'}{\dot {\theta }}^{2}}{2}}+V'\left(\theta \right)$

where

$V'(\theta )={\frac {I_{x'}}{2}}{\frac {\left(b-a\cos \theta \right)^{2}}{\sin ^{2}\theta }}+Mgl\cos \theta .$

We can make further progress by looking at $u=\cos \theta$ which will simplify the potential. Let's write

${\dot {u}}^{2}=\sin ^{2}\theta {\dot {\theta }}^{2}=\left(1-u^{2}\right){\dot {\theta }}^{2}$

so we can write $E'$ in terms of $u$ and its time derivative

$E'={\frac {I_{x'}}{2}}{\frac {{\dot {u}}^{2}}{1-u^{2}}}+{\frac {I_{x'}}{2}}{\frac {\left(b-au\right)^{2}}{1-u^{2}}}+Mglu.$

What remains is to solve for ${\dot {u}}^{2}$ to get the boundaries of the motion

${\dot {u}}^{2}=\left(E'-Mglu\right){\frac {2}{I_{x'}}}\left(1-u^{2}\right)-\left(b-au\right)^{2}$

The Motion of a Symmetric Top

The equation for ${\dot {u}}^{2}$ is cubic so it will generally have three roots (designated at $u_{1,2,3}$ -- it must have at least one root because for large positive values of $u$ ${\dot {u}}^{2}$ is positive and for large negative values of $u$ ${\dot {u}}^{2}$ is negative. As the top moves from straight up to straight down, $u$ ranges from 1 to -1, so only the roots within this range are important. Furthermore, the equation for ${\dot {\phi }}$ has a root at $u_{0}=b/a$.

The figure shows a typical curve for ${\dot {u}}^{2}$ as a function of $\theta$. The motion is bounded between $u=0.5$ and $u=0.8$. Depending on the value of $b/a$, the motion in the $\phi$-direction can be monotonic or not. If $b/a=0.8$, the motion will have a cusp at the upper end of the range.

Motion of the top with b/a=0.75

Motion of the top with b/a=0.8

Motion of the top with b/a=0.85

The figures to the right show the motion of the $z'-$axis as a function of time for several values of $b/a$ with the curve for ${\dot {u}}^{2}$ above. The bounds of the motion in all three cases is from $u=0.5$ to $u=0.8$. How can we set up a top to execute these three classes of motion? If we release the axis of the top from a stationary position, the axis will first drop (nutation) because of gravity. As it drops the angular momentum about the vertical axis must be conserved so the tip of the top starts moving sideways as well (precession) -- yielding the motion in the middle figure. On the other hand if we release the tip of the top with some motion in the direction of the precession, we will get the lower figure. Finally to get the rather freaky looking motion in the top figure, we release the top with some motion opposing the direction of the precession.

Specifically let's examine the case where we release the top from a stationary position, so $u_{0}=b/a=u_{2}$ where $u_{2}$ is the middle root of the ${\dot {u}}^{2}$ curve. Initially we have

$E'=Mgl\cos \theta =Mgl{\frac {b}{a}}$

so we can write

${\dot {u}}^{2}=\left({\frac {b}{a}}-u\right)Mgl{\frac {2}{I_{x'}}}\left(1-u^{2}\right)-\left(b-au\right)^{2}$

      $=\left(u_{0}-u\right)a^{2}\left[{\frac {2Mgl}{I_{x'}a^{2}}}\left(1-u^{2}\right)-\left(u_{0}-u\right)\right]$

The Fast Top

Unfortunately, the motion of the top in general cannot be solved in closed form, but we can make some progress in understanding the motion of the so-called fast top. If we look at the expression above, we have the ratio $2Mgl/(I_{x'}a^{2}).$ The numerator is the potential energy of the top while the denominator is the related to the kinetic energy of the top about its spin axis. The condition for a fast top is

$I_{x'}a^{2}=I_{z'}\omega _{z'}^{2}{\frac {I_{z'}}{I_{x'}}}\gg \ 2Mgl.$

In a fast top, the kinetic energy is much larger than the potential energy; consequently, the extent of the nutation ($u_{0}-u$) is small, so

$1-u^{2}\approx 1-u_{0}^{2}=\sin ^{2}\theta _{0}$

so we can write

${\dot {u}}^{2}\approx \left(u_{0}-u\right)a^{2}\left[{\frac {2Mgl}{I_{x'}a^{2}}}\sin ^{2}\theta _{0}-\left(u_{0}-u\right)\right].$

The extent of the nutation is

$u_{0}-u_{1}\approx {\frac {2Mgl}{I_{x'}a^{2}}}\sin ^{2}\theta _{0}={\frac {I_{x'}}{I_{z'}}}{\frac {2Mgl}{I_{z'}\omega _{z'}^{2}}}\sin ^{2}\theta _{0}$

and we can rewrite the equation for ${\dot {u}}^{2}$ as

${\dot {u}}^{2}=\left(u_{0}-u\right)\left(u-u_{1}\right)a^{2}.$

Let's take the time derivative of both sides to yield

$2{\dot {u}}{\ddot {u}}=a^{2}\left[\left(u_{0}-u\right){\dot {u}}-\left(u-u_{1}\right){\dot {u}}\right]$

${\ddot {u}}=a^{2}\left({\frac {u_{0}+u_{1}}{2}}-u\right).$

Let's define $y=u-(u_{0}+u_{1})/2$ so we have

${\ddot {y}}=-a^{2}y$

with the solution

$y={\frac {u_{1}-u_{0}}{2}}\cos at$ and $u={\frac {u_{1}+u_{0}}{2}}+{\frac {u_{1}-u_{0}}{2}}\cos at.$

The frequency of the nutation is

$a={\frac {I_{z'}}{I_{x'}}}\omega _{z'}$

and

${\dot {\phi }}={\frac {a\left(u_{0}-u\right)}{\sin ^{2}\theta }}={\frac {Mgl}{I_{z'}\omega _{z'}}}\left(1-\cos at\right)$

and

$\phi =\phi _{0}+{\frac {Mgl}{I_{z'}\omega {z'}}}\left(t-{\frac {1}{a}}\sin at\right).$

Sleeping Top

A sleeping top is a top that is set in motion with its spin axis vertical and the axis remains vertical. When the top is vertical it is impossible to distinguish changes in $\phi$ from changes in $\psi$ so the relationship between ${\dot {\phi }}$ and $a$ and $b$ isn't very useful. Let's look at the equation for ${\dot {u}}^{2}$. We have

${\dot {u}}^{2}=\left(E'-Mglu\right){\frac {2}{I_{x'}}}\left(1-u^{2}\right)-\left(b-au\right)^{2}.$

When the top points upwards, $u=1$. We let go of the top with ${\dot {u}}=0$ so we have

$0=\left(E'-Mgl\right){\frac {2}{I_{x'}}}\left(1-1^{2}\right)-\left(b-a\right)^{2}.$

so

$b=a$ and

${\dot {u}}^{2}=\left(1-u\right)a^{2}\left[\left(E'-Mglu\right){\frac {2}{I_{x'}a^{2}}}\left(1+u\right)-\left(1-u\right)\right].$

Using the definition of $E'$ we find

$E'=E-{\frac {I_{z'}\omega _{z'}^{2}}{2}}={\frac {I_{x'}{\dot {\theta }}^{2}}{2}}+{\frac {I_{x'}}{2}}{\frac {\left(b-a\cos \theta \right)^{2}}{\sin ^{2}\theta }}+Mgl\cos \theta =Mgl$

and

${\dot {u}}^{2}=\left(1-u\right)^{2}a^{2}\left[{\frac {2Mgl}{I_{x'}a^{2}}}\left(1+u\right)-1\right]=\left(1-u\right)^{2}a^{2}{\frac {2Mgl}{I_{x'}a^{2}}}\left[u-\left({\frac {I_{x'}a^{2}}{2Mgl}}-1\right)\right],$

so there is a double root at $u=1$ and another root at

$u={\frac {I_{x'}a^{2}}{2Mgl}}-1.$

If this root is at $u>1$, then the top will sleep, so

$I_{x'}a^{2}>4Mgl$

is the condition for a sleeping top.