Rotation about CM - Can I use Moment of Inertia?

[RoshanBBQ]

Rotation about CM -- Can I use Moment of Inertia?

Homework Statement

I have a program I coded to simulate the movement of a some point masses connected to each other rigidly. Each mass is the same. I am trying to code in the correct equations for rotation, but I am having some difficulty with a nonzero product of inertia.

The Attempt at a Solution

I have the x, y, and z coordinates to these masses where x y and z are n by 1 matrices (where I have n masses).

I calculated the center of mass as:
xcm = average(x)
ycm = average(y)
zcm = average(z)

Then I defined new coordinates, xx yy and zz:
xx = x - xcm
yy = y - ycm
zz = z- zcm

The following Wikipedia article then says, "With respect to a coordinate frame whose origin coincides with the body's center of mass, they can be expressed in matrix form as:"

and then it shows that the rotation depends only on the moments of inertia.
http://en.wikipedia.org/wiki/Newton–Euler_equations

However, I then compute the product of inertia:
sum(yy.*zz), which does not equal 0!

where .* means element-wise multiplication.

What is going on here? Am I computing something incorrectly, is Wikipedia wrong, or am I just misunderstanding something?

 

[Andrew Mason]

I calculated the center of mass as:
xcm = average(x)
ycm = average(y)
zcm = average(z)

This is the first problem you have. This is not the definition of the centre of mass. You have to multiply the mass elements by their displacements from an origin, sum them and then divide by the total mass. That gives you the displacement of the centre of mass from the origin.

The second problem you have is with the moment of inertia. The moment of inertia relates the torque to the angular acceleration produced by that torque: $\vec{\tau} = I\vec{\alpha}$

An unconstrained rigid body will rotate about an axis through the centre of mass of the body. You have to then find the moment of inertia of the body about that axis. This can be a difficult calculation and depends on the geometry of the object. See, for example: http://hyperphysics.phy-astr.gsu.edu/hbase/mi.html#mi

You will have to have a solid understanding of rotational motion before you can write your program.

AM

 

[RoshanBBQ]

This is the first problem you have. This is not the definition of the centre of mass. You have to multiply the mass elements by their displacements from an origin, sum them and then divide by the total mass. That gives you the displacement of the centre of mass from the origin.

The second problem you have is with the moment of inertia. The moment of inertia relates the torque to the angular acceleration produced by that torque: $\vec{\tau} = I\vec{\alpha}$

An unconstrained rigid body will rotate about an axis through the centre of mass of the body. You have to then find the moment of inertia of the body about that axis. This can be a difficult calculation and depends on the geometry of the object. See, for example: http://hyperphysics.phy-astr.gsu.edu/hbase/mi.html#mi

You will have to have a solid understanding of rotational motion before you can write your program.

AM

Since all masses are the same, the center of mass is the average of the coordinates.

Let us assume we have n point masses of mass m located on the x-axis with coordinates x_1, x_2, ..., x_n.
$$M = \sum_{k=1}^n m=m\sum_{k=1}^n 1=mn$$
$$x_{cm} = \frac{1}{M}\sum_{k=1}^n x_n m=\frac{m}{M}\sum_{k=1}^n x_n=\frac{m}{nm}\sum_{k=1}^n x_n=\frac{1}{n}\sum_{k=1}^n x_n$$

The result is known as the arithmetic average of the coordinates.

Also, the moment of inertia calculation is not difficult as you claimed it can be since I am dealing with point masses. It is simply the sum of squared coordinates scaled by m.

 

[Andrew Mason]

Since all masses are the same, the center of mass is the average of the coordinates.

Let us assume we have n point masses of mass m located on the x-axis with coordinates x_1, x_2, ..., x_n.
$$M = \sum_{k=1}^n m=m\sum_{k=1}^n 1=mn$$
$$x_{cm} = \frac{1}{M}\sum_{k=1}^n x_n m=\frac{m}{M}\sum_{k=1}^n x_n=\frac{m}{nm}\sum_{k=1}^n x_n=\frac{1}{n}\sum_{k=1}^n x_n$$

