Mechanics of Continua and Structures
Dynamic Euler Equation: Lie-Group Integration Scheme
The dydnamic Euler equation is \(\begin{align} \bs{J}\dot{\bs{\omega}} + \bs{\Omega}\bs{J}\omega - \bs{T} = 0, \end{align}\) where $\bs{\omega}$ and $\bs{J}$ are the rigid body angular velocity vector and inertia tensor with respect to the center of mass, expressed in the reference frame attached to the body; $\bs{T}$ is the external torque related to the body center of mass, expressed in the same reference frame; and \(\begin{align} \bs{\Omega} = \dot{\bs{R}}\bs{R}^{\text{T}} \end{align}\) is the angular velocity tensor, which is a skew symmetric tensor with an axial vector, i.e., the angular velocity vector $\bs{\omega}$.
We use the integration scheme presented by Terze et al. [1] to compute the instant rotation tensor $\bs{R}$ and angular velocity vector $\bs{\omega}$ at time $t$ with the initial conditions $\bs{R}(0)=\bs{R}_0$ and $\bs{\omega}(0)=\bs{\omega}_0$.
Suppose we know the rotation tensor $\bs{R}$ and angular velocity vector $\bs{\omega}$ at $n_{th}$ time step, denoted as $\bs{R}_n$ and $\bs{\omega}_n$. Then we can update the rotation tensor and velocity vector as follows:
\(\begin{align} \bs{Y}_{n+\frac{1}{2}} &= \exp\left( -\frac{h}{2} \bs{\Omega}_n \right) \left( \bs{J} \bs{\omega}_n + \frac{h}{2}\bs{T}_n \right)\\ \bs{\omega}_{n+\frac{1}{2}} &= \bs{J}^{-1} \bs{Y}_{n+\frac{1}{2}} \\ \bs{R}_{n+1} &= \bs{R}_n \exp\left(h \bs{\Omega}_{n+\frac{1}{2}}\right) \\ \bs{Y}_{n+1} &= \exp\left( -\frac{h}{2} \bs{\Omega}_{n+\frac{1}{2}} \right) \left[\exp\left( -\frac{h}{2} \bs{\Omega}_{n+\frac{1}{2}} \right) \left( \bs{J} \bs{\omega}_n + \frac{h}{2}\bs{T}_n \right) + \frac{h}{2} \exp\left( \frac{h}{2} \bs{\Omega}_{n+\frac{1}{2}} \right) \bs{T}_{n+1} \right]\\ \bs{\omega}_{n+1} &= \bs{J}^{-1} \bs{Y}_{n+1} \end{align}\) If the body in motion is free of external torque, then $\bs{T} = 0$. The $\bs{\Omega}_n$ is the skew symmetric tensor with its axial vector given by $\bs{\omega}_n$. The $\exp(\cdot)$ is the exponential map of a skew-symmetric tensor [2]. For example, given a skew symmetric tensor $\bs{B}$ with matrix representation \(\begin{align} [\bs{B}] = \begin{bmatrix} 0 & -c & b \\ c & 0 & -a \\ -b & a & 0 \end{bmatrix}, \end{align}\) we have \(\begin{align} \exp(\bs{B}) = \bs{I} + \frac{\sin\alpha}{\alpha} \bs{B} + \frac{1-\cos\alpha}{\alpha^2}\bs{B}^2, \end{align}\) where $\alpha = \sqrt{a^2+b^2+c^2}$.
Numerical Example
-
A Tossed Book in the air
Mathematica Code is available here
-
A Stumbling T-handle in Space
-
Anti-Gravity Wheel (Gyroscopic Precession)
References:
[1]. Zdravko Terze, Andreas Müller and Dario Zlatar, An Angular Momentum and Energy Conserving Lie-Group Integration Scheme for Rigid Body Rotational Dynamics Originating From Stormer–Verlet Algorithm, J. Comput. Nonlinear Dynam 10(5):051005, 2015.
[2]. Jean Gallier and Dianna Xu, Computing Exponentials of Skew Symmetric Matrices And Logarithms of Orthogonal Matrices, International Journal of Robotics and Automation, 17(4):1-11, 2002.