Volume and moment of inertia of a torus

P8.2.3Volume and moment of inertia of a torus

Consider a torus of uniform density, unit mass, average radius $$R$$ and cross sectional radius $$r$$ . The volume and moments of inertia of such a torus may be evaluated analytically and give the results:

\begin{align*} V &= 2\pi^2 Rr^2,\\ I_z &= R^2 + \textstyle \frac{3}{4}r^2,\\ I_x = I_y &= \textstyle \frac{1}{2}R^2 + \frac{5}{8}r^2, \end{align*}

where the centre of mass of the torus is at the origin and the $$z$$ axis is taken to be its symmetry axis.

Here we take a numerical approach. In cylindrical coordinates, $(\rho, \theta, z)$ , it may be shown that:

\begin{align*} V &= 2\int_0^{2\pi} \int_{R-r}^{R+r} \int_0^{\sqrt{r^2 - (\rho-R)^2}} \rho \;\mathrm{d}z\,\mathrm{d}\rho\,\mathrm{d}\theta,\\ I_z &= \frac{2}{V}\int_0^{2\pi} \int_{R-r}^{R+r} \int_0^{\sqrt{r^2 - (\rho-R)^2}} \rho^3 \;\mathrm{d}z\,\mathrm{d}\rho\,\mathrm{d}\theta,\\ I_x = I_y &= \frac{2}{V}\int_0^{2\pi} \int_{R-r}^{R+r} \int_0^{\sqrt{r^2 - (\rho-R)^2}} (\rho^2\sin^2\theta + z^2)\rho \;\mathrm{d}z\,\mathrm{d}\rho\,\mathrm{d}\theta. \end{align*}

Evaluate these integrals for the torus with dimensions $$R=4$$ , $$r=1$$ and compare with the exact values.

Show solution

Here is one solution.

torusmoi

import numpy as np
from scipy.integrate import tplquad

# Torus radius, cross sectional radius
R, r = 4, 1  # R, r

# Integral limits defining the torus in cylindrical coordinates
# theta limits
a, b = 0, 2 * np.pi
# r limits
gfun, hfun = lambda theta: R - r, lambda theta: R + r
# z limits
qfun = lambda theta, rho: 0
rfun = lambda theta, rho: np.sqrt(r**2 - (rho - R) ** 2)

lims = (a, b, gfun, hfun, qfun, rfun)

# Torus volume
V, _ = tplquad(lambda z, rho, theta: rho, *lims)
V *= 2
Vexact = 2 * np.pi**2 * R * r**2
print(f"V = {V} (exact: {Vexact})")

# Integrand for Ix
f = lambda z, rho, theta: ((rho * np.sin(theta)) ** 2 + z**2) * rho
Ix, _ = tplquad(f, *lims)
Ix *= 2 / V
Ix_exact = 0.5 * R**2 + 5 / 8 * r**2
print(f"Ix = {Ix} (exact: {Ix_exact})")

Iz, _ = tplquad(lambda z, rho, theta: rho**3, *lims)
Iz *= 2 / V
Iz_exact = R**2 + 0.75 * r**2
print(f"Iz = {Iz} (exact: {Iz_exact})")

Output:

V = 78.95683520871499 (exact: 78.95683520871486)
Ix = 8.624999999999678 (exact: 8.625)
Iz = 16.749999999999808 (exact: 16.75)