第5章：接触力学模型

章节大纲

5.1 库仑摩擦与LuGre模型

5.1.1 经典库仑摩擦定律
5.1.2 静摩擦与动摩擦的转换
5.1.3 LuGre动态摩擦模型
5.1.4 摩擦参数的物理意义

5.2 赫兹接触理论

5.2.1 弹性接触的基本假设
5.2.2 接触面积与压力分布
5.2.3 接触刚度的计算
5.2.4 非赫兹接触情况

5.3 柔性接触模型

5.3.1 罚函数方法
5.3.2 Hunt-Crossley模型
5.3.3 阻尼与能量耗散
5.3.4 接触参数的选择

5.4 粘滑转换与摩擦锥

5.4.1 粘滑现象的物理机制
5.4.2 摩擦锥的几何表示
5.4.3 最大耗散原理
5.4.4 各向异性摩擦

开篇段落

接触力学是物理仿真中最具挑战性的领域之一。当两个物体接触时，会产生复杂的力学相互作用，包括法向力、切向摩擦力以及可能的粘附力。这些力的准确建模对于实现逼真的物理仿真至关重要，特别是在机器人抓取、操作和运动控制等任务中。本章将深入探讨各种接触力学模型，从经典的库仑摩擦定律到现代的动态摩擦模型，从理想的赫兹接触理论到实用的柔性接触方法。我们将特别关注这些模型在强化学习环境中的应用，分析不同模型对训练稳定性和策略迁移的影响。

学习目标

完成本章学习后，您将能够：

理解并实现各种摩擦模型，包括静态和动态摩擦
掌握赫兹接触理论的数学基础及其适用范围
设计适合特定任务的柔性接触模型
处理粘滑转换等接触中的非线性现象
为强化学习任务选择合适的接触参数
理解接触模型对仿真稳定性和真实性的影响

5.1 库仑摩擦与LuGre模型

5.1.1 经典库仑摩擦定律

库仑摩擦定律是接触力学中最基本的模型，由Charles-Augustin de Coulomb在1785年通过大量实验总结得出。尽管该模型做了许多简化假设，但其简洁性和实用性使其成为物理仿真中最常用的摩擦模型。

库仑摩擦的核心思想是将摩擦力分为静摩擦和动摩擦两种状态：

$$\mathbf{f}_t = \begin{cases} -\mathbf{f}_{\text{applied}} & \text{if } |\mathbf{f}_{\text{applied}}| \leq \mu_s f_n \text{ (静摩擦)} \\ -\mu_k f_n \frac{\mathbf{v}_t}{|\mathbf{v}_t|} & \text{if } |\mathbf{v}_t| > 0 \text{ (动摩擦)} \end{cases}$$ 其中：

$\mathbf{f}_t$ 是切向摩擦力
$f_n$ 是法向力大小
$\mu_s$ 是静摩擦系数
$\mu_k$ 是动摩擦系数（通常 $\mu_k < \mu_s$）
$\mathbf{v}_t$ 是相对切向速度

这个看似简单的模型在数值实现时会遇到几个挑战：

不连续性问题：在静动转换点，摩擦力会发生突变，导致数值积分不稳定
零速度奇异性：当 $|\mathbf{v}_t| \to 0$ 时，摩擦力方向变得不确定
静态不确定性：在静摩擦状态下，摩擦力可以取范围内的任意值

5.1.2 静摩擦与动摩擦的转换

为了解决库仑摩擦的不连续性问题，实际仿真中常采用正则化方法。最常见的是速度正则化： $$\mathbf{f}_t = -\mu(v) f_n \frac{\mathbf{v}_t}{\max(|\mathbf{v}_t|, v_{\text{min}})}$$ 其中摩擦系数是速度的函数： $$\mu(v) = \mu_k + (\mu_s - \mu_k) e^{-\alpha |v|}$$ 这种指数衰减模型能够平滑地过渡静动摩擦，参数 $\alpha$ 控制过渡的陡峭程度。在强化学习应用中，这种平滑过渡对于梯度计算尤为重要。

