刚体惯量与转动方程

设想将两个质量相同的砝码安装在可绕固定轴转动的轻杆上：当把砝码沿杆向外推移后再施加同样大小的驱动力矩时，轻杆角速度的增长会明显变慢。这是因为增大砝码到转轴的距离，不仅提高了同一角速度下砝码的切向线速度，也成倍增加了达到该转速所需的动能。从这一基本物理图景出发，我们可以直接导出转动惯量与轴距平方之间的定量关系。

第二章已经通过旋转矩阵求导给出了刚体上任一点的速度向量公式；为了看清标量转动惯量的几何起源，这里先直接从圆周轨迹的弧长与转角关系出发推导切向速率，再代入动能表达式（角速度的正方向依然遵循第二章的右手螺旋定则）。

考虑一个质量为 m 的质点绕空间中的固定直线轴转动，设质点到转轴的垂直距离为 d。从质点向转轴作垂线，随着刚体转动，质点在空间中描绘出以垂足为圆心、以 d 为半径的圆周。我们只需先由这段圆周轨迹求出质点的瞬时速率，即可立即写出它的动能。

在弧度制下，转角 \theta 定义为圆周上的有向弧长 s 与半径 d 的比值，即 s=d\theta。由于转轴固定且刚体内部距离不变，半径 d 为常数；对时间求导便得到质点的切向速度分量与速率：

式中角速度 \omega=\dot\theta 的正负号区分正反两个转向，速率 v 则取其绝对值。这一正比关系直观地说明：在相同时间内扫过同一转角的两个质点，离转轴越远，走过的圆弧越长，线速率也就越大。

例如，在图中取轴距 d=0.8\,\mathrm m、角速度 \omega=1\,\mathrm{rad/s}、时间间隔 \Delta t=0.6\,\mathrm s，质点走过的圆弧长为 0.48\,\mathrm m，沿圆弧的平均速率为 0.8\,\mathrm{m/s}；当逐步缩短观察时间间隔 \Delta t 时，端点弦长与时间之比也迅速收敛到同一数值。将所得切向速率 v=d|\omega| 代入质点的动能公式，可得：

式中由质量与轴距平方构成的组合量 md^2，称为该质点关于这根转轴的转动惯量（moment of inertia）。当多个质点刚性固连在一起、以相同的角速度 \omega 绕同一根固定轴转动时，将各质点的动能逐项相加，即可得到质点系的总转动动能：

在本节中，转轴在空间中的位置已经固定，下标 n 用于标识转轴的单位正方向。转动惯量 I_n 的国际单位为 \mathrm{kg\,m^2}：每个质点对转动惯量的贡献与其到转轴垂直距离的平方成正比，而质点沿转轴方向的位置移动则不会改变这一贡献。

在动力学上，绕固定轴转动的刚体满足定轴转动方程 \tau_n=I_n\dot\omega，其中 \tau_n 是作用在刚体上的合外力矩沿转轴正方向的投影分量，该力矩做功的瞬时功率为 P=\tau_n\omega。由于转轴固定在刚体上，各质点到转轴的距离 d_i 保持恒定，因此标量惯量 I_n 是一个与时间无关的常数。需要注意的是，固定转轴的轴承支撑还会提供垂直于转轴的约束反力矩以维持轴向不变，而沿转轴方向的分量 \tau_n 则单独决定了转速 \omega 的快慢变化。当刚体在三维空间中自由运动时，转轴方向本身也会随时间改变，这就需要把三个空间方向的质量分布与力矩耦合统一纳入矩阵方程。

在空间中自由运动的刚体同时包含整体平动与绕某一点的转动。如果选取刚体的质心（center of mass）作为基准参考点，不仅能用质心一个点的运动完整代表整个刚体的平动，还能使刚体各部分质量相对参考点的位置偏差相互平衡。这一平衡性质将在接下来的动能分解与平行轴定理推导中直接消去平动与转动的交叉耦合项。

设刚体由离散质点组成，第 i 个质点的质量为 m_i、世界坐标位置为 \boldsymbol x_i。刚体的总质量 m 为各质点质量之和，质心位置 \boldsymbol c 则定义为各质点位置的质量加权平均：

由定义可知，质量越大的质点对质心位置的影响权重越大。例如，将质量分别为 1\,\mathrm{kg} 和 3\,\mathrm{kg} 的两个质点放在一维数轴的 0\,\mathrm m 和 2\,\mathrm m 处，它们的质心便偏向较重一侧，位于 (1\times0+3\times2)/4=1.5\,\mathrm m 处。

记第 i 个质点相对于质心的位置向量为 \boldsymbol r_i=\boldsymbol x_i-\boldsymbol c。将质心定义式（4.2）改写为相对位置的形式，即可得到一个在刚体力学中反复使用的基本恒等式：

连续质量分布

