0
$\begingroup$

I'm a mechanical engineer who's trying to implement a physics engine for a 3D game simulation, so I apologize for being incorrect or simply ignorant of some aspects of computation. I'm implementing rigid body kinematics and I'm at the point where I need to integrate the equations of motion for a rigid body. I have a doubt on how to do it for the angular momentum. I know that the derivative of the total moment of angular momentum (excuse me if the term it's wrong, I'm a non native english speaker) is (sorry I don't know how to write equations here):

$$ \frac{d\vec{K_{G}}}{{d}t} = I \frac{d\vec{w}}{dt} + \vec{w} \times I\vec{w} $$

The inertia tensor is calculated with respect to a fixed frame with origin in the center of mass $G$, so it's constant, and the moment of angular momentum is calculated with respect to the center of mass. Total moment of angular momentum equals the sum of all external torques with respect to the center of mass, so:

$$ I \frac{d\vec{w}}{dt} + \vec{w} \times I\vec{w} = \vec{M_{G}} $$

So far so good.

I need the first derivative of $\vec{w}$, so I solve for it (again, sorry for the notation):

$$ \frac{d\vec{w}}{dt} = {I}^{-1} (\vec{M_{G}} - \vec{w} \times I\vec{w}) $$

Since I'm using explicit Euler to integrate the equation, I calculate the angular velocity vector for the next frame using the angular velocity from the current frame and the angular acceleration just calculated:

$$ \vec{w}_{i+1} = \vec{w}_{i} + I^{-1} (\vec{M}_{G} - \vec{w} \times I\vec{w}) * \text{frametime} $$

It works fine, the simulation is good as you can see here

But I understand that I can use the impulse-momentum formula as well:

$$ \vec{K}_{{G}_{1}} - \vec{K}_{{G}_{0}} = \int_{a}^{b}\vec{M}_{G}dt $$

with:

$$ \vec{K}_{{G}_{1}} - \vec{K}_{{G}_{0}} = I(\vec{w}_{1} - \vec{w}_{0}) $$

so the integration formula becomes:

$$ \vec{w}_{i+1} = \vec{w}_{i} + I^{-1} \vec{M_{G}} * \text{frametime} $$

where I assumed $\vec{M_G}$ constant for the frame duration.

The simulation runs fine as well, so my question is: I can't understand why I'm using two different expressions for the previous value of $\vec{w}$ and it works fine in both cases (basically I'm not evaluating the change of $\vec{w}$ due to the rotation of the frame): which is the correct approach? and as a bonus question: why I can't see any substantial difference in behaviour?

EDIT I used some latex for the equations

$\endgroup$
  • $\begingroup$ I don't quite understand your formulas, but LaTeX formatting is available. Surround your equations with $ or $$, like so $ \frac{dKg}{dt} = I * \frac{dw}{dt} + w * Iw $ $\endgroup$ – jpaugh Oct 15 '18 at 17:57
  • $\begingroup$ @jpaugh Ok, thanks, I edited the question and added some LaTeX $\endgroup$ – Luca Oct 15 '18 at 18:15
  • $\begingroup$ It would be good to know more about how you reckon the solution is satisfactory, even when the frame rotation term is omitted. Are the trajectories identical? I have a recollection that, even if you omit that term, things can still look more-or-less OK, including conservation of energy, despite the fact that the equations are not being solved correctly. $\endgroup$ – user197851 Oct 15 '18 at 19:29
  • $\begingroup$ The difficulty with this is what the variables in the equations actually mean. If you are in a fixed coordinate system, for an arbitrary shaped body the numerical values of the elements of $I$ depend on the position of the body as it rotates. On the other hand, if you work in a coordinate system fixed to the body, you have to refer the external forces to that coordinate system, which is rotating, but (in the general case) not with constant angular velocity. The "best" way out of this dilemma involves Lie groups, Lie algebra, and quaternions - which is not math that most engineers have learned. $\endgroup$ – alephzero Oct 15 '18 at 21:10
  • $\begingroup$ @alephzero, nope, I don't know what Lie groups or Lie algebra are. I do know about quaternions, though. I know that it makes sense when I project the equation on one or another frame. What I do not understand is how I can write the impulse momentum equation for the rotation in the inertial reference. $\endgroup$ – Luca Oct 15 '18 at 21:30
1
$\begingroup$

EDIT: the OP has provided the answer in a comment, see the bottom of this post. Most of my answer (until the dividing line) is irrelevant to resolving the problem, but I'm leaving it in place as it provides background information and advice on checking the numerical method.

I'm going to try and interpret your equations and offer an explanation. Your notation is different from what I am used to, so I can't guarantee that we'll agree on everything. I'm working from H Goldstein "Classical Mechanics".

Your first two equations together express the relation between angular momentum $\vec{K}$ and applied torque $\vec{M}$ in the stationary, non-rotating, space-fixed frame. So, you seem to be using the subscript $G$ to denote this frame, or perhaps it is to indicate moments taken about the centre of mass. To avoid ambiguity, I am going to use $S$ for this space-fixed frame. $$ \frac{d \vec{K}_S}{d t} = \vec{M}_S $$ This is the basis of your impulse-momentum formula, which applies in the stationary frame. It is possible to integrate this equation forward $$ \vec{K}_S(t+\Delta t) \approx \vec{K}_S(t) + \vec{M}_S\Delta t $$ Some people do tackle the rotational equations of motion this way. But as I think you have noticed, obtaining $\vec{\omega}_S(t+\Delta t)$ involves taking account of the rotation of the inertia tensor (in space-fixed coordinates), over the timestep. You haven't indicated how you are handling this part of the equations of motion. Certainly one can write $$ \vec{\omega}_S(t)=I_S^{-1}(t)\vec{K}_S(t), \qquad \vec{\omega}_S(t+\Delta t) =I_S^{-1}(t+\Delta t) \vec{K}_S(t+\Delta t) $$ if one has $I_S$ at both times. When you write $$ \vec{K}_{G_1} - \vec{K}_{G_0} = I(\vec{\omega}_1-\vec{\omega}_0) $$ it is not clear how you are transforming the body-fixed quantities on the right, with $I$ assumed constant, into the space-fixed form needed on the left, or indeed what $G_0$ and $G_1$ are. So I have some doubts about these equations.

The time derivative of the angular momentum in rotating, body-fixed coordinates may be written (following Goldstein, in your notation, but using the subscript $B$ to indicate vectors resolved in the body-fixed frame) $$ \frac{d \vec{K}_B}{d t} + \vec{\omega_B}\times\vec{K}_B = \vec{M}_B $$ For these quantities we can write $\vec{K}_B=I_B\vec{\omega}_B$ with $I_B$ independent of time. This is your second equation. It is possible to solve the equations of motion in this form, rearranging them as you have done $$ \frac{d \vec{\omega}_B}{d t} = I_B^{-1}\left[\vec{M}_B -\vec{\omega_B}\times I_B\vec{\omega}_B\right] $$ I believe that these equations are correct.

That leads me on to the question: how would you test whether it is working properly? I'm not sure how important this is for your application. In quantitative physics applications, to test the solutions, it is helpful to check some things that are meant to be conserved, rather than just relying on the visual appearance of the dynamics. Suppose we consider one rotating body, being acted on by a torque $\vec{M}$ which is derived from a potential energy function. Then the equations of motion should conserve the total energy (kinetic + potential, or KE+PE). This can be seen because $$ \frac{d\, \text{KE}}{dt} = \vec{\omega}_B\cdot\vec{M}_B =\vec{\omega}_S\cdot\vec{M}_S $$ and it can be shown that the right hand side is the negative of the rate of change of PE, because of the way $\vec{M}$, the torque, is defined. However, if the body-fixed equations of motion are solved without the $\vec{\omega_B}\times I_B\vec{\omega}_B$ term, they will still conserve energy, because this term is perpendicular to the angular velocity. It may be worth checking how well your method conserves kinetic energy in the special case of zero torque.

However, I believe that the incorrect equations (missing that term) will not satisfy conservation of total angular momentum $\vec{K}_S$ in the space-fixed frame. It may be possible to test this with zero external torque $\vec{M}_S$ on the body. For a general inertia tensor $I$, so $\vec{\omega}$ and $\vec{K}$ are not parallel, and $\vec{M}_S=0$ you should see the body wobble around, so $\vec{\omega}$ changes noticeably in the space-fixed frame, but $\vec{K}_S$ should be constant. I recommend doing that test, as well as the energy conservation test.

Just to recap: I'm not sure that your second method based on the impulse-momentum formula, is just the body-fixed equation missing the $\vec{\omega_B}\times I_B\vec{\omega}_B$ term; it may be something more complicated than that. If so, my comments about the consequences of omitting that term may be irrelevant. It is certainly possible to solve the rotational equations, by applying the impulse-momentum form in space-fixed coordinates, and taking care with the rotation of the inertia tensor. But I do have my doubts about the way this is implemented.

And it is also possible to solve the equations in body-fixed form, as you seem to be doing, correctly as far as I can judge. In this case, my only worries are about whether the numerical method is accurate enough, conserving energy and momentum. For this form of the equations (with $\vec{\omega}$ appearing on the right) it may be that predictor-corrector methods are more accurate than explicit Euler. And then, I have no idea how you are tackling the associated orientational equations of motion, which actually rotate your rigid body. Of course these details may not be critical in your application.


EDIT following OP comment.

In the event that the applied torque takes the form of an instantaneous impulse $\vec{M}=\vec{Q}\delta(t)$, where $\vec{Q}$ is a constant vector, and $\delta(t)$ is a Dirac delta function, then the apparent discrepancy disappears. The OP has correctly resolved the problem in this case. In the space-fixed frame, time integrating over the infinitesimal interval during which the impulse is applied gives $$ \vec{K}_S' - \vec{K}_S = \vec{Q}_S \qquad \Rightarrow \qquad \vec{\omega}_S' - \vec{\omega}_S = I_S^{-1}\vec{Q}_S $$ where $'$ denotes the values after the impulse. The space-fixed inertia tensor does not change during this infinitesimal time interval.

In the body-fixed frame the equation for $d\vec{K}_B/dt$ includes an extra term $\vec{\omega_B}\times\vec{K}_B$. However, during the infinitesimal interval of the impulse, $\vec{\omega_B}$ and $\vec{K}_B$ will have only step function discontinuities, no delta-functions, and so integrating them over this interval will give zero. Therefore $$ \vec{K}_B' - \vec{K}_B = \vec{Q}_B \qquad \Rightarrow \qquad \vec{\omega}_B' - \vec{\omega}_B = I_B^{-1}\vec{Q}_B $$ where everything is expressed in body-fixed axes. So, for a delta-function impulse, no error arises from using this equation.

It is easy to check angular momentum conservation (in the space-fixed frame) by converting these body-fixed equations to that frame, using relations such as $\vec{Q}_B= R\,\vec{Q}_S$ and $I_B=R\, I_S\, R^T$, where $R$ is the rotation matrix describing the transformation.

$\endgroup$
  • $\begingroup$ I sorted out the problem: when I am writing $$ K_{G} $$ I'm considering the total moment of momentum with respect to the center of mass G (I use G since I need a point with velocity parallel to G, or fixed in the inertial frame. Since I'll consider the body's velocities relative to G it make sense to use G as the pole for the moments). What I saw in my books was the following formula: $$ \triangle K_{G} = \int M_{G}dt $$ The book considers an impulsive moment acting for an infinitely small instant (dirac delta), so I consider the inertia tensor I constant, hence $$ \triangle w = I^{-1}M_{G} $$ $\endgroup$ – Luca Oct 27 '18 at 14:37
  • $\begingroup$ I believe that this condition (delta-function impulsive torque) does indeed resolve the problem, and you have indeed sorted it out. It's a pity that my answer turned out to be mostly irrelevant; nonetheless I've added a couple of edits at the top and bottom of it, to clarify. I'm deleting my other comments here to tidy up. $\endgroup$ – user197851 Oct 28 '18 at 11:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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