Tumbling Through Space: A Dive into Rotational Dynamic Simulations

Have you ever spun your phone in the air? If you have you probably noticed that when you tossed it in a specific orientation it spun better than others. In this article we will go into why that is and how you can predict the motion of irregular objects when they rotate.

When someone says the word tumbling, what do you think of? I think of the people chasing a wheel of cheese down a hill but we can talk about that in a different article. Tumbling is a complex, sometimes chaotic motion that is a combination of external forces and the rotation of an object, but what if the object in question if floating through space? This is a very important concept if you are designing a satellite. You want to ensure that you can point it in the desired direction for communication or science payloads. And these dynamics are important for designing a control system.

While we aren’t going to go that far in this article, I’d like to introduce the topic with a classic example that is best explained in zero-gravity, The Rotation of a T-handle. First, we will go through the math that describes the rotation of an object. Then we will review the results of a simulation for a solid T-shaped object. Finally, we will walk through the entire process of creating our object in CAD software, writing the code for our simulation, and exporting the simulation data into animation software to create videos.

The Tennis Racket Theorem

This is also called the intermediate axis theorem. It describes the motion of a rigid body with 3 principle moments of inertia. There are two axes that have stable rotations, these happen to be the axes with the largest and smallest moments of inertia. If a rotation is applied around the axis with the moment of inertia that is between the other two it exhibits unique oscillatory behavior. This effect can be easily observed in objects such as wingnuts, T-handles, tennis rackets, cereal boxes, and cellphones.

We can see the strange behaviors of these spinning objects in videos taken in microgravity. Here are a few good examples.

The Rotation Equation

Euler’s Rotation Equation describes how the inertia of an object is related to it’s angular velocity and angular momentum. The angular momentum of an axis is the product of the angular velocity and inertia. This set of equations can be used to determine how much torque you will need to produce an angular acceleration.
Euler’s Rotation Equation is:
$$M=I\dot{\omega}+\omega \times (I\omega)$$ Where \(M\) is the angular momentum, \(\omega\) is the angular velocity, and \(I\) is the moment of inertia.
It can be rewritten as a series of differential equations that relate the angular momentum around each of the 3 principal axes with the moment of inertia and angular velocity of each of the axes of rotation.
$$
\begin{matrix}
I_1\dot{\omega_1}+(I_3-I_2)\omega_2\omega_3=M_1 \\
I_2\dot{\omega_1}+(I_1-I_3)\omega_3\omega_1=M_2 \\
I_3\dot{\omega_1}+(I_2-I_1)\omega_1\omega_2=M_3 \\
\end{matrix}$$

Simulating the T-Handle Dynamics

The T-Handle is a shape that has some interesting dynamic responses to intermediate axis rotations that we can simulate and visualize. We can use a model to easily represent the differential equations that describe the motion of our object. For our simulation model we will use an Inertia tensor computed for the dimensions of our object. We will also simulate the handle in zero gravity as this is where the dynamic behavior is most evident. We can also compare our simulation to actual footage of these rotational dynamics in a microgravity environment.

The first thing that we need to do is use Euler’s equations for rigid body dynamics with zero external torque
$$
\begin{matrix}
I_1\dot{\omega_1}+(I_3-I_2)\omega_2\omega_3=0 \\
I_2\dot{\omega_2}+(I_1-I_3)\omega_3\omega_1=0 \\
I_3\dot{\omega_3}+(I_2-I_1)\omega_1\omega_2=0 \\
\end{matrix}$$
These can be re-written to isolate the terms that represent the changes in angular velocity \(\omega\) around each principal axis.
$$
\begin{matrix}
\dot{\omega_1}=\frac{(I_2-I_3)}{I_1}\omega_2\omega_3 \\
\dot{\omega_2}=\frac{(I_3-I_1)}{I_2}\omega_3\omega_1 \\
\dot{\omega_3}=\frac{(I_1-I_2)}{I_3}\omega_1\omega_2 \\
\end{matrix}$$
These are the differential equations that describe the changes in our state variables, the angular velocity, of each of the 3 axes of rotation.

We can then use the numeric integrator of our choice to propagate the states through each timestep in our simulation.

