Mass Moment of Inertia of Arbitrary Triangle

In this post I try to explain how the moment of inertia of an arbitrary triangle may be derived (full disclosure: this is based entirely on this excellent video here).

Also it is recommended that before reading this post, the post here which (again, based on the same Youtube channel) covers how the centroid of an arbitrary triangle may be calculated, since the techniques used are essentially the same.

Triangle consisting of 3 points: a, b, and c

Consider a triangle consisting of points a, b and c, as follows. The goal is to calculate the moment of inertia around the axis passing through vertex a (which is also assumed to be the origin), perpendicular to the plane, .


Start with the stanard formula for the mmoi.

\[I=\int r^2\,dm=\rho\int r^2\,dA\]

where:

r is the vector pointing to dA from the origin (= vertex a), and ρ is the density of our "lamina".


r can of course be expressed as being some units (scalar u) along axis u and some units (scalar v) along axis v.

\[\vec r=u\vec u+v\vec v\]

where:

\[0\leqslant u\leqslant 1\]

\[0\leqslant v\leqslant 1\]

As in the calculation of the centroid, we make use of the following relationship.

\[dA=2A\,du\,dv\]

where A is the total area of the triangle.


We can express the mmoi equation as follows.

\[\begin{equation}\begin{split}I &= 2\rho A\int \left| u\vec u+ v\vec v \right|^2du\,dv\\&=2\rho A\int \left(u^2 \left|\vec u  \right|^2 +v^2 \left|\vec v  \right|^2+2uv \vec u \cdot\vec v\right)du \,dv\end{split}\end{equation}\]

As with the centroid calculation, we make use of the Green Theorem to convert the above integral over an area to an integral over a line, to make solving the above equation easier.

Green's Therem

\[{\displaystyle \oint _{C}(L\,dx+M\,dy)=\iint _{D}\left({\frac {\partial M}{\partial x}}-{\frac {\partial L}{\partial y}}\right)dx\,dy}\] 

We will consider the elements in converting the area integral into a line integral.

Taking the first element in the brackets, we can come up with the following by Green Theorem equivalent (assume that M is zero).

\[\begin{equation}\begin{split} I_u =\oint \frac{1}{3}u^3\left| \vec u \right|^2 dv= \int \int u^2 \left|\vec u  \right|^2 du \,dv\end{split}\end{equation}\]

We need to perform the above integral over the line = outline of the triangle anti-clockwise, ie vertex a -> vertex b -> vertex c -> vertex a.

For a->b, v = 0 and dv = 0, hence the line integral is zero.

For c->a, u = 0, and du = 0, hence the line integral is zero.

So we only are left with vertex b = vertex c. The line from vertex b to vertex c can be expressed as the following (just think of the equation of a line joining (1,0) and (0,1) in a normal cartesian coordinate system).

\[u+v=1\]

It follows that \(du = -dv\). Hence:

\[\begin{equation}\begin{split} I_u &= \frac{1}{3}\left| \vec u \right|^2 \int _0 ^1 (1-v)^3 dv\\&=\frac{1}{3}\left| \vec u \right|^2\left[ \frac{-(1-v)^4}{4} \right]_0 ^1&=\frac{1}{12}\left| \vec u \right|^2\end{split}\end{equation}\]

Similarly, the second element would be:

\[\begin{equation}\begin{split} I_v &= \frac{1}{12}\left| \vec v \right|^2\end{split}\end{equation}\]

For the third element below, we turn to the Green's theorem to convert the area integral into a line integral.

\[\begin{equation}\begin{split} I_{uv} = \int \int 2uv \vec u \cdot \vec v\, du \,dv\end{split}\end{equation}\]

Referring back to the Green's theorem, assume that L = 0 and partial derivative of M against u is the 2uv, arrive at the following:

\[\begin{equation}\begin{split} \oint u^2v \,\vec u \cdot \vec v\,du= \int \int 2uv\, \vec u \cdot \vec v\, du \,dv\end{split}\end{equation}\]

Now we need to perform the line integral against the outline of the triangle. Similar to the other elements, paths a->b and c->a turns out to be zero, so we only need to worry about path b->c (which is expressed as u+v = 1).

