📚 课程目录第4章姿态动力学建模与仿真

4. 姿态动力学建模与仿真

姿态动力学回答的是“力矩如何改变航天器的转动状态”。在仿真中，如果只建立运动学模型，我们只能描述姿态如何随时间变化，却无法解释变化为何发生。动力学模型把执行机构力矩、环境干扰力矩与星体角速度联系起来，是姿态控制软件闭环仿真的核心。本节围绕刚体转动惯量、角动量定理、欧拉动力学方程、陀螺力矩以及数值仿真实现展开。

一、刚体转动惯量与惯量矩阵

刚体转动惯量描述质量相对转轴的分布对转动阻碍的程度。对连续体，绕某轴的转动惯量为 I = ∫ r² dm，其中 r 是质量微元到转轴的垂直距离。航天器并非质点，绕不同轴的转动惯量通常不同，因此完整描述需要 3×3 惯量矩阵：

I = [ Ixx  -Ixy  -Ixz
     -Ixy  Iyy  -Iyz
     -Ixz  -Iyz  Izz ]

对角线元素是绕本体坐标轴的主转动惯量，非对角元素是惯性积，反映质量分布的不对称。若本体坐标系与惯性主轴重合，惯性积为零，惯量矩阵变为对角阵，欧拉方程会显著简化。工程上常通过质量特性计算或试验测量获得惯量矩阵，并在仿真初始化时以常量或随质量变化（如帆板展开、燃料消耗）的时变参数输入。

需要特别说明坐标系约定：本文所有角速度、力矩和惯量均在本体坐标系 O_bx_by_bz_b 下描述，本体坐标系固连于星体，原点取质心。角速度矢量记为 ω = [ω_x, ω_y, ω_z]^T，外力矩记为 M = [M_x, M_y, M_z]^T。后续公式与图示均遵循这一约定。

二、角动量定理与欧拉动力学方程

刚体绕质心的角动量定义为 H = Iω。角动量定理指出，本体坐标系下角动量的变化率与外力矩满足：

dH/dt + ω × H = M

其中 ω × H 是由于本体坐标系旋转产生的附加项。将 H = Iω 代入并展开，得到欧拉动力学方程。若惯量矩阵为对角阵，即 I_xx、I_yy、I_zz 为主惯量，则三分量形式为：

Ixx·dωx/dt = (Iyy - Izz)·ωy·ωz + Mx
Iyy·dωy/dt = (Izz - Ixx)·ωz·ωx + My
Izz·dωz/dt = (Ixx - Iyy)·ωx·ωy + Mz

这三个方程是姿态动力学仿真的核心。它们表明：即使外力矩为零，只要三轴主惯量不相等且存在角速度耦合项，角速度仍可能发生复杂变化。例如绕中间轴旋转的刚体可能出现翻滚，这正是欧拉方程非线性耦合的体现。仿真中通常把上述方程写成状态方程 ω̇ = f(ω, M)，再与姿态运动学方程联立积分。

三、陀螺力矩与含飞轮/动量轮的情形

当航天器装有高速旋转的飞轮、动量轮或控制力矩陀螺时，总角动量由星体角动量与转子角动量共同组成。设转子相对本体的角动量为 h_w，则总角动量 H = Iω + h_w。代入角动量定理可得：

I·ω̇ + ω × (I·ω + h_w) = M_ext - ḣ_w

其中 ω × h_w 就是陀螺力矩项。它的物理意义是：当星体存在角速度时，旋转转子的角动量方向被本体坐标系带动，产生一个与 ω 和 h_w 都垂直的力矩。陀螺力矩既可能成为干扰，也可能被控制力矩陀螺用来产生控制力矩。仿真时必须把转子转速、转子惯量和安装方向一并建模，否则无法复现陀螺效应。

在仿真初始化时，建议先检查惯量矩阵是否正定、主惯量是否满足三角不等式（Ixx + Iyy ≥ Izz 等），否则刚体转动惯量在物理上不可实现，积分结果会出现异常。

四、仿真实现与结果判读

下面给出一个最小可运行的欧拉方程数值仿真示例，采用四阶龙格-库塔法积分，时间步长 0.01 s，仿真 10 s。示例参数为：Ixx = 100 kg·m²，Iyy = 120 kg·m²，Izz = 80 kg·m²，初始角速度 ω = [0.05, 0.02, 0.01] rad/s，外力矩为零，用于观察自由转动下的角速度耦合。

import numpy as np

I = np.diag([100.0, 120.0, 80.0])
Ixx, Iyy, Izz = np.diag(I)

def dynamics(w, M):
    wx, wy, wz = w
    dwx = ((Iyy - Izz) * wy * wz + M[0]) / Ixx
    dwy = ((Izz - Ixx) * wz * wx + M[1]) / Iyy
    dwz = ((Ixx - Iyy) * wx * wy + M[2]) / Izz
    return np.array([dwx, dwy, dwz])

def rk4_step(w, M, dt):
    k1 = dynamics(w, M)
    k2 = dynamics(w + 0.5*dt*k1, M)
    k3 = dynamics(w + 0.5*dt*k2, M)
    k4 = dynamics(w + dt*k3, M)
    return w + dt*(k1 + 2*k2 + 2*k3 + k4)/6.0

w = np.array([0.05, 0.02, 0.01])
M = np.zeros(3)
dt, T = 0.01, 10.0
for _ in range(int(T/dt)):
    w = rk4_step(w, M, dt)
print("final omega =", w)

运行后可以观察三轴角速度的周期性交换：由于 Iyy 最大、Izz 最小，ω_x 与 ω_z 之间通过耦合项相互激励，而 ω_y 变化相对平缓。若把初始角速度改为绕中间轴 ω = [0, 0.05, 0]，则会看到 ω_x 和 ω_z 被逐渐激发，这正是中间轴不稳定的典型现象。判读结果时，应关注角速度是否出现非物理发散、总角动量模长是否守恒（无外力矩时 |H| 应保持常数）、以及能量是否在数值积分中缓慢漂移。若 |H| 明显漂移，应减小步长或改用辛积分方法。

把上述动力学模块与姿态运动学模块串联，就得到完整的姿态仿真链路：控制力矩 → 欧拉方程 → 角速度 → 运动学积分 → 姿态四元数/欧拉角 → 敏感器测量 → 控制律。实际工程中还需考虑帆板挠性、燃料晃动、执行机构饱和与延迟等因素，但刚体欧拉方程始终是这条链路的基准模型。建议在仿真报告中同时给出角速度曲线、角动量模长曲线和能量曲线，三者互相印证，才能判断动力学积分是否可信。