# Solving Rigid Body ODEs in rotating frame

I'm trying to solve numerically the motion of a rigid body in a fluid. To this end I have to solve the fluid PDEs, and at every time step the 6 ODEs of the rigid motion degrees of freedom:

$$(m\mathcal{I} + \mathcal{A}) \frac{d\mathbf{U}}{dt} + \boldsymbol{\Omega} \times ((m\mathcal{I}+\mathcal{A}) \mathbf{U}) = m_1 \mathbf{g} + \mathbf{F}\\ (\mathcal{J}+\mathcal{D}) \frac{d\Omega}{dt} + \boldsymbol{\Omega} \times ((\mathcal{J}+\mathcal{D})\boldsymbol{\Omega}) + \mathbf{U}\times(\mathcal{A}\mathbf{U}) = \mathbf{M}$$

where $$\mathcal{A}, \mathcal{D}$$ model the added mass effects of the body.

These ODEs are expressed in a system of axes that rotates over time the same way as the body axes, but has a fixed origin. Therefore I need the components of the gravity vector $$\mathbf{g}$$ in this relative system. My first attempt is to say that since $$\mathbf{0} = \frac{d}{dt}\mathbf{g} = \left( \frac{d \mathbf{g}}{dt}\right)_r + \boldsymbol{\Omega}\times \mathbf{g}$$ is the relation between the time derivatives in the fixed and rotating frame, I can simply add to the 6 ODEs of the body above the other 3 ODEs $$\dot{\mathbf{g}}_r = - \boldsymbol{\Omega}\times \mathbf{g}.$$

Therefore, every time step, I update the values of the force $$\mathbf{F}$$ and torque $$\mathbf{M}$$ and advance the 9 ODE system with a Runge Kutta method, and I get the $$\mathbf{g}$$ components in the rotated reference for free with the new $$\mathbf{\Omega}$$ and $$\mathbf{U}$$.

Is it correct though? What is the relation that gives the transformation of the components of $$\mathbf{g}$$ from the laboratory reference to the rotated one at every time step? What are the angles that I'm using here? Because I know for example that Euler angles are generally used in these situations, but I haven't introduced any convention here, and I'm not worrying about the order in which I do the three successive rotations.

Finally, is it better to use quaternions in such problems, to avoid singularities?

• The time derivative of g is wrong, the components of g must be transferred with the rotation matrix R, additionally you get out of the rotation matrix equation where $\Omega=f(\varphi,\dot \varphi)$ where $\varphi$ are 3 Euler angle
– Eli
Commented May 9, 2020 at 6:38
• @Eli thanks, do you have a reference or know of any examples I could use to understand better the steps I have to do?
– gc11
Commented May 9, 2020 at 10:06
– Eli
Commented May 9, 2020 at 13:19

How to simulate

$$(m\mathcal{I} + \mathcal{A}) \frac{d\mathbf{U}}{dt} + \boldsymbol{\Omega} \times ((m\mathcal{I}+\mathcal{A}) \mathbf{U}) = m_1 \mathbf{g} + \mathbf{F}\\ (\mathcal{J}+\mathcal{D}) \frac{d\Omega}{dt} + \boldsymbol{\Omega} \times ((\mathcal{J}+\mathcal{D})\boldsymbol{\Omega}) + \mathbf{U}\times(\mathcal{A}\mathbf{U}) = \mathbf{M}$$

we want to simulate those equations in body fixed system (B-System), thus all vectors components must be given in B-system.

the transformation matrix $$R$$ between B-system and inertial system (I-System) can be build out of three matrices.

$$R_x(\phi)= \left[ \begin {array}{ccc} 1&0&0\\ 0&\cos \left( \phi \right) &-\sin \left( \phi \right) \\ 0&\sin \left( \phi \right) &\cos \left( \phi \right) \end {array} \right]$$ $$R_y(\theta)=\left[ \begin {array}{ccc} \cos \left( \theta \right) &0&\sin \left( \theta \right) \\ 0&1&0\\ -\sin \left( \theta \right) &0&\cos \left( \theta \right) \end {array} \right]$$ and $$R_z(\psi)=\left[ \begin {array}{ccc} \cos \left( \psi \right) &-\sin \left( \psi \right) &0\\ \sin \left( \psi \right) &\cos \left( \psi \right) &0\\0&0&1\end {array} \right]$$ where $$\phi\,,\theta\,,\psi$$ are the Euler angles