另一种更精确的方法是使用隐式积分器配合互补性条件： $$\begin{aligned} 0 \leq \lambda_n &\perp g(\mathbf{q}) \geq 0 \quad \text{(法向)} \\ |\mathbf{f}_t| &\leq \mu \lambda_n \quad \text{(摩擦锥约束)} \\ \mathbf{f}_t &= -\mu \lambda_n \frac{\mathbf{v}_t + \Delta\mathbf{v}_t}{|\mathbf{v}_t + \Delta\mathbf{v}_t|} \quad \text{if sliding} \end{aligned}$$

5.1.3 LuGre动态摩擦模型

LuGre模型由Canudas de Wit等人在1995年提出，它通过引入微观鬃毛变形的概念，将摩擦建模为动态系统： $$\begin{aligned} \dot{z} &= v - \sigma_0 \frac{|v|}{g(v)} z \\ g(v) &= \mu_c + (\mu_s - \mu_c) e^{-(v/v_s)^2} \\ f &= \sigma_0 z + \sigma_1 \dot{z} + \sigma_2 v \end{aligned}$$ 其中：

$z$ 是平均鬃毛变形（内部状态变量）
$\sigma_0$ 是鬃毛刚度
$\sigma_1$ 是微观阻尼系数
$\sigma_2$ 是粘性摩擦系数
$v_s$ 是Stribeck速度

LuGre模型的优势在于：

能够捕捉预滑移现象
自然地处理粘滑转换
包含速度依赖的摩擦特性（Stribeck效应）

        摩擦力
          ^
          |     
    μs fn |______
          |      \___
    μk fn |          \_______ Stribeck曲线
          |                    
          |____________________> 速度
                   vs

5.1.4 摩擦参数的物理意义

理解摩擦参数的物理意义对于参数辨识和仿真调试至关重要：

静摩擦系数 $\mu_s$：

物理意义：表面微观互锁的强度
典型值：钢-钢(干)：0.6-0.8；橡胶-混凝土：0.7-1.0
RL影响：过大会导致机器人"粘"在地面，过小则抓取困难

动摩擦系数 $\mu_k$：

物理意义：滑动时的能量耗散率
典型值：通常为静摩擦的60-80%
RL影响：影响滑动控制的稳定性

Stribeck速度 $v_s$：

物理意义：边界润滑到流体润滑的转换速度
典型值：0.001-0.1 m/s
RL影响：影响低速运动的精度

在强化学习训练中，摩擦参数的域随机化策略： $$\begin{aligned} \mu_s &\sim \mathcal{U}[0.5\mu_{s,\text{nom}}, 1.5\mu_{s,\text{nom}}] \\ \mu_k &\sim \mathcal{U}[0.6\mu_s, 0.9\mu_s] \\ v_s &\sim \log\mathcal{U}[10^{-4}, 10^{-1}] \end{aligned}$$

5.2 赫兹接触理论

5.2.1 弹性接触的基本假设

赫兹接触理论由Heinrich Hertz在1882年建立，是弹性接触力学的基石。该理论基于以下假设：

小变形假设：接触区域远小于物体的特征尺寸
弹性假设：材料在弹性极限内，遵循胡克定律
光滑表面：忽略表面粗糙度的影响
无摩擦：只考虑法向接触，不考虑切向力
半空间假设：物体可视为弹性半空间

在这些假设下，赫兹推导出了著名的接触力-变形关系： $$F = \frac{4}{3} E^* \sqrt{R} \delta^{3/2}$$ 其中：

$E^*$ 是等效弹性模量：$\frac{1}{E^*} = \frac{1-\nu_1^2}{E_1} + \frac{1-\nu_2^2}{E_2}$
$R$ 是等效曲率半径：$\frac{1}{R} = \frac{1}{R_1} + \frac{1}{R_2}$
$\delta$ 是法向压入深度
$\nu_i$ 是泊松比

这个3/2次幂关系是赫兹理论的标志性特征，与线性弹簧模型形成鲜明对比。

5.2.2 接触面积与压力分布

