Rigid bodies: inertia tensor, Euler's equations and gyroscopes

PHYS-130 treated rotation about a fixed axis, where angular momentum and angular velocity point the same way and the moment of inertia is a single number. Real bodies rotating freely in three dimensions are more interesting. A thrown phone or book tumbles when spun about one of its axes but not the others. A spinning top does not fall over but circles slowly. The Earth's axis traces a cone in the sky over 26,000 years. To understand these motions we need the inertia tensor, which relates angular momentum to angular velocity when they do not point the same way, and Euler's equations, which describe how a body's spin evolves in its own frame.

Angular momentum and the inertia tensor

A rigid body rotating about a fixed point (or about its centre of mass) with angular velocity \(\boldsymbol\omega\) has each particle moving at \(\mathbf v_i = \boldsymbol\omega\times\mathbf r_i\). Its angular momentum is

using the vector identity \(\mathbf a\times(\mathbf b\times\mathbf c) = \mathbf b(\mathbf a\cdot\mathbf c) - \mathbf c(\mathbf a\cdot\mathbf b)\). Each component of \(\mathbf L\) is a linear combination of the components of \(\boldsymbol\omega\), so we can write \(\mathbf L = \mathsf I\boldsymbol\omega\) with a symmetric 3×3 matrix, the inertia tensor:

The diagonal elements are the familiar moments of inertia about the three coordinate axes. The off-diagonal elements, the products of inertia, measure how lopsided the mass distribution is. When they are not zero, \(\mathbf L\) is not parallel to \(\boldsymbol\omega\). A badly balanced car wheel is an example: spinning about its axle, its angular momentum points slightly off the axle and rotates with the wheel, so the bearings must supply a rotating torque, felt as vibration. Balancing the wheel with small weights removes the products of inertia. The rotational kinetic energy is \(T = \tfrac12\boldsymbol\omega\cdot\mathbf L = \tfrac12\boldsymbol\omega^{\mathsf T}\mathsf I\boldsymbol\omega\).

Principal axes

Because \(\mathsf I\) is real and symmetric, linear algebra guarantees three mutually perpendicular directions, its eigenvectors, along which \(\mathsf I\boldsymbol\omega\) is parallel to \(\boldsymbol\omega\). These are the principal axes, and the eigenvalues \(I_1 \le I_2 \le I_3\) are the principal moments. In principal-axis coordinates the tensor is diagonal: \(\mathbf L = (I_1\omega_1, I_2\omega_2, I_3\omega_3)\) and \(T = \tfrac12(I_1\omega_1^2 + I_2\omega_2^2 + I_3\omega_3^2)\). Any symmetry of the body helps find them: an axis of rotational symmetry is always principal, and so is any axis perpendicular to a plane of mirror symmetry.

The parallel-axis theorem generalizes to tensors. If \(\mathsf I^{\text{cm}}\) is the inertia tensor about the centre of mass, the tensor about a point displaced by \(\mathbf a\) is \(I_{jk} = I^{\text{cm}}_{jk} + M(a^2\delta_{jk} - a_ja_k)\). The program uses the corner-based tensor of a uniform box, finds its principal moments with Jacobi's method (repeated rotations that zero the off-diagonal elements one at a time), and compares them with the principal moments about the centre.

import math

# 1. Inertia tensor of a uniform box (a phone: 0.15 x 0.07 x 0.008 m, 0.20 kg) about a corner,
#    then its principal moments by Jacobi rotations
def box_inertia_corner(m, a, b, c):
    Ixx, Iyy, Izz = m * (b * b + c * c) / 3, m * (a * a + c * c) / 3, m * (a * a + b * b) / 3
    return [[Ixx, -m * a * b / 4, -m * a * c / 4], [-m * a * b / 4, Iyy, -m * b * c / 4], [-m * a * c / 4, -m * b * c / 4, Izz]]

def jacobi_eigs(A, sweeps=50):
    A = [row[:] for row in A]
    n = len(A)
    for _ in range(sweeps):
        off = sum(A[i][j] ** 2 for i in range(n) for j in range(n) if i != j)
        if off < 1e-30: break
        for p in range(n):
            for q in range(p + 1, n):
                if abs(A[p][q]) < 1e-30: continue
                th = 0.5 * math.atan2(2 * A[p][q], A[q][q] - A[p][p])
                c, s = math.cos(th), math.sin(th)
                for k in range(n):                      # rotate columns p, q
                    akp, akq = A[k][p], A[k][q]
                    A[k][p], A[k][q] = c * akp - s * akq, s * akp + c * akq
                for k in range(n):                      # rotate rows p, q
                    apk, aqk = A[p][k], A[q][k]
                    A[p][k], A[q][k] = c * apk - s * aqk, s * apk + c * aqk
    return sorted(A[i][i] for i in range(n))

