面面接触算法初探

接触是有限元非线性分析中非常重要、同时也相对复杂的一类问题。与普通的材料非线性或几何非线性不同，接触问题不仅需要求解结构本身的平衡方程，还需要在计算过程中不断判断两个物体是否发生接触、接触发生在什么位置、接触面的法向和切向如何定义，以及接触约束应该以什么方式加入整体有限元方程。因此，一个完整的接触算法实际上涉及几何搜索、接触运动学、约束条件、接触本构以及非线性求解等多个环节。

本文主要讨论一种基于 slave-side Gauss integration 和 pointwise projection 的 Surface-to-Surface 接触方法[^1]，目的是让读者建立一张完整的接触计算“地图”。其核心问题可以概括为：给定 slave 面上的一个积分点，如何在 master 面上找到与之对应的投影点，并根据两者之间的相对位置构造法向间隙（normal gap）和切向滑移（tangential slip）；随后，再根据接触状态建立相应的接触约束，并将产生的接触力以及一致切线刚度组装到整体有限元方程中。所以，本文将沿着接触几何 → 运动学变量 → 接触本构 → 接触虚功 → 有限元残差 → 接触刚度这一条主线展开，介绍接触算法。

下面先固定本文使用的符号约定。

• \(s\)：secondary，从面；
• \(m\)：master，主面；
• \(g\)：代码中的法向 gap；
• \(g\le 0\)：穿透/闭合；
• \(g>0\)：分离；
• \(\mathbf n\)：主面法向；
• \(\mathbf t_1,\mathbf t_2\)：主面切向基。

0. 点面接触 vs 面面接触

在有限元接触算法中，常见的接触离散方式包括点—面接触（Node-to-Surface, N2S ）和面—面接触（Surface-to-Surface, S2S）。

点—面接触以 slave 节点作为基本接触对象。算法将 slave 节点投影到 master 面上，根据节点与投影点之间的相对位置判断接触状态并计算接触力。其实现比较直接，但接触约束集中在离散节点上，因此通常对 slave/master 的选择以及两侧网格密度更加敏感。经典 Node-to-Segment 接触算法及其离散特性可参考 Zavarise 和 De Lorenzis [^3] 。

面—面接触则从两个表面之间的分布式接触作用出发。有限元离散后，通常在 slave 面上设置积分点，将 slave 高斯点投影到 master 面，在积分点处计算接触状态和接触力，再通过数值积分将这些贡献传递到两侧节点。

本文主要讨论第二种情况。后面的推导将从一个slave 高斯点如何投影到 master 面开始，逐步介绍接触运动学、接触本构、虚功、有限元残差以及接触刚度。

1. 接触运动学 ( Kinematics )

下面采用有限元接触力学中常见的 slave/master surface 参数化和投影描述 [^1, 6] 。

1.1 Slave 与 Master 表面的离散

设一个从面单元有 \(n_s\) 个节点，一个主面单元有 \(n_m\) 个节点。

从面上的点由形函数插值：

\[\mathbf x_s(\boldsymbol\xi_s) = \sum_{I=1}^{n_s}N_I^s(\boldsymbol\xi_s)\mathbf x_I^s \]

主面投影点为：

\[\mathbf x_m(\boldsymbol\xi_m) = \sum_{J=1}^{n_m}N_J^m(\boldsymbol\xi_m)\mathbf x_J^m \]

对于 S2S contact，我们通常从 slave 面上的一个积分点 \(\mathbf x_s\) 出发，在 master 面上寻找与之对应的投影点 \(\mathbf x_m\)。

因此，接触运动学首先要解决的问题是：

给定一个 slave Gauss point，它对应 master surface 上的哪个位置？

1.2 法向 gap

设 master 投影点处的法向为 \(\mathbf n\)，法向 gap 定义为：

\[g = \mathbf n\cdot(\mathbf x_s-\mathbf x_m) \]

因此：

• 两个表面分离时，\(g>0\)；
• 从面向主面靠近时，\(g\) 减小；
• 恰好接触时，\(g=0\)；
• 从面穿过主面时 ( penetration ) ，\(g<0\)。

1.3 Master surface 上的投影

