有限元分析中的接触和摩擦模拟（二）

接触问题求解的效率和稳定性，不仅依赖于对接触界面的合理离散化，更依赖于算法的质量。目前已经发展了一系列的接触算法，分别具有不同的适用范围，对不同的问题各有优缺点。其中以罚函数法、拉氏乘子法和直接约束法最为通用，已经实现于多种商用有限元分析软件中。

拉氏乘子技术是将数学约束条件引入一个系统的最完美的数学描述。使用拉氏乘子引入施加接触体必须满足的非穿透约束条件，用以描述约束极值问题，可以使界面约束条件得到充分的满足，确保法向不可贯入性。其不足之处是增加了方程的自由度数，并使系统矩阵主对角线元素为零。这就需要在数值方案中处理非正定系统，数学上将发生困难，需实施额外的操作才能保证计算精度，从而使计算费用增加。另外，由于拉格朗日乘子与质量无关，使得包含惯性项的接触问题的有限元方程和显式时间积分的求解格式不协调，导致这种由拉格朗日乘子描述的接触算法难以直接用于显式动力撞击问题分析。

拉格朗日乘子技术还经常用于采用特殊的界面单元（GAP单元）描述接触的问题分析。该方法需要预先知道接触发生的确切部位，以便施加界面单元。这样的额外要求对于撞击、压力加工等通常事先并不知道准确接触区域所在的问题是难于满足的。

罚函数法是另一种施加接触约束的数值方法。其原理是一旦接触区域发生穿透，罚函数便夸大这种误差的影响，从而使系统的求解(满足力的平衡和位移的协调)无法正常实现。换言之，只有在约束条件满足之后，才能求解出有实际物理意义的结果。

用罚函数法施加接触约束的方法可以类比成在物体之间施加非线性弹簧所起的作用。该方法不增加未知量数从而增加系统矩阵带宽，数值上实施比较容易，在显示动力分析中被广泛应用。罚函数的缺点是，约束条件只能被近似的满足，界面上有少量的相互贯入产生，由于接触压力是假定与贯入量成正比的，因此算得的接触压力分布通常是振荡的。理论上讲，增大罚参数的数值可令计算精度提高，但在实际分析中罚函数的大小受到严格的限制，过大的罚参数将对系统的数值求解造成不良影响。一是罚参数愈大，显式解法时间步长的临界值降低得愈多，二是过大的罚参数可能使接触体的相对运动发生虚假的反向，从而使求解过程不稳定。

罚函数法的数值实现相对简洁，因此得到了最广泛的应用。如果罚参数的取值合理，多数情况下罚函数法能够提供健壮稳定的求解。但是对于由变形驱动的接触过程，拉氏乘子法通常能给出更精确更稳定的数值结果。

直接约束法也是常用的一种接触算法。用直接约束法处理接触问题是追踪物体的运动轨迹，一旦探测出发生接触，便将接触所需的运动约束（即法向无相对运动，切线可滑动）和节点力（法向压力和切向摩擦力）作为边界条件直接施加在产生接触的节点上，两接触体的运动约束转化成了节点自由度的约束和节点力的约束。这种方法对接触的描述精度高，具有普遍适应性。不需要增加特殊的界面单元，也不涉及复杂的接触条件变化。该方法不增加系统自由度数（由于接触关系的变化会增加系统矩阵带宽）。但是，界面约束的有效引入强烈依赖于界面上的离散情况。对于变形体与刚性体接触的情况，直接约束法效果明显，但是对于两个变形体接触的情况，如果两者在界面上的网格不匹配，该方法无法直接使用，需要进行特殊处理。

摄动拉格朗日法，是罚函数法和拉氏乘子法的混合模式。该方法使用的压力分布插值函数通常比位移分布插值函数低一阶。该方法能够克服纯粹的拉氏乘子法所导致的一系列数值困难。

增广拉格朗日法可视为罚函数法与拉氏乘子法之间的折中方案，该方法综合了两者的优点。该方法可用于分析无摩擦接触问题，也可用于大变形摩擦接触问题。增广拉格朗日技术可以与Uszawa算法结合在一起使用，在计算流程中使用嵌套的双重循环，内循环用以处理接触约束条件，外循环用以更新拉格朗日乘子。嵌套迭代方式增加了总的迭代次数，但是使算法的数值实现变得简洁。