The result is known as the arithmetic average of the coordinates.

Ok.

Also, the moment of inertia calculation is not difficult as you claimed it can be since I am dealing with point masses. It is simply the sum of squared coordinates scaled by m.

So would the moment of inertia about the axis of rotation (z axis) not be:

$$m\sum_{i=1}^{n} |x_i-x_{cm}|^2 + |y_i-y_{cm}|^2$$?

AM

 

[D H]

The following Wikipedia article then says, "With respect to a coordinate frame whose origin coincides with the body's center of mass, they can be expressed in matrix form as:"

and then it shows that the rotation depends only on the moments of inertia.
http://en.wikipedia.org/wiki/Newton–Euler_equations

However, I then compute the product of inertia:
sum(yy.*zz), which does not equal 0!

where .* means element-wise multiplication.

What is going on here? Am I computing something incorrectly, is Wikipedia wrong, or am I just misunderstanding something?

That wikipedia article is poorly written, which results in you misunderstanding something.

That first equation in that wikipedia article essentially says that translational and rotational dynamics are decoupled in the special case in which the origin of the frame describing the rigid body has its origin at the body's center of mass. It does not say that "rotation depends only on the moments of inertia." The J_c in that wiki article's equation is the inertia tensor, not the moments of inertia.

To understand this equation, you need to know a bit about what I'll call body frames and structure frames. A structure frame has its origin at some pre-determined point on the rigid body and has some pre-determined set of orthogonal axes that are based on the structure. A body frame is also based on the rigid body, but has its origin at the center of mass of the body.

For example, an aircraft structure frame often has the origin in front of and well below the nose of the aircraft. The x-hat axis points toward the rear of the aircraft along the aircraft's longitudinal axis, the the z-hat axis points up (away from the wheels), and the y-hat axis completes a right-handed coordinate system. This makes the x and z coordinates of any point inside the aircraft positive. A typical aircraft body frame has its origin at the aircraft center of mass, with x-hat and z-hat pointing in the opposite directions as the structural x-hat and z-hat axes. Thus the body frame's x-hat axis points toward the nose, y-hat points out the right wing, and z-hat points down.

These body axes almost certainly are not the principal axes of the aircraft. The principal axes are convenient for simplifying the equations of motion but are quite useless in reality. Think of a plane with sumo wrestlers seated on the left side, a bunch of school kids seated on the right, and a somewhat unbalanced cargo below. The principal axes vary from flight to flight, and even within flight as the plane burns fuel.

You can define a structure frame and a body frame for your set of rigidly connected point masses. Just pick some orthogonal axes that make sense. Working with the center of mass (a body frame) is a lot easier than some other point as the origin (a structure frame) because translational and rotational dynamics decouple with this choice. Translational motion is described by good old $\vec F_{\text{ext}} = d{\vec p}/dt = m \ddot{\vec x}_{\text{cm}}$. Rotational dynamics is a bit hairier, but a lot less hairy than using something other than the center of mass as your origin.For a number of reasons, the angular velocity is typically expressed as angular velocity of the body frame with respect to inertial but expressed in body frame coordinates. One reason is tradition: This goes way back to Euler. Another is that the inertial tensor is a time-varying beast when expressed in inertial coordinates but it is constant (assuming a true constant mass rigid body) in the body frame. Yet another is that sensor readings are in some frame whose orientation is fixed with respect to structure / body.

This representation of angular velocity in body frame coordinates means that the rotational equations of motion are not the simple rotational equivalent of F=dp/dt=ma, $\vec {\tau}_{\text{ext}} =\mathbf J \dot{\vec{\omega}}$, where J is the inertia tensor. The relation F=dp/dt does have a rotational equivalent; it is
$$\vec {\tau}_{\text{ext}} = \left(\frac{d\vec L}{dt}\right)_I$$
where that subscript "I" on the derivative means "from the perspective of an inertial observer."