For this example our T-Handle will have the following moment of inertia tensor
$$\begin{bmatrix}
.299 & 0 & 0 \\
0 & .408 & 0\\
0 & 0 & .116
\end{bmatrix}kg/ cm^2$$
We can simulate the motion of the handle for a rotation around each principle axis by using the initial conditions of 4 rad/s around the principle axis that we would like to investigate, and 0.1 rad/s around the other two.

Axis with the Smallest Moment of Inertia

Starting with the axis with the smallest moment of inertia (the Z axis) we see that the rotation is stable. Looking at the plot of angular velocities we see that the z axis rotates while the x and y axes stay mostly the same. This is confirmed in the animation generated from this simulation.

Z-Axis Rotation

Axis with the Largest Moment of Inertia

Next we look at the axis with the largest moment of inertia (the Y axis). From the graph we can see that the rotation is also stable with the x and z axes oscillating around our initial perturbations. We can visually confirm this behavior by observing the animation of this rotation.

Y-Axis Rotation

Axis with the Intermediate Moment of Inertia

Finally, for the axis with the intermediate moment of inertia we see a very strange behavior. By looking at the graph of angular velocity there seems to be some tumbling behavior where the x-axis starts out rotating the fastest, then it slows, changes direction, and continues this oscillatory behavior. We see that the y axis always has a positive rotation and exhibits peak oscillations in sync with the peak oscillations of the z axis. The Z axis alternates between forward and backward rotations in line with the peaks on the y-axis. When observing the animation for this behvior we can see that our handle flips back and forth in the same way that we saw from the videos at the beginning of this article.

X-Axis Rotation

Using CAD software to find the Moment of Inertia

An easy way to determine the moment of inertia for an irregular object is to use commercial software CAD tools. We will create a model of our T-handle using two cylinders and then have the Autodesk Inventor CAD software calculate the moments of inertia that we can then use in a simulation. First we will open up a new part in Inventor. You should see the following screen.

We need to change the default units to mm and kg so we will click on the tool bar in the upper menu and click document settings.

This will bring up the settings for the document. We want to click on the units tab and change the length unit to mm and the mass unit to kg. Hit OK.

We will start by making a 50mm x 10mm cylinder in the x direction. To do this we will click Start 2D Sketch and then select the yz plane.

Then we make a 10mm circle centered on the origin.

Click Finish Sketch. Now we will extrude the surface by selecting the extrude command from the 3D model tab, selecting our sketch, changing the distance to 50mm, and clicking Ok.

We now need to make the other part of the T-handle. To do this we create a new sketch and select the XY Plane from the origin folder.

We can then draw a second 10mm circle centered on the origin and click finish sketch.

As with our previous extrusion we select our sketch and this time we extrude 100mm and select symmetric. This means that the total extrusion will be 50mm on each side of the sketch plane.

Now we have the completed model for our T-handle. Finally we will need to select a material, build the moment of inertia tensor, and export the model for future animations.
First click File -> iProperties to open up the properties of the part. In the physical tab Change the material to Titanium, And select the Inertial properties from the center of gravity. This window shows us the lower-triangular part of the moment of inertia tensor.

The moment of inertia tensor is symmetrical around the principle axes so our final inertia tensor is
$$\begin{bmatrix}
.299 & 0 & 0 \\
0 & .408 & 0\\
0 & 0 & .116
\end{bmatrix}kg/cm^2$$

Finally, we can export or model as an STL file by selecting File -> Export -> CAD Format, and then saving it in a location of your choosing.

Coding a Simulation

Once again, we will be using the Python language to write our simulation. This is beacauseP Python has easy syntax and many useful libraries. We will start with the imports. The operator libary will be used for mapping functions onto lists, the matplotlib library will be used for generating and saving graphs, and the numpy library will be used for list operatons. The csv and datetime libraries will be used for saving the simulation into timestamped output files.

Our simulation will use 4 input variables, the inertia tensor, the initial angular velocity, the integration timestep, and the total time of the simulation. These are passed into the main simulation function called rotation_sim.