如果采用一个给定的法向场 \(\mathbf n(\boldsymbol\xi)\)，slave 点和 master 投影点之间可以建立关系：

\[\boxed{ \mathbf x_s = \mathbf x_m(\boldsymbol\xi) + g\,\mathbf n(\boldsymbol\xi) } \]

其中

\[\mathbf x_m(\boldsymbol\xi) = \sum_A N_A(\boldsymbol\xi)\mathbf x_A\\ \]

如果法向由 master 节点法向插值得到，则：

\[\mathbf n(\boldsymbol\xi) = \operatorname{normalize} \left( \sum_A N_A(\boldsymbol\xi)\mathbf n_A \right) \]

这里未知量是：\(\boxed{\boldsymbol\xi,\quad g}\)

这个投影问题可以写成一个三维非线性方程：

\[\boxed{ \mathbf r(\xi^1,\xi^2,h) = \mathbf x_s - \mathbf x_m(\xi^1,\xi^2) - g\,\mathbf n(\xi^1,\xi^2) = \mathbf 0 } \]

对于一般曲面，\(\mathbf x_m\) 和 \(\mathbf n\) 都依赖于 \(\boldsymbol\xi\)，因此这是一个非线性方程组，可以采用 Newton–Raphson 方法求解。

令

\[\mathbf q = \begin{bmatrix} \xi^1\\ \xi^2\\ g \end{bmatrix} \]

则 Newton 迭代写成：

\[\mathbf J^{(k)} \Delta\mathbf q = -\mathbf r^{(k)} \]

以及

\[\mathbf q^{(k+1)} = \mathbf q^{(k)} + \Delta\mathbf q \]

直到投影残差满足给定容差。

1.4 Closest-point projection 与节点法向投影

需要注意，接触投影并不只有一种定义。

经典的 closest-point projection 要求 slave-master 连线垂直于 master surface 的局部切平面。

定义 master surface 的两个协变切向量：

\[\mathbf a_\alpha = \frac{\partial\mathbf x_m} {\partial\xi^\alpha}, \qquad \alpha=1,2 \]

其中：

\[\mathbf a_\alpha = \sum_A N_{A,\alpha}\mathbf x_A \]

closest-point projection 满足：

\[(\mathbf x_s-\mathbf x_m) \cdot \mathbf a_1 = 0 \]

以及：

\[(\mathbf x_s-\mathbf x_m) \cdot \mathbf a_2 = 0 \]

求得投影坐标 \(\bar{\boldsymbol\xi}\) 后，可以根据 master surface 的真实局部几何计算法向：

\[\boxed{ \mathbf n_{\mathrm{local}} = \frac{ \mathbf a_1\times\mathbf a_2 }{ \|\mathbf a_1\times\mathbf a_2\| } } \]

而另一种实现方式是首先在 master 节点处构造节点法向 \(\mathbf n_A\)，然后通过形函数插值得到：

\[\mathbf n_{\mathrm{interp}} = \operatorname{normalize} \left( \sum_A N_A\mathbf n_A \right) \]

再利用这个平滑法向场进行投影。

1.5 切向运动

在接触点处建立局部正交坐标系：

\[\{ \mathbf n, \mathbf t_1, \mathbf t_2 \} \]

其中 \(\mathbf t_1\) 和 \(\mathbf t_2\) 位于 master surface 的局部切平面。

为了描述摩擦，需要进一步定义 slave/master 之间的相对切向运动。需要注意的是，在有限滑移问题中，切向 slip 一般不能简单理解为当前 slave 点和 master 投影点连线的切向分量。特别是在 closest-point projection 下，\(\mathbf x_s-\mathbf x_m\) 本身沿法向，因此其切向投影理论上为零。

摩擦真正关心的是 slave 接触点相对于 master surface 的切向运动历史。因此本文把切向滑移作为历史变量，使用：

\[s_\alpha \]

表示局部切向滑移变量，而其增量记为：

\[\Delta s_\alpha \]

实际有限滑移算法中，\(\Delta s_\alpha\) 需要根据前后两个构型中接触投影点的变化以及局部切向基进行更新。

到这里，接触几何已经被转化为一组局部运动学变量：

\[\boxed{ g,\quad s_1,\quad s_2 } \]