The time derivative of any vector quantity $\vec u$ from the perspective of an inertial (non-rotating) observer and a rotating observer are related via
$$\left(\frac{d\vec q}{dt}\right)_I = \left(\frac{d\vec q}{dt}\right)_R + \vec{\omega}\times \vec q$$
where the subscript "R" on the derivative on the right hand side means "from the perspective of a rotating observer."

The angular momentum as expressed in body frame coordinates is $\vec L = \mathbf J \vec{\omega}$. Using the rotational equivalent of F=dp/dt and the expression that relates inertial and rotating time derivatives, the rotational equations of motion become
$$\vec {\tau}_{\text{ext,body}} = \mathbf J \frac{d\vec{\omega}}{dt} + \vec{\omega} \times \left(\mathbf J \vec{\omega}\right)$$
Here $\vec {\tau}_{\text{ext,body}}$ are the real external torques on the body expressed in body frame coordinates (rather than inertial).One final note: The above gives the angular velocity in body frame coordinates. It does not give the orientation of the body frame with respect to inertial. You cannot simply integrate angular velocity. Do you need to know the orientation of your body, and if so, do you know how to compute that?

 

[RoshanBBQ]

Ok.
So would the moment of inertia about the axis of rotation (z axis) not be:

$$m\sum_{i=1}^{n} |x_i-x_{cm}|^2 + |y_i-y_{cm}|^2$$?

AM

Yes, that is how I computed it. I used xx yy and zz, which are relative to CM.

That wikipedia article is poorly written, which results in you misunderstanding something.

That first equation in that wikipedia article essentially says that translational and rotational dynamics are decoupled in the special case in which the origin of the frame describing the rigid body has its origin at the body's center of mass. It does not say that "rotation depends only on the moments of inertia." The J_c in that wiki article's equation is the inertia tensor, not the moments of inertia.

To understand this equation, you need to know a bit about what I'll call body frames and structure frames. A structure frame has its origin at some pre-determined point on the rigid body and has some pre-determined set of orthogonal axes that are based on the structure. A body frame is also based on the rigid body, but has its origin at the center of mass of the body.

For example, an aircraft structure frame often has the origin in front of and well below the nose of the aircraft. The x-hat axis points toward the rear of the aircraft along the aircraft's longitudinal axis, the the z-hat axis points up (away from the wheels), and the y-hat axis completes a right-handed coordinate system. This makes the x and z coordinates of any point inside the aircraft positive. A typical aircraft body frame has its origin at the aircraft center of mass, with x-hat and z-hat pointing in the opposite directions as the structural x-hat and z-hat axes. Thus the body frame's x-hat axis points toward the nose, y-hat points out the right wing, and z-hat points down.

These body axes almost certainly are not the principal axes of the aircraft. The principal axes are convenient for simplifying the equations of motion but are quite useless in reality. Think of a plane with sumo wrestlers seated on the left side, a bunch of school kids seated on the right, and a somewhat unbalanced cargo below. The principal axes vary from flight to flight, and even within flight as the plane burns fuel.

You can define a structure frame and a body frame for your set of rigidly connected point masses. Just pick some orthogonal axes that make sense. Working with the center of mass (a body frame) is a lot easier than some other point as the origin (a structure frame) because translational and rotational dynamics decouple with this choice. Translational motion is described by good old $\vec F_{\text{ext}} = d{\vec p}/dt = m \ddot{\vec x}_{\text{cm}}$. Rotational dynamics is a bit hairier, but a lot less hairy than using something other than the center of mass as your origin.For a number of reasons, the angular velocity is typically expressed as angular velocity of the body frame with respect to inertial but expressed in body frame coordinates. One reason is tradition: This goes way back to Euler. Another is that the inertial tensor is a time-varying beast when expressed in inertial coordinates but it is constant (assuming a true constant mass rigid body) in the body frame. Yet another is that sensor readings are in some frame whose orientation is fixed with respect to structure / body.

