3. 多体系统动力学方程建立

多体系统动力学的核心任务，是把由多个刚体或柔性体通过铰链、滑移副、接触副连接而成的机械系统，写成一组可数值求解的微分方程或微分-代数方程。对于接触碰撞仿真而言，方程建立方式直接决定了能否稳定地处理碰撞瞬间的冲击、约束突变与能量交换。本节围绕拉格朗日方程、虚功原理、质量矩阵、广义力与约束反力五个关键概念，说明它们如何共同构成多体系统的动力学方程。

一、从虚功原理到拉格朗日方程

虚功原理是建立动力学方程最本质的出发点。它指出：系统在任意虚位移下，惯性力、外力与约束力所做虚功之和为零。设系统用广义坐标 q 描述，则虚位移 δq 满足约束相容条件。对每个刚体，惯性力虚功可写为 δW_i = -δr_i^T m_i r̈_i，其中 r_i 为质心位置。将所有刚体求和并代入 r_i = r_i(q)，利用链式法则 ṙ_i = J_i q̇，即可把虚功原理转化为广义坐标下的形式。

拉格朗日方程正是虚功原理在完整约束系统下的等价表达。定义动能 T = ½ q̇^T M(q) q̇，势能 V(q)，拉格朗日函数 L = T - V，则第二类拉格朗日方程为：

d/dt(∂L/∂q̇) - ∂L/∂q = Q_nc

其中 Q_nc 为非保守广义力，包括阻尼力、驱动力、接触力等。展开后可得标准形式：

M(q) q̈ + C(q, q̇) q̇ + Φ_q^T λ = Q_e

这里 M(q) 为质量矩阵，C(q, q̇)q̇ 为科氏力与离心力项，Φ_q^T λ 为约束反力项，Q_e 为外加广义力。该式即为多体系统动力学方程的一般形式。

二、质量矩阵的构造与物理意义

质量矩阵 M(q) 是动能二次型的系数矩阵，对 n 个自由度的系统为 n×n 对称正定矩阵。对单个刚体，若广义坐标为质心位置与姿态角，则质量矩阵由平动质量 m 与转动惯量张量 I 组成。当多个刚体通过铰链连接时，质量矩阵通过雅可比矩阵 J_i 组装：

M(q) = Σ J_i^T M_i J_i

其中 M_i = diag(m_i, m_i, m_i, I_ix, I_iy, I_iz)。质量矩阵随位形变化，这是多体系统与线性振动系统的本质区别。在接触碰撞中，质量矩阵的时变性会导致碰撞瞬间加速度突变，因此数值积分必须采用能处理刚性问题的隐式方法，如 Newmark-β 或广义-α 法。

三、广义力与约束反力

广义力 Q 是把所有外力、驱动力、阻尼力、接触力投影到广义坐标方向的结果。对作用在点 p 的外力 F，其广义力为 Q = J_p^T F。接触碰撞中，法向接触力常用罚函数模型 F_n = k δ + c δ̇，其中 δ 为穿透量，k 为接触刚度，c 为阻尼系数。切向力则用库仑摩擦模型 F_t = μ F_n 近似。

约束反力由拉格朗日乘子 λ 给出，Φ_q^T λ 表示约束方程 Φ(q,t)=0 对应的约束力。对理想约束，约束反力不做虚功，因此不出现在虚功原理的右端。但在接触碰撞中，接触约束是单边约束，需用互补条件 Φ ≥ 0, λ ≥ 0, Φ^T λ = 0 描述，这使方程从纯微分方程变为微分-代数方程（DAE）。

四、接触碰撞中的方程特性与求解要点

接触碰撞使多体系统动力学方程呈现三个显著特征。第一，质量矩阵虽仍对称正定，但接触刚度 k 通常取 10^6~10^8 N/m 量级，导致方程刚性极强，显式积分步长需小于 2/ω_max，实际中往往小到 10^-7 s 以下，因此工程上更常用隐式或半隐式方法。第二，约束反力 λ 在碰撞瞬间发生阶跃，若用拉格朗日乘子法直接求解，需同时满足位移约束与速度约束，容易产生数值振荡；罚函数法通过允许微小穿透 δ 把约束转化为力，回避了乘子求解，但需控制 k 与 c 的匹配，通常取 c = 2ζ√(k m_eff)，ζ 取 0.1~0.3。第三，广义力中的接触力依赖穿透量与穿透速度，使方程右端成为状态量的非线性函数，需在每个时间步内迭代求解。

以一个典型示例说明：两刚体质量分别为 m₁ = 2 kg、m₂ = 3 kg，初始相对速度 v_rel = 1 m/s，接触刚度 k = 1×10^7 N/m，阻尼 c = 200 N·s/m。取有效质量 m_eff = m₁m₂/(m₁+m₂) = 1.2 kg，则接触频率 ω_n = √(k/m_eff) ≈ 2887 rad/s，对应周期约 2.18 ms。若采用显式积分，步长应小于 0.1 ms；若采用 Newmark-β 隐式格式，步长可放宽到 0.5 ms 左右，但需在每个步内更新 M(q) 与 Φ_q。这些参数关系说明：质量矩阵、广义力与约束反力并非独立概念，而是共同决定方程刚性与求解策略的整体。

五、方程建立中的常见误区

• 忽略质量矩阵的位形相关性：把 M 当作常数矩阵，会导致大范围转动时动能计算错误，碰撞响应失真。
• 混淆约束反力与外力：约束反力由乘子 λ 产生，不能直接并入 Q_e，否则约束方程与动力学方程耦合关系被破坏。
• 接触刚度取值过大或过小：过大导致刚性过强、步长过小；过小导致穿透明显、物理失真。建议先用 k = 10^6~10^7 N/m 试算，再根据穿透量调整。
• 阻尼系数与刚度不匹配：阻尼过大会使碰撞能量被过度耗散，过小则碰撞后持续振荡。按 c = 2ζ√(k m_eff) 估算较为稳妥。

综上，多体系统动力学方程的建立，本质上是把虚功原理通过广义坐标展开为拉格朗日方程，再显式写出质量矩阵、广义力与约束反力三部分。对接触碰撞问题，还需引入单边约束与互补条件，使方程成为微分-代数形式。理解这一结构，才能在后续数值求解与仿真调试中准确判断问题来源。