计算凹二维多边形相对于其原点的惯性矩

Question

我想计算（2D）凹多边形的转动惯量。我在互联网上找到了这个。但我不太确定如何解释公式...

公式 http://img101.imageshack.us/img101/8141/92175941c14cadeeb956d8f.gif

1) 这个公式正确吗？
2）如果是这样，我到 C++ 的转换是否正确？

float sum (0);
for (int i = 0; i < N; i++) // N = number of vertices
{
    int j = (i + 1) % N;
    sum += (p[j].y - p[i].y) * (p[j].x + p[i].x) * (pow(p[j].x, 2) + pow(p[i].x, 2)) - (p[j].x - p[i].x) * (p[j].y + p[i].y) * (pow(p[j].y, 2) + pow(p[i].y, 2));
}
float inertia = (1.f / 12.f * sum) * density;

马丁

Answer 1

我认为除了将公式转换为代码之外，您还有更多工作要做。您需要准确理解这个公式的含义。

当您有一个 2D 多边形时，您可以相对于给定坐标系计算三个惯性矩：关于 x 的力矩、关于 y 的力矩和极惯性矩。平行轴定理允许您从一个坐标系转换到另一个坐标系。

你知道这个公式适用于哪个力矩和坐标系吗？

下面是一些可能对您有帮助的代码，以及证明其有效的 JUnit 测试：

import java.awt.geom.Point2D;

/**
 * PolygonInertiaCalculator
 * User: Michael
 * Date: Jul 25, 2010
 * Time: 9:51:47 AM
 */
public class PolygonInertiaCalculator
{
    private static final int MIN_POINTS = 2;

    public static double dot(Point2D u, Point2D v)
    {
        return u.getX()*v.getX() + u.getY()*v.getY();
    }

    public static double cross(Point2D u, Point2D v)
    {
        return u.getX()*v.getY() - u.getY()*v.getX();
    }

    /**
     * Calculate moment of inertia about x-axis
     * @param poly of 2D points defining a closed polygon
     * @return moment of inertia about x-axis
     */
    public static double ix(Point2D [] poly)
    {
        double ix = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double sum = 0.0;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY();
                sum += (poly[n].getY()*poly[n].getY() + poly[n].getY()*poly[n+1].getY() + poly[n+1].getY()*poly[n+1].getY())*twiceArea;
            }

            ix = sum/12.0;
        }

        return ix;
    }

    /**
     * Calculate moment of inertia about y-axis
     * @param poly of 2D points defining a closed polygon
     * @return moment of inertia about y-axis
     * @link http://en.wikipedia.org/wiki/Second_moment_of_area
     */
    public static double iy(Point2D [] poly)
    {
        double iy = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double sum = 0.0;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY();
                sum += (poly[n].getX()*poly[n].getX() + poly[n].getX()*poly[n+1].getX() + poly[n+1].getX()*poly[n+1].getX())*twiceArea;
            }

            iy = sum/12.0;
        }

        return iy;
    }

    /**
     * Calculate polar moment of inertia xy
     * @param poly of 2D points defining a closed polygon
     * @return polar moment of inertia xy
     * @link http://en.wikipedia.org/wiki/Second_moment_of_area
     */
    public static double ixy(Point2D [] poly)
    {
        double ixy = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double sum = 0.0;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY();
                sum += (poly[n].getX()*poly[n+1].getY() + 2.0*poly[n].getX()*poly[n].getY() + 2.0*poly[n+1].getX()*poly[n+1].getY() + poly[n+1].getX()*poly[n].getY())*twiceArea;
            }

            ixy = sum/24.0;
        }

        return ixy;
    }

    /**
     * Calculate the moment of inertia of a 2D concave polygon
     * @param poly array of 2D points defining the perimeter of the polygon
     * @return moment of inertia
     * @link http://www.physicsforums.com/showthread.php?t=43071
     * @link http://www.physicsforums.com/showthread.php?t=25293
     * @link http://stackoverflow.com/questions/3329383/help-me-with-converting-latex-formula-to-code
     */
    public static double inertia(Point2D[] poly)
    {
        double inertia = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double numer = 0.0;
            double denom = 0.0;
            double scale;
            double mag;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                scale = dot(poly[n + 1], poly[n + 1]) + dot(poly[n + 1], poly[n]) + dot(poly[n], poly[n]);
                mag = Math.sqrt(cross(poly[n], poly[n+1]));
                numer += mag * scale;
                denom += mag;
            }

            inertia = numer / denom / 6.0;
        }

        return inertia;
    }
}