This representation of angular velocity in body frame coordinates means that the rotational equations of motion are not the simple rotational equivalent of F=dp/dt=ma, $\vec {\tau}_{\text{ext}} =\mathbf J \dot{\vec{\omega}}$, where J is the inertia tensor. The relation F=dp/dt does have a rotational equivalent; it is
$$\vec {\tau}_{\text{ext}} = \left(\frac{d\vec L}{dt}\right)_I$$
where that subscript "I" on the derivative means "from the perspective of an inertial observer."

The time derivative of any vector quantity $\vec u$ from the perspective of an inertial (non-rotating) observer and a rotating observer are related via
$$\left(\frac{d\vec q}{dt}\right)_I = \left(\frac{d\vec q}{dt}\right)_R + \vec{\omega}\times \vec q$$
where the subscript "R" on the derivative on the right hand side means "from the perspective of a rotating observer."

The angular momentum as expressed in body frame coordinates is $\vec L = \mathbf J \vec{\omega}$. Using the rotational equivalent of F=dp/dt and the expression that relates inertial and rotating time derivatives, the rotational equations of motion become
$$\vec {\tau}_{\text{ext,body}} = \mathbf J \frac{d\vec{\omega}}{dt} + \vec{\omega} \times \left(\mathbf J \vec{\omega}\right)$$
Here $\vec {\tau}_{\text{ext,body}}$ are the real external torques on the body expressed in body frame coordinates (rather than inertial).One final note: The above gives the angular velocity in body frame coordinates. It does not give the orientation of the body frame with respect to inertial. You cannot simply integrate angular velocity. Do you need to know the orientation of your body, and if so, do you know how to compute that?

