1 研究背景

接触分析中的一个关键问题是如何对接触界面进行离散，常见的方法有点-点接触（Node to Node，NTN），点-面接触（Nodes to Surface，NTS）和面-面接触（Surface to Surface，STS）。其中点-点接触要求接触的节点之间必须重合，且在变形过程中不发生显著的相对位移，无法处理大滑移问题；而点-面接触虽然不要求网格重合，也可以处理大滑移问题，但其不满足接触分片试验，计算精度较低。上述两类方案均无法很好的处理三维大滑移接触问题。

面-面接触算法的早期形式是高斯点-面（Gauss Points to Surface，GPTS），即在从面面片的高斯积分点上施加接触约束条件，该方法在二维问题中有一定优势，但推广到三维问题时精度不佳，且无法达到最佳收敛效率。Mortar方法是一种用于连接不同网格密度求解域的数值方法，最早广泛应用于区域分解方法中，可以保证离散的最佳收敛效率。上世纪初，Mortar方法开始应用于接触问题，经过二十余年发展，从二维到三维，小滑移到大滑移，无摩擦到有摩擦，目前已经成为面-面接触算法中最成熟的界面离散方案，广泛应用于各类商用有限元分析软件中。

2 Mortar接触界面离散

Mortar方法中，接触界面的接触力虚功表达式为：

其中是从节点A处的正压力，采用罚函数法时，定义为：

是从节点A处的平均外法线。定义Mortar积分：

代入虚功表达式得到：

分别定义从节点和主节点上的MORTAR接触内力：

与传统的面-面接触离散方法不同，Mortar接触离散的核心是如何精确计算Mortar积分和。由于常规有限元离散中单元面片之间只满足连续性，因此无法直接在主面或者从面上进行数值积分，需要采用投影-裁切（Projection-Clipping Algorithm）构造积分域。

图1 投影裁切算法

如图所示，构造Mortar积分的计算区域方法如下：

（a） 计算从面片的平均外法线，构造一个平均面

（b） 将主面片投影到从面片的平均面上

（c） 调用切分算法，生成主从面相交的凸多边形

（d） 对凸多边形进行三角剖分，生成Mortar积分三角片。

由此可见，Mortar界面离散方法采用了一定的平均化和光滑化手段，避免了离散下的积分区域导致的性能退化问题。

3 Mortar接触算法的实现

Mortar接触算法的实现类似于经典的NTS算法，其接触变量（正压力，摩擦力，滑移量等）均存储在从节点上。采用隐式方法求解时，需要对Mortar接触内力进行线性化。以从节点接触力为例：

线性化结果为：

从上式可以看出，线性化结果需要分别计算从节点间隙的线性化，从节点平均外法线的线性化以及Mortar积分的线性化。由于Mortar积分的积分区域和变形过程相关，因此上述变量的线性化过程极其复杂，这里不做具体展开。

得到每个接触力的线性化结果后，最终的Mortar接触力线性化方程为：

维度为，和分别为从节点和主节点的个数，的具体形式为：

4 数值算例

SiPESC平台目前已经实现了基于Mortar的面-面接触算法，采用罚函数方法和PDASS拉格朗日乘子法，支持低阶面片（三角形，四边形）。如无特殊说明，本章节的算例的默认材料属性为弹性模量2E11，泊松比0.3，对标商用软件ANSYS，单元类型SOLID45（线性四面体，六面体），接触设置为SurfaceToSurface，LagrangeMultiplier，OnNodesSurfaceProjection，收敛精度为相对位移1E-6，后文不再赘述。

4.1 算例1 简单接触分片试验

接触分片试验（Contact Patch Test）用于检验接触算法能否正确传递均匀的接触变形和接触压力，大部分基于点-面离散的接触算法无法通过接触分片试验。本算例用于验证SiPESC平台实现的Mortar面-面接触算法能否通过接触分片试验。

图2 算例1有限元模型

如图2所示，两个尺寸为1x1x1的立方块对其摆放，中间预留0.01初始间隙。X和Y方向设置对称边界条件，Z方向底部约束，顶部施加-0.0101大小的指定位移，使之产生0.0001的法向接触变形。上方块采用线性四面体单元划分网格，定义为从面，下方块采用线性六面体单元划分网格，定义为主面，摩擦系数为0。

(1) 最终变形状态对比

图3 算例1最终状态变形（左ANSYS，右SiPESC）

(2) 最终应力状态对比

由于结构的所有边界条件均对称，因此只产生均匀的Z向正应力。

图4 算例1最终Z向正应力状态（左ANSYS，右SiPESC）

从应力结果的对比可以看出，SiPESC平台实现的Mortar面-面接触算法可以通过简单的接触分片试验。

4.2 算例2 半圆挤压

图5 算例2有限元模型

如图5所示，一个尺寸为1x0.5x0.1的方块在中心移除一块半径为0.25的半圆块，将该半圆块向上平移0.02作为初始状态。加载时，方块下表面全约束，半圆块上表面X，Z方向全约束，Y方向施加一个-0.029的指定位移用于模拟半圆块重新压入方块的过程。该问题的难点在于如何保持接触计算结果的对称性和光滑性，是一个经典的Benchmark算例。本算例采用20步加载，从面为半圆块，主面为方块，摩擦系数为0。

(1) 最终状态结果对比

图6 算例2最终状态X方向位移（上ANSYS，下SiPESC）

图7 算例2最终状态Mises应力云图（上ANSYS，下SiPESC）

由结果可见，SiPESC平台的Mortar接触算法与商用软件ANSYS对比结果良好，在保持结果的对称性以及应力的光滑性方面与ANSYS基本保持一致。

（2） 迭代收敛性对比

图8 SiPESC整体迭代残差历史

4.3 算例3 刚性球穿越柔性管道

图9 算例3有限元模型

如图9所示，一个直径为500的刚性球穿过一个外径相同，但内径略小的柔性管道。加载时，柔性管道的右端全约束，刚性球沿柔性管道轴线方向移动400。该问题的难点在于接触面是曲面，主从面的刚度差异较大，且在运动中的接触状态变化大，导致收敛困难。本算例采用100步加载，摩擦系数为0。

(1) 最终状态结果对比

图10 算例3最终状态X方向位移（左ANSYS，右SiPESC）

图11 算例3最终状态Mises应力（左ANSYS，右SiPESC）

(2) 迭代收敛性对比

图12 SiPESC整体迭代残差历史

由迭代收敛性的对比可见，SiPESC平台的Mortar接触算法采用了PDASS技术，相对于采用传统震荡条件的拉格朗日乘子算法，大幅度提高了复杂大滑移接触问题的收敛速率。

(3) 变形动画

4.4 总结

综上，通过三个典型的三维面-面接触算例的对比可以看出，SiPESC目前实现的MORTAR面-面接触算法和商用软件的同类型算法相比，精度相当，且收敛速度更快，具备一致性收敛速率。

5 结论

接触分析，特别是大变形有限滑移问题是工程数值模拟中的重要的分析内容，也是当前技术中的难点问题。SiPESC平台目前已经实现了高效的三维接触分析功能，支持点-面、面-面（MORTAR）等接触面离散方式，支持罚函数、拉格朗日乘子和增广拉格朗日等约束施加方法，多种单元类型（四面体/六面体）等，与商用软件精度对比良好。后续在此基础上将进一步开展摩擦接触，热力耦合接触等分析功能的研发。

文中ANSYS测试算例由项目合作方提供。