赫兹理论不仅给出了力-变形关系，还预测了接触区域的几何形状和压力分布。对于球形接触，接触区域是半径为 $a$ 的圆形： $$a = \sqrt{R\delta}$$ 接触压力呈半椭球分布： $$p(r) = p_0 \sqrt{1 - \left(\frac{r}{a}\right)^2}$$ 其中最大接触压力为： $$p_0 = \frac{3F}{2\pi a^2} = \frac{3}{2} \left(\frac{E^*}{\pi}\right)^{2/3} \left(\frac{F}{R^2}\right)^{1/3}$$ 这种非均匀的压力分布对于理解接触应力集中和疲劳失效至关重要。在机器人抓取任务中，了解压力分布有助于：

预测接触损伤
优化抓取力分配
设计柔顺控制策略

5.2.3 接触刚度的计算

接触刚度是仿真中的关键参数，它决定了接触的"硬度"。赫兹接触的刚度是非线性的： $$k = \frac{dF}{d\delta} = 2E^* \sqrt{R\delta}$$ 注意刚度随压入深度增加，这与线性弹簧截然不同。对于小振动分析，可以在工作点线性化： $$k_{\text{lin}} = 2E^* \sqrt{R\delta_0}$$ 在数值仿真中，接触刚度的选择需要平衡：

过大的刚度：导致数值刚性，需要极小的时间步长
过小的刚度：产生过度穿透，失去物理真实性

经验公式用于估算合适的接触刚度： $$k_{\text{sim}} = \alpha \cdot \frac{m \cdot g}{d_{\text{tol}}}$$ 其中 $\alpha \approx 10-100$，$d_{\text{tol}}$ 是允许的穿透深度。

5.2.4 非赫兹接触情况

实际接触常常偏离赫兹假设，需要修正或替代模型：

粘附接触（JKR模型）

当表面能不可忽略时，Johnson-Kendall-Roberts模型考虑粘附力： $$F = \frac{4E^*a^3}{3R} - \sqrt{8\pi\gamma E^* a^3}$$ 其中 $\gamma$ 是表面能。这在微纳尺度或软材料接触中很重要。

粗糙表面接触（Greenwood-Williamson模型）

真实表面具有多尺度粗糙度，GW模型将接触视为多个微凸体的集合： $$F = \frac{4}{3} n E^* \int_d^{\infty} \sqrt{R_a} (z-d)^{3/2} \phi(z) dz$$ 其中 $n$ 是微凸体密度，$\phi(z)$ 是高度分布函数。

塑性接触

当接触应力超过屈服极限时，需要考虑塑性变形： $$F = \begin{cases} \frac{4}{3} E^* \sqrt{R} \delta^{3/2} & \delta < \delta_y \\ F_y + H \cdot A_{\text{plastic}} & \delta \geq \delta_y \end{cases}$$ 其中 $H$ 是材料硬度，$\delta_y$ 是屈服点压入深度。

在强化学习环境中，选择合适的接触模型需要考虑：

计算效率：赫兹模型计算快速，适合大规模并行仿真
物理真实性：复杂模型更准确但计算昂贵
任务相关性：抓取软物体需要考虑粘附，金属接触需要考虑塑性

5.3 柔性接触模型

5.3.1 罚函数方法

罚函数方法是物理仿真中最广泛使用的接触处理方法。其核心思想是允许小的穿透，并施加与穿透深度成比例的恢复力： $$\mathbf{f}_n = \begin{cases} k_n \delta \mathbf{n} + c_n \dot{\delta} \mathbf{n} & \delta > 0 \\ 0 & \delta \leq 0 \end{cases}$$ 其中：

$k_n$ 是法向刚度
$c_n$ 是法向阻尼
$\delta$ 是穿透深度
$\mathbf{n}$ 是接触法向

这种方法的优势在于：

实现简单：不需要求解约束方程
连续性好：力是位置和速度的连续函数
并行友好：每个接触可以独立计算

但也存在固有缺陷：

允许穿透违背了刚体假设
需要仔细调节参数以平衡穿透和稳定性
过大的刚度导致数值刚性问题

参数选择策略：