m, a, b, c = 0.20, 0.15, 0.07, 0.008
I_corner = box_inertia_corner(m, a, b, c)
print("principal moments about the corner (g m^2):", ", ".join(f"{x * 1000:.4f}" for x in jacobi_eigs(I_corner)))
I_cm = [m * (b * b + c * c) / 12, m * (a * a + c * c) / 12, m * (a * a + b * b) / 12]
print("principal moments about the CM    (g m^2):", ", ".join(f"{x * 1000:.4f}" for x in sorted(I_cm)))

# 2. Torque-free Euler equations: spin about the intermediate axis is unstable (the tennis racket flip)
I1, I2, I3 = sorted(I_cm)
def deriv(w):
    w1, w2, w3 = w
    return [(I2 - I3) / I1 * w2 * w3, (I3 - I1) / I2 * w3 * w1, (I1 - I2) / I3 * w1 * w2]
def rk4(w, dt):
    k1 = deriv(w); k2 = deriv([x + dt / 2 * y for x, y in zip(w, k1)])
    k3 = deriv([x + dt / 2 * y for x, y in zip(w, k2)]); k4 = deriv([x + dt * y for x, y in zip(w, k3)])
    return [x + dt / 6 * (p + 2 * q + 2 * r + s) for x, p, q, r, s in zip(w, k1, k2, k3, k4)]
for axis in (0, 1, 2):
    w = [0.01, 0.01, 0.01]; w[axis] = 20.0           # spin at 20 rad/s plus a tiny wobble
    w0 = w[:]; L0 = math.sqrt(sum((I * x) ** 2 for I, x in zip((I1, I2, I3), w)))
    t, dt, flips, sign = 0.0, 1e-4, [], 1
    while t < 3.0:
        w = rk4(w, dt); t += dt
        if w[axis] * sign < 0: flips.append(round(t, 3)); sign = -sign
    L1 = math.sqrt(sum((I * x) ** 2 for I, x in zip((I1, I2, I3), w)))
    print(f"spin about axis {axis + 1}: reversals in 3 s at {flips[:4]}, |L| drift {abs(L1 - L0) / L0:.1e}")
lam = 20.0 * math.sqrt((I2 - I1) * (I3 - I2) / (I1 * I3))
print(f"predicted growth rate of the wobble about axis 2: {lam:.2f} per s")

# 3. Gyroscope precession: Omega = m g d / (I w_s)
mg, d, I_s, ws = 0.15 * 9.81, 0.04, 1.2e-4, 2 * math.pi * 50
Om = mg * d / (I_s * ws)
print(f"toy gyroscope at 50 rev/s: precession {Om:.3f} rad/s, one turn every {2 * math.pi / Om:.1f} s")

principal moments about the corner (g m^2): 0.1278, 1.7047, 1.8294
principal moments about the CM    (g m^2): 0.0827, 0.3761, 0.4567
spin about axis 1: reversals in 3 s at [], |L| drift 1.8e-15
spin about axis 2: reversals in 3 s at [0.66, 1.836], |L| drift 9.2e-16
spin about axis 3: reversals in 3 s at [], |L| drift 7.8e-15
predicted growth rate of the wobble about axis 2: 15.82 per s
toy gyroscope at 50 rev/s: precession 1.561 rad/s, one turn every 4.0 s

About its centre of mass, the phone-shaped box has its smallest moment, 0.083 g m², about the axis along its length, an intermediate moment of 0.376 g m² about the axis across its width, and its largest moment, 0.457 g m², about the axis through its face. About a corner the principal moments are larger, and the principal axes are tilted away from the edges, because the products of inertia are no longer zero.

Euler's equations

In an inertial frame the law of rotation is \(d\mathbf L/dt = \boldsymbol\tau\). But \(\mathsf I\) is constant only in a frame fixed in the body, which rotates at \(\boldsymbol\omega\). Applying the rotating-frame rule of lesson 10 to \(\mathbf L\), with the body frame as the rotating frame,

In principal-axis coordinates, \((d\mathbf L/dt)_{\text{body}} = (I_1\dot\omega_1, I_2\dot\omega_2, I_3\dot\omega_3)\), and writing out the cross product gives Euler's equations, published by Leonhard Euler in 1765:

Even with no torque, the angular velocity of a body changes, because \(\boldsymbol\omega\) and \(\mathbf L\) are not parallel in general and the body turns under its own fixed angular momentum. Two quantities are conserved in torque-free motion: the kinetic energy \(T = \tfrac12\sum I_i\omega_i^2\) and the magnitude \(L^2 = \sum I_i^2\omega_i^2\). Viewed in the body frame, the tip of \(\mathbf L\) must stay on both a sphere (fixed \(L\)) and an ellipsoid (fixed \(T\)); it moves along their intersection curves, called polhodes.