对于质量连续分布的实体刚体，可以将其划分为无穷多个微小体积元，每个微元各自具有位置 \boldsymbol x 与质量微元 \mathrm dm。当划分无限加细时，离散求和自然过渡为体积分。设刚体在位置 \boldsymbol x 处的体密度为 \rho(\boldsymbol x)，体积微元为 \mathrm dV，则质量微元为 \mathrm dm=\rho(\boldsymbol x)\mathrm dV，总质量、质心位置以及相对位置的质量加权平衡关系分别写为：

上述积分遍历刚体所占据的全部空间区域，向量积分则按三个坐标分量分别积分（例如横坐标分量为 \int x\,\mathrm dm）。一般而言，体密度 \rho(\boldsymbol x) 可以随空间位置变化；只有当刚体材质均匀（\rho 为常数）时，才能将密度提到积分号外。对于密度均匀的对称几何体（如均匀长方体、均匀球体），由于任一质量微元关于几何中心都存在质量相等、相对位置相反的镜像微元，其质心必然与几何中心重合。

作为一维连续分布的简单例子，考察一根沿 x 轴从 0 延伸到 L、总质量为 m 的均匀细杆。其线密度为常数 \lambda=m/L，长度微元 \mathrm dx 对应的质量微元为 \mathrm dm=\lambda\,\mathrm dx，直接积分即可求得其质心位于中点 c=L/2。同理，二维均匀薄板使用面密度乘面积微元，三维实体使用体密度乘体积微元，三者在物理本质上都是对质量微元的加权求和。

质心的平动方程

当刚体上各点受到外力作用时，记所有外力的矢量和为合外力 \boldsymbol F。对质心定义求二阶时间导数可知，刚体质心的运动满足 m\ddot{\boldsymbol c}=\boldsymbol F；也就是说，整个刚体的整体平动完全等同于一个集中了全部质量 m 并承受全部合外力 \boldsymbol F 的质点。

质心平动定理将合外力与刚体的整体平移直接联系起来，余下的运动便是各质点围绕质心的相对转动。接下来我们将各质点的速度分解式代入总动能，继续利用式（4.3）的质量加权平衡性质分离平动动能与转动动能。

下方实验将三个离散质点置于同一平面内，便于同时观察质心位置与轴距的变化。尝试调整第三个质点的质量或位置，观察质心如何向重质点偏移，并对比绕质心竖轴与绕固定原点竖轴的转动惯量读数。

可以看到，增大某个质点的质量会将质心拉向该质点；而在保持各质点质量与位置不变的前提下切换转轴参考点（质心或原点），改变的则是各质点到转轴的距离及相应的切向速度分布。下一节我们将把绕任意方向转轴的这种距离平方加权关系整理为统一的惯量矩阵。

设刚体的质心线速度为 \boldsymbol v_C，世界系下的角速度为 \boldsymbol\omega。根据第二章的刚体点速度公式，第 i 个质点的世界速度可写为质心平动速度与绕质心转动速度之和：\boldsymbol v_i=\boldsymbol v_C+\boldsymbol\omega\times\boldsymbol r_i。将其代入全体质点的总动能表达式，质心条件 \sum_i m_i\boldsymbol r_i=\boldsymbol0 会使平动与转动的交叉项恒等于零，从而将总动能干净地分解为质心平动动能与绕质心转动动能两项：

当刚体绕某一固定方向的转轴转动时，我们可以按第 1 节的方法先算出该轴的标量惯量 I_n，再乘以角速度平方的一半。然而在一般的三维运动中，角速度向量 \boldsymbol\omega 的方向随时间连续变化，若每次改变方向都重新对全体质点求和既繁琐又低效。注意到式（4.4）第二项中角速度 \boldsymbol\omega 的各个分量总是以二次齐次形式出现，我们可以把仅与质量分布相关的几何系数提取出来，预先汇总为一个 3\times3 实对称矩阵，称为刚体关于质心的惯量矩阵（inertia matrix）或惯性张量（inertia tensor）：

矩阵分量与质量分布

将质量微元相对于参考点的直角坐标写为 \boldsymbol r=(x,y,z)^{\mathsf T}，展开 \|\boldsymbol r\|^2E-\boldsymbol r\boldsymbol r^{\mathsf T}，可得单个质量微元 \mathrm dm 对惯量矩阵的贡献：

矩阵的三个对角元（称为转动惯量，moments of inertia）分别是质量微元到 x,y,z 三根坐标轴的垂直距离平方（如 y^2+z^2 为到 x 轴的距离平方）乘以质量微元后的积分。六个非对角元（其相反数或本身常称为惯性积，products of inertia）则刻画了质量分布在不同坐标轴方向之间的耦合不对称性：例如 I_{12}=-\int xy\,\mathrm dm，位于第一、三象限（x,y 同号）的质量微元对其产生负贡献，位于第二、四象限（x,y 异号）的质量微元产生正贡献。若刚体的质量分布关于坐标平面具有镜像对称性，这些正负贡献便会成对抵消，使非对角元归零。