That makes much more sense. Wikipedia calls J "moment of inertia about the center of mass." They should call it the inertia tensor if it is the inertia tensor. I was previously using an inertia tensor in my program when I was using a structure frame (But at the time, I didn't know that meant a coupling between the two types of motion. When I figured it out, I promptly switched to a body frame and threw away my tensor stuff because of this Wikipedia article). You should really edit this Wikipedia article's definition of J if you have the knowledge to do it confidently.

As for your question, I do not think I need to know the angle with respect to the inertial observer. I am modeling a boat with 10 states in my state space: x y z (relative to starting point), their derivatives, roll and yaw, and their derivatives. I am assuming pitch is always zero to simplify the controls problem of directing the boat's motion.

So basically, my code outputs this:

a = tensor\torque.'; % inv(tensor)*torque.'

where a is the instantaneous angular acceleration of pitch, roll, and yaw about the CM. "torque" is a 1 by 3 vector (which is why I transpose it during its interaction with the inertia tensor) defined as
$$\vec {\tau}_{\text{ext,body}} - \vec{\omega} \times \left(\mathbf J \vec{\omega}\right)$$I then have the output go into 2 integrators in a Simulink model. Roll/Yaw and dRoll/Dyaw are fed back into the function for use in computing buoyancy forces etc. as well as torque in the above equation.I have one question for now: Let us assume I have a 3x3 inertia tensor, T, relative to CM for a certain orientation I call "0 degrees from x, y, and z". Is there a way to transform T into T' where T' is the inertia tensor for that same object rotated a, b, and c degrees from x, y, and z respectively? That would save a lot of computation. As of now, I simply recompute the inertia tensor each time (and I have 100,000s of point masses). The code is called 100s of times, so it would be nice to cut off a few ms. As of now, the MATLAB code runs .5 seconds the first time (when it computes a lot of stuff, most burdensome being the x,y, and z coordinates of the hull). Subsequent calls to the function use persistent variables to keep many results, taking only .1 seconds. To take out the tensor calculation would probably cut .03 seconds or so since it is 100,000s of multiplications, sums, and squares. The simulink finally compiles the MATLAB code into C++ code, so it is hard to tell how fast that executes.

PS
Before I made this thread, I tried to do some serious self-studying on the topic. Take a look at that misleading image I uploaded from a book I was using. It seemed all sources agreed you could use only moments of inertia to analyze motion about the CM (wikipedia, books they cite, the book i quote, etc.), but failed to stress you would still need the principal axes for that frame of reference. The book I quote does, however, mention the idea about principal axes. I was left to assume one of two things: I could analyze that motion with that equation about the CM GIVEN I HAVE CHOSEN AXES THAT WERE PRINCIPAL or I could analyze motion about the CM BECAUSE THOSE AXES ARE ALWAYS PRINCIPAL AXES. That is just not a nice ambiguity. You cleared it up well though -- apparently it was the former. I am doing this for an engineering introduction to controls class. I am not well-versed even in 1-d physics, so I have had to do a lot of thinking and self-studying to formulate my 3-d model.

 

[D H]

RoshanBBQ,

I am a bit troubled by this last post of yours. What exactly are you trying to do here? What problem are you trying to solve? You talked in the initial post about rigid bodies, but now you are talking about thousands (hundreds of thousands) of point masses! So which is it?

In a rigid body formalism, you have essentially twelve independent variables that describe each body. With hundreds of thousands of point masses, each with six degrees of freedom (position and velocity), you have on the order of a million independent variables! Whoa! A mesh of hundreds of thousands points would stress even the best finite element analysis package. You might want to rethink what you are doing.

Compared to your million or so degree of freedom system, a rigid body has but twelve degrees of freedom:
- Three for the position of some selected point on the rigid body,
- Three for the time derivative (velocity) of that point,
- Three for the orientation of the body or structure frame, and
- Three for the rotation rate (angular velocity) of the body frame.

There are reasons to use a finite elements approach. Real bodies are deformable rather than rigid. With a rigid body formalism you are forgoing looking at what happens inside the body. This is obviously a disadvantage, but the advantage of a rigid body formalism is that the number of independent variables shrinks by orders and orders of magnitude. There's no reason to use a finite elements approach if the object of interest is fairly rigid and you don't care about behaviors within the body.

The issues that need to be addressed in a rigid body formalism are
- How to physically describe the body,
- How to represent the body's dynamical state,
- How to formulate the equations of motion, and
- How to integrate those equations of motion with respect to time.How to physically describe the body
A point mass has but one physical characteristic, its mass. A rigid body has mass plus an inertia tensor. The inertia tensor of a rigid body is constant in a frame fixed with respect to the body. Assuming the body truly is rigid, you need to calculate the inertia tensor but once.

Here's where you can make use of your hundreds of thousands of point masses. Calculate the inertia tensor for this system of masses. You'll need to transform to the body frame if the frame in which the locations of the point masses are represented is not the body frame. Suppose you calculate the inertia tensor in some frame A, but you want it in frame B. Let T be the transformation matrix that transforms a vector as represented in frame A to the representation of that same vector in frame B. The inertia tensor transforms according to
$$\mathbf I_B = \mathbf T\, \mathbf I_A\,\mathbf T^T$$How to represent state
Rigid bodies have a translational and rotational state. The translational and rotational equations of motion decouple if the translational state is represented as the position and velocity of the center of mass of the rigid body. There are reasons to a structure frame whose origin is not the center of mass, but I recommend that you start with a center of mass body frame. The equations of motion for the translational state such a frame are fairly simple: It's just F_ext=m_bodya. It is rotational state that is the tough nut to crack.

One way to represent rotational state is via the transformation matrix. This is in a sense the canonical form. After all, the mathematics of rotation in Euclidean three space is the group SO(3), where SO is short for special orthogonal (matrices). There is an issue with this representation scheme: Too many numbers. A 3x3 matrix has nine elements, but rotation has but three degrees of freedom. This will lead to challenges when it comes to integrating the equations of motion. There are ways around this. The problem is not intractable.

F=ma is a second order differential equation. The rotational equations of motion are also second order. Some rotational analog of velocity is needed here. The obvious choice is angular velocity. One advantage of using angular velocity is that it is not over-specified. Angular velocity comprises three parameters, and there are three degrees of freedom. For a number of reasons, angular velocity is typically represented as angular velocity of the rotating body with respect to inertial but represented in body frame coordinates.How to formulate the equations of motion
The translational equations of motion are simply F=ma for a center of mass representation. The tricky issue is rotational state. The equations of motion for angular velocity is given in post #5. That leaves rotation itself. If you represent the orientation of your rigid body by the transformation matrix T from inertial to body coordinates, the time derivative of this matrix is
$$\dot{\mathbf T} = -\mathbf X(\vec{\omega}) \mathbf T$$
where $\vec{\omega}$ is the angular velocity as represented in body frame coordinates and $\mathbf X(\vec{\omega})$ is the cross product matrix generated from the vector $\vec{\omega}$. If you instead use the transformation from body to inertial, this becomes
$$\dot{\mathbf T} = \mathbf T \mathbf X(\vec{\omega})$$How to integrate those equations of motion with respect to time
Position, velocity, and angular velocity are conceptually easy. Just use some numerical integration technique. Once again it is rotation that is tricky. Use a typical numerical integration technique on either of the above derivative equations and you will get some matrix that is not orthogonal. One way to get around this is to force the matrix to be orthogonal. Pick some row and normalize (make that row a unit vector). Next, remove the projection of that vector from some second row and normalize. Finally, remove the projections of both of these unit vectors from the third row and normalize. If you use this approach it is best to alternate your choices of this first, second, and third row with each integration step.

Another approach is to use a better integration technique. Google the term "geometric integration". Yet another approach is to use a better representation of orientation. There are many different ways to represent SO(3) that are better suited to numerical integration. Some of these alternative approaches are unit quaternions, Rogrigues parameters, and modified Rogrigues parameters.
With regard to some specific comments in your last post,

That makes much more sense. Wikipedia calls J "moment of inertia about the center of mass." They should call it the inertia tensor if it is the inertia tensor. I was previously using an inertia tensor in my program when I was using a structure frame (But at the time, I didn't know that meant a coupling between the two types of motion. When I figured it out, I promptly switched to a body frame and threw away my tensor stuff because of this Wikipedia article). You should really edit this Wikipedia article's definition of J if you have the knowledge to do it confidently.

Thanks, but no thanks. Wikipedia is a good starting point (sometimes), but too often it is not that good. Too many cooks spoil the broth, and wikipedia has far too many cooks already.

In their defense, they did use a bold J, which indicates a matrix or tensor, and even texts often use moment of inertia when it would be technically more correct to say moment of inertia tensor. Personally, I use a bold J to indicate the moment of inertia tensor about something other than the center of mass and a bold I to indicate the moment of inertia tensor about something the center of mass. What about the identity matrix, which is also denoted by a bold I? I always use a blackboard bold 1 (the LaTeX here isn't up to snuff to show this).

As for your question, I do not think I need to know the angle with respect to the inertial observer. I am modeling a boat with 10 states in my state space: x y z (relative to starting point), their derivatives, roll and yaw, and their derivatives. I am assuming pitch is always zero to simplify the controls problem of directing the boat's motion.

You do need to know this if you want to avoid the million or so independent variables problem.

PS
Before I made this thread, I tried to do some serious self-studying on the topic. Take a look at that misleading image I uploaded from a book I was using. It seemed all sources agreed you could use only moments of inertia to analyze motion about the CM (wikipedia, books they cite, the book i quote, etc.), but failed to stress you would still need the principal axes for that frame of reference. The book I quote does, however, mention the idea about principal axes. I was left to assume one of two things: I could analyze that motion with that equation about the CM GIVEN I HAVE CHOSEN AXES THAT WERE PRINCIPAL or I could analyze motion about the CM BECAUSE THOSE AXES ARE ALWAYS PRINCIPAL AXES. That is just not a nice ambiguity. You cleared it up well though -- apparently it was the former. I am doing this for an engineering introduction to controls class. I am not well-versed even in 1-d physics, so I have had to do a lot of thinking and self-studying to formulate my 3-d model.

This is more of an issue of you trying to teach yourself -- especially from wikipedia. There's something to be said for taking a class, or set of classes. You essentially are trying to learn the field of mechanical engineering on your own. Lacking that, look to some texts. A good sophomore/junior level classical mechanics text is good start, as are texts from an intermediate mechanical engineering course.