Skip to content

Instantly share code, notes, and snippets.

@Memresable
Last active September 13, 2024 20:29
Show Gist options
  • Save Memresable/2f8935ec8d61406e8cc758cb60ef8da8 to your computer and use it in GitHub Desktop.
Save Memresable/2f8935ec8d61406e8cc758cb60ef8da8 to your computer and use it in GitHub Desktop.
The Fundamentals Of Rotational Dynamics In 3D Rigid Body Simulation

The Fundamentals Of Rotational Dynamics In 3D Rigid Body Simulation

part_nil_minus_1

In this article, I hopefully clear the most important parts in rotational dynamics related to 3D rigid body simulation and introduce them in an intuitive way, this is mainly for non-physicists (I'm not even a physicist myself, so that's that...), so if you're a physicist or already familiar with moments of inertia well, then I don't think you're going to get anything new out of this, otherwise then this is probably for you!

It's mainly focused on engineers and programmers, so there's not a whole lot of walls of math text and giggles, mostly focused on the very fundamentals and applications of them

The article is also mainly to make sure I have gotten certain concepts right for the most part, so there may be some mistakes here and there, be warned


Now the way to think of a regular rigid body is not as a totally solid shape per se, but as tons of point masses tied together to form a solid shape (box, sphere, cylinder etc.), and they are all constrained to each other almost perfectly that they literally become a solid shape altogether, this mindset is helpful since we want to approximate the orientation of the body more precisely that way, and apply rotations from specific points on a rigid body, we will derive things out of this idea alone so keep that in mind

You may also realize that this means if you want to take a bunch of point masses and constraint their distance to be hard individually towards each other, then you've still gotten a rigid body that way and you can still approximate the orientation of the body in some way, so really there are other ways you could represent a rigid body, however taking a known static shape and deriving things analytically is much faster which is also the common (preferred) way on most physics engines

This is also exactly the reason why most physics engines has the concept of contact points when detecting collisions between bodies, they are helpful in terms of determining which exact point mass on the rigid body got overlapped with the other body so we apply a rotational force from that specific point mass relative to the nearby distribution of the body mass


In physics, cross products don't work the way you would normally think, so things start to get interesting from here

Say you have two vectors A and B, crossing them together (A x B) would give you a vector that's perpendicular to the oriented plane with it's length equal to the parallelogram area between the two vectors, it's orientation would depend on the product order as going from one vector to another, it's not associative, switching the product order to B x A would give an opposite orientation and the normal vector would be on the opposite side of the plane for instance

part_nil_0

This is AB oriented plane and it's orientation order is counter-clockwise like in the image (the orientation goes from A towards B), with it's correspondence purple vector that's normal to the plane, you could also think of this oriented plane as something similar to Euler's 2D rotation formula, where it's really just an exponential complex number e_^ix, the plane AB is the imaginary number and it's area is the real number, combined in a single exponential complex number

part_nil_1

So if you extend this to quaternions for instance (which would be in 3D in this case) this would be like how much you want to orient in three planes instead of one (in which case the imaginary parts are three planes, XY, XZ, YZ, and the real part is saying how much you want to orient on all of these planes simultaneously)

For more info about Euler's formula, you could look up this link, it's excellent: http://www.songho.ca/math/euler/euler.html

So when it comes to physics (or even in general), cross products are just the resultant of an oriented plane spanned by the two vectors and it's area parallelogram being the magnitude orientation on that plane, it's also just an exponential complex number


Angular Velocity is similar to the usual velocity, but it's angular, the usual one would be called linear velocity in our case, angular velocity is the amount of angle in polar form and it's unit is rad/s, it's defined as omega_w = r x v, r is the distance from the "center of mass" of the body to the point mass, v is the point mass velocity, x here is the cross product, I will describe center of mass later

angular_velocity

(Quick note, the omega vector symbol is going to be represented as w for the rest of the article)

Angular Acceleration accelerates angular velocity, and it's unit is rad/s^2, and is defined as alpha_a = r x a, a is point mass acceleration

