mi [ω 2 ri ’ (ω · r i )2 ] = ω · L,

2

K=

2 2

i

using the expression (15.49) for L.

Let us summarize the steps to obtain the equations for simulation of the

angular motion of a rigid body:

(i) Determine the center of mass r cm and the total mass M of the body.

(ii) Determine the moment of inertia of the body and its principal axes. If

not obvious from the symmetry of the body, then ¬rst obtain I in an

arbitrary body-¬xed coordinate system (see (15.52)) and diagonalize

it subsequently.

(iii) Assume the initial (cartesian) coordinates and velocities of the con-

stituent particles are known in a given space-¬xed coordinate system.

From these determine the initial rotation matrix R and angular ve-

locity ω.

(iv) Determine the forces F i on the constituent particles, and the total

force F tot on the body.

(v) Separate the c.o.m. motion: F tot /M is the acceleration of the c.o.m.,

and thus of all particles. Subtract mi F tot /M from every force F i .

(vi) Compute the total torque T in the body-¬xed principal axes coordi-

nate system, best by ¬rst transforming the particle forces (after c.o.m.

correction) with R to the body-¬xed frame, to obtain F i . Then apply

F i — ri .

T= (15.70)

i

(vii) Determine the time derivative of ω from (15.67) and the time deriva-

tive of R from (15.61). Integrate both equations simultaneously with

an appropriate algorithm.

The di¬erential equation for the rotation matrix (last step) can be cast in

several forms, depending on the de¬nition of the variables used to describe

the matrix. In the next section we consider three di¬erent possibilities.

412 Lagrangian and Hamiltonian mechanics

Z

Z´ ‘ X´

ψ

Y

•

K

X

Figure 15.2 Euler angles de¬ning the rotation of coordinate axes XY Z to X Y Z .

For explanation see text.

15.7.2 Unit vectors

A straightforward solution of the equation of motion for R is obtained by

applying (15.61) directly to its nine components, i.e., to the body-¬xed basis

vectors a, b, c. There are two possibilities: either one rotates ω to the space-

¬xed axes and then applies (15.61):

ω = RT ω , (15.71)

a = ω—a

™ (15.72)

(and likewise for b and c), or one applies (15.63) term for term. Both

methods preserve the orthonormality of the three basis vectors, but in a

numerical simulation the matrix may drift slowly away from orthonormality

by integration errors. Provisions must be made to correct such drift.

15.7.3 Euler angles

The traditional description of rotational motion is in Euler angles, of which

there are three and no problem with drifting constraints occur. The Euler

angles do have other severe problems related to the singularity of the Euler

equations when the second Euler angle ‘ has the value zero (see below). For

that reason Euler angles are not popular for simulations and we shall only

give the de¬nitions and equations of motion for the sake of completeness.

The Euler angles •, ‘, ψ, de¬ned as follows (see Fig. 15.2). We rotate

XY Z in three consecutive rotations to X Y Z . First locate the line of

15.7 Rigid body motion 413

nodes K which is the line of intersection of the XY - and the X Y -planes

(there are two directions for K; the choice is arbitrary). Then:

(i) rotate XY Z around Z over an angle • until X coincides with K;

(ii) rotate around K over an angle ‘ until Z coincides with Z ;

(iii) rotate around Z over an angle ψ until K coincides with X .

Rotations are positive in the sense of a right-hand screw (make a ¬st and

point the thumb of your right hand in the direction of the rotation axis;

your ¬ngers bend in the direction of positive rotation). The rotation matrix

is

⎛ ⎞

cos • cos ψ+ sin • cos ψ+ + sin ‘ sin ψ

⎜ ’ sin • cos ‘ sin ψ cos • cos ‘ sin ψ ⎟

⎜ ⎟

⎜ ⎟

R = ⎜ ’ cos • sin ψ+ ’ sin • sin ψ+ + sin ‘ cos ψ ⎟ . (15.73)

⎜ ⎟

⎝ ’ sin • cos ‘ cos ψ cos • cos ‘ cos ψ ⎠

’ cos • sin ‘

sin • sin ‘ cos ‘

The equations of motion, relating the angular derivatives to the body-¬xed

angular velocities, can be derived from (15.61). They are

⎛ ⎞ ⎛ sin • cos ‘ cos • cos ‘ ⎞⎛ ⎞

•™ ωx

1

sin ‘ sin ‘

⎝ ‘ ⎠ = ⎝ cos • ⎠⎝ ωy ⎠ .

™ (15.74)

sin •

™ sin • cos •

’ sin ‘ ωz

ψ sin ‘

15.7.4 Quaternions

Quaternions12 [q0 , q1 , q2 , q3 ] = q0 + q1 i + q2 j + q3 k are hypercomplex num-

bers with four real components that can be viewed as an extension of the

complex numbers a + bi. They were invented by Hamilton (1844). The

2

normalized quaternions, with i qi = 1, can be conveniently used to de-

scribe 3D rotations; these are known as the Euler“Rodrigues parameters,

described a few years before Hamilton™s quaternions by Rodrigues (1840).

Subsequently the quaternions have received almost no attention in the me-

chanics of molecules,13 until they were revived by Evans (1977). Because

equations of motion using quaternions do not su¬er from singularities as

12 See Altmann (1986) for an extensive review of quaternions and their use in expressing 3D

rotations, including a survey of the historical development. Another general reference is Kyrala

(1967).

13 In the dynamics of macroscopic multi-body systems the Euler“Rodrigues parameters are well-

known, but not by the name “quaternion”. See for example Shabana (1989). They are com-

monly named Euler parameters; they di¬er from the three Rodrigues parameters, which are

de¬ned as the component of the unit vector N multiplied by tan( 1 φ) (see (15.79)).

2

414 Lagrangian and Hamiltonian mechanics

those with Euler angles do, and do not involve the computation of gonio-

metric functions, they have become popular in simulations involving rigid

bodies (Allen and Tildesley, 1987; Fincham, 1992).

The unit quaternions14 1 = [1, 0, 0, 0], i = [0, 1, 0, 0], j = [0, 0, 1, 0] and

k = [0, 0, 0, 1] obey the following multiplication rules ([q] is any quaternion):

1[q] = [q]1 = [q], (15.75)

i2 = j 2 = k 2 = ’1, (15.76)

ij = ’ji = k and cyclic permutations. (15.77)

A general quaternion can be considered as the combination [q0 , Q] of a scalar

q0 and a vector Q with components (q1 , q2 , q3 ). The multiplication rules then

imply (as the reader is invited to check) that

[a, A][b, B] = [ab ’ A · B, aB + bA + A — B]. (15.78)

According to Euler™s famous theorem,15 any 3D rotation with ¬xed origin