Rotor in a Reflective 3D Box

Introduction

The goal of this animation is to provide a rotational analog to the usual hard sphere collision problem. The simplest extended rotating model is that of the symmetric rigid rotor with massive spheres at each end of a thin rigid rod. Just as in the hard sphere collision problem, we will make the assumption that the forces of the collisions are always directed at the centers of the end spheres. This is equivalent to saying that, during the collision, there is no tangential friction at the surface of the end sphere. As a result of this assumption, spin about the axis of the rod cannot be induced by collisions since there is no torque moment about this axis. That is a great simplification, since then we do not have to solve the Euler equations that usually couple the spins about the three axes of the rotor even under torque-free conditions.

1. Response of a Rotor to an Impulse Applied to One of Its Ends

Before we start, I should say that the bold character symbols in both text and equations denote vectors or matrices. This avoids the need for putting arrows over the vectors and matrices.

The basic driver involved in hard object collisions is an impulse. An impulse is the product of a force times a very small time increment which, of course, leads to a change in momentum like mδv where m is the mass and δv is the change in the velocity. The time increment is small enough that there will be no significant rotation of the rotor or displacement of its center of mass within the duration of the impulse. The significant changes in angle and position will occur after the impulse as results of changes in angular spin rate or linear velocity.

The specific diagram for this problem is shown below. To make the problem primitive and simple, initial center of mass speed and rotation rate are zero. The final motion will be a combination of a center of mass velocity, δv_x, and a final rotation, at rate δω, about the center of mass. To start, the impulse is along the -x direction as shown.

Figure 0: Showing the rotor being hit by a mallet at an oblique angle with impulse P_x. The impulse will cause a combination of rotation about the Roll (R) and Yaw (Y) axes. The rotation about the pitch (P) axis will not be affected because there is no torque moment about that axis.

The change in x momentum is a combination of the change in center of mass velocity and rotation speed δω_z. Here P_x is the impulse F_xdt where F_x is the force in the -x direction; Note that P_x is negative.:

We have to take into account that there are 2 body axes about which the force impulse can have an effect. We will arbitrarily call these the roll axis and the yaw axis while the pitch axis will be along the length of the rod connecting the two end spheres. Initially the pitch axis will be the x axis while the roll axis will be along the y direction and the yaw will be along the z direction (out of and into of the screen). The program will track the orientations of these three axes and resolve the applied torque into its roll and yaw components.

The coordinates perpendicular to ρ can be resolved into parts, one of which is along the roll axis and one along the yaw axis. The pitch axis is not important since it can't receive any torque.

$\begin{array}{l} R = r_{x} x + r_{y} y + r_{z} z \\ Y = y_{x} x + y_{y} y + y_{z} z \end{array}$

(1.1)

While the body axes can be important when we have body-induced forces like on-board thrusters, for the cases analyzed here, all of the forces will come from external sources like impulses or by being hit by other bodies so we will just compute the torques and reaction changes in rotation rates and vectors in the laboratory (x,y,z) coordinate system.

Let's take the case where the impulse is along an arbitrary unit vector, c, not aligned with any of the xyz axes. Then we have the following equations:

$\begin{array}{l} δ v • \hat{c} + (\hat{c} \times ρ) • δ ω = \frac{P}{m} \\ δ ω • \hat{c} = 0 \\ δ ω • ρ = 0 \end{array}$

(1.2)

where δv is the velocity of the center of mass and ρ is the vector from the center of the rotor to one of its end spheres.

The first of these equations just says that the change in the component of velocity along the c direction is equal, by Newton's law, to P/m, where m is the total mass of the rotor. The second and third equation say that the change in the vector rotation rate, δω, is perpendicular to both c and ρ.

We can re-write equation (1.2) in terms of its xyz components:

$\begin{array}{l} δ v_{x} c_{x} + δ v_{y} c_{y} + δ v_{z} c_{z} + (_{c \times ρ) x} δ ω_{x} + (_{c \times ρ) y} δ ω_{y} + (_{c \times ρ) z} δ ω_{z} = \frac{P}{m} \\ c_{x} δ ω_{x} + c_{y} δ ω_{y} + c_{z} δ ω_{z} = 0 \\ ρ_{x} δ ω_{x} + ρ_{y} δ ω_{y} + ρ_{z} δ ω_{z} = 0 \end{array}$

(1.3)

