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?

 

I think what is happening here is under Unity's configuration of PhysX (perhaps this logic is built-in to PhysX itself?), it follows that velocity and angular velocity are conserved from frame-to-frame unless acted upon by a force. In the real world, momentum and angular momentum are conserved instead.

Using a more realistic configuration, any applied force is simply added to the rigidbody's momentum vector, and likewise any applied torque is simply added to the rigidbody's angular momentum vector. Unity on the other hand takes the approach of only considering mass and inertia when a force or torque needs to be converted to linear or angular acceleration.

This is why angular velocity doesn't change from frame-to-frame when you have a free-spinning object regardless of what shape it takes. It doesn't get updated unless acted upon by an event that produces acceleration.

I believe what happens when you apply an impule at a position is:

1. Divide the force vector by the object's mass, and then add it to the rigidbody's velocity.
2. Calculate the torque vector (Vector3.Cross(pointForceIsAppliedTo-centerOfMass, forceVector)), rotate the torque vector into the inertia tensor's local space (Quaternion.Inverse(rigidbody.rotation * rigidbody.inertiaTensorRotation) * torqueVector), scale the torque vector component-wise by the inertia tensor, unrotate the torque vector back into world space, and then add the resulting vector to the rigidbody's angular velocity.


I'm guessing this was done to improve stability and performance, and perhaps conserving angular velocity results in more intuitive behaviour. Also, if angular velocity changes as it rotates, the evolution of a rigidbody's rotation is heavily dependent on what time steps you use.

 

Hey Joe, thanks for the inspiring concept! I really struggled with the problems of unities rotation physics and upon trying yours it seemed more realistic.

However I noticed that this way to preserve inertia seems to not consider energy conservation within the object. When testing I saw that only the rotation along the axis with lowest inertia was stable - every other initial rotation including around highest inertia would always drift towards the lowest inertia axis.
To check it further I spawned unit mass objects in the pattern of a block and calculated their kinetic energy while rotating around the centre of mass. And this value kept growing towards a maximum when the rotation ends up along the smallest inertia axis.



Now I have the centrifugal forces of the masses available but I'm not sure how to use it to correct for the drift...