In the first part of our simulation method we use the datetime library to get the current timestamp. This must be formatted into a string that can be used as a valid filename, therefore we use the strftime() command to define the format as Ymd_HMS, where Y is the year, m is the moth, d is the day, H is the hour, M is the minute, and S is the second.

Next we save the inputs to the simulation using the save_config() function. This function opens a new text file with the timestamp as the prefix of the filename. It then writes the Inertia tensor, angular velocity, timestep, and end time to the config file on separate lines, and closes the file.

Going back to the rotation_sim function we need to set up our data structures to store the outputs of our simulation. These outputs are the angular velocities and quaternion rotation components.

Next we initialize our rotation. Quaternions are a 4 parameter set that describes a rotation around an euler axis. They represent the vector equation
$$q = q_0+iq_1+jq_2+kq_3$$
where \(i\), \(j\), and \(k\), are scaled unit vectors in a 3D cartesian coordinate system. They are more computationally efficient and numerically stable compared to their Euler Angle counterparts as they don’t exhibit the gimbal lock phenomena when two or more of the axes align with each other.

Because the magnitude of a quaternion is always 1, we need to initialize our attitude to the components (1,0,0,0). After that, we set up our simulation time using the np.arrange() method. This method gives a linear mapping of values from a start to an end with a specified timestep. We will use 0 as our start time, end as our end time, and t as our timestep.

Moving along we get to our main simulation loop. Our state variable is going to be a 3D vector ofthe angular velocity around each axis. We then loop through all of the elements in our time series T, integrate our equations of motion, use the resulting angular velocity to update the orientation quaternion, and then save the angular velocity, and quaternions into our output structures.

Looking at the first part of the loop, we can see how the state variable is updated with the current state using the forward Euler method. This method is simple to implement and understand with the downside of being inaccurate and possibly unusable for specific classes of problems.

Looking at the forward_euler method we see that it takes a function \(f\), a state variable \(x\), a timestep \(dt\) and an inertia tensor \(I\). The inertia tensor is needed because we pass it into our differential equation for rotation along with our state variable. We have two temporary variables \(a\) and \(b\). \(a\) stores the value of our differential function evaluated at the current state \(x\) using the inertia tensor of our object. This is then multiplied elementwise by the timestep to obtain a vector for the derivative of our angular velocity. This is then added to the existing state vector to get our new state vector at the new timestep.

Next we can look at the function t_handle_derivs() for the implementation of Euler’s Rotation Equations. As discussed above, these are the re-written equations for the angular momentum. The only inputs that we have are the current angular velocity and the inertial tensor of our object, and our return value is the vector of the derivative of the angular velocity, which would be acceleration.

Going back to our simulation loop we save the outputs of our new state and then go on to calculating rotation.

The body angular rates can be described by the following equation which relates the body angular velocities to quaternions using the following matrix:$$\begin{bmatrix}
\dot{q_1} \\
\dot{q_2} \\
\dot{q_3} \\
\dot{q_4}
\end{bmatrix}=
\frac{1}{2}
\begin{bmatrix}
-q_2 & -q_3 & -q_4 \\
q_2 & -q_4 & q_3 \\
q_4 & q_1 & -q_2 \\
-q_3 & q_2 & q_1
\end{bmatrix}
\begin{bmatrix}
\omega_1 \\ \omega_2 \\ \omega_3
\end{bmatrix}$$

Back to our main simulation loop, we save the rotations at each timestep in our output structures.

Next we can look at the plot_results() function. This function plots all 3 angular velocity components on a graph, adds a title and a legend to the graph, saves the graph as a timestamped jpg file and displays the graph as a new figure.

Next we save the angular velocities in a csv file using the save_csv() function. This function takes our angular velocity components and a filename, creates a new file using that filename, and saves the angular velocity components at each timestep. This is where the writerow() method of the csv library comes in handy. Once we are done writing the csv data we close the file.

Back at our main simulation we then plot the rotation quaternion for each timestep using the plot_quaternions() method. It’s very difficult to imagine exactly what a rotation would look like given these time series rotations but it is a good test to make sure that you are getting the output that you are expecting. This function also saves the graph as a jpg file.

