计算二维惯性张量

Question

我正在研究如何找到二维形状的惯性。这个形状的轮廓由几个点网格化，每个点的x和y坐标是已知的。

我知道Ixx、Iyy和Ixy的表达式，但是body没有质量。我该如何进行？

Answer 1

不管你有什么形状，你都需要把它分成三角形，然后分别处理每个三角形。然后最后使用以下规则组合结果

总体

% Combined total area of all triangles
total_area = SUM( area(i), i=1:n )
total_mass = SUM( mass(i), i=1:n )
% Combined centroid (center of mass) coordinates
combined_centroid_x = SUM( mass(i)*centroid_x(i), i=1:n)/total_mass
combined_centroid_y = SUM( mass(i)*centroid_y(i), i=1:n)/total_mass
% Each distance to triangle (squared)
centroid_distance_sq(i) = centroid_x(i)*centroid_x(i)+centroid_y(i)*centroid_y(i)
% Combined mass moment of inertia
combined_mmoi = SUM(mmoi(i)+mass(i)*centroid_distance_sq(i), i=1:n)

现在每个三角形。

考虑向量坐标的三个角顶点，点A、B和C

a=[ax,ay]
b=[bx,by]
c=[cx,cy]

以及以下点积和叉积（标量）组合

a·a = ax*ax+ay*ay
b·b = bx*bx+by*by
c·c = cx*cx+cy*cy
a·b = ax*bx+ay*by
b·c = bx*cx+by*cy
c·a = cx*ax+cy*ay
a×b = ax*by-ay*bx
b×c = bx*cy-by*cx
c×a = cx*ay-cy*ax

三角形的属性是（t(i)厚度和rho质量密度）

area(i) = 1/2*ABS( a×b + b×c + c×a )
mass(i) = rho*t(i)*area(i)
centroid_x(i) = 1/3*(ax + bx + cx)
centroid_y(i) = 1/3*(ay + by + cy)
mmoi(i) = 1/6*mass(i)*( a·a + b·b + c·c + a·b + b·c + c·a )

按组件以上是

area(i) = 1/2*ABS( ax*(by-cy)+ay*(cx-bx)+bx*cy-by*cx)
mmoi(i) = mass(i)/6*(ax^2+ax*(bx+cx)+bx^2+bx*cx+cx^2+ay^2+ay*(by+cy)+by^2+by*cy+cy^2)

附录

这里有一点理论。使用

找到每个三角形的面积

Area = 1/2 * || (b-a) × (c-b) ||

其中×是向量叉积，|| .. ||是向量范数（长度函数）。

三角形由两个变量 t 和 s 参数化，使得二重积分 A = INT(INT(1,dx),dy) 给出总面积

% position r(s,t) = [x,y]
[x,y] = [ax,ay] + t*[bx-ax, by-zy] + t*s*[cx-bx,cy-by]

% gradient directions along s and t
(dr/dt) = [bx-ax,by-ay] + s*[cx-bx,cy-by]
(dr/ds) = t*[cx-bx,cy-by]

% Integration area element
dA = || (dr/ds)×(dr/dt) || = (2*A*t)*ds*dt
%
%   where A = 1/2*||(b-a)×(c-b)||

% Check that the integral returns the area
Area = INT( INT( 2*A*t,s=0..1), t=0..1) = 2*A*(1/2) = A

% Mass moment of inertia components

         /  /  /  | y^2+z^2  -x*y    -x*z   |
I = 2*m*|  |  | t*|  -x*y   x^2+z^2  -y*z   | dz ds dt
        /  /  /   |  -x*z    -y*z   x^2+y^2 |

% where [x,y] are defined from the parametrization

Answer 2

我想稍微纠正一下优秀的 John Alexiou 回答：

三角形 mmoi 算法期望角是相对于三角形质心（质心）定义的。所以在计算 mmoi
Shape mmoi 算法期望每个三角形的质心相对于形状质心（组合质心）定义。所以在计算 mmoi 之前从每个三角形质心减去组合质心。

因此结果代码如下所示：

public static float CalculateMMOI(Triangle[] triangles, float thickness, float density)
{
    float[] mass = new float[triangles.Length];
    float[] mmoi = new float[triangles.Length];
    Vector2[] centroid = new Vector2[triangles.Length];
    
    float combinedMass = 0;
    float combinedMMOI = 0;
    Vector2 combinedCentroid = new Vector2(0, 0);
    
    for (var i = 0; i < triangles.Length; ++i) 
    {
        mass[i] = triangles[i].CalculateMass(thickness, density);
        mmoi[i] = triangles[i].CalculateMMOI(mass[i]);
        centroid[i] = triangles[i].CalculateCentroid();

        combinedMass += mass[i];
        combinedCentroid += mass[i] * centroid[i];
    }
    combinedCentroid /= combinedMass;


    for (var i = 0; i < triangles.Length; ++i) {
        combinedMMOI += mmoi[i] + mass[i] * (centroid[i] - combinedCentroid).LengthSquared();
    }

    return combinedMMOI;
}


public struct Triangle
{
    public Vector2 A, B, C;
    
    public float CalculateMass(float thickness, float density)
    {
        var area = CalculateArea();
        return area * thickness * density;
    }
    
    public float CalculateArea()
    {
        return 0.5f * Math.Abs(Vector2.Cross(A - B, A - C));
    }

    public float CalculateMMOI(float mass)
    {
        var centroid = CalculateCentroid()
        var a = A - centroid;
        var b = B - centroid;
        var c = C - centroid;
        
        var aa = Vector2.Dot(a, a);
        var bb = Vector2.Dot(b, b);
        var cc = Vector2.Dot(c, c);
        var ab = Vector2.Dot(a, b);
        var bc = Vector2.Dot(b, c);
        var ca = Vector2.Dot(c, a);
        return (aa + bb + cc + ab + bc + ca) * mass / 6f;
    }

    public Vector2 CalculateCentroid()
    {
        return (A + B + C) / 3f;
    }
}


public struct Vector2
{
    public float X, Y;
    
    public float LengthSquared()
    {
        return X * X + Y * Y;
    }
    
    public static float Dot(Vector2 a, Vector2 b)
    {
        return a.X * b.X + a.Y * b.Y;
    }

    public static float Cross(Vector2 a, Vector2 b)
    {
        return a.X * b.Y - a.Y * b.X;
    }
}

计算二维惯性张量

Computing tensor of Inertia in 2D

c++

physics

scientific-computing