📚 课程目录第4章转动惯量张量的计算与数值实现

4. 转动惯量张量的计算与数值实现

在刚体动力学中，转动惯量并不是一个标量，而是一个二阶张量。只有当刚体绕某个特定轴旋转时，我们才把它投影成一个标量值。真正驱动欧拉方程、决定姿态演化的是完整的 3×3 惯量张量 I。本节围绕“如何从质量分布得到惯量张量、如何用数值方法实现、如何找到惯性主轴”这条主线展开，并给出一个可复现的离散积分案例。

一、从连续积分到离散求和

设刚体占据区域 Ω，密度为 ρ(r)，则惯量张量定义为：

I = ∫ ρ(r) [ (r·r) E − r rᵀ ] dV

其中 E 为 3×3 单位阵，r = (x, y, z)ᵀ。展开后对角元为 Ixx = ∫ρ(y²+z²)dV，Iyy = ∫ρ(x²+z²)dV，Izz = ∫ρ(x²+y²)dV；非对角元为惯性积，例如 Ixy = −∫ρxy dV。注意惯性积带负号，这是很多实现错误的来源。张量对称：Ixy = Iyx，因此只需存储 6 个独立分量。

数值实现时，把刚体离散为 N 个质点（或体素、有限元单元），每个质点质量 mᵢ、位置 rᵢ，则：

I = Σ mᵢ [ (rᵢ·rᵢ) E − rᵢ rᵢᵀ ]

对体素法，mᵢ = ρ·ΔV；对三角网格，可用四面体分解或面片积分公式。离散化的关键是保证质量总和与质心位置正确，否则惯量会整体偏移。

二、平行轴定理：先移到质心

惯量张量必须在质心系下使用才有简洁形式。若已知绕某参考点 O 的张量 I_O，质心相对 O 的位置为 d，总质量 M，则质心系张量为：

I_C = I_O − M [ (d·d) E − d dᵀ ]

这就是张量形式的平行轴定理。实践中常见流程是：先在建模坐标系（如 CAD 原点）下累加得到 I_O，再求质心 d = Σmᵢrᵢ / M，最后用上式平移到质心。若跳过这一步直接用 I_O 代入欧拉方程，角速度演化会出现虚假的耦合项。

检验技巧：平移后应满足 I_C 的三个对角元均不大于 I_O 对应对角元（因为减去的项半正定），可作为数值实现的快速自检。

三、惯性主轴与主转动惯量

惯量张量是实对称矩阵，必可正交对角化：I_C = R Λ Rᵀ，其中 Λ = diag(I₁, I₂, I₃) 为主转动惯量，R 的列向量即惯性主轴在质心系中的方向。数值上就是对 3×3 对称矩阵做特征分解。三个主惯量按大小排序后，可判断刚体类型：三者互异为不对称陀螺，两个相等为对称陀螺，三个相等为球形陀螺。

主轴的意义在于：绕主轴旋转时角动量与角速度共线，欧拉方程退化为无耦合形式，数值积分更稳定。因此工程上常先把模型旋转到主轴坐标系再仿真，避免大惯量积带来的刚性耦合。

四、数值积分案例：长方体惯量张量

取一个教学示例：均质长方体，尺寸 0.4 m × 0.2 m × 0.1 m，质量 2 kg，密度 ρ = 2/(0.4×0.2×0.1) = 250 kg/m³，几何中心在原点。解析解为 Ixx = M(b²+c²)/12 = 2×(0.04+0.01)/12 ≈ 0.00833 kg·m²，Iyy = 2×(0.16+0.01)/12 ≈ 0.02833，Izz = 2×(0.16+0.04)/12 ≈ 0.03333，惯性积全为零（因对称）。下面用体素离散验证：把每个方向分成 n 段，共 n³ 个体素，每个体素质量 m = M/n³，中心坐标取体素几何中心。

import numpy as np

def inertia_tensor_voxel(L, M, n):
    # L: 三边长 (a,b,c)，M: 总质量，n: 每方向分段数
    a, b, c = L
    m = M / n**3
    I = np.zeros((3, 3))
    # 体素中心坐标
    xs = (np.arange(n) + 0.5) * a / n - a / 2
    ys = (np.arange(n) + 0.5) * b / n - b / 2
    zs = (np.arange(n) + 0.5) * c / n - c / 2
    for x in xs:
        for y in ys:
            for z in zs:
                r = np.array([x, y, z])
                I += m * (np.dot(r, r) * np.eye(3) - np.outer(r, r))
    return I

I = inertia_tensor_voxel((0.4, 0.2, 0.1), 2.0, 20)
print(np.round(I, 5))

取 n = 20 时，数值结果与解析解的对角元误差在 0.1% 量级，非对角元因对称性相互抵消接近 0。若把长方体沿 x 方向平移 0.3 m（质心移到 (0.3, 0, 0)），累加得到的 I_O 中 Iyy、Izz 会各增加 M·d² = 2×0.09 = 0.18，而 Ixx 不变。此时用平行轴定理减回质心，应恢复原值，这正是验证平移公式是否写对的直接手段。

五、实现要点与常见陷阱

质量守恒检查：累加前先确认 Σmᵢ 与设计质量一致，体素法尤其容易因边界体素被重复或遗漏而偏差。
质心先行：先算质心再算张量，或算完张量后立即平移，避免在非质心系下做后续动力学。
对称性利用：只累加上三角 6 个分量，最后镜像填充，可减少一半计算量并天然保证对称。
特征值排序：np.linalg.eigh 返回升序特征值，主轴方向需与特征值一一对应，不要错位。
单位统一：CAD 常以 mm 建模，惯量对长度是二次量纲，换算时是 10⁻⁶ 而非 10⁻³，这是最常见的量级错误。

若数值得到的惯量张量出现明显非对角元而模型本应对称，优先检查体素中心坐标是否相对几何中心正确偏移，以及是否误把惯性积的负号写成了正号。

完成上述步骤后，得到的 I_C 与主轴 R 即可直接送入欧拉方程或刚体仿真器。下一阶段的所有角速度、角动量计算，都应在这个质心主轴系或经 R 变换后的坐标系中进行，以保证数值稳定与物理一致。