📚 课程目录第3章刚体动力学基本方程

3. 刚体动力学基本方程

上一节我们把刚体转动描述为“绕轴转”，但真实刚体很少只绕一条固定轴旋转。卫星在太空中翻滚、无人机做横滚—俯仰—偏航耦合机动、机械臂末端执行器快速换向，都会让角速度方向不断改变。此时仅用一个标量转动惯量已经不够，必须引入转动惯量张量，并建立角动量与力矩之间的矢量关系。本节围绕这条主线展开：先明确张量的物理含义，再推导欧拉方程，最后讨论数值求解时最先遇到的几个问题。

一、转动惯量张量：从标量到 3×3 矩阵

在刚体上建立固连坐标系（随体坐标系），原点取在质心或指定参考点。刚体绕任意方向单位矢量 n 的转动惯量可写成二次型 I_n = n^TIn，其中

        ⎡ Ixx  Ixy  Ixz ⎤
I   =   ⎢ Iyx  Iyy  Iyz ⎥
        ⎣ Izx  Izy  Izz ⎦

Ixx = ∫(y²+z²)dm,  Iyy = ∫(x²+z²)dm,  Izz = ∫(x²+y²)dm
Ixy = Iyx = -∫xy dm,  Ixz = Izx = -∫xz dm,  Iyz = Izy = -∫yz dm

对角元是绕三个坐标轴的转动惯量，非对角元称为惯性积，带负号是工程力学的通用约定。张量是对称矩阵，因此总可以找到一个正交变换使其对角化，对角化后的三个轴称为主轴，对应的对角元为主转动惯量 I₁、I₂、I₃。数值模拟中，若几何体规则或可近似为长方体、圆柱、椭球，通常直接在主轴上建模，令惯性积为零，可显著简化方程。

需要强调：转动惯量张量是随体量。刚体转动时，它在惯性系中的分量会变化，但在固连系中保持不变。这一性质是后面推导欧拉方程的关键。

二、角动量与力矩的关系

刚体对参考点的角动量为 L = Iω。注意 L 与 ω 一般不共线，只有当 ω 沿主轴方向时才共线。这一点常被初学者忽略，也是陀螺效应、进动等复杂运动的根源。

力矩与角动量的关系为 M = dL/dt。在惯性系中直接对 L 求导，需要跟踪张量的姿态变化；更实用的做法是转到固连系。设固连系相对惯性系的角速度为 ω，则任意矢量在惯性系中的时间导数与固连系中的时间导数满足

(dL/dt)_惯性 = (dL/dt)_固连 + ω × L

代入 L = Iω，并利用固连系中 I 为常量，得到欧拉方程：

I·(dω/dt) + ω × (I·ω) = M

写成分量形式（取主轴坐标系，惯性积为零）：

I1·dω1/dt + (I3 - I2)·ω2·ω3 = M1
I2·dω2/dt + (I1 - I3)·ω3·ω1 = M2
I3·dω3/dt + (I2 - I1)·ω1·ω2 = M3

这三个方程就是刚体转动数值模拟的核心。它们是一阶非线性常微分方程组，耦合项 (I3−I2)ω2ω3 等正是陀螺力矩的体现。当 I₁=I₂=I₃ 时耦合项消失，退化为三个独立的一阶方程，对应球对称刚体。

三、数值求解的初步考虑

欧拉方程给出角速度的演化，但要得到刚体姿态，还需补充姿态运动学方程。常用四元数 q=(q0,q1,q2,q3) 描述姿态，其导数为

dq/dt = 0.5 · q ⊗ (0, ω1, ω2, ω3)

其中 ⊗ 表示四元数乘法。这样得到 7 个状态量（3 个角速度 + 4 个四元数）的联立方程组，可用四阶龙格—库塔法（RK4）推进。下面给出一个可直接参考的 RK4 单步实现框架（Python 风格伪代码，参数为示例假设）：

def f(state, I, M):
    w = state[0:3]; q = state[3:7]
    dw = inv(I) @ (M - cross(w, I @ w))
    dq = 0.5 * quat_mul(q, [0, w[0], w[1], w[2]])
    return concat(dw, dq)

# 示例参数：I = diag(0.02, 0.03, 0.05) kg·m²，M = [0.1, 0, 0] N·m
# 步长 h = 1e-3 s，初值 w = [0.5, 0.1, 0.2] rad/s，q = [1,0,0,0]
k1 = f(s, I, M)
k2 = f(s + 0.5*h*k1, I, M)
k3 = f(s + 0.5*h*k2, I, M)
k4 = f(s + h*k3, I, M)
s_next = s + h/6*(k1 + 2*k2 + 2*k3 + k4)
s_next[3:7] /= norm(s_next[3:7])   # 四元数归一化

数值实现中有三个要点必须注意。第一，四元数归一化：每步积分后都要除以模长，否则舍入误差会让姿态漂移。第二，步长选择：欧拉方程含 ω² 量级的耦合项，步长应满足 h·|ω|_max 远小于 1，示例中 1e-3 s 对应约 0.5 rad/s 的角速度是合理的，若角速度升到 50 rad/s 则需缩小到 1e-5 s 量级。第三，能量与角动量守恒检验：无外力矩时，转动动能 T = ½ω^TIω 和惯性系角动量 |L| 应保持不变，可作为判断积分器是否可靠的依据。

四、一个可观察的示例：非对称刚体的自由翻滚

取示例参数 I₁=0.02、I₂=0.03、I₃=0.05 kg·m²，初始角速度 ω₀=(0.5, 0.1, 0.2) rad/s，外力矩为零，步长 1e-3 s，积分 10 s。由于三个主转动惯量互不相等，角速度分量会周期性交换能量，表现为翻滚运动；若把初值改为绕最小惯量轴（I₁）或最大惯量轴（I₃）附近的小扰动，则运动相对稳定，而绕中间惯量轴（I₂）的扰动会迅速放大——这就是经典的“网球拍定理”。数值上可通过观察 |L| 与 T 的相对漂移判断积分质量：若 10 s 内相对漂移超过 1e-4，应减小步长或改用辛积分器。

建模时先判断能否取主轴坐标系：若刚体有对称面或近似规则，令惯性积为零可把 6 个独立张量分量减为 3 个主惯量，既减少输入错误，也让欧拉方程的耦合结构一目了然。

不要用标量转动惯量直接替代张量做三轴耦合仿真。只有当角速度始终沿同一主轴时，标量模型才成立；一旦出现换轴或外力矩方向变化，忽略惯性积会带来明显的姿态误差。

小结本节逻辑：转动惯量张量刻画刚体质量分布的各向异性，角动量与角速度通过张量联系且一般不共线，力矩等于惯性系中角动量的变化率，转到固连系后得到含陀螺耦合项的欧拉方程；数值上以角速度和四元数为状态量，用 RK4 推进并每步归一化，再用动能与角动量守恒检验结果可信度。这条链路是后续所有刚体转动仿真的共同基础。