\[\begin{equation}\begin{split} I_{uv}&=\oint u^2v\vec u\cdot\vec v\, dv\\&=\vec u\cdot\vec v\int_0^1v(1-v)^2\,dv\\&=\vec u\cdot\vec v\int_0^1v(1+v^2-2v)\,dv\\&=\vec u\cdot\vec v\left[\frac{1}{2}v^2+\frac{1}{4}v^4-\frac{2}{3}v^3  \right]_0^1\\&=\frac{1}{12}\vec u\cdot\vec v\end{split}\end{equation}\]

Now we can add the elements together to arrive at the equation for the mmoi.

\[\begin{equation}\begin{split} I&=2\rho A(\frac{1}{12}\left| \vec u \right|^2+\frac{1}{12}\left|\vec v \right|^2+\frac{1}{12}\vec u \cdot \vec v)\\&=\frac{1}{6}M(\left| \vec u \right|^2+\left|\vec v \right|^2 +\vec u \cdot \vec v)\\&=\frac{1}{6}M(\vec u \cdot \vec u +\vec v \cdot \vec v + \vec u \cdot \vec v)\end{split}\end{equation}\]

where M is the product of the density and the area, ie mass of the triangular lamina.


Useful References

https://youtu.be/uFDTxRv7O8M 

https://physics.stackexchange.com/questions/670439/moment-of-inertia-of-a-triangular-lamina-with-vertices-at-the-origin-p-and-q 


Mass Momentum of Inertia of Capsule

To calculate the mmoi of a capsule, we consider it as the sum of a rectangle and two semi-circles. Then we can add the mmoi of each of the components and add them together.

We already know the mmoi of a rectangle, around the z-axis through its centre of mass.