由于惯量矩阵的每一个元素都只取决于刚体的几何形状与质量分布，在刚体自身的局部坐标系下只需预先计算一次。一旦有了惯量矩阵 I，对于任意给定的单位转轴方向 \boldsymbol n，绕该轴的标量转动惯量只需通过二次型直接取出：

这一二次型在数值上仍然严格等于各质点到转轴 \boldsymbol n 的垂直距离平方与质量的乘积之和，因此恒满足非负性 I_n\ge0。只有当刚体的全部质量都恰好退化在某一条直线轴上时，绕该轴的转动惯量才会为零；只要任意直线轴外都分布有正质量，惯量矩阵 I 就是严格对称正定矩阵。实际物理仿真中具有有限三维体积的正质量刚体（以及平面内不共线的多质点系统）均满足这一正定条件。

回到上一节的实验 02，保持默认参数不变，将转轴参考点从“质心”切换为“固定原点”，即可直接对比两种转轴下的惯量与动能读数。

当整个系统改绕过原点的固定竖轴以同样的角速度 2\,\mathrm{rad/s} 转动时，直接按原点的坐标计算竖轴惯量得 I_{O,33}=1(0.8)^2+1(0.8)^2+2(0.3^2+1^2)=3.46\,\mathrm{kg\,m^2}，相应的总动能增至 \tfrac12(3.46)(2^2)=6.92\,\mathrm J。之所以比绕质心转动多出 6.92-4.74=2.18\,\mathrm J 的动能，是因为此时质心本身也在以半径 \|\boldsymbol c\|=\sqrt{0.15^2+0.5^2}\approx0.5220\,\mathrm m 绕原点做圆周平动，其质心平动动能恰好等于 \tfrac12 m(\|\boldsymbol c\|\omega)^2=\tfrac12(4)(0.2725)(4)=2.18\,\mathrm J。第 5 节将把这一差值推广为任意参考点之间的平行轴定理。

考察一个质量为 m、三条棱长分别为 a,b,c 的均匀长方体。将坐标原点设在长方体质心，令三条坐标轴分别平行于长方体的三条棱边，则长方体的常体密度为 \rho=m/(abc)，质量微元为 \mathrm dm=\rho\,\mathrm dx\,\mathrm dy\,\mathrm dz，三个坐标的积分区间分别为关于原点对称的 [-a/2,a/2]、[-b/2,b/2] 和 [-c/2,c/2]。计算绕横轴（x 轴）的转动惯量时，每个质量微元的轴距平方权重为 y^2+z^2，而横坐标 x 只决定沿轴向叠加的厚度。对三个坐标方向分别积分，即可得到均匀长方体关于质心的对角惯量矩阵：

从式（4.6）可以看出，在保持总质量 m 不变的前提下，拉长横向边长 a 会同时增大绕纵轴和绕竖轴的转动惯量，但绕横轴自身的转动惯量 I_{11}=m(b^2+c^2)/12 却只取决于截面边长 b 和 c。这与第 1 节“质点沿转轴方向平移不改变轴距”的结论完全吻合。至此，离散质点求和与连续几何体积分都归结为了同一种结构的惯量矩阵。

角动量

为了进一步分析外力矩如何改变刚体的转动状态，需要引入与平动动量相对应的转动物理量。在质点平动中，线动量定义为质量与速度的乘积；在刚体转动中，对每个质量微元计算相对质心位置向量 \boldsymbol r_i 与其线动量 m_i\boldsymbol v_i 的叉乘，再对整个刚体求和，便得到刚体关于质心的角动量（angular momentum）\boldsymbol\ell。将刚体各点的速度分解式代入角动量定义，利用质心平衡条件与二重叉乘恒等式，可以直接将角动量表示为惯量矩阵与角速度向量的乘积：

角动量 \boldsymbol\ell 的国际单位为 \mathrm{kg\,m^2/s}（或 \mathrm{N\,m\,s}）。值得注意的是，由于刚体在不同方向上的惯量一般不同，经矩阵 I 线性变换后，角动量向量 \boldsymbol\ell 的方向往往会偏离角速度向量 \boldsymbol\omega 的方向。

主轴与主惯量

如果刚体绕某个特定的单位轴方向 \boldsymbol n 转动时，角动量恰好与角速度保持共线，即满足特征方程 I\boldsymbol n=I_n\boldsymbol n，则称该方向为刚体的惯量主轴（principal axis of inertia），对应的比例系数 I_n 称为主惯量（principal moment of inertia）。从线性代数看，惯量主轴与主惯量正是实对称矩阵 I 的特征向量与特征值。

