Boxes: a benchmark for free-floating rigid bodies

This test simulates free-floating rigid bodies with and without gravity and verifies that their response matches the predicted analytical solutions.

The relevant physical principles are momentum conservation and the rotational stability of rigid bodies, sometimes referred to as the tennis racket theorem or the Dzhanibekov effect. Depending on the shape of an object, it may have unstable rotations about a particular axis, which causes tumbling (excellent video here).

Two scenarios are chosen in this benchmark to represent simple and complex trajectories. The "simple" scenario has no gravity and stable rotations, and the "complex" scenario has gravity and unstable rotations.

Simple scenario

Without gravity, the linear momentum of the rigid body remains constant. The expected trajectory of the center of mass is to follow a straight line with constant velocity. Additionally, the initial angular velocity is aligned with a stable axis of rotation. In addition to constant angular momentum, the angular velocity should be constant, and the cumulative rotation angle will vary linearly with time.

Complex scenario

With gravity, the linear momentum changes linearly with time. The expected trajectory of the center of mass is parabolic. Additionally, the initial angular velocity is slightly mis-aligned from the unstable axis of rotation. Since no external torques act, the angular momentum will remain constant, though that manifests in a tumbling angular velocity.

Definition of coordinate frames and variables

Consider an inertial frame $O$ and a rigid body with coordinate frame $c$ attached at the center of mass:

• The position of the center of gravity (cg) in $O$ is given by $\textbf{c}$.
• The orientation of $c$ with respect to $O$ is given by the quaternion $\textbf{q}$.
• A rotation matrix from $O$ to $c$: $\textbf{R}(\textbf{q})$
• Angular velocity in frame $c$: $\boldsymbol{\omega}$

The rigid body has the following inertial parameters:

• mass $m$
• inertia matrix $\textbf{I}$ expressed in $c$

With uniform gravity vector $\textbf{g}$, momentum and energy of the rigid body are defined as follows:

• Linear momentum in frame $O$: $\textbf{p}(t) = m \dot{\textbf{c}}(t)$
• Angular momentum with respect to $c$, expressed in $O$: $\textbf{L}(t) = \textbf{R}^T(\textbf{q}(t))\textbf{I}\boldsymbol{\omega}(t)$
• Kinetic energy: $T(t) = \frac{1}{2} m \dot{\textbf{c}}^T(t) \dot{\textbf{c}}(t) + \frac{1}{2} \boldsymbol{\omega}^T(t) \textbf{I} \boldsymbol{\omega}(t)$
• Potential energy: $V(t) = -m \textbf{g}^T \textbf{c}(t)$
• Total energy: $E(t) = T(t) + V(t)$

Initial conditions

For time $t$, denote the initial conditions:

• $\textbf{c}(t=t_0) = \textbf{c}_0$,
• $\dot{\textbf{c}}(t=t_0) = \dot{\textbf{c}}_0$,
• $\textbf{q}(t=t_0) = \textbf{q}_0$,
• $\boldsymbol{\omega}(t=t_0) = \boldsymbol{\omega}_0$,
• $\textbf{p}(t=t_0) = \textbf{p}_0$,
• $\textbf{L}(t=t_0) = \textbf{L}_0$,
• $T(t=t_0) = T_0$,
• $V(t=t_0) = V_0$,
• $E(t=t_0) = E_0$.

Inertial parameters for box

These tests use a box of uniform density with relative dimensions 1x4x9 in local coordinate frame $xyz$.

The inertial properties are computed below:

The box's $x$ axis has the largest inertia, and the $z$ axis has the smallest inertia. The $y$ axis is thus the unstable rotational axis.

Analytical solution

With gravity as the only external force, Newton 2nd law implies that the time derivative of linear momentum is equal to the gravitational force:

• $\dot{\textbf{p}}(t) = m \textbf{g}$

The center of mass acceleration is constant

• $\ddot{\textbf{c}}(t) = \textbf{g}$,

the velocity is linear in time

• $\dot{\textbf{c}}(t) = \dot{\textbf{c}}_0 + \textbf{g} (t-t_0)$,

and the position is parabolic in time

• $\textbf{c}(t) = \textbf{c}_0 + \dot{\textbf{c}}_0(t-t_0) + \frac{1}{2} \textbf{g} (t-t_0)^2$.

With no external torques acting, the angular momentum is conserved:

• $\textbf{L}(t) = \textbf{L}_0$,
• $\textbf{R}^T(\textbf{q(t)})\textbf{I}\boldsymbol{\omega}(t) = \textbf{R}^T(\textbf{q}_0)\textbf{I}\boldsymbol{\omega}_0$

This implies that the angular momentum vector will be constant in a fixed frame.

Additionally, gravity is a conservative force, so the system is conservative and energy is constant.

• $E(t) = E_0$

Analytical solution: simple scenario

There are some simplifications of the analytical solution for the simple scenario. Since gravity is not acting ($\textbf{g} = \textbf{0}$), the linear velocity reduces to a constant and the position to a linear function of time. The potential energy is zero, so energy conservation implies that the kinetic energy is constant:

• $T(t) = T_0$

The specific initial conditions are given as follows:

• $\textbf{c}_0 = [0, 0, 0]^T$
• $\textbf{q}_0$ is identity
• $\dot{\textbf{c}}_0 = [-0.9, 0.4, 0.1]^T$
• $\boldsymbol{\omega}_0 = [0.5, 0, 0]^T$

Note that the angular velocity is aligned with the $x$ axis, which has the largest moment of inertia and thus rotational stability.

The initial angular momentum and energy are computed as:

The scenario is simulated over of a time interval of 10 seconds.

Analytical solution: complex scenario

The gravity and specific initial conditions for the complex scenario are given as follows:

• $\textbf{g} = [0, 0, -9.8]^T$
• $\textbf{c}_0 = [0, 0, 0]^T$
• $\textbf{q}_0$ is identity
• $\dot{\textbf{c}}_0 = [-2.0, 2.0, 8.0]^T$
• $\boldsymbol{\omega}_0 = [0.1, 5.0, 0.1]^T$

Note that the angular velocity is slightly mis-aligned from the unstable $y$ axis, which has neither the largest nor smallest moment of inertia value and thus rotational instability.