How to approximate a triple integral over a rotated grid?

조회 수: 13(최근 30일)
Oscar Engelmark
Oscar Engelmark 2022년 1월 7일
댓글: Oscar Engelmark 2022년 1월 9일
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!
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);

채택된 답변

David Goodmanson
David Goodmanson 2022년 1월 9일
편집: David Goodmanson 2022년 1월 9일
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
  댓글 수: 1
Oscar Engelmark
Oscar Engelmark 2022년 1월 9일
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!

댓글을 달려면 로그인하십시오.

추가 답변(0개)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by