根据实对称矩阵的正交对角化定理，任意刚体都存在三根相互正交的单位惯量主轴。将这三根主轴按右手定则排列为正交矩阵 U 的三列，则惯量矩阵可对角化为 I=U\operatorname{diag}(I_1,I_2,I_3)U^{\mathsf T}；前文均匀长方体平行于棱边的三根几何对称轴，就自然构成了这样一组主轴坐标系。对于任意单位转轴方向 \boldsymbol n，其标量轴惯量 I_n=\boldsymbol n^{\mathsf T}I\boldsymbol n 等于三个主惯量的非负加权平均（权重为 \boldsymbol n 在三根主轴上方向余弦的平方），因此必有上下界 \min_j I_j\le I_n\le\max_j I_j，且仅当转轴与相应主轴重合时取到极值。在实际计算中，若直接将刚体局部坐标系选在惯量主轴上，惯量矩阵便简化为对角阵，角动量的三个分量只需由各方向角速度分量分别乘以对应的主惯量即可求得。

在下方实验中，将转轴方向 \boldsymbol n=(\cos\beta,\sin\beta,0)^{\mathsf T} 从横轴（\beta=0^\circ）逐渐转到纵轴（\beta=90^\circ），可以先在三维视图中观察角速度向量与角动量向量之间的张角变化，再切换到曲线视图核对标量轴惯量 I_n=I_1\cos^2\beta+I_2\sin^2\beta 的过渡规律。

可以看到，只有当转轴严格对齐横轴（\beta=0^\circ）或纵轴（\beta=90^\circ）这两根主轴时，角速度与角动量的方向箭头才完全重合；在所有中间角度上，角动量箭头均偏向主惯量较大的纵轴一侧，这直观地体现了非主轴转动下的特征分量缩放差异。与此同时，轴惯量曲线从 I_1 平滑上升至 I_2，并在 \beta=45^\circ 处取到两者的算术平均值；若转轴进一步向垂直于该平面的竖轴偏转，轴惯量还会继续向最大主惯量 I_3 增长。

在物理引擎中，刚体的惯量参数通常作为常量保存在随体局部坐标系中，而外力、力矩与角速度则经常在世界坐标系下表达；此外，多体关节动力学还常将参考点从质心移至铰接支点或子树质心。这两类操作分别对应坐标架旋转（换轴）与参考点平移（换点）：前者改变同一惯性张量在不同基底下的矩阵分量，后者则改变了各质量微元到参考点的几何距离。

首先保持参考点在质心不变，考察坐标架旋转的影响。设刚体在随体坐标系下的质心惯量矩阵为常量 I_b，刚体当前从物体系到世界系的姿态矩阵为 R。由于同一个角速度向量和角动量向量在两个坐标系下的分量均通过左乘 R 转换，世界坐标系下的质心惯量矩阵 I 满足如下正交合同变换：

平行轴定理

接下来保持坐标轴的方向不变，将参考点从质心 C 平移到另一点 O。记从新参考点 O 指向质心 C 的位移向量为 \boldsymbol d=\boldsymbol c-\boldsymbol x_O，则原先相对于质心位置为 \boldsymbol r 的质量微元，相对于新参考点 O 的位置变为 \boldsymbol r+\boldsymbol d。将其代入惯量矩阵定义，即可导出著名的平行轴定理（parallel axis theorem）：

在下方实验中，可以先将参考点的横、纵坐标均设为零（与质心重合），单独改变长方体的偏航角：此时世界系惯量矩阵的非对角元随姿态角周期变化，但矩阵的三个特征值（主惯量）保持恒定。随后固定姿态角、将参考点 O 向外平移，对比过参考点与过质心的两条平行竖轴，新增的竖轴惯量恰好等于 m(x_O^2+y_O^2)。

由于图中刚体仅绕竖轴改变姿态角，质心处的竖轴惯量 I_{C,33} 始终保持为 I_3，而水平面内的分量 I_{11},I_{22},I_{12} 则按式（4.8）随姿态发生混合；当拖动参考点偏离质心时，平行轴修正项按水平距离的平方叠加到竖轴惯量上，将参考点移回质心时又恢复为 I_3。

若进一步以刚体上（或与刚体固连的延伸空间中）异于质心的任意固连点 O 作为运动参考点，记点 O 的世界线速度为 \boldsymbol v_O，则各质量微元的世界速度可写为 \boldsymbol v=\boldsymbol v_O+\boldsymbol\omega\times(\boldsymbol r+\boldsymbol d)。将其代入总动能积分，平动与转动之间便在包含各自平方项的同时，还会多出一项依赖于质心偏移 m\boldsymbol d 的交叉耦合项：

式（4.9a）中的中间项将参考点线速度 \boldsymbol v_O 与刚体角速度 \boldsymbol\omega 耦合在一起：只有当参考点选在质心（\boldsymbol d=\boldsymbol0）或参考点在空间中固定不动（\boldsymbol v_O=\boldsymbol0，此时仅剩绕定点的纯转动能 \tfrac12\boldsymbol\omega^{\mathsf T}I_O\boldsymbol\omega）时，该交叉项才会消失；对于一般运动的关节铰接点，这三项必须同时保留。