The final step in our simulation is to save the actual rotations in a csv file so that we can export them to an animation program, or simply save the data for later analysis. This is done with the save_quat_csv() function. Just like the function for saving the angular velocities, it writes each of the vector components to a column and saves the vector at each timestep. It’s always good practice to close an external file when you’re done writing to it.

Using Blender to Animate our Simulation

Blender is a powerful open-source 3D animation software tool that we will use for visualizing our simulation data in a way that is more intuitive than time-series graphs. The reason that I am using it is that It uses a custom API for python that can be used to write automation scripts.

We will now animate our simulation. We have the 3d model that we created as an stl file, and we have a set of quaternion rotations from the simulation. We are going to import both of them into the animation software Blender to create a high-quality physics animation.

When we first open up blender we get a window that looks like this, it has a defulat camera, light source, and a cube object. We can right-click and delete the cube as we will replace it with the STL file that we saved earlier.

Next Click File -> Import -> (.stl) and navigate to the directory where you save the .stl file.

Cllick Import STL. In the scene collections menu at the top right, rename the object to “T-Handle”

We then need to scale the model down to fit in the camera field of view. We can do this by clicking ‘s’ and then moving the mouse to scale down the model. Once we have scaled it to an appropriate size we are all set to import the simulation data.

Click on scripting to open the scripting perspective.

Click on the New button to start a new script. We then need to import the blender API and the csv module for reading the simulation data.

Next we select our handle object inside the context of the blender scene

Just to make sure that we selected the object properly we can add a translation line to move the object on the positive Z axis.

When we run this script we can see that our handle moves in the z direction and we can also see that the z coordinate of our object has changed relative to the origin.

Now that we know we can manipulate the object, The next thing that we need to do is load the csv file that contains the quaternion rotations of the handle from our simulation. We can use the following python function to load the rotations into 4 separate arrays.

Our simulation has a time resolution of 1ms, but our animation is going to be running at 24 fps. Because our simulation timesteps are 1ms we need to take the value at every 41 timesteps and add that into our animation as a keyframe. We can do this in our animation function.

Now that we have arrays with each of our quaternion components we can add the rotations and save the keyframes for every frame of our animation.

We then call these two functions to import our simulation data.

Finally, before running the script we need to change the rotation mode to quaternion from the default Euler.

Click the Play button in the menu to run your script. We can then switch the animation mode by selecting it from the top menu.

In the animation mode we select the T-Handle object and we can see all of our keyframes in the timeline at the bottom. By pressing play we can see the result of our animation.

Now we can render our animation into a movie file. In the scene properties window, set the output file format to FFmpeg, and the encoding container to MPEG-4. This tells Blender to output our animation as an .mp4 movie file in the tmp directory. In the top menu bar click Render Animation under the Render tab to render the animation.

Once the rendering has completed, we can navigate to the tmp directory and view our newly generated movie.

By changing the imported file, we can generate animations for rotations around each of the 3 principle axes. This method was used to generate the animations that we looked at at the beginning of the article.

Conclusion

This was a bit of a long post but I hope that you stuck through till the end! The code for this project and the generated graphs and simulation data are available on GitHub if you have the urge to experiment with it and I hope that this inspires someone to make their own simulations.

Sources

Class 6 – Quadrotor Dynamics, 2019. https://www.youtube.com/watch?v=UC8W3SfKGmg.

COMSOL. “Why Do Tennis Rackets Tumble? The Dzhanibekov Effect Explained….” Accessed January 13, 2023. https://www.comsol.com/blogs/why-do-tennis-rackets-tumble-the-dzhanibekov-effect-explained/.

“Gyro Noise and Allan Deviation + IMU Example,” May 9, 2021. https://mwrona.com/posts/gyro-noise-analysis/.

“Intermediate Axis Spin Simulation (Tennis Racket / Dzhanibekov Effect) – YouTube.” Accessed October 19, 2022. https://www.youtube.com/watch?v=21scyr5_Spg.

“Quaternions.” Accessed February 9, 2023. https://danceswithcode.net/engineeringnotes/quaternions/quaternions.html.

ThatsMaths. “The Intermediate Axis Theorem,” December 12, 2019. https://thatsmaths.com/2019/12/12/the-intermediate-axis-theorem/.

Keep In Touch

Twitter Linkedin-in