for example

$$R=R_z(\phi)\,R_y(\theta)\,R_z(\psi)\tag 1$$

or any other combination, each rotation matrix has singularity in one of the Euler angle.

The vectors components in B-system:

$$\vec{g}\mapsto R^T\,\vec{g}$$ $$\vec{F}\mapsto R^T\,\vec{F}$$ $$\vec{M}\mapsto R^T\,\vec{M}$$

with:

$$\dot{R}=R\,\left[ \begin {array}{ccc} 0&-\Omega_{{z}}&\Omega_{{y}} \\ \Omega_{{z}}&0&-\Omega_{{x}}\\ -\Omega_{{y}}&\Omega_{{x}}&0\end {array} \right]$$

thus:

$$\vec{\Omega}=\underbrace{\left[ \begin {array}{ccc} -\sin \left( \theta \right) \cos \left( \psi \right) &\sin \left( \psi \right) &0\\ \sin \left( \theta \right) \sin \left( \psi \right) &\cos \left( \psi \right) &0\\ \cos \left( \theta \right) &0&1 \end {array} \right]}_{J_R} \,\underbrace{\begin{bmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi}\\ \end{bmatrix}}_{\vec{\dot{\varphi}}}\tag 2$$

to see where the singularity is, you invert the matrix $$J_R$$

$$J_R^{-1}=\left[ \begin {array}{ccc} -{\frac {\cos \left( \psi \right) }{\sin \left( \theta \right) }}&{\frac {\sin \left( \psi \right) }{\sin \left( \theta \right) }}&0\\\sin \left( \psi \right) &\cos \left( \psi \right) &0\\ {\frac {\cos \left( \theta \right) \cos \left( \psi \right) }{\sin \left( \theta \right) }}&-{\frac {\cos \left( \theta \right) \sin \left( \psi \right) }{\sin \left( \theta \right) }}&1\end {array} \right]$$

so the singularity in this case is for $$\theta=0$$.

with equation (2) you obtain:

$$\vec{\dot{\Omega}}=J_R\vec{\ddot{\varphi}}+\dot{J}_R\,\vec{\dot{\varphi}}\tag 3$$

put equation (2) and (3) in your ODE's , and multiply the second equation by $$J_R^T$$ you get 6 differential equations

$$\frac{d\vec U}{dt}=\ldots$$ $$\frac{d^2\vec{\varphi}}{dt^2}=\ldots$$

to do the numerical simulation, you must transfer those ODE's to first order differential equations $$\vec{\dot{y}}=\vec{f}(\vec{y})$$

remarks:

if you want other singularity, you can change the combination of your transformation matrix $$R$$, for example $$R=R_x(\phi)\,R_y(\theta)\,R_z(\psi)$$

Edit

other why to simulate :

from equation (2) you get:

$$\vec{\dot{\varphi}}=J_R^{-1}\,\vec{\Omega}$$

now all your ODE's are first order

$$\vec{\dot{U}}=\ldots$$ $$\vec{\dot{\Omega}}=\ldots$$ $$\vec{\dot{\varphi}}=\ldots$$

• Thanks just one thing: if I take the transpose of the expression for $\dot{R}$ that you wrote, and I multiply it for the vector g in the body reference, don't I get exactly $\dot{g} = -\Omega \times R^T g$, which was the expression for the derivative I had. If that is the case, why can't I add those 3 equations and solve for the vector $\Omega$ directly instead of introducing the three Euler angles?? In other words, regarding the original approach I posted in the opener, I would like to understand at what point it is that it fails.
– gc11
Commented May 9, 2020 at 20:19
• Is $(\dot R)^Tg=\dot g$ I don’t think so
– Eli
Commented May 9, 2020 at 20:27
• my interpretation was that in $g = R^T g$, at the LHS the components of g are changing in time, because they are in the relative system, on the RHS instead only $R^T$ changes in time, so if I take time derivatives on left and right I have $\dot{g} = \left(\dot{R}\right)^T g$
– gc11
Commented May 9, 2020 at 20:42
• $R=R(\vec{\varphi})$, $\dot{U}=f(\varphi)$ ,$\dot{\Omega}=f(\varphi)$, where is your differential equation to solve $\vec{\varphi}$ ?
– Eli
Commented May 9, 2020 at 20:52