Rigid Body Dynamics

Chapter 2 Rigid Body Dynamics 2.1 Coordinates of a Rigid Body A set of N particles forms a rigid body if the distance between any 2 particles is fixed: rij ≡ |ri − rj | = cij = constant. (2.1) Given these constraints, how many generalized coordinates are there? If we know 3 non-collinear points in the body, the remaing points are fully determined by triangulation. The first point has 3 coordinates for translation in 3 dimensions. The second point has 2 coordinates for spherical rotation about the first point, as r12 is fixed. The third point has one coordinate for circular rotation about the axis of r12 , as r13 and r23 are fixed. Hence, there are 6 independent coordinates, as represented in Fig. 2.1. This result is independent of N , so this also applies to a continuous body (in the limit of N → ∞). Figure 2.1: 3 non-collinear points can be fully determined by using only 6 coordinates. Since the distances between any two other points are fixed in the rigid body, any other point of the body is fully determined by the distance to these 3 points. The translations of the body require three spatial coordinates. These translations can be taken from any fixed point in the body. Typically the fixed point is the center of mass (CM), defined as: 1 X R= mi ri , (2.2) M i where mi is the mass of the i-th particle and ri the position of that particle with respect to a fixed origin and set of axes (which will notationally be unprimed) as in Fig. 2.2. In the case of a continuous body, this definition generalizes as: 1 R= M Z rρ(r) dV, (2.3) V where ρ(r) is the mass density at position r and we integrate over the volume V. Figure 2.2: The three translational coordinates correspond to the position of the Center of Mass, and the three rotational coordinates correspond to the three angles necessary to define the orientation of the axis fixed with the body. Rotations of the body are taken by fixing axes with respect to the body (we will denote these body fixed axes with primes) and describing their orientation with respect to the unprimed axes by 3 angles (φ, θ, ψ). A particularly useful choice of angles are called Euler angles. The angle φ is taken as a rotation about the z-axis, forming new x̃- and ỹ-axes while leaving the z-axis unchanged, as shown in Fig. 2.3. The angle θ is then taken as a rotation about the x̃-axis, forming new ỹ 0 - and z 0 -axes while leaving the x̃-axis unchanged, as shown in Fig. 2.4. Finally, the angle ψ is taken as a rotation about the z 0 -axis, forming new x0 - and y 0 -axes while leaving the z 0 -axis unchanged, as shown in Fig. 2.5. (The x̃-axis is called the line of nodes, as it is the intersection of the xy- and x̃ỹ-planes.) Figure 2.3: First rotation is by φ around the original z axis. Figure 2.4: Second rotation is by θ around the intermediate x̃ axis. Figure 2.5: Final rotation is by ψ around the final z 0 axis. Rotations can be described by 3 × 3 matrices U . This means each rotation step can be described as a matrix multiplication. Where r = (x, y, z), then    cos(φ) sin(φ) 0 x r̃ = Uφ r = − sin(φ) cos(φ) 0 y  . (2.4) 0 0 1 z Similar transformations can be written for the other terms: r̃0 = Uθ r̃ , r0 = Uψ r̃0 = Uψ Uθ r̃ = Uψ Uθ Uφ r. Defining the total transformation as U , it can be written as: U ≡ Uψ Uθ Uφ ⇒ r0 = U r. (2.5) Care is required with the order of the terms since the matrices don’t commute. Writing U out explicitly:     cos(ψ) sin(ψ) 0 1 0 0 cos(φ) sin(φ) 0 (2.6) U = − sin(ψ) cos(ψ) 0 0 cos(θ) sin(θ)  − sin(φ) cos(φ) 0 . 0 0 1 0 − sin(θ) cos(θ) 0 0 1 All rotation matrices, including Uφ , Uθ , Uψ , and U are orthogonal. Orthogonal matrices W satisfy W > W = W W > = 1 ⇔ W > = W −1 , (2.7) where 1 refers to the identity matrix and > to the transpose. This ensures that the length of a vector is invariant under rotations: r02 = r> (W > W )r = r2 . (2.8) Orthogonal matrices W have 9 entries but need to fulfill 6 conditions from orthogonality, leaving only 3 free parameters, corresponding to the 3 angles necessary to determine the rotation. We can also view r0 = U r as a transformation from the vector r to the vector r0 in the same coordinate system. This is an active transformation, as opposed to the previous perspective which was a passive transformation. Finally, note that inversions like     −1 0 0 1 0 0 U =  0 −1 0  or U = 0 1 0  (2.9) 0 0 −1 0 0 −1 are not rotations. These have det(U ) = −1, so they can be forbidden by demanding that det(U ) = 1. All orthogonal matrices have det(W ) = ±1 because det(W > W ) = (det(W ))2 = 1. In the language of group theory, the restriction to det(W ) = 1 gives the special orthogonal group SO(3) as opposed to simply O(3), the orthogonal group. We disregard the det(U ) = −1 subset of transformations because it is impossible for the system to undergo these transformations continuously without the distance between the particles changing in the process, so it would no longer be a rigid body. Intuitively, we could rotate the coordinates (x, y, z) directly into the coordinates (x0 , y 0 , z 0 ) by picking the right axis of rotation. In fact, the Euler theorem states that a general displacement of a rigid body with one point fixed is a rotation about some axis. This theorem will be true if a general rotation U leaves some axis fixed, which is satisfied by Ur = r (2.10) for any point r on this axis. This is an eigenvalue equation for U with eigenvalue 1. To better understand this, we need to develop a little linear algebra. Although the notion of an eigenvalue equation generally holds for linear operators, for now the discussion will be restricted to orthogonal rotation matrices U . The eigenvalue equation is U ξ = λξ, (2.11) where ξ is an eigenvector and λ is the associated eigenvalue. Rewriting this as (U − λ1)ξ = 0 (2.12) requires that det(U − λ1) = 0, so that U − λ1 is not invertible and the solution can be non-trivial, ξ = 6 0. det(U − λ1) = 0 is a cubic equation in λ, which has 3 solutions, which are the eigenvalues λα for α ∈ {1, 2, 3}. The associated eigenvectors are ξ (α) and satisfy U ξ (α) = λα ξ (α) , (2.13) where no implicit sum over repeated indices is taken. Forming a matrix from the resulting eigenvectors as columns:   ↑ ↑ ↑ X = ξ (1) ξ (2) ξ (3)  (2.14) ↓ ↓ ↓ then we can rewrite Eq.(2.13) as U X = X · diag(λ1 , λ2 , λ3 ) ⇒ X −1 U X = diag(λ1 , λ2 , λ3 ) (2.15) This means X diagonalizes U . Since U is orthogonal, the matrix X is unitary (X † X = XX † = 1). Note that > indicates transposition whereas † indicates Hermitian conjugation (i.e. complex conjugation ? combined with transposition > ). Next we note that since det(U ) = 1, then λ1 λ2 λ3 = 1. Additionally, |λα |2 = 1 for each α because: U ξ = λξ ⇒ ξ † U > = λ? ξ † ⇒ ξ † ξ = ξ † U > U ξ = |λ|2 ξ † ξ . (2.16) Finally, if λ is an eigenvalue, then so is λ? : U ξ = λξ ⇒ U ξ ? = λ? ξ ? (2.17) where ξ ? is still a column vector but with its elements undergoing complex conjugation with respect to ξ. Without loss of generality, let us say for a rotation matrix U that λ2 = λ?3 . Then 1 = λ1 |λ2 |2 = λ1 , so one of the eigenvalues is 1, giving Eq.(2.10), and thus proving Euler’s Theorem. The associated eigenvector ξ (1) to the eigenvalue λ1 = 1 is the rotation axis, and if λ2 = λ?3 = eiΦ then Φ is the rotation angle about that axis. In fact, we can make good use of our analysis of Euler’s theorem. Together the rotation axis and rotation angle can be used to define the instantaneous angular velocity ω(t) such that: |ω| = Φ˙ and ω k ξ (1) . (2.18) The angular velocity will play an important role in our discussion of time dependence with rotating coordinates in the next section. If we consider several consecutive displacements of ˙ so ω changes at each the rigid body, then each can have its own axis ξ (1) and its own Φ, instance of time, and hence ω = ω(t) (for the entire rigid body). 2.2 Time Evolution with Rotating Coordinates Lets use unprimed axes (x, y, z) for the fixed (inertial) axes, with fixed basis vectors ei . We will also use primed axes (x0 , y 0 , z 0 ) for the body axes with basis vectors e0i . If we consider any vector then it can be decomposed with either set of basis vectors: X X b= bi e i = b0i e0i . (2.19) i i For fixed axes basis vectors by definition ėi = 0, while for those in the body frame, ė0i = ω(t) × e0i (2.20) meaning vectors of fixed length undergo a rotation at a time t. The derivation of this result is shown in the figure on the right, by considering the change to the vector after an infinitesimal time interval dt. Summing over repeated indices, this means: b˙ = b˙ i ei = b˙ 0i e0i + ω(t) × (b0i e0i ) = b˙ 0i e0i + ω(t) × b dR d the time evolution in as the time evolution in the fixed (F) frame and dt dt the rotating/body (R) frame, then vectors evolve in time according to Defining db dR b = + ω × b. dt dt (2.21) As a mnemonic we have the operation “(d/dt) = dR /dt + ω×” which can act on any vector. Let us apply this to the position r of a particle of mass m, which gives dr dR r +ω×r = dt dt ⇔ vF = vR + ω × r. (2.22) Taking another time derivative gives us the analog for acceleration: F dvF dR v F = = + ω × vF m dt dt dR v R dR ω dR r = + ×r+ω× + ω × vR + ω × (ω × r) . dt dt dt Rr As ddt = vR is the velocity within the rotating body frame and within the body frame, then dR vR dt maR = F − mω × (ω × r) − 2mω × vR − m (2.23) = aR is the acceleration dR ω ×r dt (2.24) gives the acceleration in the body frame with respect to the forces that seem to be present in that frame. The terms −mω × (ω × r) and −2mω × vR are, respectively, the centrifugal and Coriolis ficticious forces respectively, while the last term −m dRdtω × r is a ficticious force that arises from non-uniform rotational motion, so that there is angular acceleration within the body frame. The same result could also have been obtained with the Euler-Lagrange equations for L in the rotating coordinates: L= m (ṙ + ω × r)2 − V , 2 (2.25) and you will explore this on a problem set. Note that the centrifugal term is radially outward and perpendicular to the rotation axis. To see this, decompose r into components parallel and perpendicular to ω, r = rk + r⊥ , then ω × r = ω × r⊥ , so −ω × (ω × r⊥ ) = w2 r⊥ . This term is present for any rotating body. On the other hand, the Coriolis force is nonzero when there is a nonzero velocity in the rotating/body frame: vR = 6 0. Example: Consider the impact of the Coriolis force on projectile motion on the rotating 2π ≈ 7.3 × 10−5 s−1 . We work out the Earth, where the angular velocity is ωEarth = 24×3600 s cross-product −ω × vr as shown in Fig. 2.6 for a particle in the northern hemisphere, where ω points to the north pole. Thus a projectile in the northern/southern hemisphere would be perturbed to the right/left relative to its velocity direction vr . Figure 2.6: For a projectile, in the Northern Hemisphere, the Coriolis pushes it to its right, relative to its direction of motion. Example: Consider a Foucault pendulum which hangs from a rigid rod, but is free to oscillate in two angular directions, as shown in Fig. 2.2. For θ  1 and working to first order in the small ω, the result derived from the Coriolis force gives φ˙ ≈ ωEarth sin(λ). Here λ is the latitude angle measured from equator. The precession is clockwise in the northern hemisphere, and is maximized at the north pole where λ = 90◦ . (This is proven as a homework problem.) Example: Consider the Coriolis deflection of a freely falling body on Earth in the northern hemisphere. We use the coordinate system shown below, where z is perpendicular to the surface of the earth and y is parallel to the earth’s surface and points towards the north pole. Working to first order in the small ω gives us maR = mv̇R = −mgẑ − 2mω × v, (2.26) where the centrifugal terms of order O(ω 2 ) are dropped. As an initial condition we take v(t = 0) = v0 ẑ. The term −ω × v points along x̂, so: z̈ = −g + O(ω 2 ) ⇒ vz = v0 − gt Moreover implementing the boundary condition that ẋ(t = 0) = 0:  g  ẍ = −2(ω × v)x = −2ω sin(θ)vz (t) ⇒ ẋ = −2ω sin(θ) v0 t − t2 . 2 (2.27) (2.28) Taking also x(t = 0) = 0, and integrating one further time, the motion in the x direction is: v g  0 2 t − t3 . (2.29) x(t) = −2ω sin(θ) 2 6 Lets consider this effect for a couple simple cases. If the mass m is dropped from a height z(t = 0) = hmax with zero velocity, v0 = 0, then: g z = hmax − t2 2 (2.30) and the mass reaches the floor at time s t1 = 2hmax . g (2.31) From Eq.(2.28) we see that ẋ > 0 for all t, and that: x(t = t1 ) = 8ω sin(θ)h3max > 0. 3g 2 However, if the mass m is thrown up with an initial ż(t = 0) = v0 > 0 from the ground (z = 0), then : g z(t) = v0 t − t2 > 0. 2 (2.32) Here the particle rises to a maximum height z = v02 /(2g) at time t = v0 /g, and then falls back to earth. Using Eq.(2.28) we see that ẋ < 0 for all t. If t1 is the time it reaches the ground again (t1 = 2vg0 ), then: x(t = t1 ) = − 2.3 4ω sin(θ)v03 < 0. 3g 2 (2.33) Kinetic Energy, Angular Momentum, and the Moment of Inertia Tensor for Rigid Bodies Returning to rigid bodies, consider one built out of N fixed particles The kinetic energy is best expressed using CM coordinates, where R is the CM and we here take ri to be the displacement of particle i relative to the CM. Once again making sums over repeated subscripts as implicit, the kinetic energy (T ) of the system is given by: 1 ˙2 1 T = MR + mi ṙ2i . 2 2 (2.34) As the body is rigid, then points cannot translate relative to the body but can only rotate so that ṙi = ω × ri . The rotational kinetic energy is then  1 1 1  TR = mi ṙ2i = mi (ω × ri )2 = mi ω 2 r2i − (ω · ri )2 . 2 2 2 (2.35) Labeling Cartesian indices with a and b to reserve i and j for particle indices, then we can write out this result making the indicies all explict as 1X (2.36) TR = mi (δab r2i − ria rib )ωa ωb . 2 i,a,b It is convenient to separate out the parts in this formula that depend on the shape and distributions of masses in the body by defining the moment of inertia tensor Iˆ for the discrete body as X Iˆab ≡ mi (δab ri2 − ria rib ) . (2.37) i The analog for a continuous body of mass density ρ(r) is: Z ˆ Iab ≡ (r2 δab − ra rb )ρ(r) dV. (2.38) V In terms of the moment of inertia tensor, the kinetic energy from rotation can now be written as: 1Xˆ 1 TR = (2.39) Iab ωa ωb = ω · Iˆ · ω , 2 a,b 2 where in the last step we adopt a convenient matrix multiplication notation. The moment of inertia tensor can be written with its components as a matrix in the form  2  yi + zi2 −xi yi −xi zi X Iˆ = mi  −xi yi x2i + zi2 −yi zi  , (2.40) 2 2 i −xi zi −yi zi xi + yi where the diagonal terms are the “moments of inertia” and the off-diagonal terms are the “products of inertia”. Note also that Iˆ is symmetric in any basis, so Iˆab = Iˆba . Special case: if the rotation happens about only one axis which can be defined as the z-axis for convenience so that ω = (0, 0, ω), then TR = 12 Iˆzz ω 2 which reproduces the simpler and more familiar scalar form of the moment of inertia. Lets now let ri be measured from a stationary point in the rigid body, which need not necessarily be the CM. The angular momentum can be calculated about this fixed point. Since vi = ω × ri , we can write the angular momentum as:   L = mi ri × vi = mi ri × (ω × ri ) = mi ri2 ω − (ω · ri )ri . (2.41) Writing out the components X X X La = mi (r2i ωa − (ω · ri )ria ) = ωb mi (δab r2i − ria rib ) = Iˆab ωb , i i,b b (2.42) which translates to the matrix equation: L = Iˆ · ω . (2.43) This allows us to write the corresponding rotational kinetic energy as: 1 TR = ω · L. (2.44) 2 Note that in general, L is not parallel to ω. We will see an explicit example of this below. Also note that the formula that we used for Iˆ in this case is the same as we gave above. We use these formulas whether or not ri is taken with respect to the CM. It is useful to pause to see what precisely the calculation of Iˆ depends on. Since it involves components of the vectors ri it depends on the choice of the origin for the rotation. Furthermore the entries of the matrix Iˆab obviously depend on the orientation of the axes used to define the components labeled by a and b. Given this, it is natural to ask whether given the result for Iˆab with one choice of axes and orientation, whether we can determine an Iˆ0 a0 b0 for a different origin and axes orientation. This is always possible with the help of a couple of theorems. The parallel axis theorem: Given IˆCM about the CM, it is simple to find IˆQ about a different point Q with the same orientation for the axes. Referring to the figure below, we define r0i as the coordinate of a particle i in the rigid body with respect to point Q and ri to be the coordinate of that particle with respect to the CM, so that: r0i = R + ri . By definition of the CM: X i mi ri = 0 and we let (2.45) M= X i mi . (2.46) The tensor of inertia around the new axis is then: Q 0 0 Iˆab = mi (δab r02 i − ria rib ) = mi (δab (r2i + 2ri · R + R2 − ria rib − ria Rb − Ra rib − Ra Rb ) , (2.47) (2.48) where the cross terms involving a single ri or single component ria sum up to zero by Eq.(2.46). The terms quadratic in r are recognized as giving the moment of inertia tensor about the CM. This gives the parallel axis theorem for translating the origin: Q CM Iˆab = M (δab R2 − Ra Rb ) + Iˆab , (2.49) If we wish to carry out a translation between P and Q, neither of which is the CM, then we can simply use this formula twice. Another formula can be obtained by projecting the parallel axis onto a specific axis n̂ where n̂2 = 1 (giving a result that may be familiar from an earlier classical mechanics course): n̂ · IˆQ · n̂ = M (R2 − (n̂ · R)2 ) + n̂ · IˆCM · n̂ = M R2 [1 − cos2 (θ)] + n̂ · IˆCM · n̂ = M R2 sin2 (θ) + n̂ · IˆCM · n̂ (2.50) where n̂ · R ≡ R cos(θ). Example: Lets consider an example of the calculation of Iˆ for a situation where L is not parallel to ω. Consider a dumbbell made of 2 identical point passes m attached by a massless rigid rod (but with different separations r1 and r2 from the axis of rotation), spinning so that ω = ωẑ and so that the rod makes an angle α with the axis of rotation, as shown We define body axes where the masses lie in the yz-plane. Here, r1 = (0, r1 sin α, r1 cos α) and r2 = (0, −r2 sin α, −r2 cos α). (2.51) Then using the definition of the moment inertia tensor: Izz Ixx Iyy Ixy Ixz Iyz = m(x12 + y12 ) + m(x22 + y22 ) = m(r12 + r22 ) sin2 α = m(y12 + z12 ) + m(y22 + z22 ) = m(r12 + r22 ) = m(x21 + z12 ) + m(x22 + z22 ) = m(r12 + r22 ) cos2 α = Iyx = m(−x1 y1 − x2 y2 ) = 0 = Izx = m(−x1 z1 − x2 z2 ) = 0 = Izy = m(−y1 z1 − y2 z2 ) = −m(r12 + r22 ) sin α cos α (2.52) Plugging these into L = Iˆ · ω, recalling that only ωz is non-zero, this gives L = (0, Iyz ω, Izz ω). (2.53) Thus in this example L is not parallel to ω. Next, instead of translating the axes in a parallel manner, let us keep the origin fixed and rotate the axes according to an orthogonal rotation matrix U satisfying U > U = U U > = 1. Vectors are rotated as L0 = U L , ω0 = U ω and therefore ω = U > ω 0 . (2.54) Putting these together ˆ >) · ω0 L0 = U Iˆ · ω = (U IU ⇒ ˆ >, Iˆ0 = U IU (2.55) where Iˆ0 is the new moment of inertia tensor. (The fact that it transforms this way defines it as a tensor.) This allows us to calculate the new moment of inertia tensor after a rotation. ˆ there always exists a rotation from an orthogonal matrix For a real symmetric tensor I, U that diagonalizes Iˆ giving a diagonal matrix Iˆ0 :   I1 0 0 IˆD =  0 I2 0  . (2.56) 0 0 I3 The entries of the diagonal moment of inertia tensor, Iα , are real and positive. This is just a special case of saying a Hermitian matrix can always be diagonalized by a unitary transformation (which is often derived in a Quantum Mechanics course as part of showing that a Hermitian matrix has real eigenvalues and orthogonal eigenvectors). The positivity of diagonal matrix follows immediately from the definition of the moment of inertia tensor for the situation with zero off-diagonal terms. The axes that make Iˆ diagonal are called the principal axes and the components Iα are the principal moments of inertia. We find them by solving the eigenvalue problem Iˆ · ξ = λ ξ, (2.57) where the 3 eigenvalues λ give the principal moments of inertia Iα , and are obtained from (α) solving det(Iˆ − λ1) = 0. The corresponding are the  (1)  3 real orthogonal eigenvectors ξ > (2) (3) principal axes. Here U = ξ ξ ξ , where the eigenvectors vectors fill out the columns. Then, without summing over repeated indices: 1X (2.58) Iα ωα2 , Lα = Iα ωα and T = 2 α where Lα and ωα are the components of L and ω, respectively, evaluated along the principal axes. To summarize, for any choice of origin for any rigid body, there is a choice of axes that ˆ For T to separate into translational and rotational parts, we must pick the diagonalizes I. origin to be the CM. Often, the principal axes can be identified by a symmetry of the body. Example: for a thin rectangle lying in the yz-plane with one edge coinciding with the z-axis, and the origin chosen as shown below, then Iyz = 0 as the body is symmetric under z ↔ −z, while Ixz = Ixy = 0 as the body lies entirely within x = 0. Hence these are principal axes. Sometimes, symmetry allows multiple choices for the principal axes. Example: for a sphere, any orthogonal axes through the origin are principal axes. Example: for a cylinder whose central axis is aligned along the z-axis, because of rotational symmetry any choice of the x- and y-axes gives principal axes. Example: Lets consider an example where the principal axes may not be apparent, which we can solve through the eigenvalue problem. Consider a uniform cube with sides of length a, mass m, and having the origin at one corner, as shown below. By symmetry we have Ixx = Iyy = Izz Ixy = Iyz = Ixz Z Z Z m a a a 2 2 = 3 (x + y 2 ) dx dy dz = ma2 , a 0 0 0 3 Z aZ aZ a m 1 −xz dx dy dz = − ma2 . = 3 4 a 0 0 0 (2.59) Thus the matrix is  + 23 − 14 − 14    Iˆ = ma2 − 14 + 23 − 14  . − 14 − 14 (2.60) + 23 The principal moments of inertia are found from  2   11 1 2 2 det(Iˆ − λ1) = ma − λ ma − λ = 0 . 12 6 (2.61) This gives I1 = λ1 = 16 ma2 . Solving (Iˆ − λ1 1)ξ (1) = 0 we find ξ (1) = (1, 1, 1). (2.62) The remaining eigenvalues are degenerate: I2 = I3 = λ2 = λ3 = 1 ma2 12 (2.63) ˆ 2 1)ξ (2,3) = so there is some freedom in determining the corresponding principal axes from (I−λ (1) 0, though they still should be orthogonal to each other (and ξ ). One example of a solution is: ξ (2) = (1, −1, 0) and ξ (3) = (1, 1, −2) (2.64) Using these principal axes and the same origin, the moment of inertia tensor becomes   1 0 0 ma2  1  IˆD = (2.65) 0 2 0  . 6 1 0 0 2 In contrast, if we had chosen the origin as the center of the cube, then one choice for the principal axes would have the same orientation, but with IˆCM = 16 ma2 1. This result could be obtained from Eq. (2.65) using the parallel axis theorem. 2.4 Euler Equations Consider the rotational motion of a rigid body about a fixed point (which could be the CM but could also be another point). We aim to describe the motion of this rigid body by exploiting properties of the body frame. To simplify things as much as possible, for this fixed point, we choose the principal axes fixed in the body frame indexed by α ∈ {1, 2, 3}. Using the relation between time derivatives in the inertial and rotating frames, the torque is then given by: dL dR L τ = +ω×L (2.66) = dt dt where ω = ω(t). For example: τ1 = dR L1 + ω2 L3 − ω3 L2 . dt (2.67) Not summing over repeated indices and using the formula for angular momentum along the principal axes gives Lα = Iα ωα . Since we have fixed moments of inertia within the body we have dR Iα /dt = 0. Note that dω/dt = dR ω/dt + ω × ω = dR ω/dt, so its rotating and inertial time derivatives are the same, and we can write ω̇α without possible cause of confusion. Thus dR Lα /dt = Iα ω̇α . This yields the Euler equations: I1 ω̇1 − (I2 − I3 )ω2 ω3 = τ1 I2 ω̇2 − (I3 − I1 )ω3 ω1 = τ2 I3 ω̇3 − (I1 − I2 )ω1 ω2 = τ3 (2.68) where in all of these ω and τ are calculated in the rotating/body frame. This can also be written as τα = Iα ω̇α + αlk ωl ωk Ik , (2.69) with α fixed but a sum implied over the repeated l and k indicies. Here abc is the fully antisymmetric Levi-Civita symbol. Solving these equations gives ωα (t). Since the result is expressed in the body frame, rather than the inertial frame of the observer, this solution for ω(t) may not always make the physical motion transparent. To fix this we can connect our solution to the Euler angles using the relations ω1 = ωx0 = φ˙ sin θ sin ψ + θ˙ cos ψ , ω2 = ωy0 = φ˙ sin θ cos ψ − θ̇ sin ψ , ˙ ω3 = ωz0 = φ˙ cos θ + ψ. (2.70) These results should be derived as exercise for the student. Example: let us consider the stability of rigid-body free rotations (τ = 0). Is a rotation ω = ω1 e1 about the principal axis e1 stable? Perturbations can be expressed by taking ω = ω1 e1 + κ2 e2 + κ3 e3 , where κ2 and κ3 are small and can be treated to 1st order. The Euler equations are: ω̇1 = (I2 − I3 ) κ2 κ3 = O(κ2 ) ≈ 0, I1 (2.71) so ω1 is constant at this order, and κ̇2 = (I3 − I1 ) ω1 κ3 I2 and κ̇3 = (I1 − I2 ) ω1 κ2 . I3 Combining these two equations yields   (I3 − I1 )(I1 − I2 )ω12 κ̈2 = κ2 . I2 I3 (2.72) (2.73) The terms in the square bracket are all constant, and is either negative = −w2 with an oscillating solution κ2 ∝ cos(wt + φ), or is positive = α2 with exponential solutions κ2 ∝ aeαt + be−αt . If I1 < I2,3 or I2,3 < I1 then the constant prefactor is negative, yielding stable oscillatory solutions. If instead I2 < I1 < I3 or I3 < I1 < I2 then the constant prefactor is positive, yielding an unstable exponentially growing component to their solution! This behavior can be demonstrated by spinning almost any object that has three distinct principal moments of inertia. 2.5 Symmetric Top with One Point Fixed This section is devoted to a detailed analysis of a particular example that appears in many situations, the symmetric top with one point fixed, acted upon by a constant force. Labeling the body axes as (x, y, z) and the fixed axes as (xI , yI , zI ), as depicted in the right, symmetry implies that I1 = I2 , and we will assume that I1,2 = 6 I3 . The Euler angles are as usual (φ, θ, ψ). From the figure we see that ψ˙ is the rotation rate of the top about the (body) z-axis, φ˙ is the precession rate about the zI fixed inertial axis, and θ˙ is the nutation rate by which the top may move away or towards the zI axis. The Euler equations in this case are I1 ω̇1 − (I2 − I3 )ω2 ω3 = τ1 , I1 ω̇2 − (I3 − I1 )ω3 ω1 = τ2 , I3 ω̇3 = 0 = τ3 . (2.74) Since the CM coordinate R is aligned along the z-axis there is no torque along z, τ3 = 0, leading to a constant ω3 . There are two main cases that we will consider. Case: τ = 0 and θ˙ = 0 The first case we will consider is when τ = 0 (so there is no gravity) and θ˙ = 0 (so there is no nutation). Then dL = τ = 0 ⇒ L = constant (2.75) dt Let us define the constant: I3 − I1 Ω≡ ω3 . (2.76) I1 Then the Euler equations for this situation reduce to: ω̇1 + Ωω2 = 0 and ω̇2 − Ωω1 = 0 . (2.77) The simplest solution correspond to ω1 (t) = ω2 (t) = 0, where we just have a rotation about the z-axis. Here: L = L3 e3 where L3 = I3 ω3 ω1 = ω2 = 0 ⇒ θ˙ = φ˙ = 0 and ψ˙ = ω3 . (2.78) In this case L k ω. A more general situation is when L and ω are not necessarily parallel, and ω1 and ω2 do not vanish. In this case Eq. (2.77) is solved by: ω1 = C sin(Ωt + D) and ω2 = −C cos(Ωt + D). (2.79) The simple case corresponds to C = 0, so now we take C > 0 (since a sign can be accounted for by the constant phase D). This solution means ω precesses about the body z-axis at the rate Ω, as pictured on the right. Since ω12 + ω22 is constant, the full ω = |ω| is constant, and is given by ω 2 = C 2 + ω32 . The energy here is just rotational kinetic energy TR = 1 ω · L which is constant too, since both ω and L are 2 constant. Thus ω also precesses about L. We can picture this motion by thinking about a body cone that rolls around a cone in the fixed coordinate system, where in the case pictured with a larger cone about L we have I1 = I2 > I3 . To obtain more explicit results for the motion we can relate Eq.(2.79) to Euler angles. Since θ̇ = 0, we take θ = θ0 to be constant. The other Euler angles come from:     sin(θ0 ) sin(ψ)φ˙ C sin(Ωt + D) (2.80) ω = −C cos(Ωt + D) = sin(θ0 ) cos(ψ)φ˙  . ˙ ˙ ω3 cos(θ0 )φ + ψ Adding the squares of the 1st and 2nd components gives C 2 = sin2 (θ0 ) φ˙ 2 . (2.81) To be definite, take the positive square root of this equation to give φ̇ = C sin(θ0 ) ⇒ φ= C t + φ0 . sin(θ0 ) (2.82) The first two equations in Eq. (2.80) are then fully solved by taking ψ = π − Ωt − D, so we find that both φ and ψ have linear dependence on time. Finally the third equation gives a relation between various constants ω3 = C cot(θ0 ) − Ω. (2.83) Thus, we see that the solution has φ˙ and ψ˙ are constants with θ˙ = 0. If we had picked the opposite sign when solving Eq. (2.81) then we would have found similar results: φ̇ = − C sin(θ0 ) ⇒ ψ = −Ωt − D and ω3 = −C cot(θ0 ) − Ω . (2.84) Case: τ = 6 0 and θ˙ = 6 0 Now we consider the general case where τ 6= 0 and θ˙ 6= 0. It is now more convenient to use the Lagrangian than the Euler equations directly. Since I1 = I2 , using   sin(θ) sin(ψ)φ˙ + cos(ψ)θ˙
1   T = I1 (ω12 + ω22 ) + I3 ω32 (2.85) and ω =  sin(θ) cos(ψ)φ˙ − sin(ψ)θ˙  , 2 cos(θ)φ˙ + ψ˙ gives us the kinetic energy T =  I  2 I1  ˙2 3 θ + sin2 θ φ˙ 2 + ψ˙ + cos θ φ˙ . 2 2 (2.86) Moreover, V = mgR cos(θ), so in the Lagrangian L = T − V the variables φ and ψ are cyclic. This means that the momenta  ∂L  (2.87) = I1 sin2 (θ) + I3 cos2 (θ) φ˙ + I3 cos(θ)ψ˙ pφ = ∂φ˙ ∂L ˙ = I3 ω3 pψ = = I3 (ψ˙ + cos(θ)φ) (2.88) ∂ψ˙ are conserved (constant). Here pψ is same as the angular momentum L3 discussed in the case above. The torque is along the line of nodes, and pφ and pψ correspond to two projections of L that are perpendicular to this torque (i.e. along ẑI and ẑ). Additionally, the energy is given by  I  2 I1  ˙2 3 2 2 ˙ ˙ ˙ E =T +V = θ + sin (θ)φ + ψ + cos(θ)φ + mgR cos(θ) (2.89) 2 2 and is also conserved. Solving the momentum equations, Eq. (2.87), for φ˙ and ψ˙ gives pφ − pψ cos(θ) φ˙ = I1 sin2 (θ) pψ (pφ − pψ cos(θ)) cos(θ) ψ˙ = − . I3 I1 sin2 (θ) (2.90) Note that once we have a solution for θ(t) that these two equations then allow us to immediately obtain solutions for φ(t) and ψ(t) by integration. Eq. (2.90) can be plugged into the energy formula to give p2ψ I1 2 (pφ − pψ cos(θ))2 + E = θ̇ + + mgR cos(θ), 2 2I3 2I1 sin2 (θ) (2.91) which is a (nonlinear) differential equation for θ, since all other quantities that appear are simply constants. To simplify this result take u = cos(θ) so that: 1 − u2 = sin2 (θ) , u̇ = − sin(θ)θ˙ , θ˙2 = u̇2 . 1 − u2 Putting all this together gives:    
1 pφ − pψ u 2 2EI3 − p2ψ mgR u̇2 2 = − u 1−u − ≡ −Veff (u) , 2 2I1 I3 I1 2 I1 (2.92) (2.93) which is a cubic polynomial that we’ve defined to be the effective potential Veff (u). The solution to this from du (2.94) dt = ±p −2Veff (u) yields a complicated elliptic function, from which it is hard to get intuition for the motion. Veff u1 −1 θ =180 ◦ 0 u2 u 1◦ θ =0 Allowed region Figure 2.7: Allowed region for solutions for the top’s nutation angle θ that solve Eq. (2.95). Instead, we can look at the form of Veff (u), because 1 2 u̇ + Veff (u) = 0 2 (2.95) is the equation for the energy of a particle of unit mass m = 1, kinetic energy u̇2 /2, a potential Veff (u), and with vanishing total energy. The cubic equation will have in general three roots where Veff (u) = 0. Since the kinetic energy is always positive or zero, the potential energy must be negative or zero in the physical region, and hence the particle can not pass through any of the roots. The roots therefore serve as turning points. Furthermore, physical solutions are bounded by −1 ≤ (u = cos θ) ≤ 1. While the precise values for the roots will depend on the initial conditions or values of E, pψ , and pφ , we can still describe the solutions in a fairly generic manner. Consider two roots u1 and u2 (corresponding respectively to some angles θ1 and θ2 as u = cos(θ)) satisfying Veff (u1 ) = Veff (u2 ) = 0, where Veff (u) < 0 for u1 < u < u2 ; as shown in Fig. 2.7. We see that u1 and u2 correspond to the turning points of the motion. The region u1 < u < u2 corresponds to the region where the motion of our top lives and gives rise to a periodic nutation, where the solution bounces between the two turning points. Depending on the precise value of the various constants that appear in this Vef f this gives rise to different qualitative motions, with examples shown in Figs. 2.8–2.11. Recalling that φ˙ = (pφ − pψ u)/[I1 (1 − u2 )], we see that the possible signs for φ˙ will depend on pφ and pψ . In Fig. 2.8 the top nutates between θ1 and θ2 while always precessing in the same direction Figure 2.8: φ˙ > 0 Figure 2.9: φ˙ has both Figure 2.10: at θ2 we Figure 2.11: No nuhave φ˙ = 0, θ˙ = 0 signs tation with φ˙ > 0, whereas in Fig. 2.9 the precession is also in the backward direction, φ˙ < 0, for part of the range of motion. In Fig. 2.10 the top has φ˙ = 0 at θ2 , before falling back down in the potential and gaining φ˙ > 0 again. This figure also captures the case where we let go of a top at θ = θ2 ≥ 0 that initially has ψ˙ > 0 but φ˙ = 0. Finally in Fig. 2.11 we have the situation where there is no nutation oscillation because the two angles coincide, θ1 = θ2 .