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

摩擦数值算法简介

在过去的三十年间，计算接触力学领域发展了多种用于求解摩擦接触问题的算法。在对摩擦接触问题进行数值模拟时，一个主要困难是摩擦力与切向滑动量之间的本构关系是非光滑的，本构函数在某些点上不可微分，从而造成数值计算中迭代收敛困难，这个问题可以通过对摩擦本构关系的规则化来解决。

目前已有多种迭代方案用于带摩擦的接触分析，可大致分为以下几类：试探-校核算法、基于塑性理论中的弹塑性类比方法、基于优化理论的数学规划方法。后两类方法以严密的数学理论为基础，其可靠性高于试探-校核算法。

试探-校核算法通常是应用于小变形情况。对于摩擦接触问题，解的唯一性和存在性均不能保证，因此试探-校核算法在很多情况下并不可靠。但试探-校核算法仍然得到了成功的应用，在显式有限元分析中能够获得合理的结果。

近几十年，Coulomb定律和其他摩擦本构关系被纳入塑性理论的研究范畴，基于弹塑性理论的返回映射方案已成功应用于有限元摩擦接触分析，该方案使算法的收敛行为和可靠性产生了本质的提高。返回映射方案最初用于材料非线性问题，用以积分弹塑性本构关系。将摩擦定律类比为弹塑性本构关系，就可以将返回映射方案应用于带摩擦的接触分析。由于摩擦本构关系的非关联性，返回映射方案所得到的切线刚度矩阵通常不对称，增加了数值求解的难度。除返回映射方案之外，其他几种来源于塑性理论的方法，例如屈服极限拉氏乘子法，也已用于摩擦接触问题的数值分析中。

数学规划方案在摩擦的模拟中也有较为成功的应用。此外，其他的一些途径，例如内点法、屏障法、连续近似法等，也被尝试用于摩擦接触分析中。

摩擦接触分析中，需要区分滑动或粘结两种状态。对于粘结状态，通常可以使用一个参数（罚值）εT将其引入到罚函数或者拉格朗日格式，εT相当于切向滑动的刚度。但是这种处理使得静摩擦力的大小依赖εT的数值。对于过小的εT，可能导致物理上有问题的计算结果，即接触界面出现类似粘滑效应的行为，与实际情况不符。虽然可以采用较大的εT来解决这个问题，但是εT取值过大将导致方程组出现病态或者迭代流程失去健壮性。

8

摩擦定律与摩擦定律的正则化

 8.1 库伦摩擦定律

两接触表面的切向力一旦超过了某一个极限值，两表面不再保持粘结状态，而是发生相对的滑动。需要使用适当摩擦模型来描述这种切向行为。在实际工程中，库伦定律因其简单和适用性被广泛采用。库仑定律的表达式如下：

其中μ为摩擦系数，在经典库仑定律中被视为常数，由两接触面的材料种类决定；pN为法向接触压力；ġT为相对滑动速度；tT为切向摩擦力。

库仑定律表明，当摩擦力小于μpN时，无相对滑动，处于粘结状态；当有相对滑动时，摩擦力必然达到极限值μpN。

摩擦系数通常受多种因素的影响，例如表面粗糙度、相对滑动速度ġT、法向接触压力pN、环境温度θ等。如果考虑这些影响，可以得到经典库仑定律的改进形式：

 8.2 库仑定律的规则化

经典库伦摩擦定律是不可微的，在相对滑动速度ġT=0处发生由粘接至滑动的突变，当ġT反向时，摩擦力tT也立即发生反转。这种突变将导致数值计算中迭代的收敛困难。因此提出了一些规则化（即光滑化）的模型来代替经典库伦摩擦模型，这些规则化模型可以表达为如下形式

式中，函数φ(ġT)用来描述由粘结状态至滑动状态的光滑过渡。不同的规则化模型中，φ(ġT)可采用多种形式，分别对应于不同的规则化模型。

• 平方根规则化模型

• 双曲正切规则化模型

• 反正切规则化模型
• 分段直线模型

在以上模型中，标量参数ε为规则化参数，当ε→0即为经典库伦摩擦定律。ε不宜取值过大，否则规则化模型将与实际的接触界面粘结-滑动规律有较大出入。

上述规则化模型的前三种能够提供光滑的摩擦力-滑动量曲线，但第四种模型类似于理想弹塑性应力-应变关系，便于采用弹塑性力学中的相关方案处理。

