MatLab simulation of the intermediate axis theorem, also called the Dzhanibekov effect. I have explained the physics behind the simulation at the end of the page.
Rigid body free rotating dynamics simulation conditions:
- Cubic body (mimicing a CUBESAT). Body parameters found in
body_params.m - Free rotation → no external forces applied
- The body has initial angular velocity in the intermediate axis of inertia and a really small angular velocity in one of the other axis (hypotetical perturbation that triggers this strange behavior). Initial conditions are provided in
mian.m
NOTE: the code provided is compatible with the free rotation dynamics of any rigid body having provided its body parameters. Besides that, adding external forces can be easily done changing rot_dynamics.m.
To run the simulation, first it is necessary to derive the rotational dynamics equations. According to elementary mechanics:
To make things easier, we will assume the rotating body is a rigid body which has a main implication: the inertia matrix is constant in the body frame since there are not relative motion between the particles of the body.
The inertia matrix and the angular velocity in the body frame are expressed as follows.
Now, given that we are taking magnitudes in non-inertial frames, we will expand the right term of the 2nd Newton Law for rotating bodies in the body frame to take advantage of the rigid body conclusion.
As a result, the equation defining rotational dynamics becomes:
where the inertia matrix and the angular velocity are expressed in body frame coordinates. Besides that, since the simulation shows the free rotation of a rigid body, external forces momentum is zero.
Next, it is necessary to clear the angular velocity derivative terms in order to apply numerical methods for solving the system of equations. For simplicity, we will take the skew symetric operator to show the result in a matrix form.
Let's expand the whole system of equations to view the components of the final expression:
Now, we could solve the system and get the angular velocity in the body frame. However, to plot the "real" movement of the body we have to plot how it moves in the inertial frame. To do so, we just have to transform the angular velocity in the body frame into the angular velocity in the inertial frame. This can be accomplished thank to the following expressions with state the relationship between the angular velocities of the rotating frame and the inertial frame:
We can now clear the angular velocities in order to apply a numerical ODE solver
Finally, putting together both systems of equations we can write the whole system of equations that we have to solve as follows.
Which has been easily solved applying Runge-Kutta 4 technique
