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 mdv where m is the mass and dv 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, dv_x, and a final rotation, at rate dw, about the center of mass.    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 dw_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 r 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.

(1.1)

Now 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:

(1.2)

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

                                                                       

(1.3)

Then the matrix and vector expression becomes:

(1.4)

(1.5)

Solving equation (1.5) we obtain:

             (1.6)

We need to insert the values of dw into an energy conservation equation.

(1.7)

and then solve the resulting quadratic equation for dv_c.

 

(1.8)

                                                                       

                                                         (1.9)

The solutions for dv and dw are then:

                                                                       

(1.10)

(1.11)

(1.12)

The new spin rate vector is w+dw.

 

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

(1.13)

where

(1.14)

and dt is the time increment.

In order to keep the body axes of the rotor updated, we need to multiply both the rotor end sphere coordinates and its R and Y body axes by M_R.

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.

            (2.1)

where -  is equal to the c component of the velocity due to the rotational motion of the rotor. 

 

Also note that, by definition, r and w are always perpendicular,

(2.2)

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

(2.3)

All of these equations are valid regardless of the coordinate system in which c, w and r are described.  We simply choose to work in the xyz system and to convert 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:

 

                                                                       

                                                                       

(2.4)

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

(2.5)

 

                                                                                                                                               

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

                                                                       

(2.6)

where

(2.7)

and the new rotation axis unit vector is:

(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:

(2.9)

where w_R' and  w_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

(2.10)

Then the new value of the CM speed is

(2.11)

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

 

(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_w, 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:

 

(2.14)

while the angular momentum of the rotor about its axis is:

(2.15)

and the linear momentum of the 2 masses is:

(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 dv, then, from equation (2.16), the Sphere speed change will be -(m_r/m_s)dv.  The expression for the conservation of angular momentum of is a bit more elaborate.  The equation for the final angular momentum is:

(2.17)

while the initial angular momentum is:

(2.18)

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

 

(2.19)

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

(2.20)

(2.21)

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

(2.22)

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

(2.23)

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

             (2.24)

(2.25)

        (2.26)

The result for dv is:

(2.27)

Then dw becomes

(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.  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 one rotor's end Sphere to the center of the one of the end spheres of the other rotor.  

 

 

            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_w, 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:

 

(3.1)

The final angular momentum after the collision can be written:

(3.2)

and the initial angular momentum before the collision can be written

(3.3)

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

 

(3.4)

First let's rewrite  in terms of r₁ and r₂.  Note that the vector expression relating  is

(3.5)

Then solving for  we have

(3.6)

so we can rewrite equation (3.4)

(3.7)

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

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

(3.8)

Then re-name dw₁ to dw and compute P:

(3.9)

And then we might use our result for P in the dw₂ equation:

(3.10)

The previous equations would be fine except that we can't divide a vector into a scalar or another vector. 

Instead we define the following 3x3 matrices:

(3.11)

                                                                       

Then we can write the matrix analog of equation (3.10):

(3.12)

Equation (3.12) can be solved for dw₂ in terms of dw.  Then equation (3.7) can be solved as:

(3.13)

so that the solution for dw is:

(3.14)

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

 

The determinant of the matrix to be inverted in equation (3.12) is

(3.15)

Then naming the coefficients of the matrix product in equation (3.12)

(3.16)

The values of the elements of M₂₁ are:

(3.17)

(3.18)

(3.19)

(3.20)

(3.21)

(3.22)

(3.23)

(3.24)

(3.25)

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

(3.26)

Simplify notation for equation (3.14)

(3.27)

(3.28)

Recalling equation (3.14) for dw we have:

(3.29)

Relation Between Rotor Angular Orientation and its Spin Axis Vector, w, in Free Rotation

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

(4.1)

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