Then the matrix and vector expression becomes:

$(\begin{matrix} c_{y} ρ_{z} - c_{z} ρ_{y} & - c_{x} ρ_{z} + c_{z} ρ_{x} & c_{x} ρ_{y} - c_{y} ρ_{x} \\ c_{x} & c_{y} & c_{z} \\ ρ_{x} & ρ_{y} & ρ_{z} \end{matrix}) (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) = (\frac{P}{m} - δ v) (\begin{matrix} 1 \\ 0 \\ 0 \end{matrix})$

(1.4)

$(\begin{matrix} {(c \times ρ)}_{x} & {(c \times ρ)}_{y} & {(c \times ρ)}_{z} \\ c_{x} & c_{y} & c_{z} \\ ρ_{x} & ρ_{y} & ρ_{z} \end{matrix}) (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) = (\frac{P}{m} - δ v) (\begin{matrix} 1 \\ 0 \\ 0 \end{matrix})$

(1.5)

Solving equation (1.5) we obtain:

$\begin{array}{l} (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) = \frac{(\frac{P}{m} - δ v)}{\det} (\begin{matrix} {(c \times ρ)}_{x} \\ {(c \times ρ)}_{y} \\ {(c \times ρ)}_{z} \end{matrix}) \\ \det = - (_{c \times ρ) x} ρ_{y} c_{z} + (_{c \times ρ) x} ρ_{z} c_{y} + (_{c \times ρ) y} r ρ_{x} c_{z} - (_{c \times ρ) z} ρ_{x} c_{y} - (_{c \times ρ) y} ρ_{z} c_{x} + (_{c \times ρ) z} ρ_{y} c_{x} \\ \det = (_{c \times ρ}^{)} + (_{c \times ρ}^{)} + (_{c \times ρ}^{)} \\ (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) = \frac{(\frac{P}{m} - δ v)}{(c \times ρ) • (c \times ρ)} (\begin{matrix} {(c \times ρ)}_{x} \\ {(c \times ρ)}_{y} \\ {(c \times ρ)}_{z} \end{matrix}) \\ δ ω = \frac{(\frac{P}{m} - δ v)}{(c \times ρ) • (c \times ρ)} (c \times ρ) \end{array}$ (1.6)

We need to insert the values of δω² into an energy conservation equation.

$\frac{1}{2} (m δ v_{c}^{2} + m l^{2} δ ω^{2}) = P \frac{δ v + (\hat{c} \times ρ) • δ ω}{2}$

(1.7)

and then solve the resulting quadratic equation for δv_c. Note that:

$\begin{array}{l} δ ω^{2} = \frac{{(\frac{P}{m} - δ v)}^{2}}{{(c \times ρ)}^{2}} \\ (c \times ρ) • δ ω = \frac{P}{m} - δ v \end{array}$

(1.8)

$\frac{1}{2} (m δ v^{2} + m l^{2} \frac{{(\frac{P}{m} - δ v)}^{2}}{(c \times ρ) • (c \times ρ)}) = P \frac{δ v + (\frac{P}{m} - δ v)}{2}$ (1.9)

The solutions for δv and δω are then:

$δ v = \frac{l^{2} - | c \times ρ |^{2}}{l^{2} + | c \times ρ |^{2}} \frac{P}{m}$

(1.10)

$\frac{P}{m} - δ v = 2 \frac{P}{m} \frac{| c \times ρ |^{2}}{l^{2} + | c \times ρ |^{2}}$

(1.11)

$\begin{array}{l} (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) = \frac{2 P}{m (l^{2} + {(c \times ρ)}^{2}]} (\begin{matrix} {(c \times ρ)}_{x} \\ {(c \times ρ)}_{y} \\ {(c \times ρ)}_{z} \end{matrix}) \\ δ ω = \frac{2 P}{m (l^{2} + {(c \times ρ)}^{2}]} (c \times ρ) \end{array}$

(1.12)

The new spin rate vector is ω+δω.

To incrementally change the orientation of the rotor about this new spin axis, we apply the Rodrigues matrix below:

$M_{R} = (\begin{matrix} c + u_{x}^{2} (1 - c) & u_{x} u_{y} (1 - c) - u_{z} s & u_{x} u_{z} (1 - c) + u_{y} s \\ u_{y} u_{x} (1 - c) + u_{z} s & c + u_{y}^{2} (1 - c) & u_{y} u_{z} (1 - c) - u_{x} s \\ u_{z} u_{x} (1 - c) - u_{y} s & u_{z} u_{y} (1 - c) + u_{x} s & c + u_{z}^{2} (1 - c) \end{matrix})$