第五章引入六维空间代数时，将把线速度 \boldsymbol v_O 与角速度 \boldsymbol\omega 拼成六维空间速度，届时 6\times6 空间惯量矩阵之所以在偏离质心时出现非对角耦合块，其物理根源正是式（4.9a）中的中间交叉项。

当外力作用于刚体上的不同点时，除了推动质心产生平动加速度，还会对质心产生力矩。对于质量恒定的刚体，将所有外力关于质心的力矩相加，记为合外力矩 \boldsymbol\tau。前两节已经建立了由角速度求角动量的代数关系 \boldsymbol\ell=I\boldsymbol\omega；动力学接下来要回答的，正是合外力矩 \boldsymbol\tau 如何驱动角动量 \boldsymbol\ell 随时间变化。在惯性世界坐标系下，刚体关于质心的角动量定理（angular momentum theorem）写为：

向量分量与坐标轴的共同变化

在固连于刚体的局部坐标系中，惯量矩阵 I_b 是不随时间改变的常量，这为计算角动量提供了极大便利。但要注意的是，随体坐标系的基向量本身正随刚体一起转动，因此在对向量求时间导数时，必须同时计入体坐标分量的变化与随体坐标轴自身的转动。为了直观理解这两部分的配合，可以先考察一个在世界系中恒定不变的常向量 \boldsymbol y=\boldsymbol e_x：设刚体绕世界竖轴以常角速度 \omega 匀速旋转，偏航角为 \varphi=\omega t，则该常向量在随体坐标系中的投影分量及其时间导数分别为：

在初始时刻 \varphi=0，向量 \boldsymbol y 指向随体横轴正向，但它的体坐标分量变化率 \dot{\boldsymbol y}_b 却直接指向随体负纵轴 (0,-\omega,0)^{\mathsf T}。原因在于：随体坐标轴正向逆时针旋转时，空间中固定不动的向量相对于随体坐标轴表现为顺时针后退。只有把体坐标分量的表观变化率 R\dot{\boldsymbol y}_b 与坐标轴牵连转动带来的变化率 \boldsymbol\omega\times\boldsymbol y 相加，才能还原出世界系下的真实零导数。

将这一几何关系推广到任意随时间变化的向量 \boldsymbol y(t)。设其世界系表达与体坐标系表达通过姿态矩阵联系为 \boldsymbol y=R\boldsymbol y_b，对时间求导即可得到转动坐标系求导公式（transport theorem）：

在固连于刚体的随体坐标系中，关于质心的惯量矩阵 I_b 为常数矩阵。将体坐标角动量 \boldsymbol\ell_b=I_b\boldsymbol\omega_b 代入转动坐标系求导公式（4.11）与角动量定理（4.10），即可得到用随体分量表示的刚体转动动力学方程：

式（4.12）称为Euler 转动方程（Euler equations of rigid-body rotation）。方程右侧第一项 I_b\dot{\boldsymbol\omega}_b 描述体坐标角速度分量的瞬时变化，第二项 \boldsymbol\omega_b\times(I_b\boldsymbol\omega_b) 则称为陀螺项（gyroscopic term）或离心力矩项。从几何上看，当角动量 I_b\boldsymbol\omega_b 与角速度 \boldsymbol\omega_b 不共线时，哪怕刚体在随体坐标系下保持恒定的角速度分量（\dot{\boldsymbol\omega}_b=\boldsymbol0），角动量向量在世界空间中也会随刚体一起旋转，维持这种旋转所需施加的外力矩正是陀螺项 \boldsymbol\omega_b\times(I_b\boldsymbol\omega_b)；反之，当合外力矩为零（\boldsymbol\tau_b=\boldsymbol0）时，体坐标角速度就必须随时间发生变化，以使 I_b\dot{\boldsymbol\omega}_b 抵消陀螺项。

若刚体恰好绕某一根惯量主轴转动，则角动量 I_b\boldsymbol\omega_b 与角速度 \boldsymbol\omega_b 共线，陀螺叉乘项自动归零；此时若外力矩也始终沿该主轴方向，刚体就将一直保持绕该主轴的定轴转动，式（4.12）直接退化为第 1 节的标量方程 \tau_n=I_n\dot\omega。而在一般的三维转动中，角速度同时具有多个主轴分量，三个方向的方程便通过陀螺项非线性地耦合在一起。

世界坐标中的同一转动方程

如果将所有向量都改用世界坐标系表达，记世界角加速度为 \boldsymbol\alpha=\dot{\boldsymbol\omega}=R\dot{\boldsymbol\omega}_b，将式（4.12）两侧左乘 R，可得世界系下的 Euler 转动方程：

在世界系下，惯量矩阵 I=RI_bR^{\mathsf T} 随刚体姿态 R(t) 实时变化；直接对世界系角动量 \boldsymbol\ell=I\boldsymbol\omega 求导时，惯量矩阵的时间导数项 \dot I\boldsymbol\omega 恰好生成了同一个陀螺项 \boldsymbol\omega\times(I\boldsymbol\omega)。因此，体坐标表达与世界坐标表达提供了两条等价的计算路径：前者使用恒定的 I_b 配合转动系求导，后者则使用随姿态更新的 I(t) 直接在世界系中组装方程。

