How to approximate a triple integral over a rotated grid?

I'm trying to numerically compute the moment of inertia of a three dimensional solid cylinder of height h and radius , having an arbitrary orientation but with its center of mass at the origin. It basically involves the evaluation of integrals like

over the entire volume enclosed by the cylinder. I've successfully done this using both trapz and a riemann sum for the special case when the cylinder is centered around the z-axis such that I could use cylindrical coordinates and evaluate over the block , , . But I can't seem to figure out a proper method to do it once I have rotated the entire grid. I believe the difficulty arises due to the points no longer being spaced correctly, even in cylindrical coordinates. I'm thinking a riemann sum might be the way to go (and I'm more concerned with speed than accuracy) but I can't figure out how to properly determine the volume elements.

Below is the code I've written (with the cylinder centered around the z-axis). Any ideas on what might be an efficient solution would be appreciated!

clear

h = 2; % cylinder height

r0 = .5; % cylinder radius

values = 100;

[XR,YR,ZR] = cylPts(r0,h,[0 0 1],values); % Rotated grid

F = (XR.^2 + YR.^2).^(3/2); % Function to integrate

% Using a riemann sum

dt = 2*pi/(values+1);

dr = r0/(values+1);

dz = h/(values+1);

dV = dr*dt*dz;

I1 = sum(F,'all')*dV

% Using trapz

r = linspace(0,r0,values);

t = linspace(0,2*pi,values);

z = linspace(-.5*h,.5*h,values);

I2 = trapz(z,trapz(t,trapz(r,F,1),2),3)

function [XR,YR,ZR] = cylPts(R0,L,n,vals)

% Creates a cylinder consisting of vals^3 points with its main axis

% in the direction of n

N = norm(n);

ta = pi/2 + atan2(n(2),n(1)); % azimuth

ti = acos(n(3)/N); % inclination

% rotation matrices

rotZ = [cos(ta) -sin(ta) 0; sin(ta) cos(ta) 0; 0 0 1];

rotX = [1 0 0; 0 cos(ti) -sin(ti); 0 sin(ti) cos(ti)];

% Create arrays

t = linspace(0,2*pi,vals);

z = linspace(-.5*L,.5*L,vals);

r = linspace(0,R0,vals);

% create grids

[R,T,Z] = ndgrid(r,t,z);

X = R.*cos(T);

Y = R.*sin(T);

% rotate grids

tempGrid = [X(:),Y(:),Z(:)]*rotZ*rotX.';

sz = size(X);

% reshape rotated grids

XR = reshape(tempGrid(:,1),sz);

YR = reshape(tempGrid(:,2),sz);

ZR = reshape(tempGrid(:,3),sz);

end

Answer 1

[David Goodmanson]

Hi Oscar,

Before the cyinder is tilted, you are describing the volume with cylindrical coordinates and finding the moment of inertia (moi) with the integrand F = r*r^2, where of course r is the distance of points from the z axis. When you tip the cylinder and keep the rotation axis as still the z axis then F remains the same, nice and simple, but for the tilted cylinder the discription of the boundary surface gets far more complicated. An equivalent and I believe better approach is to leave the cylinder where it its, perfectly described by the coordinates, and go to a tilted axis of rotation. Then F is more complicated than before, but only somewhat.

What's needed is the distance of each point to the new axis of rotation, which is still assumed to go through the point (0,0) the center of mass. A unit vector along the new axis of rotation can be defined as

u = a i + b j + c k

where i j k are unit vectors along in the x y z directions, a,b,c are direction cosines, and a^2 + b^2 + c^2 = 1.

[ the rotation axis could be defined by using spherical coordinates with z as the north pole and angles mu and nu,

u = sin(mu)cos(nu) i + sin(mu)sin(nu) j + cos(mu) k

and a b c are established in this way ].

In cylindrical coordinates, the points in the volume are

p = rcos(theta) i + rsin(theta) j + z k 

For a point p and an axis defined by unit vector u (both based from the origin) the perpendicular displacement from the axis is

v = p - u(p.u)     (dot product). 

d^2, the squared distance, equals v.v, so

d^2 = (p - u(p.u)).(p - u(p.u)) = p^2 - (p.u)^2 

p^2 is of course r^2 + z^2, so

d^2 = r^2 + z^2 - (a rcos(theta) + b rsin(theta) + c z)^2 

The usual way to define the moi is by integrating d^2 times the mass density. The mass density is rho = M / (pi r0^2 h) and the volume element dV is r dr dtheta dz, so the integration is

Integral  rho d^2 r dr dtheta dz

You can use the unrotated ndgrid to define a grid on r, theta,z as you have already done with [R,T,Z] = ndgrid(r,t,z) and use the trapz method from the first part of your code for the integral. Going from 100 to 400 points for increased accuracy, and with M = 1,

for a = 0, b = 0, c = 1 (untilted cylinder)

I2 = 0.1250

which differs from the exact answer (1/2)Mr0^2 by about one part in 10^6.

for a = 1, b = 0, c = 0 (rotation about the x axis)

I2 = 0.3958

which differs from the exact answer (1/4)*M*r0^2 + (1/12)*M*h^2 by about one part in 10^5

[Oscar Engelmark]

This is ingenious, much simpler than what I was trying to do! I will be using this to animate a rotating cylinder in an educational mechanics video I'm working on. Will surely give you credit for helping me out. Thanks a lot!