(1.13)

where

$\begin{array}{l} c = \cos ω δ t \\ s = s i n ω δ t \\ u = \frac{ω}{| ω |} \end{array}$

(1.14)

and δt is the time increment.

2. Collision with a flat rigid surface

It is expected that the rigid surface collision will result in negation of the normal component of v, the total velocity due to both rotation and linear motion, but the energy of the rotor will not be changed. For this simulation it will be assumed that there is no friction between the wall and the rotor end. That is, the wall applies no impulse parallel to its surface.

Figure 1: Illustration of rotor hitting surface x=constant with rotational angle a and center of mass velocity v. The distance between the centers of the end Spheres is 2l and the radius of the end Spheres is b. The Roll (R) and Yaw (Y) axes are shown

In the case of a collision with a rigid wall, the impulse discussed in section 1 is just twice the negative of the x component of the incident momentum. That means that the normal component of the velocity of the end sphere hitting the x=+constant wall gets negated. For generality, let x=c where c is the unit vector along the outward normal to any of the 6 walls.

$v' • c + (ρ \times c) • ω' = - [v • x + (ρ \times c) • ω]$ (2.1)

where - $(c \times ρ) • ω$ is equal to the c component of the velocity due to the rotational motion of the rotor.

Also note that, by definition, ρ and ω are always perpendicular,

$ρ • ω' = ρ • ω = ρ • δ ω = 0$

(2.2)

Further, the component of ω parallel to c does not result in a change of ω.

$c • δ ω = 0$

(2.3)

All of these equations are valid regardless of the coordinate system in which c, ω and ρ are described. We simply choose to work in the xyz system and to convert, if needed, to body rotations afterward. Rewriting equations (2.1), (2.2) and (2.3) in matrix form and setting c=x for the general case we have:

$(\begin{matrix} {(ρ \times c)}_{x} & {(ρ \times c)}_{y} & {(ρ \times c)}_{z} \\ c_{x} & c_{y} & c_{z} \\ ρ_{x} & ρ_{y} & ρ_{z} \end{matrix}) (\begin{matrix} δ ω_{x} \\ δ ω_{z} \\ δ ω_{z} \end{matrix}) = (- v_{c} - v_{c}' - 2 (ρ \times c) • ω) (\begin{matrix} 1 \\ 0 \\ 0 \end{matrix})$

(2.4)

After taking the inverse of the matrix on the left, the solution for δω in equation (2.4) is:

$δ ω = (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) = \frac{(- v_{c} - v_{c}' - 2 (ρ \times c) • ω)}{(ρ \times c) • (ρ \times c)} (\begin{matrix} {(ρ \times c)}_{x} \\ {(ρ \times c)}_{y} \\ {(ρ \times c)}_{z} \end{matrix}) = \frac{(- v_{c} - v_{c}' - 2 (ρ \times c) • ω)}{(ρ \times c) • (ρ \times c)} (ρ \times c)$

(2.5)

The magnitude of the rotational rate after the collision is given by:

$ω = \sqrt{ω' • ω'}$

(2.6)

where

$ω' = ω + δ ω$

(2.7)

and the new rotation axis unit vector is:

$u_{ω} = \frac{ω'}{| ω' |} = \frac{ω'}{\sqrt{| (ω' • ω') |}}$

(2.8)

which will be the unit vector used in the Rodrigues matrix of the previous section.

We have one more expression from which to obtain the value of v_c'. The conservation of energy equation is:

$v_{c}^{' 2} - v_{c}^{2} + l^{2} [ω'^{2} - ω^{2}] = v_{c}^{' 2} - v_{c}^{2} + l^{2} [ω'_{R}^{2} + ω'_{Y}^{2} - ω_{R}^{2} - ω_{Y}^{2}] = 0$

(2.9)

where ω_R' and ω_Y' are linear functions of v'_c. Once we've inserted their v'_c dependence into equation (2.9), we can solve what will be a quadratic equation for v'_c.

The result is

$v_{c}' = - \frac{v_{c} (l^{2} - | ρ \times c |^{2}) + 2 l^{2} [(ρ \times c) • ω]}{l^{2} + | ρ \times c |^{2}}$

(2.10)

Then the new value of the CM speed is

$v'_{C M} = v_{C M} + (v_{c}' - v_{c}) c$

(2.11)

The value of δω after inserting the new value v'_CM into equation (2.5) is:

$\begin{array}{l} δ ω = \frac{(- v_{c} - 2 (ρ \times c) • ω + \frac{v_{c} (l^{2} - | ρ \times c |^{2}) + 2 l^{2} [(ρ \times c) • ω]}{l^{2} + | ρ \times c |^{2}})}{| ρ \times c |^{2}} (ρ \times c) = \\ \frac{([- v_{c} - 2 (ρ \times c) • ω] [l^{2} + | ρ \times c |^{2}] + v_{c} (l^{2} - | ρ \times c |^{2}) + 2 l^{2} [(ρ \times c) • ω])}{| ρ \times c |^{2} [l^{2} + | ρ \times c |^{2}]} (ρ \times c) = \\ - 2 \frac{(ρ \times c) • ω | ρ \times c |^{2} + v_{c} | ρ \times c |^{2}}{| ρ \times c |^{2} [l^{2} + | ρ \times c |^{2}]} (ρ \times c) = - 2 \frac{(ρ \times c) • ω + v_{c}}{l^{2} + | ρ \times c |^{2}} (ρ \times c) \end{array}$

(2.12)

(2.13)

Rotor-Sphere Collisions

Figure 2: Illustration of rotor and a free sphere. The rotor rotational angle is a and center of mass velocity v The distance between the centers of the rotor end spheres is 2l₀ and the radius of the end Spheres is b_r while the single Sphere radius is b_s Thus a collision is computed when the distance between Sphere center and either rotor end center is less than b_r+b_s

For the case of a Rotor and a free Sphere, the linear momentum associated with the Sphere and rotor center of mass (CM) velocity must be conserved and the angular momentum, L_CM, associated with their linear momentum about their CM as well as the simple rotational angular momentum, L_ω, of the Rotor must be conserved. Having used the linear and angular momentum conservation equations to compute the new values of Rotor rotation rate, we must then use the conservation of energy equation to compute the new values of their relative velocity with respect to their previous velocities. These results will always center around their c unit vector which is the vector between the center of the Sphere and the center of the rotor end with which it collides. Recall from the previous 2D Rotor calculations that the angular momentum associated with the linear motion of the 2 masses is:

$\begin{array}{l} L_{C M} = M_{r} (r_{r} - r_{s}) \times (v_{r} - v_{s}) \\ M_{r} = \frac{m_{r} m_{s}}{m_{r} + m_{s}} \end{array}$

(2.14)

where the subscripts r and s refer to the values for the rotor and free sphere and the ms are their masses, respectively.

The angular momentum of the rotor about its axis is:

$L_{ω} = m_{r} l^{2} ω$

(2.15)

and the linear momentum of the 2 masses is:

$P_{C M} = m_{r} v_{r} + m_{s} v_{s}$

(2.16)

The velocities of the two masses will change along the c vector. If we name the as yet unknown magnitude of the rotor CM speed change δv, then, from equation (2.16), the Sphere speed change will be -(m_r/m_s)δv. The expression for the conservation of angular momentum of is a bit more elaborate. The equation for the final angular momentum is:

$L' = M_{r} (r_{r} - r_{s}) \times [v_{r} + δ v c - (v_{s} - δ v c \frac{m_{r}}{m_{s}})] + m_{r} l^{2} ω'$

(2.17)

while the initial angular momentum is:

$L = M_{r} (r_{r} - r_{s}) \times [v_{r} - v_{s})] + m_{r} l^{2} ω$

(2.18)

To compute the final rotor rotational angular momentum, we need to set equations (2.17) and (2.18) equal to each other:

$\begin{array}{l} M_{r} (r_{r} - r_{s}) \times [v_{r} + δ v c - (v_{s} - δ v c \frac{m_{r}}{m_{s}})] + m_{r} l^{2} ω' = \\ M_{r} (r_{r} - r_{s}) \times [v_{r} - v_{s})] + m_{r} l^{2} ω \end{array}$

(2.19)

Solving equation (2.19) for the change in ω we obtain:

$m_{r} l^{2} (ω' - ω) = - M_{r} δ v (1 + \frac{m_{r}}{m_{s}}) [(r_{r} - r_{s}) \times c]$

(2.20)

$\begin{array}{l} m_{r} l^{2} (\begin{matrix} δ ω_{x} \\ δ ω_{z} \\ δ ω_{z} \end{matrix}) = - M_{r} δ v (1 + \frac{m_{r}}{m_{s}}) (\begin{matrix} ({[(r_{r} - r_{s}) \times c]}_{x} \\ {[(r_{r} - r_{s}) \times c]}_{y} \\ {[(r_{r} - r_{s}) \times c]}_{z} \end{matrix}) = - m_{r} δ v (\begin{matrix} {[(r_{r} - r_{s}) \times c]}_{x} \\ {[(r_{r} - r_{s}) \times c]}_{y} \\ {[(r_{r} - r_{s}) \times c]}_{z} \end{matrix}) \\ δ ω = - \frac{m_{r}}{l^{2}} δ v (r_{r} - r_{s}) \times c \end{array}$

(2.21)

With some vector algebra we can rewrite the right hand side of equation (2.21) in terms of ρ, the vector that points from the rotor axis to the center of its end sphere:

$\begin{array}{l} ρ = r_{r} - r_{s} + (b_{r} + b_{s}) c \\ (r_{r} - r_{s}) \times c = ρ \times c \end{array}$

(2.22)

which guarantees that the condition for constancy of spin angular momentum will exist:

$δ ω • ρ = - \frac{m_{r}}{l^{2}} δ v [(r_{r} - r_{s}) \times c] • ρ = - \frac{m_{r}}{l^{2}} δ v (ρ \times c) • ρ = 0$

(2.23)

The conservation of energy equation dictates that the following equation must be true:

$\begin{array}{l} m_{r} (v_{r} + δ v c) • (v_{r} + δ v c) + m_{r} l^{2} (ω + δ ω) • (ω + δ ω) + m_{s} (v_{s} - \frac{m_{r}}{m_{s}} δ v c) • (v_{s} - \frac{m_{r}}{m_{s}} δ v c) = \\ m_{r} v_{r}^{2} + m_{r} l^{2} ω^{2} + m_{s} v_{s}^{2} \end{array}$ (2.24)

$m_{r} {[2 (v_{r} • c) - 2 (v_{s} • c)] δ v + δ v^{2} (1 + {(\frac{m_{r}}{m_{s}})}^{2})} + m_{r} l^{2} [(2 ω • δ ω + δ ω^{2}) = 0$

(2.25)

$m_{r} {[2 (v_{r} • c) - 2 (v_{d} • c)] δ v + δ v^{2} (1 + {(\frac{m_{r}}{m_{s}})}^{2})} + m_{r} l^{2} [(- 2 ω • \frac{[(r_{r} - r_{s}) \times c]}{l} δ v + {(\frac{[(r_{r} - r_{s}) \times c]}{l} δ v)}^{2}) = 0$ (2.26)

The result for δv is:

$δ v = - 2 \frac{(v_{r} - v_{s}) • c - ω • [(r_{r} - r_{s}) \times c]}{1 + \frac{m_{r}}{m_{s}} + \frac{{[(r_{r} - r_{s}) \times c]}^{2}}{l^{2}}}$

(2.27)

Then δω becomes

$δ ω = - \frac{δ v}{l^{2}} (r_{r} - r_{s}) \times c = \frac{2}{l^{2}} {\frac{(v_{r} - v_{s}) • c - ω • [(r_{r} - r_{s}) \times c]}{1 + \frac{m_{r}}{m_{s}} + \frac{{[(r_{r} - r_{s}) \times c]}^{2}}{l^{2}}}} [(r_{r} - r_{s}) \times c]$

(2.28)

Rotor-Rotor Collisions

Figure 3: Illustration of 2 rotors colliding . The distance between the centers of the rotor end Spheres is 2l and the radius of the end Sphere is b_r Thus a collision is computed when the distance between Sphere center and either rotor end center is less than 2b_r. ρ is the vector from the axis of the rotor to the center of its colliding end Sphere, and c is a unit vector along the line from the center of rotor 2s end Sphere to the center of the end spheres of rotor 1.

For collisions of two Rotors, the linear momentum associated with their center of mass (CM) velocity must be conserved and the angular momentum, L_CM, associated with their linear momentum about their CM and their simple rotational angular momentum, L_ω, of the Rotors must be conserved. Having used the linear and angular momentum conservation equations to compute the new values of Rotors' rotation rates, we must then use the conservation of energy equation to compute the new values of their relative velocity with respect to their previous velocities. These results will always center around their c unit vector which is the vector between the center of one rotor end sphere and the center of the other rotor end with which it collides. Recall from the previous 2D Rotor calculations that the angular momentum associated with the linear motion of the 2 Rotor masses is:

$\begin{array}{l} L_{C M} = M_{r} (r_{1} - r_{2}) \times (v_{1} - v_{2}) \\ M_{r} = \frac{m_{1} m_{2}}{m_{1} + m_{2}} \end{array}$

(3.1)

The final angular momentum after the collision can be written:

$L' = M_{r} (r_{1} - r_{2}) \times [v_{1} + δ v c - (v_{2} - δ v c \frac{m_{1}}{m_{2}})] + m_{1} l^{2} ω_{1}' + m_{2} l^{2} ω_{2}'$

(3.2)

and the initial angular momentum before the collision can be written

$L = M_{r} (r_{1} - r_{2}) \times [v_{1} - v_{2})] + m_{1} l^{2} ω_{1} + m_{2} l^{2} ω_{2}$

(3.3)

We may set equations (3.2) and (3.3) equal getting the following equation:

$\begin{array}{l} M_{r} [(r_{1} - r_{2}) \times c] δ v (1 + \frac{m_{1}}{m_{2}}) + m_{1} l^{2} (ω_{1}' - ω_{1}) + m_{2} l^{2} (ω_{2}' - ω_{2}) = 0 \\ m_{1} [(r_{1} - r_{2}) \times c] δ v + m_{1} l^{2} (ω_{1}' - ω_{1}) + m_{2} l^{2} (ω_{2}' - ω_{2}) = 0 \end{array}$

(3.4)

First let's rewrite $(r_{1} - r_{2}) \times c$ in terms of ρ₁ and ρ₂. Note that the vector expression relating $r_{1}, r_{2} a n d c$ is

$r_{1} + ρ_{1} = r_{2} + ρ_{2} - 2 b_{r} c$

(3.5)

Then solving for $r_{1} - r_{2}$ we have

$\begin{array}{l} r_{1} - r_{2} = ρ_{2} - ρ_{1} - 2 b_{r} c \\ (r_{1} - r_{2}) \times c = (ρ_{2} - ρ_{1}) \times c \end{array}$

(3.6)

so we can rewrite equation (3.4)

$m_{1} [(ρ_{1} - ρ_{2}) \times c] δ v + m_{1} l^{2} (ω_{1}' - ω_{1}) + m_{2} l^{2} (ω_{2}' - ω_{2}) = 0$

(3.7)

Since the same (unknown) impulse, P, causes both ω₁ and ω₂ to change, there is a specific relationship between their changes.

Let both torque impulses be described in terms of the unknown value of P

$\begin{array}{l} δ τ_{1} = P (ρ_{1} \times c) \\ δ τ_{2} = - P (ρ_{2} \times c) \\ ω_{1}' - ω_{1} \equiv δ ω_{1} = \frac{δ τ_{1}}{m_{1} l^{2}} = \frac{P (ρ_{1} \times c)}{m_{1} l^{2}} \\ ω_{2}' - ω_{2} \equiv δ ω_{2} = \frac{δ τ_{2}}{m_{2} l^{2}} = - \frac{P (ρ_{2} \times c)}{m_{2} l^{2}} \end{array}$

(3.8)

Dot multiply both sides of the equation for $δ ω_{1}$ by $ρ_{1} \times c$ and we have:

$(ρ_{1} \times c) • δ ω_{1} = \frac{P (ρ_{1} \times c) • (ρ_{1} \times c)}{m_{1} l^{2}}$

(3.9)

Then the solution for P is gotten by solving equation (3.9)

$P = \frac{m_{1} l^{2} (ρ_{1} \times c) • δ ω}{(ρ_{1} \times c) • (ρ_{1} \times c)}$

(3.10)

And then we can use our result for P in the δω₂ equation:

$δ ω_{2} = - P \frac{(ρ_{2} \times c)}{m_{2} l^{2}} = - \frac{m_{1} (ρ_{2} \times c) • (ρ_{1} \times c) • δ ω}{m_{2} (ρ_{1} \times c) • (ρ_{1} \times c)}$

(3.11)

We may rewrite equation (3.11) as a matrix multiplied by δω. First, for brevity, rename some of the quantities in equation (3.11):

$\begin{array}{l} a_{1} = ρ_{1} \times c \\ a_{2} = \frac{m_{1}}{m_{2}} \frac{ρ_{2} \times c}{(ρ_{1} \times c) • (ρ_{1} \times c)} \\ s o t h a t : \\ δ ω_{2} = - a_{2} (a_{1} • δ ω) \end{array}$

(3.12)

Now writing out the components of v₂ and $v_{1} • δ ω$ we find the following matrix:

$δ ω_{2} = - a_{2} (a_{1} • δ ω) = - (\begin{matrix} a_{2 x} a_{1 x} & a_{2 x} a_{1 y} & a_{2 x} a_{1 z} \\ a_{2 y} a_{1 x} & a_{2 y} v_{1 y} & a_{2 y} a_{1 z} \\ a_{2 z} a_{1 x} & a_{2 z} v_{1 y} & a_{2 z} a_{1 z} \end{matrix}) (\begin{matrix} δ ω_{x} \\ δ ω_{y} \\ δ ω_{z} \end{matrix}) \equiv - M_{21} δ ω$

(3.13)

where M₂₁ is the dyadic tensor a₂a₁.

Then equation (3.7) can be solved as:

$m_{1} l^{2} δ ω (1 - M_{21}) = - m_{1} [(ρ_{1} - ρ_{2}) \times c] δ v$

(3.14)

so that the solution for δω is:

$δ ω \equiv δ ω_{1} = - {(1 - M_{21})}^{- 1} [(ρ_{1} - ρ_{2}) \times c] \frac{δ v}{l^{2}}$

(3.15)

where the bold 1 in equation (3.15) denotes the unit matrix.

Now we can write the conservation of energy equation for the two rotors

$\begin{array}{l} m_{1} (^{v_{1} + δ v) 2} + m_{2} (v_{2} - \frac{m_{2}}{m_{2}} δ v) + m_{1} l^{2} (^{ω_{1} + δ ω) 2} + m_{1} l^{2} (^{ω_{2} + \frac{m_{2}}{m_{1}} M_{21} δ ω) 2} = \\ m_{1} v_{1}^{2} + m_{2} v_{2}^{2} + m_{1} l^{2} ω_{1}^{2} + m_{1} l^{2} ω_{2}^{2} \end{array}$

(3.16)

Simplify notation for equation (3.15)

$δ ω = - {(1 - M_{21})}^{- 1} [(ρ_{1} - ρ_{2}) \times c] \frac{δ v}{l^{2}} \equiv - b \frac{δ v}{l^{2}}$

(3.17)

where

$b = {(1 - M_{21})}^{- 1} [(ρ_{1} - ρ_{2}) \times c]$

(3.18)

where a is a vector. Then solving equation (3.16) we obtain:

$δ v = 2 \frac{m_{1} (v_{2} - v_{1}) • c + m_{1} b • ω_{1} - m_{2} (M_{21} b) • ω_{2}}{m_{1} (1 + \frac{b^{2}}{l^{2}}) + \frac{m_{1}^{2}}{m_{2}} + m_{2} {(M_{21} b)}^{2} \frac{1}{l^{2}}}$

(3.19)

Recalling equation (3.15) for δω we have:

$δ ω = - {(1 - M_{21})}^{- 1} [(ρ_{1} - ρ_{2}) \times c] \frac{δ v}{l^{2}} = - b \frac{δ v}{l^{2}}$

(3.20)

4. Relation Between Rotor Angular Orientation and its Spin Axis Vector, ω, in Torque Free Rotation

Having obtained the value of ω consistent with the value of ρ, we need to decide how to rotate the rotor. The rotor spins about the ω vector at the rate ω. In order to display this spin motion we need to use a matrix that rotates points about an axis parallel to ω in space. The Rodrigues formula http://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula is such a matrix and for our case of ρ and ω is written:

$ρ' = ρ c o s ϕ + (\hat{ω} \times ρ) \sin ϕ + \hat{ω} (\hat{ω} • ρ) (1 - \cos ϕ)$

(4.1)

where ω with the carat on top denotes a unit vector in the ω direction and ρ' is the new vector from the axis to the Rotor end sphere and $ϕ$ is the rotation angle (usually small) about the axis.