接下来需要研究：

当 slave/master 节点发生微小位移时，这些运动学变量如何变化？

2. 变分：从几何量得到 \(B\) 矩阵

首先考虑最简单的情况：暂时冻结当前接触投影关系和局部坐标系，即忽略 \(\mathbf n\)、\(\mathbf t_\alpha\) 和投影坐标 \(\bar{\boldsymbol\xi}\) 本身的变化。

此时：

\[\delta g = \mathbf n\cdot (\delta\mathbf x_s-\delta\mathbf x_m) \\ \delta s_\alpha = \mathbf t_\alpha\cdot (\delta\mathbf x_s-\delta\mathbf x_m) \]

对 slave 节点 \(I\)：

\[\delta\mathbf x_s = N_I^s\delta\mathbf u_I^s \]

因此：

\[\delta g = +N_I^s\mathbf n^T\delta\mathbf u_I^s \\ \delta s_\alpha = +N_I^s\mathbf t_\alpha^T\delta\mathbf u_I^s \]

对 master 节点 \(J\)：

\[\delta g = -N_J^m\mathbf n^T\delta\mathbf u_J^m \\ \delta s_\alpha = -N_J^m\mathbf t_\alpha^T\delta\mathbf u_J^m \]

所以可定义：

\[\mathbf e= \begin{bmatrix} h\\s_1\\s_2 \end{bmatrix}, \qquad \delta\mathbf e = B_c\,\delta\mathbf d \]

其中 \(\mathbf d\) 包含 slave/master 位移自由度，\(B_c\) 可以理解为 contact element 的运动学矩阵。

这里需要强调，目前得到的是一种 frozen-frame linearization。严格的有限滑移问题中：

\[\bar{\boldsymbol\xi} = \bar{\boldsymbol\xi}(\mathbf d) \]

同时：

\[\mathbf n = \mathbf n(\mathbf d) \]

以及：

\[\mathbf t_\alpha = \mathbf t_\alpha(\mathbf d) \]

因此完整线性化还会出现 \(\delta\bar{\boldsymbol\xi}, \delta\mathbf n, \delta\mathbf t_\alpha\) 等项。

3. 接触本构：由 \(g,s_\alpha\) 计算接触力

接触运动学告诉我们两个表面之间“发生了什么”，接触本构则回答：

给定当前的 gap 和 slip，应该产生多大的接触力？

3.1 法向接触

理想的无粘结单边接触满足 Signorini 条件：

\[g\ge0,\qquad p\ge0,\qquad pg=0 \]

严格满足该约束通常需要 Lagrange multiplier、augmented Lagrangian 或其他约束方法。

一种更简单的处理方式是 penalty method。Penalty 方法允许一个很小的数值 penetration：

\[p= -k_n g \qquad (g<0) \\ p=0 \qquad (g>0) \]

其中 \(p\) 是接触压力，\(k_n\) 是法向 penalty 刚度。

在 active contact 状态下：

\[\frac{\partial p}{\partial g}=-k_n \]

因此 penalty stiffness 越大，允许的 penetration 越小，但系统条件数通常也会变差。

3.2 切向摩擦

摩擦具有明显的路径依赖性，因此需要保存切向历史变量。

设当前增量得到的 trial tangential traction 为：

\[ \boldsymbol\tau^{\text{trial}} = k_t \left( \mathbf s_{\text{old}}+\Delta\mathbf s \right) \]

Coulomb 摩擦条件为：

\[ \|\boldsymbol\tau\|\le\mu p \]

若\(\|\boldsymbol\tau^{\text{trial}}\| \le\mu p\)，则为 sticking \(\boldsymbol\tau = \boldsymbol\tau^{\text{trial}}\)。

若超过摩擦圆，则为 slipping，执行 radial return：\(\boldsymbol\tau = \mu p\, \frac{\boldsymbol\tau^{\text{trial}}} {\|\boldsymbol\tau^{\text{trial}}\|}\)

4. 接触虚功

slave 面上的接触牵引可以写成：

\[\mathbf t_c = p\mathbf n + \tau_1\mathbf t_1 + \tau_2\mathbf t_2 \]

考虑 slave/master 上的作用力与反作用力，接触虚功可合并为：