力矩功率与转动能

沿着刚体的真实运动轨迹对转动动能 T_{\mathrm{rot}}=\tfrac12\boldsymbol\omega_b^{\mathsf T}I_b\boldsymbol\omega_b 求时间导数，可以将外力矩做功的功率与转动动能的变化率联系起来。由于随体惯量矩阵 I_b 恒定且对称，转动动能的变化率恰好等于角速度向量与合外力矩向量的内积：

由于陀螺项与角速度恒正交，在无外力矩（\boldsymbol\tau_b=\boldsymbol0）的自由转动中，尽管角速度的三个体坐标分量随时间不断交换能量，它们按主惯量加权的平方和 2T_{\mathrm{rot}}=\sum_j I_j\omega_j^2 始终严格守恒（而角速度模长 \|\boldsymbol\omega_b\| 本身则会随分量重分配而波动）。此外，若刚体存在空间固定的支点 O，由于支点速度为零、刚体动能完全由绕支点的转动能 \tfrac12\boldsymbol\omega^{\mathsf T}I_O\boldsymbol\omega 构成，采用关于固定支点 O 的惯量矩阵与力矩同样满足完全相同形式的 Euler 方程与功率关系；下一节的圆盘稳态进动分析即建立在这一固定支点模型之上。

当刚体所受合外力矩为零（\boldsymbol\tau=\boldsymbol0）时，世界系下的角动量向量 \boldsymbol\ell 保持为常向量，同时由式（4.13）知转动动能 T_{\mathrm{rot}} 也严格守恒。在随体主轴坐标系下，体坐标角速度的三个分量 (\omega_1,\omega_2,\omega_3) 必须同时满足以下两个二次守恒方程：

在以 (\omega_1,\omega_2,\omega_3) 为坐标的三维角速度空间中，这两个守恒等式分别定义了两张中心在原点的椭球面（分别称为能量椭球与角动量椭球）。由于无外力矩转动必须同时满足两项守恒定律，体坐标角速度向量 \boldsymbol\omega_b(t) 的端点只能沿着这两张椭球面的交线运动；与此同时，刚体在三维空间中的姿态继续按第三章的四元数微分方程演进。

之所以第二个守恒式能够直接用随体角速度分量写出，是因为正交姿态变换保持向量模长不变，故随体角动量的模长 \|\boldsymbol\ell_b\|=\sqrt{\sum_j I_j^2\omega_j^2} 恒等于固定世界角动量的模长 \|\boldsymbol\ell\|。记刚体的初始转动动能为 T_0>0、初始角动量模长为 L_0>0；在三个主惯量均为正数的条件下，能量椭球与角动量椭球沿第 j 根坐标轴的半轴长度分别为 a_j=\sqrt{2T_0/I_j} 和 b_j=L_0/I_j。下方实验中的“守恒曲面”视图正是利用这两组半轴长度，在同一个体角速度空间内绘制出两张椭球面及其交线轨迹。

在下方实验中，先保留默认的“中间轴微扰”初始状态，对比左侧刚体的三维翻转运动、固定不动的方向箭头（世界角动量 \boldsymbol\ell）以及右侧的体角速度分量曲线；随后切换到“守恒曲面”视图，观察当前角速度状态点如何严格沿着两张椭球面的鞍部交线从正半轴滑向负半轴。

在默认的中间轴初始条件 \boldsymbol\omega_b(0)=(0.02,2,0.01)^{\mathsf T}\,\mathrm{rad/s} 下，角速度起初几乎完全沿第二主轴正向。随着时间推进，约在 t=5.93\,\mathrm s 时，第二主轴分量 \omega_2 穿过零点，而原本微小的第一、第三分量则增长到峰值；继续演化至 t=12\,\mathrm s 附近，\omega_2 到达约 -2\,\mathrm{rad/s}，标志着刚体的第二主轴相对于空间固定的角动量方向完成了 180^\circ 翻转（这一现象常被称为贾尼别科夫效应或网球拍定理）。从守恒曲面视图中可以清楚地看到，绕最小主轴（I_1）或最大主轴（I_3）附近的椭球交线是围绕主轴顶点的闭合小环，而经过中间主轴（I_2）附近的交线则是连接正负两个鞍点的大回路。下面我们通过线性化 Euler 方程，从数学上解释为何中间轴附近的微小扰动会呈指数放大。

中间主轴附近的运动

设刚体的三个主惯量互不相等且按大小排序为 0<I_1<I_2<I_3。当初始角速度紧邻最小主惯量轴（轴 1）或最大主惯量轴（轴 3）时，微小的横向角速度扰动只会围绕零值作有界简谐振荡；而当初始角速度紧邻中间主惯量轴（轴 2）时，微小扰动则会随时间指数增长，驱动刚体离开初始主轴并进入大角度周期性翻转（在无耗散的理想刚体系统中，状态点沿闭合交线周而复始地往返翻转）。