当采用规则化的模型后，分析中不存在粘结状态，而统一的应用上面的数学描述来进行摩擦分析。这在物理上可以更真实地描述实际摩擦现象，因为两个接触面上都存在一定的不平度，因此完全粘结的接触状态是不存在的。只要有切向摩擦力的存在，总要伴随发生一定的相对滑动，而规则化的摩擦模型正好可以描述此类物理现象。

 8.3 摩擦定律的弹塑性类比

弹塑性类比方案的关键思想在于将切向滑动gT分离为弹性滑动部分gTe和塑性滑动部分gTp，如图2。

图2：切向滑移量的分离

这种分离方可视为对摩擦定律的规则化。但是它也有实际的物理含义：gTe代表接触界面上的弹性微小位移，卸载后gTe和由gTe引起的变形将恢复为0；gTp代表不可恢复的位移，在卸载后仍然存在。

对于切向摩擦的弹性部分，可采用以下各向同性线弹性模型，

式中，cT为界面弹性常数。

也可采用各向异性线弹性模型

式中，cT为弹性系数张量。

切向摩擦的塑性部分也有类似的本构方程，可根据弹塑性理论的标准概念导出。与弹塑性理论进行类比，可以建立屈服函数fs(tT)和相应的屈服准则，如下：

屈服函数fs(tT)可以简单的表达为：

上式等价于经典库伦摩擦定律。

利用最大耗散原理，得

然后可以导出塑性滑动的流动法则

此外，还可以得到加载-卸载准则

根据上式可确定参数ƛ.

9

商用有限元软件中的接触和摩擦算法

 9.1 ABAQUS中的接触算法

ABAQUS缺省使用拉氏乘子法施加接触约束，用户也可选择使用罚函数或增广拉格朗日法。在ABAQUS接触模拟中，要在各个构件上建立表面，定义可能会发生接触的一对表面（接触对），然后按照本构关系求出表面间的法向和切向相互作用。

ABAQUS可以自动确定二维和三维实体单元的自由面，并利用它们建立一个表面。对于壳单元和刚性单元，必须指明是正法向的面还是负法向的面。

ABAQUS应用单纯的主控-从属接触算法，执行单边接触检查：从属表面的节点不能穿透主控表面的任何部分，这种算法对于主控表面没有限制，它可以在从属表面的节点之间穿透从属表面。因此必须正确选择从属与主控表面，要遵循的简单原则是：

（1）从属表面应该是网格划分的更精细的表面。

（2）如果网格密度相近，从属表面应当由更柔软的材料组成。

定义接触对时，须确定两表面之间的相对滑动是小量还是有限值。默认设置是有限滑动状态。当应用小滑动公式时，从模拟开始就建立从属表面和主控表面之间的关系，确定主控表面哪个部分与从属表面的结点发生作用。这种相互关系在分析过程中保持不变。对于有限滑动情况，两个变形表面间的有限滑动模拟仅能用于二维问题，变形面和刚性面的有限滑动计算可用于三维问题。

图3为ABAQUS接触分析的流程。ABAQUS在增量步开始前检查所有接触对的状态，判断从属结点是否进入接触，对于接触的结点，判断是滑动还是粘接，然后施加相应的约束。

图3：ABAQUS接触算法流程

在检查力或力矩的平衡前，首先检查从属结点接触条件的变化，如果检测到当前迭代步的任何接触状态改变，则标识为严重不连续迭代，不进行平衡检验，而是改变接触约束进入第二次迭代。直至接触状态不再变化，才进行平衡收敛检查，如果收敛性不满足，则进行下一次迭代。

 9.2 MSC.MARC中的接触算法

MARC软件提供了基于直接约束的接触算法，是解决所有接触问题的通用方法。特别是对大面积接触，以及事先无法预知接触发生区域的接触问题，程序能根据物体的运动约束和相互作用自动探测接触区域。当发生接触时，使用边界条件直接约束运动体，两者的运动约束转化成了节点自由度的约束和节点力的约束。

MARC中的接触检查默认采用双边检查，为提高检测速度，用户也可选择单边检查。MARC不仅可以处理变形面和刚性面之间的有限滑动，对于变形面与变形面间的有限滑动也能很好的处理。

MARC软件接触算法基本流程如图4所示。从图4可见，整个过程包括：

• 定义接触体。
• 探测接触。
• 施加接触约束。
• 模拟摩擦。
• 修改接触约束。
• 检查约束的变化。
• 判断分离和穿透。

此外，对于接触区域事先已知的情况，用户还可使用间隙单元(gap单元)方案，或者基于罚函数方法，通过用户子程序USPRNG对接触界面施加非线性弹簧。

图4：MARC接触算法流程

作者简介

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