刚度选择基于时间步长和稳定性要求： $$k_n = \frac{\alpha m}{\Delta t^2}$$ 其中 $\alpha \in [10, 100]$ 是安全系数。阻尼系数通过临界阻尼比确定： $$c_n = 2\zeta\sqrt{k_n m}$$ 典型地，$\zeta \in [0.1, 0.5]$ 提供良好的能量耗散而不会过度阻尼。

5.3.2 Hunt-Crossley模型

Hunt-Crossley模型是对线性弹簧-阻尼器模型的非线性扩展，更准确地描述了材料的粘弹性行为： $$f_n = k_n \delta^n + \lambda \delta^n \dot{\delta}$$ 其中指数 $n$ 通常取1.5（对应赫兹接触）。这个模型的关键创新是阻尼力与穿透深度成比例，避免了分离时的粘附力问题。

模型的能量耗散特性由恢复系数 $e$ 决定： $$\lambda = \frac{3k_n(1-e)}{2\dot{\delta}_{\text{impact}}}$$ 其中 $\dot{\delta}_{\text{impact}}$ 是碰撞速度。这确保了正确的能量损失： $$E_{\text{loss}} = \frac{1}{2}m v^2(1-e^2)$$ 数值实现注意事项：

if (penetration > 0) {
    // 计算穿透速度
    double v_pen = dot(v_relative, normal);

    // Hunt-Crossley力
    double delta_n = pow(penetration, 1.5);
    double f_elastic = k_n * delta_n;
    double f_damping = lambda * delta_n * v_pen;

    // 确保阻尼不会反转速度
    f_damping = min(f_damping, f_elastic);

    force = (f_elastic + f_damping) * normal;
}

5.3.3 阻尼与能量耗散

接触中的能量耗散机制对于仿真的真实性至关重要。主要的耗散来源包括：

材料阻尼：内部摩擦导致的能量损失
结构阻尼：宏观变形引起的耗散
界面耗散：接触面微滑移产生的能量损失

恢复系数模型：

恢复系数 $e$ 与碰撞速度的关系可以用经验公式描述： $$e(v) = e_0 - (e_0 - e_{\infty})(1 - e^{-v/v_c})$$ 其中：

$e_0$ 是零速度恢复系数
$e_{\infty}$ 是高速恢复系数
$v_c$ 是特征速度

能量一致性检查：

在仿真中监控能量平衡是调试的重要手段： $$\Delta E_{\text{total}} = \Delta E_{\text{kinetic}} + \Delta E_{\text{potential}} + E_{\text{dissipated}}$$ 理想情况下，$\Delta E_{\text{total}} \approx 0$（考虑数值误差）。

5.3.4 接触参数的选择

选择合适的接触参数是成功仿真的关键。以下是系统化的参数选择流程：

1. 基于物理的初始估计：

对于已知材料，使用材料属性计算初始参数： $$k_n = \frac{4}{3}E^*\sqrt{R_{\text{eff}}}$$

$$c_n = 2\sqrt{mk_n}\cdot\frac{\ln(e)}{\sqrt{\pi^2 + \ln^2(e)}}$$

2. 稳定性约束：

确保数值稳定的最大刚度： $$k_{\max} = \frac{4m}{\Delta t^2}$$

3. 穿透深度约束：

限制最大穿透深度以保持视觉真实性： $$k_{\min} = \frac{F_{\text{typical}}}{\delta_{\max}}$$

4. 域随机化策略：

在强化学习训练中，参数随机化提高泛化能力： $$\begin{aligned} k_n &\sim \log\mathcal{N}(\mu_k, \sigma_k^2) \\ e &\sim \text{Beta}(\alpha_e, \beta_e) \\ \mu_f &\sim \text{TruncNormal}(\mu_{\mu}, \sigma_{\mu}^2, 0, 1) \end{aligned}$$

参数敏感性分析：

不同任务对参数的敏感性不同：

抓取任务：对摩擦系数最敏感
弹跳任务：对恢复系数最敏感
推动任务：对接触刚度较敏感

建议的参数范围（机器人仿真）：

接触刚度：$10^3 - 10^6$ N/m
阻尼比：0.1 - 0.5
恢复系数：0.0 - 0.9
摩擦系数：0.1 - 2.0