\[\delta W_c = \int_{\Gamma_s} \left( p\,\delta g+ \tau_1\,\delta s_1 + \tau_2\,\delta s_2 \right)dA \]

定义接触应力向量：

\[\boldsymbol\sigma_c = \begin{bmatrix} p\\ \tau_1\\ \tau_2 \end{bmatrix} \]

于是：

\[\delta W_c = \int_{\Gamma_s} \delta\mathbf d^T B^T\boldsymbol\sigma_c \,dA \]

到这里，可以发现 contact element 与普通 continuum element 的数学结构非常相似。

对于普通实体单元：

\[\delta\boldsymbol\varepsilon = B\,\delta\mathbf u \]

以及：

\[\delta W_{\mathrm{int}} = \int_\Omega \delta\boldsymbol\varepsilon^T \boldsymbol\sigma \,dV \]

对于接触单元：

\[\delta\mathbf e = B_c\delta\mathbf d \]

以及：

\[\delta W_c = \int_{\Gamma_c} \delta\mathbf e^T \boldsymbol\sigma_c \,dA \]

因此可以建立下面的对应关系：

\[\boxed{ \boldsymbol\varepsilon \longleftrightarrow (g,s_1,s_2) } \]

\[\boxed{ \boldsymbol\sigma \longleftrightarrow (p,\tau_1,\tau_2) } \]

\[\boxed{ B \longleftrightarrow B_c } \]

区别在于：普通实体单元在体积域上积分，而接触单元在接触界面上积分。

5. 从虚功得到有限元残差

对于 slave surface 参数坐标：

\[dA=J_s\,d\xi\,d\eta \]

采用 Gauss quadrature：

\[\int_{\Gamma_s} (\cdot)\,dA \approx \sum_{q=1}^{n_q} w_qJ_{s,q}(\cdot)_q \]

因此

\[\delta W_c \approx \sum_q \delta\mathbf d_q^T \left( J_{s,q}w_q B_{c,q}^T\boldsymbol\sigma_{c,q} \right) \]

所以单个高斯点的接触残差为：

\[\boxed{ \mathbf r_q^c = J_{s,q}w_q B_{c,q}^T\boldsymbol\sigma_{c,q} } \]

从程序实现的角度，可以将这一过程理解为：

1. 遍历 slave facets；
2. 在 slave facet 上布置 Gauss points；
3. 根据 slave shape functions 计算 \(\mathbf x_s\)；
4. 搜索候选 master facets；
5. 将 \(\mathbf x_s\) 投影到 master surface；
6. 得到 \(\bar{\boldsymbol\xi}\)、\(\mathbf x_m\) 和局部接触坐标系；
7. 计算 normal gap 和 tangential slip increment；
8. 根据接触本构计算 \(p\) 和 \(\boldsymbol\tau\)；
9. 构造 \(B_c\)；
10. 计算 \(w_qJ_{s,q}B_{c,q}^T\boldsymbol\sigma_{c,q}\)；
11. 将 slave/master 两侧贡献组装到整体残差。

这就是从几何接触关系到有限元平衡方程的完整传递过程。

6. 接触刚度

为了使用 Newton–Raphson 方法求解，需要对接触残差进行线性化。

首先考虑接触本构本身产生的切线，设 \(\delta\boldsymbol\sigma_c = D_c\,\delta\mathbf e\)，又因为 \(\delta\mathbf e=B\,\delta\mathbf d\)，所以 \(\delta\boldsymbol\sigma_c = D_cB\,\delta\mathbf d\)。

如果暂时冻结 \(B_c\)、积分 Jacobian 以及接触局部坐标系，则得到接触刚度的本构部分：

\[\boxed{ \mathbf K_{mat}^c = \sum_q J_{s,q}w_q B_{c,q}^TD_{c,q}B_q } \]

典型情况下：

• sticking：\(D_t=k_t I_2\)

• slipping：若\(\boldsymbol\tau=\tau_c\mathbf m, \mathbf m= \frac{\boldsymbol s_{\text{trial}}} {\|\boldsymbol s_{\text{trial}}\|}\)，则径向返回切线近似为\(D_t = \frac{\tau_c} {\|\boldsymbol s_{\text{trial}}\|} \left(I_2-\mathbf m\mathbf m^T\right)\)

