4. 刚体姿态动力学基础
前三章我们建立了姿态的描述方式与坐标变换关系，但还没有回答一个根本问题：是什么让姿态发生变化？答案来自角动量与力矩的关系。本章从角动量的定义出发，引入惯性张量这一描述刚体质量分布的核心量，推导出欧拉方程，并分析主轴坐标系下方程的简化形式与耦合特性。这是后续所有姿态动力学仿真模型的数学基础。
一、角动量：姿态动力学的状态量
刚体绕质心旋转时，其角动量 L 定义为各质量微元位置与动量的叉积之和。对连续体，写成积分形式：
L = ∫ r × (ω × r) dm其中 r 是微元相对质心的位置矢量，ω 是刚体角速度。展开三重叉积可得：
L = ∫ [ (r·r)ω − r(r·ω) ] dm = I · ω这里 I 就是惯性张量。注意角动量 L 与角速度 ω 一般不共线，这是刚体动力学区别于质点动力学的关键。只有当 ω 沿主轴方向时，两者才同向。
二、惯性张量：质量分布的完整描述
惯性张量是一个 3×3 对称矩阵，在体坐标系下写为：
⎡ Ixx −Ixy −Ixz ⎤
I = ⎢ −Ixy Iyy −Iyz ⎥
⎣ −Ixz −Iyz Izz ⎦对角元素是转动惯量，例如 Ixx = ∫(y²+z²)dm，表示绕 x 轴旋转的惯性；非对角元素是惯性积，如 Ixy = ∫xy dm，反映质量分布的不对称性。惯性张量有两条重要性质：
- 对称性：Ixy = Iyx，独立分量只有 6 个；
- 正定性：对任意非零 ω，ωᵀIω > 0，对应转动动能恒为正。
惯性张量依赖坐标系。同一刚体在不同体坐标系下，矩阵元素不同，但特征值（主转动惯量）不变。这为下一节的主轴坐标系提供了理论依据。
三、欧拉方程推导
角动量定理在惯性系中形式简洁：dL/dt = M，M 为外力矩。但 L = I·ω 中的 I 在体坐标系下才是常数矩阵，因此需要把导数转到体坐标系。利用矢量导数变换关系：
(dL/dt)_惯性系 = (dL/dt)_体系 + ω × L代入 L = I·ω，并令体系下 dI/dt = 0，得到矢量形式的欧拉方程：
I · dω/dt + ω × (I · ω) = M展开到主轴坐标系（此时 I 为对角阵，Ixx=Ix, Iyy=Iy, Izz=Iz），得到标量方程组：
Ix·ω̇x − (Iy − Iz)·ωy·ωz = Mx
Iy·ω̇y − (Iz − Ix)·ωz·ωx = My
Iz·ω̇z − (Ix − Iy)·ωx·ωy = Mz这就是经典的欧拉动力学方程。三个方程通过 ωy·ωz、ωz·ωx、ωx·ωy 这些交叉项相互耦合，即使外力矩为零，姿态也可能持续变化。
四、主轴坐标系下的方程与耦合特性
主轴坐标系是使惯性积全部为零的体坐标系，此时惯性张量对角化。工程上，航天器通常设计为接近主轴对称，因此主轴坐标系是姿态动力学建模的首选参考系。在主轴系下，欧拉方程呈现最简形式，但耦合项依然存在，其物理本质是陀螺力矩。
耦合特性可以从三个角度理解：
- 交叉激励：ωy 与 ωz 的乘积项会改变 ωx 的变化率，反之亦然。这意味着绕一个轴的旋转会“泄漏”到另外两个轴。
- 惯量差决定耦合强度：若 Iy ≈ Iz，则 x 轴方程中的耦合项 (Iy−Iz)ωyωz 很小，x 轴近似解耦。球对称刚体（Ix=Iy=Iz）完全无耦合。
- 无外力矩下的自由运动：M=0 时，方程仍有非平凡解。绕最大或最小惯量主轴旋转是稳定的，绕中间惯量主轴旋转不稳定——这就是著名的“网球拍定理”或“贾尼别科夫效应”。
五、一个可复现的仿真示例
设某航天器主转动惯量为 Ix=100 kg·m²，Iy=200 kg·m²，Iz=150 kg·m²，初始角速度 ω₀=(0.1, 0.05, 0.02) rad/s，外力矩为零。用四阶龙格-库塔法以步长 0.01 s 积分 50 s，可观察到：
- ωx 在 ±0.12 rad/s 范围内振荡，周期约 8 s；
- ωy 和 ωz 同步振荡，幅值分别约 0.06 和 0.03 rad/s；
- 角动量模长 |L| 保持在 12.25 kg·m²/s 附近，波动小于 10⁻⁶；
- 转动动能 T 保持在 0.6125 J 附近，验证了无外力矩下的守恒性。
若把初始角速度改为绕中间惯量轴 ω₀=(0, 0.1, 0)，仿真会显示 ωy 在数秒后突然翻转符号，而 |L| 和 T 仍然守恒。这正是中间轴不稳定的数值体现，也是姿态控制中需要避免的工况。