这一稳定性差异直接源于主惯量差值的正负号：设刚体主要绕中间轴以常角速度 W>0 转动，若此时第一分量出现微小正偏离（\omega_1>0）而第三分量为零，Euler 方程给出 \dot\omega_3=(I_1-I_2)W\omega_1/I_3<0，使 \omega_3 向负方向增长；而一旦 \omega_3<0，第一分量的方程 \dot\omega_1=(I_2-I_3)W\omega_3/I_1>0 又反过来推动 \omega_1 继续增大，两个扰动分量形成正反馈。相反，若基本转动绕最小主轴或最大主轴，两条横向线性化方程的右端系数异号，形成负反馈的恢复力。写成二阶线性微分方程，两种情形分别对应实指数增长与简谐振荡：

需要注意的是，上述线性化方程仅在扰动分量远小于基本转速 |W| 的初始阶段有效。恰好沿任意一根主轴（无任何横向扰动）的转动都是 Euler 方程的精确平衡解，但中间主轴属于不稳定鞍点：任何微小的初始扰动都会按 e^{\lambda t} 放大，待扰动增长到有限幅度后，便由完整的非线性 Euler 方程与双椭球守恒交线接管后续的翻转过程。

在实验 05 的数值实现中，程序将四元数微分方程与 Euler 转动方程联立，采用四阶 Runge–Kutta 方法（fourth-order Runge–Kutta method, RK4）推进，并在每个时间步末尾对四元数执行归一化。面板读数同时实时监测转动动能相对误差与世界角动量误差；默认积分步长为 h=0.01\,\mathrm s，读者可调大或调小步长以观察高阶积分的收敛特性。

固定支点的水平稳态进动

在无外力矩的自由转动中，世界角动量 \boldsymbol\ell 的大小与方向均固定不变；若对高速自转的刚体持续施加一个垂直于角动量水平分量的外力矩，角动量向量就会在空间中绕竖轴匀速旋转，形成陀螺进动（gyroscopic precession）。下面考察一个可精确求解的典型模型：用一根长为 d 的无质量刚性轻杆将质量为 m、半径为 r 的均匀薄圆盘中心连接到空间固定支点 O，并赋予圆盘满足水平稳态进动条件的初始角速度，使对称轴单位向量 \boldsymbol n 始终保持在水平面内。记圆盘绕自身对称轴 \boldsymbol n 的自转角速度为 \omega_s，对称轴绕世界竖轴 \boldsymbol e_z 旋转的进动角速度为 \Omega，圆盘关于固定支点 O 的轴向主惯量与横向主惯量分别为 I_s 与 I_t，则总角速度与关于支点 O 的总角动量分别为：

在重力作用下，重力关于支点 O 产生水平方向的力矩，持续改变角动量水平分量 I_s\omega_s\boldsymbol n 的指向；与此同时，轴向自转角速度 \omega_s\ne0 与竖直角动量分量 I_t\Omega\boldsymbol e_z 均保持恒定不变。由关于固定支点 O 的角动量定理 \dot{\boldsymbol\ell}_O=\boldsymbol\tau_O，可直接解得水平稳态进动的角速度为：

下方实验采用固定的 12\,\mathrm s 物理时间窗：将进度条固定在 t=2\,\mathrm s 处，对比两种不同的自转角速度——当 \omega_s=80\,\mathrm{rad/s} 时，圆盘轴在两秒内已绕竖轴扫过约 87.82^\circ；而当 \omega_s=160\,\mathrm{rad/s} 时，同样两秒内仅扫过约 43.91^\circ。

在三维视图中，重力力矩箭头 \boldsymbol\tau_O 始终沿着水平圆周的切线方向，实时指向角动量向量 \boldsymbol\ell_O 下一刻偏转的方向：圆盘自转越快，水平角动量分量 I_s\omega_s 越长，同样大小的水平力矩在单位时间内引起的偏转角 \Omega 就越小。同时可以看到，总角动量包含恒定的竖直分量 I_t\Omega\boldsymbol e_z，由于其时间导数为零，水平稳态进动关系直接由水平分量的旋转平衡唯一确定。

本节对照本地 MuJoCo 3.10.0 源码，梳理物理引擎中质量与惯量的数据结构及核心算法。相关模型字段定义在 mjmodel.h 中，空间惯量运算函数则位于 engine_util_spatial.c。

在静态模型结构 mjModel 中，body_mass 存储每个刚体的标量总质量；body_ipos 存储质心在刚体局部坐标系中的三维位置偏移；body_iquat 存储惯量主轴坐标系相对于刚体局部坐标系的姿态四元数；body_inertia 则存储该主轴坐标系下的三个主惯量 (I_1,I_2,I_3)。这种分离设计允许用户独立设定刚体几何网格的局部参考架与质心主轴架。