Stability of rotation: the tennis racket theorem

Spin a body almost exactly about principal axis 1, so \(\omega_1 = \Omega\) is large and \(\omega_2\), \(\omega_3\) are small. To first order \(\omega_1\) stays constant, and the other two equations become \(I_2\dot\omega_2 = (I_3 - I_1)\Omega\omega_3\), \(I_3\dot\omega_3 = (I_1 - I_2)\Omega\omega_2\). Differentiating the first and substituting the second,

If the coefficient is negative the small components oscillate and the spin is stable; if positive they grow exponentially and the spin is unstable. The coefficient is negative when \(I_1\) is either the largest or the smallest moment, since then the two factors have opposite signs. It is positive when \(I_1\) is the intermediate moment. Spin about the axes of largest and smallest moment is stable; spin about the intermediate axis is not. This is the intermediate axis theorem, also called the tennis racket theorem: toss a racket spinning about the axis that lies in the plane of the head and perpendicular to the handle, and it flips over every turn or so.

The program spins the phone-shaped box at 20 rad/s about each principal axis with a tiny wobble of 0.01 rad/s. About axes 1 and 3 nothing happens in three seconds. About axis 2 the spin reverses at 0.66 s and again at 1.84 s, flipping periodically, while \(|\mathbf L|\) stays constant to \(10^{-15}\). The predicted growth rate of the wobble, \(\Omega\sqrt{(I_2 - I_1)(I_3 - I_2)/I_1I_3}\), is 15.8 per second, so the wobble grows from 0.01 to the size of the main spin in about \(\ln(2000)/15.8 \approx 0.5\ \mathrm{s}\), consistent with the first flip.

In 1985 the Soviet cosmonaut Vladimir Dzhanibekov noticed aboard the Salyut 7 station that a wing nut spinning off a bolt in zero gravity flipped over at regular intervals, and the effect is often named after him. It matters in practice: the first American satellite, Explorer 1, launched in 1958 as a long thin cylinder spinning about its axis of smallest moment, began tumbling within hours because its flexible whip antennas dissipated energy. With energy loss and fixed angular momentum, a body drifts toward the lowest-energy rotation consistent with its \(L\), which is rotation about the axis of largest moment. Spin-stabilized satellites since then are designed to spin about that axis.

Free precession of a symmetric body

For a body with an axis of symmetry, \(I_1 = I_2\), and Euler's third equation gives \(\omega_3\) constant. The other two become \(\dot\omega_1 = -\Omega_p\omega_2\), \(\dot\omega_2 = \Omega_p\omega_1\) with \(\Omega_p = (I_3 - I_1)\omega_3/I_1\). The component of \(\boldsymbol\omega\) perpendicular to the symmetry axis rotates around it at \(\Omega_p\), as seen from the body. The Earth, slightly flattened with \((I_3 - I_1)/I_1 = 1/305\), should therefore wobble about its figure axis once every 305 days if its spin axis is slightly off. Euler predicted this in the eighteenth century. Seth Carlo Chandler found the wobble in 1891, with a period of about 433 days and an amplitude of a few metres at the poles. The difference from 305 days comes from the Earth's elasticity and its fluid oceans and core, which yield to the changing centrifugal forces.

Gyroscopes and precession under torque

A gyroscope is a rapidly spinning symmetric body. Suppose its spin angular momentum \(\mathbf L_s = I_s\omega_s\hat{\mathbf n}\) is horizontal, supported at one end of its axle at distance \(d\) from its centre of mass. Gravity exerts a horizontal torque of magnitude \(mgd\), perpendicular to \(\mathbf L_s\). In time \(dt\) the angular momentum changes by \(\boldsymbol\tau\,dt\), perpendicular to itself, so its direction turns rather than its length growing. The axis sweeps around horizontally at the precession rate

Faster spin means slower precession. The program's toy gyroscope, spinning at 50 revolutions per second, precesses at 1.56 rad/s, once every 4.0 s. The formula is the leading approximation for fast spin; a complete treatment includes a small extra angular momentum from the precession itself and gives nutation, a nodding of the axis superimposed on the precession, which friction quickly damps. Gyroscopes resist changes in direction, which is why they once guided ships, aircraft and spacecraft, and why a fast bicycle wheel is hard to tilt. A spinning bullet or football flies point-first for the same reason.

The Earth is a gyroscope too. The Sun and Moon pull on its equatorial bulge, exerting a torque that tries to straighten the tilted axis. The response is precession: the axis traces a cone of half-angle \(23.4^\circ\) every 25,800 years or so. This precession of the equinoxes, discovered by Hipparchus around 130 BCE from ancient star positions, means that Polaris will not always be the pole star; in about 12,000 years the north celestial pole will lie near Vega.