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?

  • 1
    $\begingroup$ 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 $\endgroup$
    – Eli
    Commented May 9, 2020 at 6:38
  • $\begingroup$ @Eli thanks, do you have a reference or know of any examples I could use to understand better the steps I have to do? $\endgroup$
    – gc11
    Commented May 9, 2020 at 10:06
  • $\begingroup$ I will answer your question,soon $\endgroup$
    – Eli
    Commented May 9, 2020 at 13:19

1 Answer 1


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}$$


$$\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] $$


$$\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})$


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)$


other why to simulate :

from equation (2) you get:


now all your ODE's are first order

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

  • $\begingroup$ 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. $\endgroup$
    – gc11
    Commented May 9, 2020 at 20:19
  • $\begingroup$ Is $(\dot R)^Tg=\dot g$ I don’t think so $\endgroup$
    – Eli
    Commented May 9, 2020 at 20:27
  • $\begingroup$ 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$ $\endgroup$
    – gc11
    Commented May 9, 2020 at 20:42
  • $\begingroup$ $R=R(\vec{\varphi})$, $\dot{U}=f(\varphi)$ ,$\dot{\Omega}=f(\varphi)$, where is your differential equation to solve $\vec{\varphi}$ ? $\endgroup$
    – Eli
    Commented May 9, 2020 at 20:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.