因为接触条件是不等式约束，我们可以导出接触问题的等效变分不等式形式，位移场的解u必须满足该不等式。有限变形接触问题的弱形式如下：

式中，集合K定义为

如果仅考虑超弹性或者线弹性材料本构关系，对于无摩擦接触问题，上面的不等式与下面的最优化问题等价

为得到变分方程，需要用到法向间隙g_N和切向滑动量g_T的变分，我们将物体A表面采用参数坐标(ξ₁，ξ₂)描述，可以得到三维情况下法向间隙g_N的变分

计算δg_N时，不仅需要考虑坐标x^A和x^B的变分，尚需考虑x^A点处外法线方向ӣ^A对于曲面坐标(ξ1，ξ2)的变分，如下式：

因为ӣ^A是法向单位矢量，所以必然有

因此可以得到以下结果

如果物体B^A是刚性的，法向间隙的变分可以简化为

切向滑动g_T的变分也可按类似的方式导出

曲面坐标ξ_α的变分可由下式给出

其中

如果接触界面已知，即在增量形式的迭代求解中已经探测出接触约束的有效集，(10)式可以写为等式形式

式中C_c代表了接触约束活性集的贡献。

对于超弹性或者线弹性材料，上面的方程等价于对两个接触体总势能的最小化，

式中，附加的泛函∏_c包含了接触约束的贡献。

在已经判断出具体接触区域的前提下，有数种技术可用于导出C_c和∏_c的格式，包括拉氏乘子法、罚函数法、直接约束法、摄动拉格朗日法以及增广拉格朗日法等。

拉氏乘子法

应用拉格朗日乘子是将数学约束引入弱形式的经典方案。使用拉氏乘子技术，可以得到∏_c的列式如下：

式中，λ_N和λ_T分别为法向和切向拉氏乘子。g_N和g_T分别为法向间隙和切向的滑动量。

对∏_c变分，可以得到C_c的列式如下

对于两接触表面发生切向滑动的情况，切向面力矢量t_T可以根据摩擦滑动本构关系确定，因此应当改写λ_Tžδg_T→t_Tžδg_T，从而得下式

罚函数法

在罚函数法中，通过罚函数项将接触约束引入能量泛函∏

此处ε_N和ε_T分别代表法向和切向罚参数。

对上式取变分，可得

类似于拉氏乘子法，需要区分接触界面的纯粘接和切向滑动两种状态。对于滑动状态，有

直接约束法

该方法将接触区域的约束作为位移边界条件直接引入，因此不会导致未知数数量的增加。在判断出具体接触区域后，对于进入接触结点，不等式约束转换为等式约束

根据上式可得

将上式作为位移约束条件直接引入系统方程组，即可对问题进行求解。

摄动拉格朗日法

该方案是将拉氏乘子法和罚函数法相结合的一种混合格式。摄动拉格朗日法将接触的贡献∏^PL表达为

积分号中的第二项和第四项是对拉氏乘子的规则化，可以看作由拉氏乘子引起的补充能量。

对上式取变分可得

右边的第一个和第三个积分项与拉氏乘子格式相关。第二项和第四项则给出接触界面的法向和切向的本构关系：

如果将上式代入式（31）积分号内的第一和第三项，即得到标准的罚函数格式；如果令ε_N和ε_T趋向于无穷大，则可得标准的拉氏乘子格式。

增广拉格朗日法

增广拉格朗日格式提供了另外一种将不可微的接触和摩擦条件规则化的方案。其主要思路是将罚函数法或界面的本构定律与C₁可微的鞍点泛函相结合。在增广拉格朗日法中，用以下泛函引入法向接触的贡献：

其中，ƛ_N=λ_N+ε_Ng_N。

上述泛函的结构形式令其不仅适用于ƛ_N≤0的情况，也适用于ƛ_N>0的情况，后者意味着两接触体间的间隙是张开的。对上式取变分，得：

