There have been some interesting things on the interwebs, lately, about rotation, gyroscopes, precession, and helicoptors (all of it brought to my attention, or literally done by, Derek Muller of Veritasium fame). It got me thinking about the modeling I’ve done with my students, using Mathematica to simulate various rotating objects.
One of the ideas I wanted to play with is described very nicely in this arXiv paper entitled “‘It has to do down a little, in order to go around’ – Following Feynman on the gyroscope.” The basic idea is that pure precession doesn’t seem to make sense because the torque is entirely in the horizontal plane (what causes the vertical angular momentum?). Feynman suggests that the axis of the gyroscope has to dip a little first to get the precession going. The arXiv paper describes this nicely and shows a beautiful experiment where they measure it.
What I wanted to do was to see if I could look at how we model a gyroscope and see what I had to turn on to get the thing to precess. I took two very different approaches to the modeling (more on the comparison at the bottom of this post). Here I’ll just describe the second approach. I used Euler angles to describe the orientation of the gyroscope, and used a Lagrangian approach to find the equations of motion for the three Euler angles. Then I numerically integrated those to make pretty pictures.
In the Euler approach, one angle (I called it phi) twists the x- and y-axes around the z-axis so that when I bent the gyroscope down from it’s original vertical alignment, it would bend in the right direction. The angle theta then describes how much it bends from the vertical (it bends towards where the new x-axis is). Then I let psi describe the rotation of the gyroscope around its own axis.
To determine the kinetic energy of the gyroscope, it’s much easier in the body axes. In those axes:
where the ‘s are the rotation rates around the body axes. This looks nice and easy, especially since, in the body axes, the moments of inertia (the I_i’s) are constant. However, the ‘s are really hard to calculate as functions of the Euler angles and their time derivatives. However, that’s still easier that recalculating the moment of inertia tensor (all nine elements) at every time step and using the tensor equivalent of the equation above. Here’s what the whole thing looks like in Mathematica:
The first three lines are the difficult-to-derive relationships among the omegas and the Euler angles (and their time derivatives). After that I set the principle moments of inertia and the distance between the pivot point and the center of mass (h). Then I set the kinetic and potential energies, and set up the Lagrangian approach to the equations of motion.
For this particular run, I wanted to see what would happen if I allowed phi to vary (because that would be the precession angle) but I fixed theta to be 90 degrees. That’s the typical setup for the bicycle wheel gyro demo. The result was a wheel that spun in a single plane without precession!
If I allow theta to also vary, the Mathematica looks like this:
When I do that, I get the precession, nutation, and rotation we expect:
Here’s a plot of theta as a function of time, showing the nutation:
Note that, even though it starts at 90 degrees, it really oscillates at angles that take it lower than that. In face, that overall drop of the axis of the gyroscope is what causes the precession!
The first approach I took was to simply model a small number (~20) of points and calculate the kinetic and potential energy of all of them, but enforce a relationship among them so that they would be part of a rigid structure. To do this, I had to think carefully about their x-, y-, and z-coordinates (in the lab frame) as functions of the Euler angles. This wasn’t that hard, just an exercise in geometry. What I like about this method is that I think it’s easier to teach. That geometry exercise is something that typical students could do with some prodding, while the omega calculation above I’ve really struggled to teach. I get the same results with both approaches, but the code is much cleaner for the approach above. Your thoughts?