Rigid Body Dynamics

The rigid body module provides comprehensive tools for 3D rotational dynamics including Euler angles, quaternions, and specialized models.

Overview

Rigid body mechanics describes rotational motion of extended objects. The module implements:

  • Euler Angles: ZYZ convention for orientation

  • Quaternions: Singularity-free rotation representation

  • Inertia Tensors: For various geometries

  • Euler’s Equations: Rotational equations of motion

  • Symmetric Top: Precession and nutation analysis

  • Gyroscope: Steady precession and applications

Theory

Euler Angles (ZYZ Convention)

Three successive rotations define orientation:

  1. Rotation by \(\phi\) about z-axis (precession)

  2. Rotation by \(\theta\) about new y-axis (nutation)

  3. Rotation by \(\psi\) about new z-axis (spin)

Warning

Euler angles have a singularity (gimbal lock) at \(\theta = 0, \pi\). Use quaternions for numerical integration.

Quaternions

A unit quaternion \(q = q_0 + q_1 i + q_2 j + q_3 k\) with \(|q| = 1\) represents rotation without singularities.

Kinematic equation:

\[\dot{q} = \frac{1}{2} \Omega(\omega) q\]

where \(\Omega(\omega)\) is a skew-symmetric matrix built from angular velocity.

Euler’s Equations

For principal axes with moments \(I_1, I_2, I_3\):

\[\begin{split}I_1 \dot{\omega}_1 &= (I_2 - I_3) \omega_2 \omega_3 + \tau_1 \\ I_2 \dot{\omega}_2 &= (I_3 - I_1) \omega_3 \omega_1 + \tau_2 \\ I_3 \dot{\omega}_3 &= (I_1 - I_2) \omega_1 \omega_2 + \tau_3\end{split}\]

Usage Examples

Rotation Matrices from Euler Angles

from mechanics_dsl.domains.classical import RigidBodyDynamics
import numpy as np

rigid = RigidBodyDynamics()

# Define Euler angles (radians)
phi = 0.5    # Precession
theta = 0.3  # Nutation
psi = 0.2    # Spin

R = rigid.euler_to_rotation_matrix(theta, phi, psi)

print("Rotation matrix:")
print(R)

Quaternion Operations

from mechanics_dsl.domains.classical import RigidBodyDynamics
import numpy as np

rigid = RigidBodyDynamics()

# Quaternion from rotation (angle, axis)
angle = np.pi / 4  # 45 degrees
axis = np.array([0, 0, 1])  # z-axis

# q = cos(θ/2) + sin(θ/2) * (axis)
q = rigid.axis_angle_to_quaternion(axis, angle)

# Convert to rotation matrix
R = rigid.quaternion_to_rotation_matrix(*q)

# Quaternion derivative from angular velocity
omega = np.array([0, 0, 1.0])  # rad/s about z
q_dot = rigid.quaternion_derivative(q, omega)

Computing Inertia Tensors

from mechanics_dsl.domains.classical import RigidBodyDynamics
from mechanics_dsl import InertiaTensor

rigid = RigidBodyDynamics()

# Solid sphere
I_sphere = rigid.compute_inertia_sphere(mass=1.0, radius=0.5)
# I = (2/5) * m * r²

# Solid cylinder (along z-axis)
I_cylinder = rigid.compute_inertia_cylinder(mass=2.0, radius=0.3, height=1.0)

# Rectangular box
I_box = rigid.compute_inertia_box(mass=5.0, length=1.0, width=0.5, height=0.2)

# Get principal moments
eigenvalues, eigenvectors = I_box.principal_axes()
print(f"Principal moments: {eigenvalues}")

Torque-Free Rotation (Euler’s Equations)

from mechanics_dsl.domains.classical import RigidBodyDynamics
from mechanics_dsl import InertiaTensor
import numpy as np
from scipy.integrate import solve_ivp

rigid = RigidBodyDynamics()

# Asymmetric top
I = InertiaTensor(I_xx=1.0, I_yy=2.0, I_zz=3.0)

# Initial angular velocity
omega0 = np.array([1.0, 0.1, 0.1])

def euler_eom(t, omega):
    return rigid.euler_equations_torque_free(omega, I)

# Integrate
sol = solve_ivp(euler_eom, [0, 10], omega0, dense_output=True)

# Angular momentum is conserved
L0 = rigid.angular_momentum(omega0, I)
L_final = rigid.angular_momentum(sol.y[:, -1], I)
print(f"L conserved: {np.allclose(L0, L_final)}")

Symmetric Top Dynamics

from mechanics_dsl.domains.classical import SymmetricTop
import numpy as np

# Create symmetric top (I1 = I2 ≠ I3)
top = SymmetricTop(I1=1.0, I3=0.5, mass=1.0, cm_height=0.3)

# Initial conditions
theta0 = 0.2      # Nutation angle
phi_dot0 = 0.5    # Precession rate
psi_dot0 = 10.0   # Spin rate

# Compute effective potential
V_eff = top.effective_potential(theta0, psi_dot0)

# Steady precession rate (for given spin)
omega_p = top.steady_precession_rate(psi_dot0, theta0)
print(f"Precession rate: {omega_p:.4f} rad/s")

Gyroscope Analysis

from mechanics_dsl.domains.classical import Gyroscope
import numpy as np

# Create gyroscope
gyro = Gyroscope(
    I_spin=0.01,       # Moment about spin axis
    I_transverse=0.02, # Moment about transverse axes
    mass=0.5,
    cm_distance=0.1    # Distance from pivot to CM
)

# Spin rate
omega_s = 100.0  # rad/s (fast spin)

# Compute precession rate
omega_p = gyro.precession_rate(omega_s, theta=np.pi/4)

# Nutation frequency (fast wobble)
omega_n = gyro.nutation_frequency(omega_s)

print(f"Precession: {omega_p:.4f} rad/s")
print(f"Nutation: {omega_n:.4f} rad/s")

API Reference

Classes

class RigidBodyDynamics

Core rigid body dynamics calculations.

euler_to_rotation_matrix(theta, phi, psi)

Convert ZYZ Euler angles to rotation matrix.

quaternion_to_rotation_matrix(q0, q1, q2, q3)

Convert quaternion to rotation matrix.

quaternion_derivative(q, omega)

Compute \(\dot{q}\) from angular velocity.

euler_equations_torque_free(omega, I)

Euler’s equations for torque-free rotation.

euler_equations_with_torque(omega, I, torque)

Euler’s equations with external torque.

angular_momentum(omega, I)

Compute \(\mathbf{L} = I \cdot \omega\).

rotational_kinetic_energy(omega, I)

Compute \(T = \frac{1}{2} \omega^T I \omega\).

class SymmetricTop(I1, I3, mass, cm_height)

Symmetric top with I₁ = I₂ ≠ I₃.

effective_potential(theta, psi_dot)

Compute effective potential for nutation.

steady_precession_rate(psi_dot, theta)

Compute rate for steady precession.

class Gyroscope(I_spin, I_transverse, mass, cm_distance)

Gyroscope model.

precession_rate(omega_spin, theta)

Compute precession angular velocity.

nutation_frequency(omega_spin)

Compute nutation oscillation frequency.

Data Classes

class InertiaTensor

Moment of inertia tensor.

Parameters:

  • I_zz (I_xx, I_yy,) – Diagonal elements

  • I_yz (I_xy, I_xz,) – Off-diagonal elements (default 0)

to_matrix()

Convert to 3x3 numpy array.

principal_axes()

Find principal moments and axes.

See Also