Angular Momentum is similar from the concept of usual momentum, but this time it's angular, it's defined as L = r x (m * v) or L = r x P, m is mass, P is momentum, this one is going to be interesting for later


Now in case you have a collision from a specific point on a rigid body, there's going to be some sort of torque to accelerate body's angular velocity

part_nil_2

From there we need to apply torque at that specific point, in which case the torque would be something like this

part_nil_3

Torque is basically a rotational force, unlike the normal force which has an origin to rotate around but at infinity which is undefined, it's so far that it essentially becomes a translational force (like say gravity that acts on a body), torque on the other hand has an origin to rotate around, so it makes sense. in which case you could also say both are equivalent forces, but it's better to distinguish them since it's not really helpful to put them both at the same time

Torque (in this case) is defined as T = r x F, where r is a vector pointing from the center of mass of the rigid body to a specific point on the body, now center of mass is quite literally the center of the mass distribution of the rigid body, the box for instance it's center of mass distribution of the body would be right at the center, and in case of non-uniform shapes however that's not always true like so

part_nil_4

Also in the context of masses, you would probably notice that force is F = m * a, now m (mass) here is where makes everything more complicated in the context of rotations

The mass here is the moment of inertia tensor, in case of 2D it's just a scalar, in 3D it's a 3x3 matrix, which seems to explode by just moving up to one higher dimension, why is that?

In 2D, the oriented plane is just XY (relative to the two dimensional orthogonal coordinates which are just x and y), so there's nothing outside of the oriented plane to take into account, which is really easy to deal with

In 3D however, there's now three oriented planes and they are XY, XZ, YZ, and you can imagine the point masses can be outside one of the three oriented planes, which complicates things even further because rotations only happen inside an oriented plane, not outside of it!

part_nil_5

Now in this case we can try to derive moment of inertia for a point mass in a rigid body using the angular momentum equation for our scenario, which is L = m * r x (w x r) (w here is angular velocity)

You might be wondering that it's not actually the case, and angular momentum is L = r x P, well, if you think about it, (w x r) is really just point's velocity direction relative to the center of mass which is constrained to the surface of a circle/sphere in our case like so

part_nil_6

So as the time goes, you would notice the point mass is still moving and constrained to the circle/sphere

The goal is to try and isolate the angular velocity from the products and derive moment of inertia, now the following derivation might look a little daunting, so you could skip to the result if you want

Since r x (w x r) is vector triple product, we can expand it to the following

part_nil_7

And turning them into matrices would look like the following

part_nil_8

part_nil_9

We then take those separated matrices and decrement them

part_nil_10

And the result is the moment of inertia for a single point mass out of many, this is also one of many ways of deriving it, I found this to be the easiest I've seen

Now I don't think this would make intuitive sense to you whatsoever at first, but I'll give my best effort to make this a little more intuitive (don't worry, it's not fully intuitive to me either or to most people in fact)

So for the easier part, in case all the points lie on the oriented plane, then the resultant of the vector triple product r x (w x r) would be in the same direction as the angular velocity, in which case the orientation on that oriented plane would be as you would expect

part_nil_11

However, in case points don't lie on the oriented plane, this is where things get interesting, the resultant vector out of the vector triple product would not give a vector that's parallel to the angular velocity anymore like so

part_nil_12

When you start rotating on that oriented plane, it wouldn't be as stable as you think, as it rotates on that plane it would also start rotating on another plane simultaneously, which is really weird but it makes sense, as you have bigger mass somewhere that's not uniform across the body, then it would start rotating on another plane that's not the one you're primarily rotating on like so

part_nil_13

The sphere here that's attached to the box is mostly the cause for the orientation on the blue part, and orange is mostly caused by combination of the box and sphere

And you could hopefully see that on the final matrix too, the purple parts are the ones right on the planes, and the orange parts are the ones that are outside of the planes that didn't get canceled out by the other sides of the planes, as you notice the orange patterns, there are equivalent parts that appear on two rows at once, as you rotate on one it would resist and start adding some acceleration on the other row

