4. 多体系统动力学建模
多体系统动力学建模的核心任务，是把由多个刚体或柔性体通过运动副连接的装备，转化为一组可被计算机求解的微分-代数方程。与单刚体相比，多体系统的难点不在单个物体的受力，而在于约束如何进入方程、惯性力如何随位形变化、外力与力矩如何映射到广义坐标，以及求解器能否在刚性系统中保持数值稳定。本节围绕这四件事展开，并以一个两连杆机械臂的简化模型作为贯穿示例。
一、从牛顿-欧拉方程到广义坐标方程
设系统由 n 个刚体组成，第 i 个刚体的质量为 mi，质心位置为 ri，惯性张量为 Ji。牛顿-欧拉方程给出每个刚体的平动与转动平衡：
m_i · a_i = F_i^ext + F_i^con
J_i · α_i + ω_i × (J_i · ω_i) = M_i^ext + M_i^con其中 Ficon、Micon 是运动副约束力与力矩。直接对每个刚体列方程会引入大量未知约束力，因此工程上更常用拉格朗日方法：选取广义坐标 q（例如两连杆机械臂取 q = [θ₁, θ₂]），将约束力消去，得到
M(q)·q̈ + C(q, q̇)·q̇ + G(q) = Q_ext
M(q) = Σ (J_i^T · m_i · J_i + J_ωi^T · J_i · J_ωi)这里 M(q) 是质量矩阵，C(q, q̇) 包含科氏力与离心力项，G(q) 是重力项，Qext 是广义外力。质量矩阵随位形变化，这正是多体系统与单刚体最本质的区别：惯性力不是常数，而是位形与速度的函数。
二、惯性力的物理含义与计算
惯性力并非真实的外界作用力，而是把加速度项移到方程一侧后产生的等效项。在 M(q)·q̈ 中，M(q) 反映当前位形下各广义坐标之间的耦合；在 C(q, q̇)·q̇ 中，科氏力与离心力反映速度引起的非线性效应。以两连杆机械臂为例，当第二关节以较高角速度旋转时，第一关节会感受到明显的离心力矩，其大小与 θ̇₂² 成正比，方向沿关节切向。若忽略该项，仿真在高速段会出现轨迹偏差。
实际建模时，惯性力的计算依赖两个关键输入：各刚体的质量、质心位置与惯性张量，以及运动副的雅可比矩阵。雅可比矩阵把关节速度映射为质心速度，从而把物理空间的惯性效应投影到广义坐标空间。因此，质量参数和运动学模型必须一致，否则惯性力计算会出现系统性误差。
三、外力与力矩的施加方式
外力与力矩的施加，本质上是把物理空间中的力旋量映射为广义力。若在刚体 i 的局部点 p 施加力 F，则对广义坐标的贡献为
Q_ext = J_p^T · F
J_p = ∂r_p / ∂q其中 Jp 是施加点的雅可比矩阵。对于力矩 M，则通过角速度雅可比 Jω 映射：Qext = JωT·M。工程中常见的外力包括重力、接触力、驱动力矩、气动力与弹簧阻尼力。以两连杆机械臂为例，若在末端施加 50 N 的竖直向下力，需要先计算末端雅可比，再将其转置后乘以力向量，得到两个关节的等效力矩。若直接在关节上施加力矩，则无需雅可比映射，但要注意力矩方向与关节正方向的一致性。
四、求解器选择与数值稳定性
多体系统方程通常写成 M(q)·q̈ = Q(q, q̇, t) 的形式，属于二阶微分方程。求解器分为显式与隐式两大类。显式方法（如 RK4、ode45）每步计算量小，但对刚性系统步长受限；隐式方法（如 BDF、ode15s、Radau）稳定性好，适合含接触、摩擦或高刚度弹簧的系统。选择依据可归纳为三点：系统刚性程度、约束类型、以及是否需要长时间稳定仿真。
数值稳定性问题主要来自三方面：一是质量矩阵病态，当某些刚体质量或惯量相差几个数量级时，M(q) 条件数恶化；二是约束漂移，位置约束在积分过程中逐渐被破坏；三是刚性项，如高刚度接触力导致方程特征值跨度大。应对措施包括：使用隐式积分器、对质量矩阵进行缩放、采用 Baumgarte 稳定化或投影法修正约束、以及合理设置相对误差容差（如 1e-6）与绝对误差容差（如 1e-8）。
五、两连杆机械臂示例：参数与结果判读
设连杆1质量 m₁ = 2.0 kg，长度 l₁ = 0.5 m，绕质心惯量 J₁ = 0.04 kg·m²；连杆2质量 m₂ = 1.0 kg，长度 l₂ = 0.4 m，J₂ = 0.012 kg·m²。初始位形 θ₁ = 0°，θ₂ = 0°，初始角速度均为 0。在末端施加 50 N 竖直向下力，仿真时长 2 s。若采用 ode45，步长自动调整，相对容差 1e-6；若采用 ode15s，则适合后续加入接触力的场景。
结果判读关注三点：一是关节角轨迹是否平滑，若出现高频振荡，说明步长或容差设置不当；二是末端位置误差是否随时间累积，若持续增大，说明约束漂移未被抑制；三是能量曲线是否单调衰减，若能量异常增长，说明数值不稳定。以本示例为例，在 0～0.5 s 内 θ₁ 从 0° 变化到约 12°，θ₂ 从 0° 变化到约 -8°，末端竖直位移约 0.06 m，与静力估算量级一致，说明外力映射与惯性力计算基本正确。
六、建模检查清单
| 检查项 | 判断依据 | 常见问题 |
|---|---|---|
| 质量与惯量参数 | 与 CAD 或实测一致 | 单位混用、质心偏移 |
| 运动副约束 | 自由度数目正确 | 过约束或欠约束 |
| 外力映射 | 雅可比矩阵维度匹配 | 方向正负号错误 |
| 求解器选择 | 刚性程度与容差匹配 | 显式求解刚性系统 |
| 能量与约束 | 能量不异常增长，约束误差收敛 | 漂移未抑制 |
多体系统动力学建模的实质，是把物理对象的质量、约束与受力，系统地翻译为广义坐标下的方程，再交给合适的求解器积分。惯性力、外力与力矩的施加决定了方程的右端项，求解器与稳定性策略决定了方程能否被可靠地解出来。三者缺一不可，任何一环的参数或约定出错，都会在仿真结果中以轨迹偏差、振荡或发散的形式表现出来。