有限元接触的数学本质：约束、KKT 与 Active-set

接触问题简介

在有限元分析中，两个物体可以相互接触，也可以彼此分离，但不能相互穿透。这看似是一个简单的几何条件，实际求解起来却并不容易。因为在计算开始之前，我们并不知道哪些位置会发生接触。 接触区域需要随着物体的变形不断调整，而接触力又反过来影响物体的变形。那么，有限元求解器究竟是如何处理这一问题的呢？

我们先从没有接触约束的情况出发。

设两个物体占据区域 \(\Omega_1,\Omega_2\)，位移分别为 \(\mathbf u_1,\mathbf u_2\)，统一记为

\[\mathbf u=(\mathbf u_1,\mathbf u_2) \]

在不考虑接触的情况下，系统的总势能可以写成

\[\boxed{ \Pi_0[\mathbf u] = \sum_{i=1}^{2} \left[ \int_{\Omega_i}W_i(\boldsymbol\varepsilon(\mathbf u_i))\,dV - \int_{\Omega_i}\mathbf b_i\cdot\mathbf u_i\,dV - \int_{\Gamma_{t,i}}\bar{\mathbf t}_i\cdot\mathbf u_i\,dA \right] } \]

其中，\(W_i\) 是单位体积的应变能密度，\(\mathbf b_i\) 是体力，\(\bar{\mathbf t}_i\) 是给定的边界面力。

如果没有接触约束，系统的平衡状态可以通过总势能的驻值条件来确定：

\[D\Pi_0[\mathbf u](\mathbf v)=0, \qquad \forall,\mathbf v\in V_0. \]

这里，\(\mathbf v\) 表示满足齐次位移边界条件的虚位移。也就是说，在所有满足位移边界条件的位移场中，寻找使总势能取驻值的位移场。

但是，当两个物体可能发生接触时，还必须考虑一个额外的几何限制：

两个物体可以彼此分离，也可以发生接触，但不能相互穿透。

如何用数学语言描述这一限制呢？

我们定义 \(h\) 为带符号的穿透量，并约定

\[h>0 \quad\Longleftrightarrow\quad \text{发生穿透}. \]

于是不可穿透条件写成

\[\boxed{ h[\mathbf u]\leq 0 \qquad \text{on }\Gamma_c } \]

其中 \(\Gamma_c\) 表示候选接触界面。

于是，接触问题就可以表述为一个带不等式约束的最优化问题：

\[\boxed{ \begin{aligned} \min_{\mathbf u}\quad&\Pi_0[\mathbf u]\\ \mathrm{s.t.}\quad& \mathbf u_i=\bar{\mathbf u}_i &&\text{on }\Gamma_{u,i},\\ &h[\mathbf u]\leq0 &&\text{on }\Gamma_c. \end{aligned} } \]

这里考虑的是保守载荷下的准静态、无摩擦、无黏着接触问题。对于稳定平衡状态，可以从势能最小化的角度理解。与普通有限元问题相比，接触问题多了一个关键要求：不仅要满足力学平衡，还必须满足不可穿透约束。

引入接触压力

如何求解这个带不等式约束的问题？

一种自然的思路是引入拉格朗日乘子，将接触约束纳入原来的平衡方程。

将接触压力 \(p\) 作为不可穿透约束对应的乘子场，构造拉格朗日泛函：

\[\mathcal L[\mathbf u,p] = \Pi_0[\mathbf u] + \int_{\Gamma_c}p\,h[\mathbf u]\,dA \]

这里的 \(p\) 并不是预先给定的外力，而是为了满足不可穿透约束而产生的约束反力。

对应的连续 KKT（Karush–Kuhn–Tucker）条件为

\[\boxed{ \begin{aligned} D\Pi_0[\mathbf u](\mathbf v) +\int_{\Gamma_c} p\,Dh[\mathbf u](\mathbf v)\,dA&=0, &&\forall\,\mathbf v\in V_0, \\ h[\mathbf u]&\leq0,&&\text{on }\Gamma_c, \\ p&\geq0,&&\text{on }\Gamma_c, \\ p\,h[\mathbf u]&=0&&\text{on }\Gamma_c. \end{aligned} } \]

这四个条件分别表示：

• 驻点条件：内力、外力以及接触力必须满足虚功平衡；
• 原始可行性：\(h\leq0\)，两个物体不能相互穿透；
• 对偶可行性：\(p\geq0\)，无黏着接触只能传递压力；
• 互补条件：\(ph=0\)，接触压力与间隙不能同时非零。

这些条件在适当的约束资格条件下是局部最优解的必要条件；对于凸优化问题，在相应条件成立时也可以成为充分条件。

其中，最值得关注的是最后一个互补条件。

当两个物体分离时，

\[h<0\quad\Longrightarrow\quad p=0. \]

此时没有接触，自然也不应该存在接触压力。

