Is this a bug in the physics engine? (Rotational stability)

There is a well-known effect in rigid-body dynamics, that an object rotating about its "middle" axis is unstable. A small perturbation will quickly result in the object tumbling, and its rotation axis periodically reversing. Weird but true; you can see this demonstrated with a deck of cards in microgravity here.

I've been trying to reproduce this effect for our High Frontier game, and so far, I've utterly failed. I have a simulation that, as far as I can tell, is a direct implementation of the cards-in-space case; check it out here. It's rotating about its middle axis, and when you poke it, you can certainly knock that axis askew. But once you stop poking it, it rotates steadily around its new axis. Under no circumstances do I see the tumbling and axis-reversal seen in real life (or in this Russian simulation for that matter).

Trying to wrap my head around this effect, it seems that it has to do with how the moments of inertia are relative to certain axes. If an object starts spinning around some other axis, then you should compute *new* moments of inertia around that axis. Unity, it seems, isn't doing this. I asked on physics.stackexchange, and there were at least some opinions that this may be a serious flaw in Unity's physics engine.

However, there are probably people here who understand physics (and Unity's physics engine) much better than I do. Can anyone reproduce this tumbling effect? Or explain why we can't?

Many thanks,
- Joe

 

I've tried writing my own little physics code, based on the fundamental relation L = I * w, where L is angular momentum (which stays constant), I is inertia (which can be in either world or local coordinates), and w is angular velocity (in the same coordinates as I).

So what we really want is to start with L, and use that to find w: w = I_inverse * L. And we have to be careful to get from local coordinates to global coordinates.

Here's the code:

Code (csharp):

1. public class CustomPhysics : MonoBehaviour {
2.    Matrix4x4   inertia;       // mass moment of inertia, in body coordinates
3.    Matrix4x4   inertiaInverse;     // cached value of inertia.inverse
4.    Vector3     angularMomentum;   // angular momentum (which remains constant, absent forces)
5.  
6.    void Start() {
7.      inertia = Matrix4x4.Scale(new Vector3(54.9f, 40.9f, 15.7f));
8.      inertiaInverse = inertia.inverse;
9.      angularMomentum = new Vector3(0, 100, 0);
10.    }
11.  
12.    void FixedUpdate() {
13.      // First we need to find our inertia-inverse in world coordinates.
14.      // Do this by applying our current orientation.
15.      Matrix4x4 Q = Matrix4x4.TRS(Vector3.zero, transform.rotation, Vector3.one);
16.      Matrix4x4 Iinv_world = Q * inertiaInverse * Q.inverse;   // I *think* this is right, but...
17.  
18.      // Now, using w = Iinv * L, we can calculate angular velocity in world coordinates.
19.      // Note that this angular velocity is a funny representation of rotation:
20.      // its direction is the angle, and its magnitude is the speed of rotation.
21.      Vector3 w = Iinv_world.MultiplyPoint(angularMomentum);
22.  
23.      // Now, apply that angle and speed (times deltaTime) to our transform.
24.      Vector3 axis = w.normalized;
25.      float speed = w.magnitude;
26.      transform.RotateAround(transform.position, axis, speed * Time.deltaTime * Mathf.Rad2Deg);
27.    }
28. }

You can try this yourself by setting up a block with a scale of (1, 4, 7) — not that this really matters, but that will make appearances match the hard-coded inertia in the code — and then attach this script. Run, and you should see the block spinning neatly around its middle axis... but tweak its rotation around X or Z just slightly (e.g. select the X rotation in the inspector and set it to 1), and it starts tumbling and inverting its axis.

If you start its initial angular momentum at (0, 0, 100), it is much more stable. You can tweak this quite a lot more before any tumbling sets in. All this is quite a bit different from the built-in physics, which I have not been able to make tumble under any circumstances.

So, I *think* this code is doing the right thing. But I'm a little uncertain about that step where I try to convert the inverse of the inertia tensor from local to world coordinates. Unity doesn't have a built-in way (that I can find) to rotate a matrix by a quaternion, so I'm doing it old-school by converting the quaternion to a rotation matrix, and applying that.

Is anybody here comfortable enough with physics to check this code over and tell me if I'm doing something wrong?

Thanks,
- Joe

 

Awesome job Joe, looks great. I'm working on something very similar for our upcoming Game Physics 101 course.

Your simulation works fine, and results in stable rotation about the x and z axis, but unstable about y just as the "intermediate axis theorem" would suggest.

I'll be re-factoring the code to make it easier for beginners to understand, but this is a great result.

 

Thanks Ben, I'm very glad it's useful (and also glad to have an independent opinion on its correctness).

Good luck with your course!

 

Hi Joe, inspired by your solution, I have refactored it a little to produce this code.

I also noticed that the physics engine can calculate approximate inertia tensors from the mass of an object, plus the colliders of child objects.

Thanks for the leg-up!

 

I find cycclone-physics and bullet3 all use this method

I = R * local_I * R^-1 (like the way you using)
τ = I α
w = w + αΔt
(when no τ , then no change of w)
https://github.com/bulletphysics/bullet3/blob/master/src/BulletDynamics/Dynamics/btRigidBody.cpp
https://github.com/idmillington/cyclone-physics/blob/master/src/body.cpp

although "τ = I α" is only correct in some specail situation.
(I find this link form Ben's Lesson)
https://www.4physics.com/phy_demo/newton/newton_rot2.htm

But in one frame , we can still concidered "I" is const ? so
τ=ΔL/Δt
τ=Δ(Iw)/Δt
τ=IΔ(w)/Δt
τ= I α
This is real confused me.

Then when I test this form in cyclone
L = L+ΔL
I = R * local_I * R^-1
L = I w

the result is too unstable

If you set your angularMomentum bigger enough. It should has the similar result.
I don't know why my cell phone won't lead to the result like that. (Maybe somebody will tell me.)
And I guess these engines are intend to avoid unstable rotation.

 

Does anyone know how Rigidbody.inertiaTensorRotation relates to this calculation?