这是随附的 JUnit 测试：

import org.junit.Test;

import java.awt.geom.Point2D;

import static org.junit.Assert.assertEquals;

/**
 * PolygonInertiaCalculatorTest
 * User: Michael
 * Date: Jul 25, 2010
 * Time: 10:16:04 AM
 */
public class PolygonInertiaCalculatorTest
{
    @Test
    public void testTriangle()
    {
        Point2D[] poly =
        {
            new Point2D.Double(0.0, 0.0),
            new Point2D.Double(1.0, 0.0),
            new Point2D.Double(0.0, 1.0)
        };


        // Moment of inertia about the y1 axis
        // http://www.efunda.com/math/areas/triangle.cfm
        double expected = 1.0/3.0;
        double actual = PolygonInertiaCalculator.inertia(poly);

        assertEquals(expected, actual, 1.0e-6);
    }

    @Test
    public void testSquare()
    {
        Point2D[] poly =
        {
            new Point2D.Double(0.0, 0.0),
            new Point2D.Double(1.0, 0.0),
            new Point2D.Double(1.0, 1.0),
            new Point2D.Double(0.0, 1.0)
        };

        // Polar moment of inertia about z axis
        // http://www.efunda.com/math/areas/Rectangle.cfm
        double expected = 2.0/3.0;
        double actual = PolygonInertiaCalculator.inertia(poly);

        assertEquals(expected, actual, 1.0e-6);
    }

    @Test
    public void testRectangle()
    {
        // This gives the moment of inertia about the y axis for a coordinate system
        // through the centroid of the rectangle
        Point2D[] poly =
        {
            new Point2D.Double(0.0, 0.0),
            new Point2D.Double(4.0, 0.0),
            new Point2D.Double(4.0, 1.0),
            new Point2D.Double(0.0, 1.0)
        };

        double expected = 5.0 + 2.0/3.0;
        double actual = PolygonInertiaCalculator.inertia(poly);

        assertEquals(expected, actual, 1.0e-6);

        double ix = PolygonInertiaCalculator.ix(poly);
        double iy = PolygonInertiaCalculator.iy(poly);
        double ixy = PolygonInertiaCalculator.ixy(poly);

        assertEquals(ix, (1.0 + 1.0/3.0), 1.0e-6);
        assertEquals(iy, (21.0 + 1.0/3.0), 1.0e-6);
        assertEquals(ixy, 4.0, 1.0e-6);
    }
}

Answer 2

作为参考，这里有一个可变的 2D org.gcs.kinetic.Vector< /code>实现和更通用、不可变的 org.jscience.mathematics.vector 实现。本文介绍计算 2D 向量的叉积也很有帮助。

Answer 3

我用曲面细分做到了。并将 MOI 全部放在一起。

Answer 4

#include <math.h> //for abs
float dot (vec a, vec b) {
   return (a.x*b.x + a.y*b.y);
}
float lengthcross (vec a, vec b) {
   return (abs(a.x*b.y - a.y*b.x));
}
...
do stuff
...
float sum1=0;
float sum2=0;
for (int n=0;n<N;++n)  { //equivalent of the Σ
   sum1 += lengthcross(P[n+1],P[n])* 
           (dot(P[n+1],P[n+1]) + dot(P[n+1],P[n]) + dot(P[n],P[n]));
   sum2 += lengthcross(P[n+1],P[n]);
}
return (m/6*sum1/sum2);

编辑：很多小的数学变化