Hopefully this all made sense... and if not, it's still fine I guess, it's important to know what it's purpose for at least

Also, mainly for games that don't rely on realism, it's common to store the moments of inertia as a 3D vector, not a 3x3 matrix

struct Inertia
{
	float yz; // resistance on the yz oriented plane
	float xz; // resistance on the xz oriented plane
	float xy; // resistance on the xy oriented plane
};

It's the common way to avoid weird behaviors caused by moments of inertia on non-uniform mass distributed shapes that are not boxes or spheres, it's less physically accurate but it works, this also saves some memory

And also at times, if you have a sphere or a perfectly uniform box extents, then you can just store the inertia as a single float

Now in case you want to rotate the body, the moment of inertia would be invalid immediately, since it's principle of axis don't match the body's world axis anymore, to fix this you have to do a special kind of transformation that's specialized to tensors which is R x I x inverse(R), R is the body's orientation matrix, I is the moment of inertia, this would update the inertia's axis to the new axis in it's local coordinates, or you could just put the body in it's local coordinates and go on without doing nasty transformations with the tensor if you wish, both ways work, also tensors have their own way of transformation rules that's not similar to what you're already familiar with, my explanation of this is quite vague (and probably not fully correct), I recommend this video series on tensors in case you're curious (and to be clear, the only reason I'm aware of tensors is, well, because of moment of inertia, since I had to)

Now the derivation for moments of inertia analytically is a little complex for known shapes, it does require calculus for integration over the area/volume for certain shapes, but there's already a list on wikipedia for various known shapes, you could also use the above resultant moment of inertia tensor and use that for every point that's inside the shape iteratively and sum them up altogether, you could imagine this could get expensive depending on how small the error is and there are methods to speed this process up

The concept of parallel axis theorem would also be handy, here's a link for the formula, it's essentially offsetting the mass distribution of the inertia tensor over it's principle axis (think of it as "translation" like in homogenous algebra, but not exactly), this helps in case you have a shape that's made up of many shapes tied together for a single body and you want to move those around and sum these tensors together in case you did derive the moments at the origin of each of the shapes axis

Now back to torque, to apply torque on the body, you would need something like this (described in C/C++)

// NOTE: this is an example, it's not the most efficient way of doing this!

struct Rigid_Body
{
	Vector3 position;
	Quaternion orientation;

	Vector3 center_of_mass;

	Vector3 linear_velocity;
	float mass;

	Matrix3x3 inertia; // moment of inertia
	Vector3 angular_velocity;
};

void apply_torque(
	Rigid_Body* body, Vector3 force, Vector3 point,
	float delta_time)
{
	Vector3 world_com = body->position + body->center_of_mass; // world space center of mass
	Vector3 r = point - world_com; // from center of mass to the point
	Vector3 torque = cross(r, force);

	Matrix3x3 orientation_matrix = get_quaternion_matrix(body->orientation);
	Matrix3x3 world_inertia = orientation_matrix * body->inertia * matrix_inverse(orientation_matrix);

	Vector3 J = (matrix_inverse(world_inertia) * torque) * delta_time;
	body->angular_velocity += J;
}

And to update body's orientation every tick, you do the following

void update_body_orientation(Rigid_Body* body, float delta_time)
{
	Vector3 oriented_plane = body->angular_velocity;
	float orientation_magnitude = length(body->angular_velocity) * delta_time;
	Quaternion theta = quaternion_angle_axis(orientation_magnitude, oriented_plane);

	body->orientation = theta * body->orientation; // orientation composition

	body->position = body->position + quaternion_rotate_vector(body->orientation, body->center_of_mass);
}

And that's all there is to it.

Since we're pretty much done here, there are various other resources I would recommend, since this article is supposed to give you a small push into rotational dynamics in general with the fundamentals, it doesn't contain enough applications and examples, but here's one that does, it's also part of the entire Pixar's physically based modeling course, and here's another one that goes much deeper into the rotational dynamics specifically in 3D, which is also part of Chris Hecker's magazine series

That's it, thank you for reading the article.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment