2. 多体系统动力学基础

在 Simpack 中建立机械系统模型，本质上是在构造一个多体系统（Multibody System, MBS）的数学模型。理解广义坐标、约束方程、动力学方程与数值求解原理，是判断模型能否正确求解、结果是否可信的基础。本节从建模对象出发，逐步说明这些概念在 Simpack 中的具体含义与相互关系。

一、广义坐标：描述系统位形的独立变量

多体系统由若干刚体（或柔性体）通过运动副连接而成。要完整描述系统在任意时刻的位形，需要一组能够确定所有物体位置和姿态的变量，这组变量就是广义坐标。在三维空间中，一个自由刚体有 6 个自由度：3 个平动（质心位置 x、y、z）和 3 个转动（姿态角，如欧拉角或卡尔丹角）。若系统包含 n 个刚体，则最多有 6n 个广义坐标。

Simpack 中每个刚体默认使用 6 个广义坐标，通常以质心位置和一组姿态角表示。姿态角的定义方式（如 3-1-3 欧拉角或 1-2-3 卡丹角）会影响转动方程的形式，建模时需与参考坐标系约定保持一致。例如，若采用 1-2-3 卡丹角，则绕 x、y、z 轴的转动依次为 ψ、θ、φ，角速度与角速度分量之间的关系不再是简单的导数关系，而需通过变换矩阵联系。

广义坐标的选择直接影响约束方程的复杂程度。对于有相对运动的物体，若直接使用绝对坐标，约束方程数量多但形式统一；若使用相对坐标（如铰链坐标），约束方程少但需处理坐标变换。Simpack 内部通常采用绝对坐标与相对坐标混合的方式，用户通过运动副类型间接指定。

二、约束方程：限制系统运动的代数关系

实际机械系统中，物体之间通过运动副连接，运动受到限制。这些限制用约束方程表示，通常写成广义坐标的代数方程：

Φ(q, t) = 0

其中 q 为广义坐标列阵，t 为时间。约束方程分为完整约束和非完整约束。完整约束只限制位置，如旋转铰限制两个物体某点重合且相对转动只有一个自由度；非完整约束限制速度但不可积分为位置关系，如车轮纯滚动。Simpack 中绝大多数运动副产生完整约束，少数如齿轮接触、轮胎模型可能涉及非完整约束。

约束方程的数量决定了系统的自由度。若系统有 n 个广义坐标、m 个独立约束方程，则系统自由度数为 n − m。例如，一个由两个刚体通过旋转铰连接的系统，n = 12，旋转铰提供 5 个约束，自由度数为 7。Simpack 在建模时会自动统计约束并检查是否过约束或欠约束。

约束方程对时间的导数给出速度约束和加速度约束：

Φ_q · q̇ = −Φ_t
Φ_q · q̈ = −(Φ_q · q̇)_q · q̇ − 2Φ_qt · q̇ − Φ_tt

其中 Φ_q 为约束雅可比矩阵。这些关系在数值求解中用于修正位置和速度，避免积分漂移。

三、动力学方程：从牛顿-欧拉到拉格朗日

多体系统的动力学方程通常由拉格朗日方程或牛顿-欧拉方法导出。对于受约束系统，引入拉格朗日乘子 λ，得到第一类拉格朗日方程：

M(q) · q̈ + Φ_q^T · λ = Q(q, q̇, t)
Φ(q, t) = 0

其中 M 为质量矩阵，Q 为广义力（含外力、弹簧阻尼力、驱动力等），Φ_q^T · λ 为约束力。这是一个微分-代数方程组（DAE），而非纯微分方程。Simpack 求解器正是处理这类 DAE 系统。

质量矩阵 M 由各刚体的质量和惯量张量组装而成。若刚体姿态用欧拉角描述，M 中会出现与姿态相关的项，使方程非线性。广义力 Q 包括：重力、弹簧力、阻尼力、接触力、驱动约束等。在 Simpack 中，力元（Force Element）和运动副（Joint）分别贡献 Q 和 Φ。

四、数值求解原理：积分与约束稳定

Simpack 对 DAE 的求解通常采用隐式积分法，如 BDF（Backward Differentiation Formula）或 Radau IIA。这类方法在每一步迭代中同时求解广义坐标和拉格朗日乘子，适合刚性系统。求解过程可概括为：

1. 根据当前 q、q̇，组装 M、Q、Φ_q；
2. 用隐式公式预测下一时刻的 q、q̇；
3. 迭代求解非线性方程组，使动力学方程和约束方程同时满足；
4. 检查残差和约束违反量，若超限则缩小步长或调整迭代。

约束稳定是 DAE 求解的关键。由于数值误差，Φ(q, t) = 0 可能逐渐被违反。Simpack 采用 Baumgarte 稳定化或投影法，将位置和速度约束的偏差反馈到方程中，抑制漂移。Baumgarte 方法将约束方程改写为：

Φ̈ + 2αΦ̇ + β²Φ = 0

其中 α、β 为稳定参数，通常取 5～20。参数过小则约束漂移明显，过大则可能引起高频振荡。Simpack 默认值经过调校，一般无需修改；若模型出现约束违反警告，可适当调整。

五、从方程到 Simpack 模型：一个简单示例

以一个单摆为例：质量 m = 1 kg 的质点，用长 L = 0.5 m 的无质量杆悬挂，重力加速度 g = 9.81 m/s²。若用绝对坐标描述，质点位置 (x, y)，约束方程为 x² + y² = L²。广义坐标 q = [x, y]^T，约束雅可比 Φ_q = [2x, 2y]。动力学方程为：

m·ẍ = λ·2x
m·ÿ = λ·2y − m·g
x² + y² = L²

在 Simpack 中，可建立一个刚体（质点）和一个旋转铰（连接大地与质点），旋转铰自动提供约束。设置重力为 −y 方向，初始角度 30°，初始角速度 0。求解器采用默认的 Radau IIA，步长 0.001 s，仿真 2 s。结果中可观察摆角随时间变化，周期约为 2π√(L/g) ≈ 1.42 s。若约束违反量保持在 1e-8 以下，说明求解稳定。

从上述关系可见，广义坐标是描述系统位形的起点，约束方程限定了运动范围，动力学方程给出运动规律，数值求解器则负责在计算机上逐步推进。四者环环相扣：广义坐标选得合适，约束方程简洁；约束方程正确，动力学方程才能反映真实受力；数值求解稳定，结果才可信。在 Simpack 中，这些环节大多由软件自动处理，但理解其原理有助于判断模型是否合理、求解是否收敛、结果是否可用。