\[I=\frac{1}{12}M(w^2+h^2\]

But what of a semi-circle? It is in fact, half the mmoi of a full circle - if this is not obvious, go back to how we derived the mmoi of a full circle in the previous circle, consider a semi-circle as consisting of very thin semi-circular rings. The "length" of a semi-circular ring is obviously half the length of a full "hoop". The rest of the derivation is the same - hence the mmoi is:

\[I_z=\frac{1}{4}MR^2\]

Then we can use the parallel axis theorem to add the mmoi of the elements together to arrive at the I of the total.

Parallel Axis Theorem

Suppose a body of mass M is rotated about an axis z passing through the body's center of mass. The body has a moment of inertia Icm with respect to this axis. The parallel axis theorem states that if the body is made to rotate instead about a new axis z′, which is parallel to the first axis and displaced from it by a distance d, then the moment of inertia I with respect to the new axis is related to Icm by:

\[I=I_{cm}+Md^2\]

where, d is the perpendicular distance between the axes z and z′.

Hence:

\[\begin{equation}\begin{split} I &=\frac{1}{12}M(w^2+h^2)+\left(\frac{1}{4}MR^2 \right)(\frac{w}{2})^2+\left(\frac{1}{4}MR^2 \right)(\frac{w}{2})^2\\&=\frac{1}{12}M(w^2+h^2)+ \frac{1}{2}MR^2 \left(\frac{w}{2}\right)^2\end{split}\end{equation}\]

Mass Moment of Inertia: Overview

In coding collision response, we have made use of the concept of impulse, which is the change in momentum, and momentum being the product of velocity and mass.

To date, we had not delved any further than the above, but before we go into the world of rotation, we need to recognise that we were actually referring to linear momentum or translational momentum. When bodies are rotating, there is a rotational equivalent of momentum which is the product of angular velocity and something called the mass moment of inertia (MMOI).

Making sure we're talking about the right thing

Before we start exploring the concept of MMOI, it is helpful to list out a number of similar founding but very much different concepts.

First Moment of Area

  • Also commonly referred to as First Moment of Inertia
  • Measure of the distribution of an area relative to an axis. Can be used to find the centroid
  • Unit: m3.

Second Moment of Area

  • Also commonly referred to as:
    • Moment of Inertia of an Area
    • Area Moment of Inertia
    • Second Area Moment
  • Measure of the resistance to torsion
  • Unit: m4

Mass Moment of Inertia

  • Also commonly referred to as:
    • Moment of Inertia (without the Mass)
    • Angular Mass,
    • Second Moment of Mass
    • Rotational Inertia,
  • Measures resistance to rotational motion and angular (rotational) acceleration due to an applied torque (moment).
  • Unit: kg m2.

It is the third one that we are dealing with in this post.


What is the Mass Moment of Inertia?

The moment of inertia I is also defined as the ratio of the net angular momentum L of a system to its angular velocity ω around a principal axis, that is:

\[I=\frac{L}{\omega}\]

Unlike mass, which depends only on amount of matter, moment of inertia is, in addition to the amount of matter, also dependent on the position of the axis of rotation and the shape of the matter. In otherwords, even when working exclusively in 2D, the moment of inertia for a certain shape would be different depending on the particular axis of rotation.

Calculating the Mass Moment of Inertia

So how do we calculate the MMOI of our physics body? Wiki explains as follows.

The quantity that determines the torque needed for a desired angular acceleration about a rotational axis, akin to how mass determines the force needed for a desired acceleration. It depends on the body's mass distribution and the axis chosen, with larger moments requiring more torque to change the body's rate of rotation.

It is an extensive (additive) property: for a point mass the moment of inertia is simply the mass times the square of the perpendicular distance to the axis of rotation. The moment of inertia of a rigid composite system is the sum of the moments of inertia of its component subsystems (all taken about the same axis). Its simplest definition is the second moment of mass with respect to distance from an axis.

For bodies constrained to rotate in a plane, only their moment of inertia about an axis perpendicular to the plane, a scalar value, matters. 

In other words, restricting our discussions to 2D, moment of inertia I, for a point of mass m at distance r from the pivot point is simply:

\[I=mr^2\]

We can extend this concept to a body with shape by the fact that the mass moment of inertia of a rigid body is the sum of the moments of inertia of all the points making up that rigid body.

More specifically, you can consider a body (in our 2D world, a lamina with a shape) as consisting of lots of small pieces that has mass dm. To find the mmoi of the entire lamina shape, we find equation for individual "elements" and integrate over the entire shape.

Usefully, the Wikipedia page below lists the formulae for various shapes.

https://en.wikipedia.org/wiki/List_of_moments_of_inertia

If you look closely, you will find that for some of the shapes, different formulae are given for different axis of rotation. It is critical to understand that unlike mass in linear motion, there is not one mmoi figure for a shape - the mmoi for a shape will be different depending on the axis of rotation.

I will delve into the calculation of the mmoi for some of the simpler shapes just to confirm my understanding.

MMOI of a hoop

We start with a thin circular loop of radius r and total mass M.


You can think of the loop as consisting of lots of points of mass dm. Every point would be at distance r, and hence the I of the individual point would be dm x squared(r). The mmoi of the loop would be the sum of all those individual points, hence intuitively you can guess that the mmoi of the loop is M x squared(r).

Let's be more rigorous about the way we calculate this.

Again, think of the "loop" as consisting of lots of points, each with mass dm, distance r from the centre. Now, assume that each "point" has a "length" of dl, and the total "length" of the hoop is L. We can say that:

\[\frac{dm}{M}=\frac{dl}{L}\]

By definition,

\[dl=r\,d\theta\]

which means that:

\[dm=M\times \frac{d\theta}{2\pi}\]

In order to find the mmoi for the loop, we need to integrate over the entire loop from limits of 0 to 2π, as follows.

\[I=\int_0^{2\pi}dm\,r^2=\int_0^{2\pi}M\frac{d\theta}{2\pi}r^2=\frac{Mr^2}{2\pi}\int_0^{2\pi}d\theta=\frac{Mr^2}{2\pi}\times 2\pi=Mr^2\]

MMOI of a "Rod"

The rod is of length L and total mass M. In order to calculate the rod's mmoi around the axis which passes through its centre, assume that the rod is made up of elements of length dl,  mass of dm.

As before, the relation between the total mass M and the total length L vs the elemental mass dm and elemental length dl is as follows:

\[\frac{M}{L}=\frac{dm}{dl}\]

It follows that:

\[md=\frac{M}{L}dl\]

Plug the above into the mmoi integral:

\[\begin{equation}\begin{split} I&=\int_{-L/2}^{L/2}l^2\frac{M}{L}\,dl\\&=\frac{M}{L}\int_{-L/2}^{L/2}l^2\,dl\\&=\frac{M}{L}\left[ \frac{l^3}{3} \right]_{-L/2}^{L/2}\\&=\frac{M}{3L}\times\frac{2L^3}{8}\\&=\frac{1}{12}ML^2\end{split}\end{equation}\]

MMOI of a Disc

In the first example, we talked about a point element with elemental mass of dm. In the second example we talked about an element of length dl and mass of dm.

In dealing with shapes with area, we introduce the concept of elemental area dA which has elemental mass of dm. The mmoi of the shape is the sum of the mmoi of the elemental area.

Thinking about the disc as concentric loops of "width" dx, distance r from the original. 

If we denote the area of the circle by A, then we can say that:

\[\frac{dm}{M}=\frac{dA}{A}=\frac{dA}{\pi R^2}\]

We can also say that dA, which is the area of the infinitesimally thin hoop, is the product of the circumference and dx.

\[dA=2\pi R\,dx\]

It follows that:

\[dm=M\frac{dA}{\pi R^2}=M\frac{2\pi xdx}{\pi R^2}=M\frac{2xdx}{R^2}\]

Now we can stick the above into the formula for MMOI.

\[\begin{equation}\begin{split}I&=\int_0^Rx^2\frac{M2x\,dx}{R^2}\\&=\frac{2M}{R^2}\int_0^Rx^3dx\\&=\frac{2M}{R^2}\left[ \frac{x^4}{4} \right]_0^R\\&=\frac{1}{2}MR^2\end{split}\end{equation}\]



Area of Arbitrary Polygon

In this post I look at various ways of calculating the area of an arbitrary polygon with n-vertices.

  1. By Integral
  2. Trapezoid formula
  3. Triangle formula
  4. Shoelace formula



1. Calculating Area by Integral

\[\int _A dA=\iint dx\,dy\]

The edges of the polygon may be expressed as follows.

\[\begin{cases}x = x_i + (x_{i+}-x_i)\,t \\ y = y_i + (y_{i+}-y_i)\,t \end{cases}\quad \mbox{with} \quad \begin{cases} i = 0,1,2,\cdots,n-1 \\ i+=i+1\mod n \end{cases}\quad \mbox{and} \quad 0 \le t \le 1\]

We can convert the area integral into a line integral by employing Green's theorem, which states that:

Let C be a positively oriented, piecewise smooth, simple closed curve in a plane, and let D be the region bounded by C. If L and M are functions of (x, y) defined on an open region containing D and have continuous partial derivatives there, then

\[\iint \left( \frac{\partial M}{\partial x} - \frac{\partial L}{\partial y} \right) dx\,dy= \oint \left( L\,dx + M\,dy \right)\]

By substitution of \(M(x,y)=x\) and \(L(x,y)=0\) we can restate the area equation as follows:

\[A = \iint dx\,dy \\= \oint x\,dy\]

We have defined the y-coordinate of each edge as:

\[y = y_i + (y_{i+}-y_i)\,t\]

Based on which we can say that:

\[\frac{dy}{dt} = (y_{i+}-y_i)\]

It follows that:

\[\begin{equation}\begin{split} A &=  \oint x\,dy \\&= \sum_{i=0}^{n-1} \int_0^1 \left[x_i + (x_{i+}-x_i)\,t\right](y_{i+}-y_i)\,dt\\&=\sum_{i=0}^{n-1}(y_{i+}-y_i)\left[x_i\left.t\right|_0^1 + (x_{i+}-x_i)\frac{1}{2}\left.t^2\right|_0^1\right]\\&=\frac{1}{2}\sum_{i=0}^{n-1}(y_{i+}-y_i) (2x_i+x_{i+}-x_i )\\& =\frac{1}{2}\sum_{i=0}^{n-1}(x_{i+}+x_i)(y_{i+}-y_i)\end{split}\end{equation}\]

The last line when expanded looks like below:

\[\begin{equation}\begin{split}y_1x_1+y_1x_0-y_0x_1-y_0x_0\\+y_2x_2+y_2x_1-y_1x_2-y_1x_1\\+y_3x_3+y_3x_2-y_2x_3-y_2x_2\\\cdot\,\cdot\,\cdot\,\cdot\,\\+y_n(x_n+x_{n-1})-y_{n-1}(x_n+x_{x-1})\end{split}\end{equation}\]

Remembering that \(x_n=x_0\) and \(y_n=y_0\), the the area of a polygon may be expressed as: 

\[A =\frac{1}{2}\sum_{i=0}^{n-1} (x_iy_{i+}-x_{i+}y_i)\]

The formula above is consistent with the following three calculation methods, which are intuitively easier to comprehend. I have simply presented the formulas, straight from the Wikipedia here, which explains in detail.

Note that, in keeping with the original Wiki post, the sequence of points are numbered from i=1 to i=n, as opposed to i=0 to i=n-1, ie \({\displaystyle P_{i}=(x_{i},y_{i}),i=1,\dots ,n}{\displaystyle P_{i}=(x_{i},y_{i}),i=1,\dots ,n}\).

Also, for the simplicity it is assumed that: \({\displaystyle P_{0}=P_{n},P_{n+1}=P_{1}}\).

2. Trapezoid formula

\[{\displaystyle {\begin{aligned}A&={\frac {1}{2}}\sum _{i=1}^{n}(y_{i}+y_{i+1})(x_{i}-x_{i+1})\\&={\frac {1}{2}}{\Big (}(y_{1}+y_{2})(x_{1}-x_{2})+\cdots +(y_{n}+y_{1})(x_{n}-x_{1}){\Big )}\end{aligned}}}\]

3. Triangle fomula

\[{\displaystyle {\begin{aligned}A&={\frac {1}{2}}\sum _{i=1}^{n}(x_{i}y_{i+1}-x_{i+1}y_{i})={\frac {1}{2}}\sum _{i=1}^{n}{\begin{vmatrix}x_{i}&x_{i+1}\\y_{i}&y_{i+1}\end{vmatrix}}={\frac {1}{2}}\sum _{i=1}^{n}{\begin{vmatrix}x_{i}&y_{i}\\x_{i+1}&y_{i+1}\end{vmatrix}}\\&={\frac {1}{2}}{\Big (}x_{1}y_{2}-x_{2}y_{1}+x_{2}y_{3}-x_{3}y_{2}+\cdots +x_{n}y_{1}-x_{1}y_{n}{\Big )}\end{aligned}}}\]

4. Shoelace formula

\[{\displaystyle {\begin{aligned}2A&=\underbrace {{\begin{vmatrix}x_{1}&x_{2}\\y_{1}&y_{2}\end{vmatrix}}+{\begin{vmatrix}x_{2}&x_{3}\\y_{2}&y_{3}\end{vmatrix}}+\cdots +{\begin{vmatrix}x_{n}&x_{1}\\y_{n}&y_{1}\end{vmatrix}}} \\&={\begin{vmatrix}x_{1}&x_{2}&x_{3}\cdots &x_{n}&x_{1}\\y_{1}&y_{2}&y_{3}\cdots &y_{n}&y_{1}\end{vmatrix}}&{\text{(horizontal shoelace scheme)}}\\&={\begin{vmatrix}x_{1}&y_{1}\\x_{2}&y_{2}\\\vdots &\vdots \\x_{n}&y_{n}\\x_{1}&y_{1}\end{vmatrix}}&{\text{(vertical form)}}\end{aligned}}}\]


Useful References

geometry - Why doesn't a simple mean give the position of a centroid in a polygon? - Mathematics Stack Exchange

Shoelace formula - Wikipedia


Centroid of Arbitrary Triangle


In the previous post we worked through the calculation of the centroid of an isosceles triangle. In this post we delve into the calculation of the centroid of an arbitrary triangle.

Not this post is pretty much a "slower" version of the content of this great video, here.

Triangle consisting of 3 points: a, b, and c


Consider a triangle consisting of points a, b and c.


Imagine that the points are defined with reference to a special frame of reference local to this triangle, with u-axis, and v-axis, which are in fact parallel with a->b and a->c.

Imagine also that the origin is a, and b is 1 unit along the u-axis + 0 along the v-axis, and c is at 0 along the u-axis and 1 along the v-axis.

Now, think of the triangle as being comprised of very small elements each of area dA.


The centroid of this triangle would be given by:
\[\vec r _o = \frac{1}{A}\int \vec r \, dA\]
where:
r is the vector pointing to dA from the origin (which we have taken to be vertex a)
A is the total area of the triangle


r can of course be expressed as being some units (scalar u) along axis u and some units (scalar v) along axis v.
\[\vec r=u\vec u+v \vec v\]
where: \[0\leqslant u \leqslant 1\] and \[0\leqslant v \leqslant 1\]

How can we define dA?


Imagine the parallelogram as follows. You will note that dA would be directly proportional to the area of the parallelogram, which is 2A (twice the area of the triangle)

Hence, we can say that: \[dA=2A\,du\,dv\]
Substituting this into the centroid formula above then gives us: \[\vec r_o = \iint  2\vec r\,du\,dv\]
Replacing r with a column vector of u and v then gives: \[\vec r_o=\iint  \begin{pmatrix} 2\vec u \\ 2\vec v \end{pmatrix} du\,dv\]
Then something called Green's theorem can be used to transform the above integral over an area (i.e. over the area of the triangle) to an integral over a line (specifically, the outline of the triangle - from vertex a to vertex b to vertex c), which is slightly easier to solve. 

Green's Therem


Wiki says: Let C be a positively oriented, piecewise smooth, simple closed curve in a plane, and let D be the region bounded by C. If L and M are functions of (x, y) defined on an open region containing D and have continuous partial derivatives there, then \[{\displaystyle \oint _{C}(L\,dx+M\,dy)=\iint _{D}\left({\frac {\partial M}{\partial x}}-{\frac {\partial L}{\partial y}}\right)dx\,dy}\] What we want to do is to convert our centroid integral from the right part of the Green's theorem to the left part of the Green's theorem (think of the x as our u, and y as our v), like so: \[\oint_C\begin{pmatrix} \vec u^2\,dv \\ \vec v^2\,du \end{pmatrix}= \int \int_D\begin{pmatrix} 2\vec u \\ 2\vec v \end{pmatrix} du\,dv\] We are now required to perform an integral along the path a->b->c. Thinking through this element by element: Path a->b: We are moving along axis u with v = 0. By definition dv is also 0. Hence the bit inside the bracket in the above integral is zero. Similarly, Path a->c (or strictly speaking, path c->a) is zero for similar reasons. So, we just need to think about the integral along path b->c. We noted earlier that any point on the triangle can be expressed as a linear function of u and v where the coefficients u+v = 1. It follows that du = -dv. Hence: \[\vec {r}_o = \oint _{C}\begin{pmatrix} \vec {u}^2\,dv \\ -\vec {v}^2\, du \end{pmatrix} =\int _{(1,0)} ^{(0,1)}\begin{pmatrix} -\vec {u}^2 \, du\\ \vec {v}_2\, dv \end{pmatrix} \] We can make this even simpler by flipping the u-axis "upside down" and getting rid of the minus in the column matrix. \[\vec {r}_o =\int _{(1,0)} ^{(0,1)}\begin{pmatrix} -\vec {u}^2 \, du\\ \vec {v}_2\, dv \end{pmatrix} =\int _{(0,0)} ^{(1,1)}\begin{pmatrix} \vec {u}^2 \, du\\ \vec {v}_2\, dv \end{pmatrix}=\begin{bmatrix} \frac{1}{3}\vec u \\ \frac{1}{3}\vec v \end{bmatrix} =\frac{\vec u + \vec v}{3} \] Geometrically, u+v is the vertex at the "other end" of the parallelogram, as below. Hence the above formula is saying that the centroid is a third of the way along the vector shown.



 By replacing u and v by the vertices, adding vertex a (which has been assumed to be the origin in the above calculation) gives us: \[\vec {r}_o = \frac{\vec u + \vec v}{3}+\vec a = \frac{(\vec b - \vec a)+(\vec c - \vec a)}{3} +\vec a =\frac{\vec a+\vec b+\vec c }{3} \]
A remarkably simple result.

Useful References


https://youtu.be/X7JMIGqXBBA

Numerical Integration: 2nd and 4th order Runge-Kutta schemes

Second-order Runge-Kutta scheme (RK2) There are several variations of the second order Runge-Kutta schemes of numerical integration. They ar...