当存在正的接触压力时，

\[p>0\quad\Longrightarrow\quad h=0. \]

此时两个物体恰好接触，不能发生穿透。

需要注意，\(h=0\) 且 \(p=0\) 也是允许的，它表示两个物体恰好接触，但没有产生接触压力。

至此，接触问题的数学结构已经清楚了。

我们不仅需要求解位移场 \(\mathbf u\) 和接触压力场 \(p\)，还需要判断哪些区域发生接触，哪些区域处于分离状态。困难在于，接触区域并不是预先已知的，而是随着物体的变形不断变化。因此，我们不能简单地在整个候选接触界面上施加 \(h=0\)，而必须在求解力学平衡的同时，确定哪些位置需要激活接触约束。

这正是接触问题区别于普通非线性有限元问题的关键：不仅要求解平衡方程，还需要确定哪些约束处于激活状态。

如果接触区域已经知道

为了理解求解过程，不妨先考虑一个简单的情况：假设我们已经知道哪些区域发生了接触。将候选接触界面划分为

\[\Gamma_c=\Gamma_A\cup\Gamma_I, \]

其中：

• \(\Gamma_A\)：当前认为处于接触状态的 active set；
• \(\Gamma_I\)：当前认为处于分离状态的 inactive set。

如果这个划分是正确的，那么 KKT 条件就可以分别处理：

\[\begin{array}{ll} \Gamma_A\text{ 上：}&h[\mathbf u]=0,\quad p\text{ 待求},\\[1mm] \Gamma_I\text{ 上：}&p=0,\quad h[\mathbf u]\text{ 由位移决定}. \end{array} \]

也就是说，在 active set 上强制满足接触约束，而在 inactive set 上不施加接触压力。这样一来，原来的不等式约束问题就暂时转化成了一个等式约束问题：

\[\boxed{ \begin{aligned} D\Pi_0[\mathbf u](\mathbf v) + \int_{\Gamma_A} p,Dh[\mathbf u](\mathbf v),dA &=0, &&\forall,\mathbf v\in V_0, \\ h[\mathbf u] &=0, &&\text{on }\Gamma_A. \end{aligned} } \]

第一行要求整个系统满足力学平衡，第二行要求当前接触区域既不分离，也不发生穿透。这里，位移 \(\mathbf u\) 和接触压力 \(p\)​ 都是未知量，需要同时求解。从这个角度看，接触压力的物理意义也就清楚了：接触压力是为了维持 \(h=0\) 这一几何约束，同时满足整体力学平衡而产生的约束反力。因此，只要接触区域已经确定，我们就可以按照等式约束问题求解位移和接触压力。但真正的问题是：我们事先并不知道这个划分是否正确。

Active-set：猜测、求解、检查、修正

既然无法提前确定接触区域，一个直接的办法就是：先猜测哪些区域发生了接触，求解以后再检查这个猜测是否正确。

这就是 Active-set 方法的基本思想：

\[\boxed{ \text{猜测接触区域} \rightarrow \text{求解位移和压力} \rightarrow \text{检查 KKT 条件} \rightarrow \text{修正接触区域} } \]

首先猜测一个 active set \(\Gamma_A\)，并按照上面的等式约束问题求解 \(\mathbf u\) 和 \(p\)。

求解完成以后，再检查这次猜测是否满足原来的 KKT 条件。这里，主要有两种情况需要修正：

• 如果某处被认为分离，但求解后：\(h[\mathbf u]>0\)。说明发生了穿透。那里不能继续保持分离，应加入接触区域。

• 如果某处被认为接触，但求解后：\(p<0\)。说明保持 \(h=0\) 需要一个拉力。而我们考虑的接触只能受压，因此那里应释放为分离状态。

由此，active-set 的基本更新规则可以概括为

\[\boxed{ \begin{aligned} \text{分离区发生穿透}&\ \Longrightarrow\ \text{加入接触区},\\ \text{接触区出现拉力}&\ \Longrightarrow\ \text{移出接触区}. \end{aligned} } \]

更新接触区域以后，再重新求解位移和接触压力，并继续检查 KKT 条件。如此反复，直到接触区域不再变化，整体平衡得到满足，并且所有位置都满足

\[h\leq0,\qquad p\geq0,\qquad ph=0. \]

至此，我们就找到了满足接触约束的平衡解。

回过头来看，Active-set 的核心并不复杂：它将接触区域未知的问题，转化为一系列接触区域暂时已知的等式约束问题。

每一次迭代都在回答两个问题：当前认为分离的区域是否发生了穿透？当前认为接触的区域是否产生了不合理的拉力？

通过不断修正接触区域，最终使整个系统同时满足力学平衡与 KKT 互补条件。

因此，KKT 条件告诉我们什么样的接触状态是正确的，而 Active-set 方法提供了一条寻找这个状态的数值求解路径。