I am trying to grasp the concepts of forces applied to a rigid body resulting in a net change of the rotation of the object.
This is not a home assignment and I read some resources on both the web and YouTube but haven't found the answer to this.
The problem
I was thinking of how I can calculate the resulting rotational matrix $R$, given a timestep $dt$, the force $F$ and the vector $r$ with the inertia matrix $I$ for the object in question.
- $dt$: The timestep
- $F$: is a general 3-dimensional force.
- $r$: is the position where this force acts relative the object's center of rotation.
- $I$: is the 3×3 inertia matrix (which also is assumed to be invertible).
We have torque $\tau = F \times r$.
We also have $\tau = I \cdot \alpha$, where $\alpha$ is the angular acceleration.
To calculate a net rotation, assuming $I$ is invertible, we have:
equation(1): $\alpha = I^{-1} \cdot \tau = I^{-1} \cdot (F \times r)$
My initial thought was to integrate over a time step $dt$ and calculate the angular velocity:
$\alpha \cdot dt$
Then from this calculate a normalized vector around where the rotation occurs and resulting rotation angle around this axis. Let's call this vector $w$ and the angle $a$
After this calculate the matrix from:
$K = \begin{bmatrix}
0 & -w_z & w_y \\
w_z & 0 & -w_x \\
-w_y & w_x & 0
\end{bmatrix}$
We then can use Rodrigues' rotation formula ($I$ in the formula below is the identity matrix):
$R = I + sin(a)\cdot K + (1-cos(a))\cdot K^2$
Question:
Given the integrated angular velocity for a small timestep $\alpha \cdot dt$, where $\alpha$ comes from equation(1) above. How can I calculate the normalized vector $w$ and the angle $a$?
Alternatively can I calculate a generalized vector $w$ where the length $|w|$ is proportional to the angle $a$?
Best Answer
Nothing in Dynamics is simple
Your question can be split up into several steps
Given a known orientation ${\rm R}$ and angular velocity $\vec{\omega}$ of a rigid body and the net torque $\vec{\tau}_C$ summed about the center of mass, find the angular acceleration of the body $\vec{\alpha}$.
Usually, the mass moment of inertia of a body is given on body-fixed coordinates $$\overline{I}_{\rm body} = \pmatrix{I_1 & & \\ & I_2 & \\ & & I_3}$$
Transform the MMOI matrix into the inertial orientation, summed about the center of mass with $$ \overline{I}_C = {\rm R}\, \overline{I}_{\rm body} {\rm R}^\intercal$$
Angular momentum vector of the body summed about the center of mass is $$ \vec{L}_C = \overline{I}_C \vec{\omega}$$
Rotational equations of motion for the body are $$ \vec{\tau}_C = \overline{I}_C \vec{\alpha} + \vec{\omega} \times \vec{L}_C $$
Rotational acceleration of the body is $$ \vec{\alpha} = \overline{I}^{-1}_C \left( \vec{\tau}_C - \vec{\omega} \times \vec{L}_C \right)$$
Integrate angular acceleration $\vec{\alpha}$ and angular velocity $\vec{\omega}$ over a time step ${\rm d}t$ to find new values for the velocity and the rotation matrix ${\rm R}$.
$$ \vec{\omega} \leftarrow \vec{\omega} + \int \vec{\alpha}\,{\rm d}t$$
To solve this problem, people usually use quaternions $q$ which encode the same information as ${\rm R}$ but with 4 values instead of 9 in the rotation matrix. In fact any quaternion is a combination of a scalar value and a vector $$q = (q_s |\, \vec{q}_v)$$ These are operations you need to work with quaternions
Construct a quaternion $q$ from a rotation axis $\hat{z}$ and angle $\theta$ $$q = ( \cos \left( \tfrac{\theta}{2} \right) | \, \hat{z}\sin \left( \tfrac{\theta}{2} \right) )$$
A sequence of two rotations is represented by the quaternion product $$ a b = ( a_s b_s - \vec{a}_v \cdot \vec{b}_v | a_s \vec{b}_v + b_s \vec{a}_v + \vec{a}_v \times \vec{b}_v)$$
Where $\cdot$ and $\times$ are the vector dot and vector cross product respectively.
The rotation matrix of quaternion $q = (q_s |\, \vec{q}_v)$ is extracted by $$ {\rm R} = {\bf 1} + 2 q_s [\vec{q}_v\times] + 2 [\vec{q}_v\times] [\vec{q}_v\times]$$ where ${\bf 1}$ is the identity matrix, and $[\vec{q}_v\times]$ is the 3×3 skew-symmetrix matrix similar to the $K$ matrix above representing the cross product as matrix multiplication. The inverse rotation matrix is calculated similarly but with a negative sign in front of $q_s$.
The rate of change of quaternion based on a body rotational velocity of $\vec{\omega}$ is evaluated with the following quaternion product $$\begin{aligned}\tfrac{{\rm d}}{{\rm d}t}(q_{s}|\,\vec{q}_{v}) & =\tfrac{1}{2}(0|\,\vec{\omega})(q_{s}|\,\vec{q}_{v})\\ & =\tfrac{1}{2}(-\vec{\omega}\cdot\vec{q}_{v}|\,q_{s}\vec{\omega}+\vec{\omega}\times\vec{q}_{v}) \end{aligned} $$
The integration of the orientation of the body is done with $$q \leftarrow q + \int \tfrac{1}{2} ( 0 |\, \vec{\omega}) \,q\,{\rm d}t$$ which surprisingly does not diverge away from a rotation very quickly. In practice a renormalization process of $$ q \leftarrow ( \frac{q_s}{\sqrt{q_s^2 + \vec{q}_v \cdot \vec{q}_v}} | \, \frac{\vec{q}_v}{\sqrt{q_s^2 + \vec{q}_v \cdot \vec{q}_v}})$$ is needed every few dozen integration steps.
Since there isn't an analytical solution to the equations above, a numerical solver for ODEs is needed to carry out the time steps.
I wanted to state that in practice, the above process is never done in terms of linear and angular velocity. What is preferred is to keep track of linear $\vec{p}$ and angular $\vec{L}_C$ momentum and perform the integrations just as follows
$$ \vec{p} \leftarrow \vec{p} + \int \vec{F} {\rm d}t $$
$$ \vec{L}_C \leftarrow \vec{L}_C + \int \vec{\tau}_C {\rm d}t $$
where $\vec{F}$ and $\vec{\tau}_C$ are the applied forces and torques.
To extract the motion from momentum use
$$ \vec{v}_C = \tfrac{1}{m} \vec{p} $$ $$ \vec{\omega} = \overline{I}_C^{-1} \vec{L}_C $$
The above process is far simpler than trying to calculate $\vec{\alpha}$ and then integrating it.
Here is a summary of each time step in a rigid body simulation