Rigid bodies mainly have 6 degrees of freedom in 3 dimensions (3 translational and 3 rotational), you can represent the particles as rigid bodies with only 3 translational degrees which is straightforward and their equations of motion aren't difficult to formalize, but what if you want to extend from particles to much more than that? Say... system of particles connected to each other to create a whole body!
First, let's suppose you have two particles, one that has an infinite mass which cannot represent motion unless an infinite force is applied, and the other one which has some finite amount of mass which is able to move around, but now let's actually link them together by a massless rigid rod or an infinitely stiff spring like so
Here because the rigid rod is infinitely stiff and the yellow particle is unable to move, then that means the green particle is unable to move away from or closer to the yellow particle because remember, the rod is infinitely stiff and it doesn't deform at all, and the link itself doesn't really have any side-effects that puts into the green particle, all of this starts to create a very interesting motion.
And the motion it creates is in the following video
rotational-motion-test.mp4
As you can see, in both the particle simulation and the rigid body simulation, the green particle was collided by the blue particle, the green particle didn't move in a straight line in both cases, but it did move in a circular line, and it's speed was exactly the speed of the blue particle throughout time, so you can think of the green particle as if it was moving in a straight line still (Which would become more clear later on!). And since it's a rotational motion, we can actually represent or approximate it's motion in polar coordinates.
Note
in the particle simulation case, the green particle was losing some bits of energy after the collision because of the way it works, it mainly uses system of constraint forces to simulate it. So the lower the step size is the better, and imagine the case of infinitely small step size where it's exactly correct for the rest of the simulation, but that's impossible with the constrained computer resources and floating point precession issues.
So suppose you have linear velocity going in a straight line, then converting it into an angular velocity is going to require us an origin to rotate around and it's distance from the origin, we pick an
Now the formula may not make sense at first, but suppose the particle from the origin is at a distance of 1, then it's speed is exactly aligned with the rotation speed, however if it's more or less than 1, then it's speed would still be assumed at a distance of 1 which is not correct! Because of this we will get speed that's faster or slower than expected. To fix this, we essentially divide the linear velocity by the magnitude of
To give a more concrete example, suppose we have three particles at radius of 1, 1.5 and 2. (The manifold sizes are not to scale here! But the angles are accurate)
And each particle is at speed of
Notice the particle that is in radius 2, it's half angle of radius 1, but their arc lengths are the same! The arc length in this context is the distance travelled over the circular manifold's surface, it's defined as
So, multiplying the angle
Within desmos, there are three points on top of the ones that are on the circular manifolds, now in the first row try to increase the
So one way to think of rotating points is that they are essentially moving in a straight line in their own space, but relative to you they are rotating, and as we've noticed their distance from their origin makes a difference for us in our familiar cartesian coordinate system. This is a very important point to digest so take a moment to think about it and keep playing around with the desmos demo.
Okay so that out of the way, there's a problem here, the angular velocity is only a magnitude in this case, but it does not have a representation for direction! We mainly need to know which direction it's rotating at as well, this is where the concept axis of rotation comes into play.
Important
In case you don't know what an axis of rotation is, I highly recommend reading the Appendix A which is at the bottom of this article, it contains a short explanation as well as pointers to existing resources.
We could do a neat little trick and use a product that does this for us, the one we're mainly interested in is called the cross product. By definition
Where
The mysterious
Notice the rotational direction, it rotates in the same order as the product
And compressing the 3D formula gives us the following
So really, all we did was to try and get the axis of rotation which also encodes a direction and then scale it by the original rotation rate magnitude, sweat!
rotational-velocity.mp4
Now that's nice and all, but how do we accelerate the particle as well as taking it's mass (which is it's resistant from external forces) into account? We can first look into the concept of momentum.
By turning vectors into scalars (for now which is more correctly defined in 2D), the linear momentum of a particle is defined as
Where
Or to be more specific
We can turn it into angular momentum by
Then we can now separate mass
Here notice that we get both
Warning
The concept of moments of inertia is usually understood in the wrong way unfortunately, the idea of "how hard it is to rotate a body" is vague and is not clear on what it actually means. The chapter on system of particles later on would get into it in more detail using collisions with the conservation laws. It's not obvious right now on what it is exactly, but we'll go over the acceleration part quickly for it.
Alright that seems handy, however we still need to accelerate the particle. We can take the time derivative of the angular momentum to get the equivalent of linear force
Here
angular-acceleration.mp4
Now getting back to the moment of inertia, suppose now you have a torque to be applied to a particle of mass 1, then
So, assuming that
Then if force
Notice that it moves at half the angle, but remember that it's only an angle! The arc length however is the same as the one of radius 1, using our previous formula on getting the arc length then gives us
Now we haven't talked about the case where the free particle doesn't hit the constrained particle at a perfect angle, what happens then? The
Ok so far we have derived pretty much the equations of motion for a single particle around an infinite mass, but what if we extend it to more than just one particle? Which then becomes system of particles tied altogether!
It's actually surprisingly straightforward, remember when we said that rods generally are rigid or infinitely stiff? Since they don't make side-effects, then we can treat multiple particles as one! (Sort of...)
And to prove that is the case, I've again modified the simulation to handle it
system-particle-rigid-comp.mp4
As shown, you would notice that the forces are being balanced exactly the same way as the rigid body equations of motion method. But notice now that it is harder to rotate than before and the blue particle now bounces off a little.
To explain all of this, we need to first look into how linear motion works in one dimension with system of particles, which should hopefully make it easier to understand the rotational case with a clever example.
Suppose you have multiple particles tied together with rigid rods, then they are essentially one whole body, which also means we can treat them as a giant particle. Let's ignore the rotational aspect for now and focus on the linear motion, we have 4 particles and each particle has a mass of 1, then the total mass of the system is 4, very simple
If the system is moving at a velocity of one, then the total momentum is
That's the same as a single particle moving at a velocity of one with mass 4.
So rather than thinking of a body as multiple particles, we could instead unify the particles into what's called center of mass and base our equations of motion on it. Essentially, center of mass is defined as
The moving center of mass is just the whole body moving as well, taking it's derivative gives us the velocity of the system
This may not sound important for now, but it will be important shortly once we dive further into rotations.
If we have a particle of mass 1 and we collide a body with it, then we get the following behavior using the particle simulation
one-vs-four-collision.mp4
The behavior shown was an elastic collision. Notice in the above system the blue particle stopped since it's entire kinetic energy was entirely transferred towards the green particle, in the bottom system however the gray particle didn't stop and it instead bounced off with less energy than before.
We could verify this using the collision force formula which makes sure the conservation of momentum and the conservation of energy laws are taken into account. If you want to understand how the collision forces are derived and used then seek Appendix C which uses the conservation laws only, along references of [12] and [13].
Now calculating the impulse for the above system is as such
Then by newton's third law
Then applying it to both particles
Which is exactly what we saw in the video! Doing it again for the underneath system gives us
Then again applying it to both particles
Which is also correct!
So now that we know how linear motions work and how they conserve momentum, let's move on to the rotational part which isn't going to be any different really.
It's important to notice that the linear motion was actually a rotation around a point at infinity with an infinite mass and with an infinite radius.
You may ask: "What? How is that even possible?". It's better to describe with an image of various polar coordinates with various center of masses colored in pink with the connected particles colored in green acting as a system
Notice in the coordinate with an infinite radius the masses are essentially approximately moving in a straight line, as the radius gets smaller the lines get curved and the center of masses get closer, eventually the connected green particles start to look like they are rotating about some origin within some small range of radius.
So what does this have to do with our linear motion consisting system of particles? Well, notice in the coordinate with an infinite radius, the two particles are within a system moving in a straight line, if we had a particle colliding with a system in a straight line then we get the behavior we just analyzed before.
To prove this, let's now compare both linear and angular motions of two systems of particles (or bodies) that get collided. We will assume the center of mass of the linear part is at infinity (even though it's in the middle of the system of particles, but that's besides the point).
Here's one with equal masses
equal-masses.mp4
And another one with un-equal masses
un-equal-masses.mp4
Notice that both linear and rotational motions behave exactly the same way in this case! So our theorem is proven.
Though why stop there with the demos, let's look at the math as well.
Suppose the following formula from Appendix C
It would seem daunting at first, but our goal is to get to the same formulation as the linear impulse assuming the radius is at infinity.
Let's assume the radius for both bodies A and B are at infinity, so that means
And that applies to
Then
Bingo! So re-arranging the original formula based on what we found gives us
Which means that it is in fact a rotational motion at infinity as well. And by the way, this also applies to rotations around some near origins as long as the
So now that we're done explaining the conservation behaviors within systems. Let's get into the concept of angular momentum with system of particles, the total angular momentum of the system is defined as
You could also take the derivative of the total angular momentum to get the total torque as well
Then we can now get the moment of inertia for the system as follows
Then multiplying the total angular momentum by the inverse of the total angular velocity gives us the total moment of inertia!
Which essentially equals
One thing to notice however is that let's suppose you have a ring of radius 1 that contains few particles, and one particle is in another ring of radius 0.5 like so
As the outstanding particle starts to get closer and closer to the center of mass, it starts to have less influence on the inertia of the system, and it makes sense, the infinite mass has no influence on rotations at all to begin with!
There's also another case where you have a very far point mass and a close one, if you collide with the close one you may realize that you need a lot of energy to rotate the whole body than if you collide with the far one. This is because the amount of distance to travel has increased significantly, hence
Also so far we have assumed we had an infinite mass particle as our center of rotation so far, but we haven't talked about the case where there isn't an infinite mass, what happens?
Our body has been rotating around an infinite mass particle, that particle defines our center of rotation which is also at the same time the center of mass since it's the heaviest mass in the system. You can sum the masses and average them out like so
This formula would give you a point that's the center of the system. To prove this, suppose one of the particles has a large mass of 1000 in comparison to the other particles that have mass of 1, then summing them up and dividing by the total mass of the system gives us the center of rotation which lies the closest to the large mass particle, which as you expect if you have an infinite mass then it directly sits within that mass.
A better example would be a line shape consisting of two particles of mass 1
The pink point that's at the center of the system is the center of mass.
The
Which essentially says the vector that takes from the center of mass
Okay so now let's talk about both linear and rotational motions combined, what kind of motions do we get in which particular situation?
Suppose in a scenario where we do a collision of one particle that goes right towards the center of mass of a body like so
collide-com.mp4
There are no rotations involved, just purely linear motions. The reason for this is that all the energy has gone straight through the center of mass which is the entire system by itself.
However that's not primarily the case when we don't collide into it like so
collide-rot.mp4
Notice that after one of the purple particles were collided, the yellow particle stopped completely! And the collided purple particle started moving in a straight line, but it was also grabbing the other one it was attached to, and since it was grabbing it, it was also being pushed towards it since there's two forces that are being extorted by the rod according to Newton's third law, this back and forth behavior creates a rotational motion as a result that keeps going on forever like so
This behavior would be important later on when we touch on 3 dimensions, especially with a term called "inertial torque" and the conservation of angular momentum of a system.
Also as you would've probably noticed from the video above, the center of mass was moving in a straight line!
Looking at the motion from the perspective of the center of mass for the rigid body shows that it is rotating and moving at a constant velocity
rotation-com-persp.mp4
Also these sorts of motion have a connection with the
Handling rotations in 2D is mostly straightforward, however in 3D it gets a lot more complicated, but the exact same principles transfers over from 2D to 3D just fine, just need a slightly different way of looking at it.
And there's really nothing new in 3 dimensions except for one interesting term that appears when taking the derivative of the 3D angular momentum and the fact that the derived moment of inertia usually isn't aligned with the orientation of the rigid body in world space.
So far all our particles have been rotating within a plane where a whole rigid body has 1 degree of rotation, in 3D however a rigid body has 3 degrees of rotations, each degree of rotation represents a rotational plane.
To extend our momentum equations to handle various oriented planes, we need to use the cross product. The angular momentum of a particle is
Which so far seems similar to the 2D formula, nothing mysterious yet...
Now the moment of inertia in 3D however is interesting, we will try to derive the total angular momentum again using the cross product like so
Where
Also be aware that it's an object and not a transformation operation, which makes it a so called "rank-2 tensor" and that raises a question: why?
Normally if you analyze the angular motion of a 3D rigid body you would notice that the total angular momentum doesn't point in the same direction as it's total angular velocity, which seems weird when you think about it, though there's another way of looking at it but that's for a bit later.
Taking the time derivative of the angular momentum is also interesting, the derivative of the moment of inertia would be quite dense to derive by hand so we will use a concept of body and world space to make things easier.
Note
There's a document that includes the derivation for the derivatives of moments of inertia by hand in 3 dimensions in reference [14], specifically section 3.0.2 called "Angular Motion".
To turn the angular momentum from body space into world space we use the rotation matrix
Taking the time derivative of the world space angular momentum then becomes
The total sum of the external torques and the inertial torque is sometimes referred to as "Euler's second law of motion". Also notice the
Also the term is more known to be called "Torque-free precession", the Torque-free part is where it's interesting.
Note
In case you're confused by the
To get a velocity of a point in space that's rotating around an axis, we can do
Which is essentially just the derivative of
The time derivative of
Which gets us the velocity in world frame. Now we could interpret the world angular velocity as a tensor
One important thing to mention is the cross product skew-symmetric matrix, it's defined as
It's a rank-two tensor! Which is essentially
Then suppose we have an angular velocity tensor that's in the body frame and we want to transform it into the world frame, the general rule for transforming the coordinates of rank-two tensors in rotating frames is
Look look up reference [18] for a quick proof or other resources on tensors. Now we could reformulate our derivative of
Notice
Caution
The inertial torque term would make no sense looking at the formula alone, you have to observe it by looking at a rigid body in motion with antisymmetric mass distribution and do analysis with the term, most resources would be completely hand-wavy about it by bringing up known effects that are caused by it, the most known being the Dzhanibekov Effect, but that is not the best way to understand it! It's actually a fairly intuitive term that causes an obvious motion which isn't difficult to digest, but was completely overlooked by ridiculous examples unfortunately.
There's a very good video on doing the analysis of the inertial torque term that's in reference [1]. I did however try my best to explain what it does in the next few paragraphs which is a little different than the conventional way of teaching it, it certainly worked for me.
Suppose you have two particles A and B that have the same mass that are also on a sphere's manifold surface separated by some angle, and they are constrained to some infinite mass particle by some distance so A and B cannot move away from it, also A and B can't move away from each other, we want to apply an impulse
Then it starts to move, but it also drags particle B into it as well, and that also in return pulls A into B by the rod's forces as per Newton's third law, and this would cause a rotational motion and a translational motion on the sphere's surface manifold, and this already sounds very familiar to the 2D case!
3d-rot.mp4
Looking at it from the perspective of the particles center of rotation shows an interesting behavior
3d-rot-persp.mp4
This looks very similar to the case we've analyzed before and it's happening again here. Now you may be wondering: "What does this have to do with the inertial torque term we just derived?"
Remember when we mentioned there are cases where the angular momentum isn't always aligned with the angular velocity direction? That is mainly dependent on the mass distribution of the system, if you apply an impulse you may not get the angular velocity direction you would expect. To give an example, we have to rely on an actual physical intuition rather than the equations for things to make sense.
As shown above, when particle A started moving, it started to "drag" particle B into it, but how exactly does it do that with only using angular velocity? By taking the conservation of angular momentum into account!
More specifically, the momentum of the system is by average rotating in a very specific motion depending on the initial conditions. when A started moving, it had to drag B into it as well, and that also means both must rotate around their principal axis of rotation since their angular velocity is not constant and is not parallel to the angular momentum! It's almost as if we are doing both translation and rotation at the same time within the spherical surface manifold.
So the angular momentum must be conserved throughout time, it's direction and length must not change if there are no external torques applied on the system (which is where the torque-free part comes from within the "Torque-free precession" term), which is a constraint! The following video shows the behavior more clearly
ang-moment-conserve.mp4
The red vector is the angular momentum, the yellow vector is the angular velocity and the green vector is the inertial torque which is pushing angular velocity in it's direction. Notice how the entire system is rotating about the angular momentum as we would expect!
And here's a comparison between the particle sim (where it approximates angular momentums and the inertial torques as vectors as shown) and the rigid body sim, you would notice their motions match and conserve momentums in the exact same way
inertla-torque-compare.mp4
To explain the behavior in more detail, notice in the first frame of rotation, particle A was rotating around particle B, but as the time goes it started to rotate around it's center of rotation and then later on stopped and particle B started to rotate around particle A instead, this cycle continues forever, this is exactly the behavior we got in the 2D case! Only now that it looks more obfuscated than before.
Now I'll be trying to shove the 3D case into the 2D case and emulate it in 2D, notice the angular velocities are right above the center of rotation
Each step is rotating the systems around the dotted angular velocities, notice their center of rotation is essentially moving in a straight line, you can try doing this by hand and show that this behavior actually works. This is almost like a workaround for the fact that you don't have a translation state to keep track of when moving in a spherical manifold.
Here's another one doing the exact same thing but in 3D
All of this can be linked back to other types of behaviors like Gyroscopic Precession or the Dzhanibekov Effect, the important thing to remember here is that all of these behaviors are caused by the body's own inertia which we saw how that happens already. In case you're interested, look up the two references [3] and [4] which are mainly about the analysis of the Dzhanibekov Effect and the Gyroscopic Precession behaviors respectively.
Now that out of the way, there's another thing to mention, if a rigid body gets rotated, then it's moments of inertia changes! Not because it deforms for our case, but mainly because the distribution of the particles in the system changes, which also affects how the energies gets distributed by the external and the internal forces. A good example is the following
Here we apply in both systems an impulse to one of the particles by the blue vectors, notice their resultant angular velocities (in yellow) aren't the same! So how do we calculate the new inertia tensor? Do we need to re-do all the work again?
Not really, since we know that points never deform in our case, then the distance from each other must be constant. We could instead orient a radius vector
Note that both the dot products and the outer products have their equivalents
A good way of interpreting
Notice the
Applying forces in either the body's frame or the world frame give the exact same results. The
Now that we're done explaining everything we needed to know, we're ready to implement the dynamics, which, for a basic one, is pretty simple.
Warning
The following implementation is not meant to be the best way of doing it, it's meant to be a good starting point which is probably good enough for most people. And be aware that there are tons of ways to implement rigid bodies, this is not the only way of doing things at all.
There are many sorts of things you need to handle properly to get a very robust dynamics, this article in particular does seem to do a very good job at covering it https://mizerski.medium.com/motion-of-spinning-bodies-b51175194512, along Erin Catto's slides on numerical methods as well https://box2d.org/files/ErinCatto_NumericalMethods_GDC2015.pdf
A typical rigid body data structure is something like
struct Rigid_Body
{
Vector3 position;
Quaternion orientation;
Vector3 center_of_mass;
Vector3 inv_mass;
Vector3 linear_velocity;
Matrix3x3 inertia;
Vector3 angular_velocity;
};Now the inertia tensor is the interesting part, as per the equation we defined earlier, you could have each individual particle's inertia calculated and summed up altogether, or if you have a pre-defined shape like a box or a capsule then you may need to pre-calculate them using a continuous integration over a volume, I've done some derivations in Appendix B already so you may wanna check that out.
We will use the box shape and calculate it's inertia like so
Matrix3x3 calculate_box_inertia(float min_width, float max_width, float min_height, float max_height, float min_depth, float max_depth, float mass)
{
float dx = max_width - min_width;
float dy = max_depth - min_depth;
float dz = max_height - min_height;
Matrix3x3 inertia = {};
inertia[0][0] = mass * (dy * dy + dz * dz) / 12.0;
inertia[1][1] = mass * (dx * dx + dz * dz) / 12.0;
inertia[2][2] = mass * (dx * dx + dy * dy) / 12.0;
return inertia;
}We will define some helper functions which would abstract some work for us to be used later. First, we may want to get the inertia tensor in the world frame, we just do the following
Matrix3x3 get_inertia_world_frame(Rigid_Body& body)
{
// world space orientation matrix
Matrix3x3 world_orient = convert_quat_to_matrix3x3(body.orientation);
// update the inertia to be in the world space configuration
Matrix3x3 world_inertia = world_orient * body.inertia * transpose(world_orient);
return world_inertia;
}
Matrix3x3 get_inverse_inertia_world_frame(Rigid_Body& body)
{
// world space orientation matrix
Matrix3x3 world_orient = convert_quat_to_matrix3x3(body.orientation);
// take the inverse and update the inertia to be in the world space configuration
Matrix3x3 world_inertia = world_orient * inverse(body.inertia) * transpose(world_orient);
return world_inertia;
}Then getting the world center of mass is like so
Vector3 get_world_center_of_mass(Rigid_Body& body)
{
return rotate_vector3_by_quat(body.orientation, body.center_of_mass) + body.position;
}Now It may depend on how you want to handle applying forces to bodies, the common way is to apply forces as impulses which is straightforward
void apply_linear_impulse(Rigid_Body& body, Vector3 impulse)
{
body.velocity += body.inv_mass * impulse;
}
void apply_angular_impulse(Rigid_Body& body, Vector3 impulse)
{
Matrix3 inv_inertia = get_inverse_inertia_world_frame(body.inertia);
body.angular_velocity += inv_inertia * impulse;
}
void apply_impulse(Rigid_Body& body, Vector3 world_point, impulse)
{
apply_linear_impulse(body, impulse);
Vector3 world_COM = get_world_center_of_mass(body);
Vector3 r = world_point - world_COM;
Vector3 H = cross_vector3(r, impulse);
apply_angular_impulse(body, H);
}Then a box rigid body constructer
Rigid_Body create_box(Vector3 positiion, Quaternion orientation, Vector3 extent, float mass)
{
Rigid_Body body = {};
// set an initial position and orientation
body.position = position;
body.orientation = orientation;
// the inv_mass is used to describe infininte masses, 0 is infinite,
// mass of 20 of a body is just it's inverse 1/20 and so on. to get back
// the normal mass you just invert it "mass = 1/inv_mass".
body.inv_mass = (mass != 0.0) ? 1/mass : 0.0;
// here we use the "extent" which is the maximum width/height/depth
// of the box altogether.
body.inertia = calculate_box_inertia(-extent.x, extent.x, -extent.y, extent.y, -extent.z, extent.z, mass);
return body;
}Now updating the rigid body's state would have one minor issue, the inertial torque cannot be calculated if we don't have a second moment of inertia in world space, the workaround is to work within the local space since the inertia tensor is constant there, hence we can just calculate it like so
void update_rigid_body_state(Rigid_Body& body, float dt)
{
body.position += body.linear_velocity * dt;
// in case the angular velocity is not a zero vector then update the orientation,
// this is mainly to make sure we don't produce bogus results out of zeros
if body.angular_velocity != get_zero_vector()
{
Matrix3x3 world_orient = convert_quat_to_matrix3(body.orientation);
// angular velocity in the body frame
Vector3 local_omega = transpose(world_orient) * body.angular_velocity;
// this is inertial torque converted to an impulse
Vector3 inertial_impulse = cross_vector3(local_omega, body.inertia * local_omega) * dt;
// calculate the jacobian of the impulse
Matrix3x3 J = body.inertia + dt * (skew_matrix3x3(local_omega) * body.inertia - skew_matrix3x3(body.inertia * local_omega));
// single iteration of the newton's method
local_omega = local_omega - inverse(J) * inertial_impulse;
// transform the angular velocity back into world space and then assign it into the rigid body's angular velocity
body.angular_velocity = world_orient * local_omega;
// get the new delta of the rotation
Quaternion theta = quaternion_angle_axis(
/* angle */ length(body.angular_velocity) * dt,
/* axis */ normalize(body.angular_velocity)
);
// make sure to normalize the quaternion so things don't diverge too much by the accumlating floating point precision error
body.orientation = normalize_quat(theta * body.orientation);
// calculate the new position relative to the oriented center of mass in case it's not within the origin of the local coordinate frame
Vector3 world_COM = get_world_center_of_mass(body);
body.position = world_COM + rotate_vector3_by_quat(theta, body.position - world_COM);
}
}Believe it or not, that's pretty much the bare minimum to accurately simulate a rigid body!
One thing to note however is within the update_rigid_body_state(...) function, the inertial torque is calculating an implicit force, to do this we needed to use newton's method, this requires multiple iterations since a single iteration is not enough from my experience and sometimes it can produce really bad results. The following is an example of how to do proper integration to make newton's method more into an effect with substepping
void update_simulation(Physics_State& state, float dt)
{
// get time difference between frames
// ...
// apply gravity to bodies
// ...
// some other stuff applied to bodies
// ...
// do multiple iterations on every rigid body with substepping to make
// the inertial torque more "stable"
float num_of_steps = 10;
float substep = dt / num_of_steps;
for int step = 0; step < num_of_steps; step += 1
{
for &body in state.rigid_bodies {
update_rigid_body_state(body, substep)
]
}
}This should make the inertial torque effect more stable since now we're doing more iterations than one, you can increase in certain cases if 10 iterations isn't enough, for me in the extreme cases (which are unlikely in the normal scenarios) I had to use between 100-225. Also adding air drag may stabilize the motion as well. For more details on this, refer to Erin's slides on numerical methods.
Now what can we do with all of this? Pretty much anything you could think of! Multibody dynamics which involves gluing multiple rigid bodies into one giant ragdoll or a vehicle or whatever using various types of constraints. Maybe also combining it with fluids like air or water as well.
All of that sounds nice, but maybe we should do something easier and more achievable, like attaching a rigid body to a spring or simulating the Dzhanibekov Effect or even collisions with the floor.
This is simply done using Hooke's law, the force of a spring is described as
Where
Implementing that in code
Rigid_Body body = /*...*/;
// point mass that's from within the rigid body, which is
// also in the rigid body's local space and NOT the world space!
Vector3 point_mass = Vector3(/* ... */);
// state varbiels
float previous_length = 0.0;
// simulation loop
while keep_simulating
{
// get the time difference between frames
float dt = /* ... */;
// apply gravity
apply_linear_impulse(body, Vector3(0,-9.8,0) * dt);
// the point mass in world space
Vector3 world_point_mass = rotate_vector3_by_quat(body.orientation, point_mass) + body.position;
// some arbitrary point in world space which is our target
Vector3 target_point = Vector3(/* ... */);
// calculate the spring length
Vector3 spring_len = (target_point - world_point_mass);
// calculate the spring force that forces the body's point mass
// to move towards the arbitrary world point
float spring_stiffness = 15.0;
Vector3 spring_force = spring_len * spring_stiffness;
// calculate the damping forces
float damping_magnitude = 1.0;
Vector3 spring_velocity = (spring_len - previous_spring_len) / dt;
Vector3 spring_damper = damping_magnitude * spring_velocity;
previous_spring_len = spring_len;
// apply forces to the body
apply_impulse(body, world_point_mass, (spring_force + spring_damper) * dt);
// ...
// update rigid body state
// ...
}Then we should get something similar to the following
spring-forces.mp4
The white point is the target point and the yellow point is the point mass that's on the box body.
This one is pretty simple, we set the rigid body at a specific initial orientation and an initial angular velocity where it makes this effect doable. For more details look up reference [3].
The idea is that the effect only happens in the intermediate axis of the inertia, we know that in 3 dimensions there are three rotating planes, so suppose that one of them represent the least moments of inertia which is
For the following, the intermediate inertial axis would be the x axis which is in between y and z axes, y is the largest moment of inertia and z is the least moment of inertia
dzhanibekov-effect.mp4
And here's another showing the angular momentum (in red) and angular velocity (in yellow) vectors (both scaled down for visibility)
dzhanibekov-velocities.mp4
And here's a comparison to a particle sim which relies on the same constraints as the rigid bodies in terms of making the effect doable
dzhanibekov-effect-particles.mp4
We could try doing some collisions with spheres, I've already covered most of the important bits in Appendix C. I think showing some code first should be self-explanatory.
The moment of inertia for a sphere is
And is calculated as follows
Matrix3x3 calculate_sphere_inertia(float radius, float mass)
{
Matrix3x3 inertia = {};
inertia[0][0] = mass * radius * radius * 2/5;
inertia[1][1] = mass * radius * radius * 2/5;
inertia[2][2] = mass * radius * radius * 2/5;
return inertia;
}And here's how to calculate the collision impulses for the sphere body
Rigid_Body sphere = /*...*/;
#define SPHERE_RADIUS 0.35
// simulation loop
while keep_simulating
{
// get the time difference between frames
// ...
// apply gravity to rigid bodies
// ...
// the collision normal of the floor
Vector3 collision_normal = Vector3(0, 1, 0);
// is the sphere interpenetrating with the floor
if dot_product(collision_normal, sphere.position) < SPHERE_RADIUS
{
// get our contact point
Vector3 contact_point = sphere.position - collision_normal * SPHERE_RADIUS;
// get absolute value of the depth
float depth = abs(contact_point.y);
// radius vector
Vector3 R = contact_point - sphere.position;
Vector3 relative_velocity = sphere.linear_velocity + cross_vector3(sphere.angular_velocity, R);
// if they are not separating already, then start applying impulses
if dot_product(relative_velocity, collision_normal) > 0
{
Matrix3x3 inv_inertia = get_inverse_inertia_world_frame(&sphere);
float linear_factor = sphere.inv_mass;
float angular_factor = dot_product(inv_inertia * cross_vector3(R, collision_normal), collision_normal);
float restitution = 1.75; // (1 + e) = (1 + 0.75) | where e = [0,1]
float normal_impulse = -restitution * dot_product(relative_velocity, collision_normal) / (linear_factor + angular_factor);
// apply normal impulse
apply_impulse(&sphere, contact_point, -normal_impulse * collision_normal);
// get the friction normal by projection the velocity into the
// collision normal plane and normalized it, only apply it
// if there's any movement within the xz plane
Vector3 tang_dir = relative_velocity - collision_normal * dot_product(collision_normal, relative_velocity);
if tang_dir != get_zero_vector()
{
tang_dir = normalize(tang_dir);
linear_factor = sphere.inv_mass;
angular_factor = dot_product(inv_inertia * cross_vector3(R, tang_dir), tang_dir);
float friction = 0.6;
float friction_impulse = -friction * dot_product(relative_velocity, tang_dir) / (linear_factor + angular_factor);
// apply static friction
apply_impulse(sphere, contact_point, friction_impulse * tang_dir);
}
}
else {
// don't do anything since they're moving away from each other!
}
// project the sphere into the floor's surface to remove the interpenetration
sphere.position += collision_normal * depth;
}
// ...
// update rigid body state
// ...
}And here's a simulation of it
sphere-collision.mp4
Now that out of the way, let's explain what's happening here. The interesting part of all of this is friction.
Without friction, the ball wouldn't rotate because the normal impulses go straight through the center of mass of the ball which just moves the entire rigid body in the normal's direction and that's to be expected. However with friction, that's not the case! The friction impulses do not go through the center of mass and we could make this more obvious with an image of how friction actually happens at the micro level
Notice within the image, in reality bodies surfaces aren't perfectly smooth and their surface usually collide with other materials at the micro level causing this behavior which also forces the ball to roll.
Note
Now you may also know that objects have two states of friction, one where they try to stick to the surface of another material, and another state where they just can't and they start slipping. Both states are called static and dynamic frictions respectively. A known model of this is called Coulomb's friction.
What we've modelled however is the static friction only for demonstrating the frictional behavior and how it relates to the rotational dynamics of the bodies, you could model the dynamic friction if you want.
Another thing to mention is that collisions on shapes that aren't spheres like boxes do cause rotations even in cases there are no frictions involved.
So far the spheres could only have one contact point with the floor, but with boxes they can have as much as four simultaneous points altogether on the floor! To demonstrate this let's show the problem in more detail with an image
The case of two contacts defines an edge collision with the floor, there will be two torques that when summed up, you get the resultant torque in circular green which rotates the entire body in one axis and hence we get
Now unfortunately though we can't simply model this in the simulation, the reason being is that if there are multiple contacts happening at the same time then they would affect each other and this creates a system of forces which is not trivial to handle. This is a very off-topic discussion and a very complicated problem to solve and does not fit this article but you can look into how people handle this problem with the so called Systems of Constraints Solvers, essentially you have a contact which could be made up of three constraints at once and you also get three other contacts which then totals the constraints to be 12 for four contacts, this cannot be solved directly with a regular system of linear equations solver however because of constraint conditions which are the inequalities (i.e. when the two bodies are already separating, then no forces are to be applied, so C >= 0 or C = 0) they eventually turn into a so called Non-Linear/Linear Complementarity Problem (usually thrown as terms like NCP and LCP respectively) also sometimes referred to as a Mixed-LCP for systems that involve different types of equality constraints altogether, I would recommend looking up three references which are [24], [25] and [26], also look into Solver2D which delves into iterative solvers that solve NCP and LCP problems in different ways for different types of constraints which are possibly not necessarily the most accurate but they are efficient for what they were designed for.
What we've derived so far throughout the article was more known as Newton-Euler Equations, It's essentially the combination of Newton's equations and Euler's equations into one, though we primarily focused on Euler's equations for the most part.
I've tried to cover the most important bits in this article and clarify what they actually are and where they came from. There are however a few things that I haven't covered like the alignment of the angular velocity and one of the inertia's principle axes (or the eigenvectors) which then causes a "stable" motion, we've covered the instability of rigid bodies which is related to that, there's also the fact that you can get a diagonalized moment of inertia tensor which doesn't involve it's off-diagonal elements which is pretty handy when doing computations on rigid bodies on paper or even simplifying computations when doing rigid body simulation, here's a resource that goes into it, you basically need to find a specific orientation by a rotation matrix
I'm looking forward to do more mechanics and such using Lagrangian Mechanics, one resource I'm aware of is called "Structure and Interpretation of Classical Mechanics" which is reference [19], another one is in reference [11] which I've seen mentioned quite a lot in various resources but I haven't delved much into it however.
Alright, that's about it. Thank you for reading the article and have a good day!
Best way to describe a rotation is by using Euler's formula
Every time you multiply the complex number by the imaginary number
The dot product of two vectors
The
Also the cross product requires a specific order, in general
In case we do
The handedness in question is a common figure, if you crawl your fingers and put your thumbs up, then your thumb is your axis of rotation and your crawled fingers are your rotation direction around the axis, it's essentially the following
Both handedness have the same axis of rotation pointing in the same direction, but the rotation direction differ like in our case (counter-clockwise for right-handed, clockwise for left-handed), because of this behavior, the vector coming out of the cross product is technically not a "vector", but instead a "pseudovector" or even a "bivector", but that's not important for now as long as we're consistent with the handedness of the coordinate system.
Now we can map
You can think of the formula as a coordinate transformation (or a linear map operation) on
And the way to see this in 3D is a rotation on an XY plane
And notice that
Alright, so to generalize this to 3D, we simply just turn the imaginary number
Note
For more proper technical reason behind the
For a geometric intuition behind it, there's a great video that goes through it which is mainly about spinors but they are "isomorphic" to quaternions by their properties https://youtu.be/ACZC_XEyg9U.
Another geometric intuition behind it is through the reflection operators from Geometric Algebra in reference [16], this one is much easier to comprehend but it doesn't really explain the reasoning behind the double cover very well and is a bit hand-wavy, it however goes on why the general rotation formula is
Covering the entirety of rotations here is going to be tedious. I'd recommend the documents & videos referenced at the bottom of this article in the Reference & Further Reading section, specifically [20] which is a friendly introduction to rotations in general, for more detailed analysis of rotations, then I highly recommend considering the other three references being [21], [22] and [23] as well which delve into the rotation groups ( i.e.
So far we've defined our inertia as the sum of the moments of every particle in a system, however it was a discrete sum and not a continuous sum, usually the bodies we deal with have their masses being continuous like a box or a sphere, so it's better to define it in terms of those and find their inertias.
The total sum for moments of inertia over a continuous volume is defined as
One important thing to notice here is the density function
Which is the mass and the volume ratio, this doesn't really mean anything on it's own to be honest, a better way to understand this is to give an example.
A very good use of the density function is the integral volume
And this is also how bodies in real world weight as much. A cube can be of a mass of 1kg relative to it's size and so on, a large volume for that mass has very low density unlike a low size which by then has a large density. So far we've been using a pure discrete sum assuming that every point mass had an infinite amount of density relative to it's infinitely small size which is definitely not realistic!
So that's pretty much all there is to it, nothing fancy about how the continuous sum is used here.
Now expanding the inertia tensor then becomes
Which is pretty verbose, but that's not important, what's important to notice here is that we have to derive the integral for every single element in the matrix, however we mainly care about the upper triangular part of the tensor and the diagonal elements, the lower triangular part is just the transpose of the upper one, so we could save some effort there.
Let's derive the moment of inertia tensor for a cube. For the first element of the tensor, the integral is
And that's for the first element of the diagonal, the rest is pretty much in the same pattern, what about the upper-triangular part? Let's do one of them
Notice that we just stopped at the first integral, the reason being is that the inner values are summed to basically zero, continuing the summations would give us zeros all the way through. and this pattern repeats for the rest of the off-diagonal elements.
Finishing off the rest of the elements, the final tensor is then
That was for the cube. We could also derive the moment of inertia of a cylinder which is more interesting.
We will be using the cylindrical coordinate system instead of the cartesian one we've used a little earlier which is defined as
Now doing our first integral
That's for the first element of the diagonal row, the second element in the row is similar, the third however has both
Now that was pretty much it, for the off-diagonal elements you could try to integrate them yourself and the values are pretty much zeros. The moment of inertia for a cylinder is then
Suppose the following
It basically says the momentum of the system for both bodies must be conserved which has to be equal before and after the collision, after the collision both were moving at the exact same speeds, so that assumes
Which is coincidentally just the velocity of the center of mass for both the particle and the body combined moving as one! Which could also be interpreted as sticking the particle into the body. Now
Now to make it applicable into both objects, we need to turn it into an impulse, we haven't defined what an impulse is.
An impulse for our case is simply an instantaneous force, it's defined as
The difference between
Turning our momentum equation into an impulse then becomes
Now by Newton's third law, if there's a force, there's an equal and opposite force as well, in our case it's an impulse or an instantaneous force
Which is great! However we haven't taken the conservation of energy so far into account, doing it that way makes the total energy non-conservative, so a way around this is to take the conservation of energy into account as well.
Kinetic energy is defined as
So the conservation of energy is
Which basically says the total energy before the collision must be the same after the collision, continuing
Now from conservation of momentum we know that
Which they would cancel in the energy formula
Re-arranging the terms then becomes
Plugging it back into the momentum formula then becomes
Which now accounts for the conservation of energy law! Turning it into an impulse then becomes
Then again from newton's third law
So we got pretty much everything we want, now however we want to generalize this concept to more than one dimension to also include rotations for 3D, we can use a so called "separating axis", essentially it's the unit/normalized axis that contains the separation between two bodies
How do we make sure they separate only along the axis of separation? It's simply just dotting the difference velocity of two bodies with the separating axis itself! This essentially transfers the problem from an n-dimensional problem into a 1-dimensional problem which makes things consistent
Now we need to re-arrange our previous formula to handle this case like so
Notice in the nominator, there's
Notice that
One other thing, if we simulate an infinitely large mass like a floor, then the other body should completely bounce off the floor, but notice that if
The workaround is to use inverse masses, if the inverse of mass of a particle is infinitely large then it cannot be moved unless there's an infinite force applied to it, essentially it would be approximately zero
So our impulse equation becomes
Try calculating impulses using this equation and see if it's equivalent to our previous one.
Now to make the elasticity variant from in-elastic to completely elastic, we can use a restitution term
Now to include rotations into the impulse equation, we need to treat rotations as linear motions to simplify things a lot, one way to achieve this is by using the following formula
Which is just the velocity of
So what does it look like in practice?
At every point mass their velocities are exactly tangent to the rotating circular manifold surface.
Remember that
Each body has a total velocity like so
And suppose we want to separate the two bodies only along the axis itself for both linear and angular combined when applying impulses
Doing some re-arrangement to our impulse formula using the total velocities then becomes
Check reference [2] for more detailed derivation of this formula if you're interested.
I've already covered the basics like Impulses and such in the Appendix C, so read that first.
Here we will derive the collision impulses based on constraints instead. Constraints for our case are represented as a function
The whole idea is that constraint forces on the system must never do any work, this is the Principle of Virtual Work which is also usually referred to as the d'Alembert's principle, I heavily recommend referring to references [8], [9] and [10] since they explain how (almost) all of this works with constraints, covering them here would be too far.
Now we will derive the impulse using a constraint condition for two bodies contained within a system and derive the constraint forces.
First, the bodies must be separate from each other, assuming that the following condition is satisfied
Where
We are mainly interested in balancing the energies within the system using the condition to separate them from each other.
We simply take the derivative of the constraint function with respect to time and we will get the following
Note that the normal is constant for our case, so it's derivative is zero (some people like to take it's derivative and is probably done for a good reason!) Here assuming the total velocity of a contact point being
Now continuing using that
Here we want to separate the velocities from the normal vectors which is our separation direction, knowing that
Then we have
The mass of the system is defined as
Now assuming we have some scalar that has a symbol as
A simpler way to interpret that is by ignoring the angular velocities for now, that gives us 6 total degrees of freedom to constraint, and we can reduce it even further to 4 degrees of freedom (2 translational degrees for body A and similarly for body B), then
This is our calculated effective mass of the system which is a scalar, assuming that
And the projected energies into the jacobian are
Now solving for
Which is our constraint impulse! It's the exact same formula we derived in Appendix C minus the restitution term.
That was for inelastic impulses, for elastic impulses you have to factor in the restitution term and it should still conserve the momentums and energies just fine.
[1] Euler's Equations of Rigid Body Dynamics Derived | Qualitative Analysis | Build Rigid Body Intuition https://youtu.be/z4-vL6lxR8U
[2] Physically Based Modeling Rigid Body Simulation https://graphics.pixar.com/pbm2001/pdf/notesg.pdf
[3] Rigid Body Motion and the Dzhanibekov Effect https://www.youtube.com/watch?v=NJLdW4DHRcA
[4] The Mystery of Gyroscopic Motion: How Does It Do That? https://www.youtube.com/watch?v=YiQVna7UTiQ
[5] Time Derivative of Rotation Matrices: A Tutorial https://arxiv.org/pdf/1609.06088
[6] Moment of inertia https://en.wikipedia.org/wiki/Moment_of_inertia
[7] Triple Integrals for Volumes of Some Classic Shapes https://sites.math.washington.edu//~aloveles/Math324Fall2013/f13m324TripleIntegralExamples.pdf
[8] A Quick Tutorial on Multibody Dynamics https://fab.cba.mit.edu/classes/865.18/design/optimization/dynamics_1.pdf
[9] Linear-Time Dynamics using Lagrange Multipliers https://www.cs.cmu.edu/~baraff/papers/sig96.pdf
[10] Physically Based Modeling Constrained Dynamics https://graphics.pixar.com/pbm2001/pdf/notesf.pdf
[11] Classical Mechanics Third Edition https://www.math.toronto.edu/khesin/biblio/GoldsteinPooleSafkoClassicalMechanics.pdf
[12] Conservation Of Momentum https://en.wikipedia.org/wiki/Conservation_of_momentum
[13] Elastic collision https://en.wikipedia.org/wiki/Elastic_collision
[14] A COMPLETE DERIVATION OF THE NEWTON-EULER EQUATIONS https://www.panikul.am/assets/Newton_Euler_Derivation.pdf
[15] A Simple And Intuitive Explanation of How Gizmos Work In 3D Space https://gist.github.com/Memresable/4b9b8921b23f86d16b23afea87ae0318
[16] Rotors: A practical introduction for 3D graphics https://jacquesheunis.com/post/rotors/
[17] Lecture L26 - 3D Rigid Body Dynamics: The Inertia Tensor https://ocw.mit.edu/courses/16-07-dynamics-fall-2009/dd277ec654440f4c2b5b07d6c286c3fd_MIT16_07F09_Lec26.pdf
[18] Transformation of axes https://www.doitpoms.ac.uk/tlplib/tensors/transformation.php
[19] Structure and Interpretation of Classical Mechanics https://tgvaughan.github.io/sicm/toc.html
[20] Lie Theory for Roboticists https://www.youtube.com/watch?v=csolG83gCV8
[21] A micro Lie Theory for state estimation in robotics https://arxiv.org/pdf/1812.01537
[22] General Lie Group/Theory resources https://www.ethaneade.com/
[23] Quaternion kinematics for the error-state Kalman filter https://arxiv.org/pdf/1711.02508
[24] SIGGRAPH'22 Course: Contact and Friction Simulation for Computer Graphics https://siggraphcontact.github.io/
[25] M3 Challenge David Baraff Lunch Talk 2016 https://youtu.be/vNTHveVpDDc?feature=shared
[26] Numerical methods for linear complementarity problems in physics-based animation https://dl.acm.org/doi/10.1145/2504435.2504443#supplementary-materials