旋转与平行轴修正

函数 mju_inertCom 负责将刚体的主惯量按式（4.8）旋转到目标坐标架，并按平行轴定理式（4.9）平移参考点：它首先计算正交合同变换 R\operatorname{diag}(I_1,I_2,I_3)R^{\mathsf T}（代码中参数 mat 为从惯量主轴架到目标坐标架的旋转矩阵），接着根据从目标参考点指向刚体质心的位移向量 dif（即 \boldsymbol d）叠加上平行轴修正项：

res[0] += mass*(dif[1]*dif[1] + dif[2]*dif[2]);
res[1] += mass*(dif[0]*dif[0] + dif[2]*dif[2]);
res[2] += mass*(dif[0]*dif[0] + dif[1]*dif[1]);
res[3] -= mass*dif[0]*dif[1];
res[4] -= mass*dif[0]*dif[2];
res[5] -= mass*dif[1]*dif[2];

在输出数组 res 中，前 6 个浮点数依次存储对称惯量矩阵的三个对角元与三个非对角元 (I_{11},I_{22},I_{33},I_{12},I_{13},I_{23})，后 4 个浮点数则依次存储一阶质量矩 m\boldsymbol d 的三个分量以及标量总质量 m。这 10 个浮点数完整编码了下一章将要介绍的六维空间惯量。

在每一步正向运动学计算中，上层函数 mj_comPos 调用 mju_inertCom 生成运行期空间惯量 cinert：其中姿态矩阵传入主轴架的世界姿态 ximat，位移向量则取刚体质心世界位置 xipos 与其所属运动学树根部子树质心 subtree_com[body_rootid] 之差——也就是说，MuJoCo 将每棵运动学树上所有刚体的空间量统一归算到该树共同的子树质心参考点上。

惯量乘积与转动方程

函数 mju_mulInertVec 负责计算 10 分量空间惯量与 6 维空间运动向量的乘积。当参考点恰好选在刚体质心（\boldsymbol d=\boldsymbol0）时，质量矩 m\boldsymbol d 的三个分量为零，此时乘积的前三个转动分量直接退化为本章的 I\boldsymbol\omega 或 I\boldsymbol\alpha。

在逆动力学与偏置力计算函数 mj_rne 中，引擎先调用 mju_mulInertVec 计算空间惯量与空间加速度的乘积，再计算空间速度与空间动量的空间叉乘（即广义陀螺项），最后将两部分相加得到刚体所需的空间合力：

mju_mulInertVec(loc_cfrc_body+6*i, d->cinert+10*i, loc_cacc+6*i);
mju_mulInertVec(tmp, d->cinert+10*i, d->cvel+6*i);
mji_crossForce(tmp1, d->cvel+6*i, tmp);
mju_addTo(loc_cfrc_body+6*i, tmp1, 6);

当参考点位于质心时，上述四行代码的角分量精确对应于本章的 Euler 转动方程右端项 I\boldsymbol\alpha+\boldsymbol\omega\times(I\boldsymbol\omega)；而在一般参考点下，六维空间形式还会自动计入平动加速度与参考点偏移带来的耦合力矩。第五章将系统建立这套六维空间分块算子的代数结构，第八章则将其直接用于多体树的递推算法。

本章配套的 inertia.py 实现了离散质点系的质心与惯量求和、均匀长方体惯量公式、坐标旋转变换、平行轴修正以及 Euler 转动方程求解。其中，瞬时角加速度通过求解如下三阶线性方程组获得：

def angular_acceleration(inertia, omega, torque):
    gyro = cross(omega, apply(inertia, omega))
    return solve3(inertia, [torque[j] - gyro[j] for j in range(3)])

调用该函数时，输入的惯量矩阵 inertia、角速度 omega 与力矩 torque 必须采用同一个坐标系与同一个参考点：若三者均取随体坐标系下的质心量，则返回体坐标角加速度 \dot{\boldsymbol\omega}_b；若三者均取世界系下的质心量，则返回世界角加速度 \boldsymbol\alpha。由于物理刚体的惯量矩阵对称正定，该三阶线性系统恒有唯一解。在进行自由转动数值积分时，程序联合调用第三章的四元数更新函数：读者可先由各质点的速度直接累加动能，验证其与矩阵二次型结果一致；接着检验坐标旋转下的动能不变性以及参考点平移下的平行轴公式；最后在无外力矩转动中同时核查转动动能与世界角动量的数值守恒。

python3 examples/check_inertia.py

inertia = box_inertia(1, [2, 1.2, 0.5])
principal = [inertia[0], inertia[4], inertia[8]]
q, w = [1, 0, 0, 0], [0.02, 2, 0.01]
for _ in range(round(12 / h)):
    q, w = step_free(q, w, principal, h)
energy = rotational_energy(inertia, w)
momentum_world = apply(to_matrix(q), apply(inertia, w))