如果 \(\tau_c=\mu p\)，还会出现切向力对法向压力的耦合项， \(\frac{\partial\boldsymbol\tau} {\partial g} \neq0\)，法向和切向接触关系发生耦合，因此最终的 contact tangent 不一定是对称矩阵。

7. 几何非线性与一致线性化

上一节得到的：

\[B_c^TD_cB_c \]

只是接触刚度中由接触本构产生的部分。

对于有限滑移 contact，接触几何本身也随节点位移变化。

例如：

\[g = \mathbf n\cdot (\mathbf x_s-\mathbf x_m) \]

其完整变分包含：

\[\delta g = \mathbf n\cdot ( \delta\mathbf x_s - \delta\mathbf x_m ) + \delta\mathbf n\cdot ( \mathbf x_s-\mathbf x_m ) \]

但实际上问题还要更加复杂，因为：

\[\mathbf x_m = \mathbf x_m(\bar{\boldsymbol\xi}) \]

而投影坐标本身也是节点自由度的函数：

\[\bar{\boldsymbol\xi} = \bar{\boldsymbol\xi}(\mathbf d) \]

因此：

\[\delta\mathbf x_m = \sum_A N_A \delta\mathbf x_A^m + \mathbf a_\alpha \delta\bar\xi^\alpha \]

类似地：

\[\delta\mathbf n \neq0 \]

\[\delta\mathbf t_\alpha \neq0 \]

同时 slave surface 的面积 Jacobian：

\[J_s \]

也可能随当前构型变化。

因此，完整的 contact residual：

\[\mathbf R_c = \int_{\Gamma_c} B_c^T \boldsymbol\sigma_c \,dA \]

线性化之后不仅产生：

\[B_c^TD_cB_c \]

还会产生由以下因素引起的几何项：

• 投影坐标 \(\bar{\boldsymbol\xi}\) 的变化；
• master normal \(\mathbf n\) 的变化；
• tangential basis \(\mathbf t_\alpha\) 的变化；
• \(B_c\) 本身的变化；
• slave/master surface geometry 的变化；
• surface integration measure 的变化。

因此可以概念性地写成：

\[\boxed{ K_c = K_{\mathrm{mat}}^c + K_{\mathrm{geo}}^c } \]

\(K_{\mathrm{geo}}^c\) 收集所有由接触几何变化产生的一致线性化项。对于大滑移、曲面接触以及高阶单元，这部分对于 Newton 收敛性尤其重要。

这部分完整推导会明显比前面复杂，本文暂时按下不表——再往下推就要加钱了。

8. 总结

最后总结一下，接触算法的本质可以概括为：

\[\boxed{ \begin{aligned} &g,s_\alpha \quad\text{定义几何约束}\\ &p,\tau_\alpha \quad\text{由接触本构计算}\\ &p\,\delta g+\tau_\alpha\delta s_\alpha \quad\text{形成接触虚功}\\ &B^T\sigma,\ B^TDB \quad\text{形成残差和切线}\\ &K\Delta u=-R \quad\text{通过 Newton 施加约束} \end{aligned} } \]

参考文献

1. G. Zavarise, P. Wriggers, A segment-to-segment contact strategy, Mathematical and Computer Modelling, 28(4–8), 497–515, 1998.
2. J. C. Simo, T. A. Laursen, An augmented Lagrangian treatment of contact problems involving friction, Computers & Structures, 42(1), 97–116, 1992.
3. G. Zavarise, L. De Lorenzis, The node-to-segment algorithm for 2D frictionless contact: Classical formulation and special cases, Computer Methods in Applied Mechanics and Engineering, 2009.
4. M. A. Puso, T. A. Laursen, A mortar segment-to-segment contact method for large deformation solid mechanics, Computer Methods in Applied Mechanics and Engineering, 193, 601–629, 2004.
5. M. A. Puso, T. A. Laursen, A mortar segment-to-segment frictional contact method for large deformations, Computer Methods in Applied Mechanics and Engineering, 193, 4891–4913, 2004.
6. P. Wriggers, Computational Contact Mechanics, Springer.