Moment of Inertia of an Ellipsoid (Arbitrary Axis)

Perform a triple integral over the ellipsoid interior, with the limits of integration dictated by the ellipsoid geometry. Determine the distance from each infinitesimal mass to the axis using the formula at this Wolfram Mathworld link . The formula for $d^2$ in terms of cross products near the end is particularly useful.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105

import math

Num = 1000

M = 1.0

a = 3.0
b = 6.0
c = 5.0

V = (4.0/3.0)*math.pi*a*b*c
rho = M/V

dx = 2.0*a/Num

#####################################

# Two points used in distance calcs

x1 = 0.0
y1 = 0.0
z1 = 0.0

x2 = 1.0
y2 = 4.0
z2 = 3.0

#####################################

I = 0.0

x = -a + dx

while x <= a - dx:

    right = 1.0 - (x/a)**2.0
    ymax = b*math.sqrt(right)

    dy = 2.0*ymax/Num

    y = -ymax + dy

    while y <= ymax - dy:

        right = 1.0 - (x/a)**2.0 - (y/b)**2.0
        zmax = c*math.sqrt(right)

        dz = 2.0*zmax/Num

        dm = rho*dx*dy*dz

        z = -zmax + dz

        while z <= zmax - dz:

            v1x = x2 - x1
            v1y = y2 - y1
            v1z = z2 - z1

            v2x = x1 - x
            v2y = y1 - y
            v2z = z1 - z

            v1sq = v1x**2.0 + v1y**2.0 + v1z**2.0

            crossx = v1y*v2z - v1z*v2y
            crossy = -(v1x*v2z - v1z*v2x)
            crossz = v1x*v2y - v1y*v2x

            dsq = (crossx**2.0 + crossy**2.0 + crossz**2.0)/v1sq

            dI = dm*dsq

            I = I + dI

            z = z + dz

        y = y + dy

    x = x + dx

#####################################

print Num
print I

#>>> 
#100
#7.52908475699

#>>> ================================ RESTART ================================
#>>> 
#200
#7.65190138196

#>>> ================================ RESTART ================================
#>>> 
#400
#7.71171460085

#>>>
#>>> 
#1000
#7.74687109529
#>>> 

Consider a point on the line along the given direction and passing through the origin:

$\vec{r}_L = t\left(\hat{i} + 4 \ \hat{j} + 3 \ \hat{k}\right)$

Any point within the cube is:

$\vec{r}_P =x \ \hat{i} + y \ \hat{j} + z \ \hat{k}$

Vector joining these two points is:

$\vec{r}_D = \vec{r}_P - \vec{r}_L$

Since we consider perpendicular distances from the axis to compute moments of inertia, the vector $\vec{r}_D$ must be perpendicular to the line:

$\implies \vec{r}_D \cdot \left(\hat{i} + 4 \ \hat{j} + 3 \ \hat{k}\right)$

This gives:

$t = \frac{x}{26}+\frac{2\,y}{13}+\frac{3\,z}{26}$

So finally, the vector $\vec{r}_D$ in terms of $x$ , $y$ and $z$ is:

$\vec{r}_D= \left(\frac{2\,y}{13}-\frac{25\,x}{26}+\frac{3\,z}{26}\right) \ \hat{i} + \left(\frac{2\,x}{13}-\frac{5\,y}{13}+\frac{6\,z}{13}\right) \ \hat{j} + \left(\frac{3\,x}{26}+\frac{6\,y}{13}-\frac{17\,z}{26}\right) \ \hat{k}$

A volume mass element is:

$dM = \frac{3}{4 \pi a b c} \ dx \ dy \ dz$ $dM = \frac{1}{120 \pi } \ dx \ dy \ dz$

The moment of inertia is therefore:

$I = \int_{V} \lvert \vec{r}_D \rvert^2 dM$ $I = \frac{1}{120 \pi }\int_{-3}^{3} \int_{-6\sqrt{1-\frac{x^2}{9}}}^{6\sqrt{1-\frac{x^2}{9}}} \int_{-5\sqrt{1-\frac{x^2}{9}-\frac{y^2}{36}}}^{5\sqrt{1-\frac{x^2}{9}-\frac{y^2}{36}}} \lvert \vec{r}_D \rvert^2 \ dx \ dy \ dz$

Integral solved numerically:

$I \approx 7.7794$

Let the position vector of a point in the ellipsoidal solid be $r = [x, y, z]^T$ , and let the unit vector of the direction of the axis be $n = \dfrac{d}{|d|}= \frac{1}{\sqrt{26}} [1, 4, 3]^T$ .

Then the required moment of inertia $I$ is defined by

$I = \rho \displaystyle \int_V | u |^2 dV$

where $\rho$ is the density of the cube, $\rho = \dfrac{m}{V}$ , and $u$ is the vector perpendicular to $n$ from the point $r$ . Using the standard projection formula, we obtain,

$u =\ (I - n n^T ) r$

Hence,

$I = \rho \displaystyle \int_V r^T (I - n n^T)(I - n n^T) r dV = \rho \int_V r^T (I - n n^T) r dV$

Further, this becomes,

$I = \rho \displaystyle \int_V r^T r - r^T n n^T r dV = \rho \int_V r^T r - n^T r r^T n dV = \rho \int_V r^T r dV - \rho n^T \int_V r r^T dV n$

Now $r^T r = x^2 + y^2 + z^2$

and

$r r^T = \begin{bmatrix} x^2 && xy && xz \\ xy && y^2 && yz \\ xz && yz && z^2 \end{bmatrix}$

It can be shown that,

$\displaystyle \int_V x^2 dV = \dfrac{1}{5} V a^2 = I_1$

$\displaystyle \int_V y^2 dV = \dfrac{1}{5} V b^2 = I_2$

$\displaystyle \int_V z^2 dV = \dfrac{1}{5} V c^2 = I_3$

For a derivation of the formula in the last three equations, the reader is referred to my solution of this problem (Courtesy of @Alak Bhattacharya ) . In addition, we have, from symmetry,

$\displaystyle \int_V xy dV = \int_V xz dV = \int_V yz dV = 0$

Thus the moment integral reduces to,

$I = \rho ( I_1 + I_2 + I_3 - n_x^2 I_1 - n_y^2 I_2 - n_z^2 I_3 ) = \frac{1}{5} m ( (1 - n_x^2) a^2 + (1 - n_y^2)b^2 + (1 - n_z^2) c^2 )$

And, therefore, by substituting the given values of $m = 1, a = 3, b = 6, c = 5$ and $n_x^2 = \frac{1}{26} , n_y^2 = \frac{16}{26} , n_z^2 = \frac{9}{26}$ , and simplifying the fraction, we get

$I = \frac{101}{13} \approx \boxed{7.769}$