用同样的方法，对于遵循经典库仑摩擦定律的切向条件，可以导出类似的增广拉格朗日格式。采用切向相对滑移增量Δg_T和切向拉氏乘子ƛ_T=λ_T+ε_TΔg_T，构造如下的泛函引入切向约束

其中μ为摩擦系数，ṗ_N(p顶部黑点代表箭头)为增广的法向接触压力。对于未接触状态(ƛ_N>0)，可采用如下的切向泛函：

在增广拉格朗日法中，拉氏乘子λ是未知的，因此需要在算法流程中采用迭代循环以更新λ的值。

接触搜索是计算接触力学中最关键的问题之一。一般说来，接触搜索可分为两个阶段：首先是空间搜索，探测出可能进入接触的物体或有限单元；然后是接触探测，确定相互接触的单元对或者结点对。下面我们对空间搜索和接触探测的算法作一简单综述。

空间搜索实际上是采用一些排序算法对固体表面的单元进行排序操作。计算接触力学中发展了以下几种方案用以执行该操作。

网格元法

该方法将一包含全部接触体的空间划分为多个形状规则的网格元。当单元均匀的分布于各网格元时，网格元法效率很高。但如果单元主要聚集于少数网格元中，该算法不具有优势。因此出现了一种自适应的网格元法，可根据有限单元空间分布的非均匀程度自动调整网格元的划分。然而，生成自适应网格将导致额外的计算花费。

二叉树法

二叉树法同样是基于一个包含多个规则网格元的空间格栅。但是二叉树结构中仅保留含有有限单元的网格元。搜索消耗的时间取决于格栅的初始构造。利用一些特殊技术，例如平衡树的分支以及最小化树的深度，可以使搜索时间进一步减少。生成二叉树的时间消耗消耗为O(NlogN)，搜索的时间消耗也为这个量级。

空间堆排序算法

在该算法中，按照沿空间轴的坐标递增的顺序对一个单元链表排序。该算法通常与基于接触体的网格元相结合使用，其搜索时间为O(NlogN)量级。空间堆排序算法的优点在于：不需要特殊的数据结构，对于对象的空间尺度不敏感，而且仅需要一个O(N)量级的数组来存储必要的数据。

在构建接触探测算法流程时，需要区分几种不同的情况：刚性体与刚性体接触、变形体与刚性体接触、变形体与变形体接触。

1. 变形体与刚性体接触

在变形体与刚性体接触的情况下，刚性体表面可用一个隐式函数表示

刚性体表面也可用于定义接触法方向，因此可将其作为主接触体或参考体。由利用隐函数f给出刚性体表面的解析描述后，判断一个从结点x_s=(x_s,y_s,z_s)是在刚性体的内部还是外部仅需要考虑f(x_s,y_s,z_s)的数值。如果f<0，x_s在刚性体内部，如果f>0则在外部。刚性体表面法方向也很容易得到：

2. 变形体与变形体接触

探测变形体与变形体的接触，远较探测变形体与刚性体的接触困难。目前已有数种有效技术用于变形体之间的接触探测。

Benson和Hallquis发展了一种可用于二维和三维问题的探测方法。该方法把一指定从结点x_s的局部接触探测分为三步：判断离x_s距离最近的主结点x_k^A、确定包含x_s在主接触面的投影ẍ^A(两点代表横杠，下同)的单元面和计算ẍ^A点的空间流动坐标。

Kane等则提出了另外一种方法，对物体进行三角剖分，可对有限元离散模型进行贯入探测。其思路是基于以下事实：当有相互贯入发生时必然出现边界段的交叉。在二维情况下，可通过计算边界段的面积来确定是否发生交叉。在三维情况下，必须考虑边界面的交叉，因此需要进行体积检查。

来源：模态空间

王朋波，清华大学力学博士，汽车结构CAE分析专家。重庆市科协成员、《计算机辅助工程》期刊审稿人、交通运输部项目评审专家。专业领域为整车疲劳耐久/NVH/碰撞安全性能开发与仿真计算，车体结构优化与轻量化，CAE分析流程自动化等。王朋波私人微信：poplewang；加微请注明：单位+姓名。