2. 多体系统动力学建模基础

多体系统动力学的建模过程，本质上是一条从“几何描述”到“方程生成”再到“数值求解”的完整链路。摩擦与间隙非线性之所以难处理，很大程度上是因为它们会同时侵入这条链路的多个环节：间隙改变了约束的拓扑与接触状态，摩擦引入了依赖法向力与相对速度的非线性力，而接触状态的切换又会让方程在连续与不连续之间反复跳变。因此，在讨论摩擦间隙之前，必须先把多体系统建模的基本框架理清楚。

广义坐标：用什么变量描述系统位形

描述一个多体系统，首先要选定一组能够唯一确定系统位形的变量，这就是广义坐标。广义坐标不一定是笛卡尔坐标，也不一定彼此独立。常见的选取方式有两类：

• 绝对坐标（笛卡尔坐标）：为每个刚体赋予 6 个坐标（三维空间中 3 个平动 + 3 个转动），系统坐标总数等于 6×刚体数。优点是约束方程形式简单、易于程序化生成；缺点是坐标数多、约束方程多。
• 相对坐标（铰坐标）：以铰链连接关系为骨架，用相邻刚体之间的相对运动量作为坐标。优点是坐标数少、约束少；缺点是约束方程和递推关系更复杂。

在本课程中，为便于统一处理间隙与摩擦，建议采用绝对坐标 + 约束方程的建模路线。设系统广义坐标向量为 q = [q₁, q₂, …, qₙ]ᵀ，其中每个刚体贡献位置向量 r 与姿态参数（欧拉角或四元数）。姿态参数的选择会影响后续约束方程的雅可比形式，若使用欧拉角需注意万向节奇异，若使用四元数则需额外引入归一化约束。

约束方程：把“连接关系”写成数学条件

约束方程描述的是系统内部必须始终满足的几何或运动条件。对多体系统而言，约束通常写成如下形式：

Φ(q, t) = 0        （完整约束，如铰链、滑轨、球铰）
Φ(q, q̇, t) = 0     （非完整约束，如纯滚动）

以平面转动铰为例，若刚体 i 与刚体 j 在点 P 处铰接，则约束可写为两个刚体上 P 点位置重合：

Φ = r_i + A_i·s_i^P − r_j − A_j·s_j^P = 0

其中 A 为姿态矩阵，s^P 为铰点在各自局部坐标系中的位置。对时间求导可得速度约束与加速度约束：

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

这里 Φ_q 是约束雅可比矩阵，γ 是加速度约束右端项。间隙的存在意味着约束方程并非固定不变：当接触点分离时，该约束暂时失效；当接触点闭合时，约束重新激活。因此间隙建模往往需要引入“接触状态变量”，在约束激活与失效之间切换。

运动学递推：从坐标到速度、加速度

在绝对坐标框架下，运动学递推的核心是：已知 q 和 q̇，如何得到各刚体的位置、姿态、速度、角速度，以及后续动力学所需的惯性量。对每个刚体 i：

• 位置：r_i = r_i(q)
• 姿态矩阵：A_i = A_i(q)
• 速度：ṙ_i = ∂r_i/∂q · q̇
• 角速度：ω_i = G_i(q)·q̇，其中 G_i 为角速度映射矩阵

若采用相对坐标，则递推沿铰链树从基座向末端逐级进行，父刚体的运动作为子刚体的输入，逐级叠加相对运动。两种路线在数学上等价，但绝对坐标更适合处理含间隙的接触约束，因为接触点位置可以直接由绝对坐标表达，无需沿树回溯。

动力学方程：从牛顿-欧拉方程到统一形式

多体系统动力学方程的标准形式为：

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

其中 M 为质量矩阵，λ 为拉格朗日乘子（对应约束力），Q 为广义力（含外力、摩擦力、接触力等）。这是一个微分-代数方程组（DAE），而非纯微分方程。摩擦与间隙非线性主要体现在 Q 中：

• 摩擦力：F_f = μ(v_rel)·F_n，其中 F_n 来自接触约束力或接触模型，v_rel 为接触点相对速度，μ 可随速度变化。
• 间隙接触力：常用罚函数模型 F_n = k·δ + c·δ̇（δ 为穿透量），或采用拉格朗日乘子/互补条件精确处理。

若采用罚函数法，约束不再严格满足，方程退化为纯 ODE，便于数值求解，但需要选择足够大的刚度 k 与合理的阻尼 c，否则会出现数值刚性问题。若采用乘子法，则需处理 DAE 的指标问题，通常将位置约束、速度约束、加速度约束联合使用，并配合 Baumgarte 稳定化或投影法抑制约束漂移。

数值求解流程：从方程到可执行的计算步骤

将上述方程交给数值求解器，通常遵循以下流程。下面给出一个可参考的实现框架，参数为示例值，用于说明流程而非实测结果。

1. 初始化：q0, q̇0, t=0, 步长 h=1e-4
2. 装配 M(q), Φ_q(q), Q(q, q̇, t)
3. 计算加速度约束右端项 γ
4. 求解线性系统：
   [ M   Φ_qᵀ ] [ q̈ ]   [ Q ]
   [ Φ_q  0   ] [ λ  ] = [ γ ]
5. 积分更新：q̇ ← q̇ + h·q̈,  q ← q + h·q̇
6. 约束修正：投影 q 到 Φ(q)=0 流形上
7. 接触状态检测：判断间隙是否闭合，更新接触力与摩擦状态
8. t ← t + h，返回步骤 2，直到 t ≥ t_end

其中第 4 步的线性系统是每个时间步的核心计算量。若系统规模较大，可利用 M 的稀疏结构或采用递归算法降低复杂度。第 6 步的约束修正非常关键：由于数值积分误差，Φ(q) 会逐渐偏离零，若不修正，约束会“漂移”，导致接触点位置失真，进而影响摩擦力方向与大小。

建模路线与摩擦间隙的衔接

把上述五个环节串起来看，摩擦间隙非线性分析的实际难点并不在单个方程，而在于它们之间的耦合：广义坐标的选择决定了约束方程的形式，约束方程决定了接触点如何表达，接触点状态决定了摩擦力是否激活，而摩擦力的非线性又反过来影响动力学方程的数值行为。因此，在后续章节中处理摩擦与间隙时，应始终回到这条链路，检查每一步的假设是否仍然成立。

上图把广义坐标、约束方程、运动学递推、动力学方程与数值求解流程组织为一条主链路，并标出摩擦间隙非线性在其中的三个侵入点：约束激活/失效、广义力中的非线性摩擦项、以及接触状态切换带来的不连续性。理解这条链路，是后续开展摩擦间隙非线性分析的基础。