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

一、接触分析的挑战性

接触是在固体力学各个领域中普遍存在的问题。对自然界中许多物理问题的描述都涉及接触现象。例如零部件装配时的配合，橡胶密封元件的防漏，轮胎与地面的相互作用，撞击问题以及压力加工行业的大量成型工艺过程等。

接触过程中两个物体在接触界面上的相互作用是复杂的力学现象，同时也是它们损伤直至失效破坏的重要原因。从力学分析角度看，接触是边界条件高度非线性的复杂问题，需要准确追踪接触前多个物体的运动以及接触发生后这些物体之间的相互作用，同时包括正确模拟接触面之间的摩擦行为和可能存在的接触间隙传热。

其中极少数的接触问题可以解析处理，绝大多数接触问题只能采用有限元、离散元、边界元等数值方法进行模拟，其中有限元法的应用最为广泛。对接触全过程进行有限元仿真，现在不仅可以实现，而且正逐步成为CAE/CAM的重要组成部分。

在实际工程中，有限元接触分析的计算结果经常用于对某些设计参数进行优化设计，例如对轮胎进行结构优化以提高安全性和寿命。如果采用基于梯度的优化算法，需要得到力学变量（位移、应力、接触反力分布状况等）相对于设计参数（材料、尺寸、形状、拓扑结构等）的变化曲线和相应的敏度（梯度）。

对于无摩擦接触情况，现有的有限元接触算法，例如拉氏乘子法、罚函数法等，能够得到足够稳定的敏度数据；但是对于有摩擦接触情况，如果不采取一些特殊的处理，则很难得到稳定的数值结果，梯度的数值通常随载荷和网格的改变而发生剧烈的振荡，不具备可用性。

虽然接触力学和相关的数值方法已经广泛应用于工程开发和科学研究，但对于接触和摩擦的物理机制，目前尚未有完全的理解。

从工程的观点来看，计算机技术和计算方法的发展，使我们能够更精确的分析接触问题以适应工程需要。对接触问题的仿真和模拟在工程设计的多个方面已经发挥了重要的作用，例如减少磨损、降低噪声和提高安全性等。

但从科研的观点来看，现有的接触数值算法流程过于复杂，需要耗费大量的内存空间和计算时间，一般经过反复的校核、修正才有可能得到符合实际情况的计算结果。迄今为止，对于带摩擦的接触分析，当前各种商用有限元软件经常不能给出精确可靠的结果。开发稳定、高效、健壮的接触算法仍然是一个亟待解决的问题。

以上挑战性不仅来自接触过程中复杂的变形和受力状况，更主要的原因是接触界面的边界条件非线性。接触界面的非线性来源于两个方面：

1、接触界面事先未知。接触界面的区域大小和相互位置以及接触状态不仅事先都是未知的，而且是随时间变化的，需要在求解过程中确定。

2、接触条件的非线性。接触条件的内容包括：两个互相接触的物体不可相互侵入；接触力的法向分量只能是压力；切向接触的摩擦条件。这些条件区别于一般约束条件，其特点是单边性的不等式约束，具有强烈的非线性。

接触界面的事先未知性和接触条件的不等式约束决定了接触分析过程需要经常插入对接触区域的搜索，需要多次迭代求解以确定接触压力和摩擦力。

另外，接触过程常常涉及材料非线性和几何非线性。例如，汽车轮胎与路面的接触是接触力学中最典型的实际工程问题，进行数值模拟时，必须考虑由于大变形引起的几何非线性，为得到可靠的计算结果，还应使用复杂的非线性材料本构关系。因此，通常要求有限元接触算法具备同时处理三种非线性（材料、几何、边界条件非线性）的能力。

二、接触问题的约束条件

1、不可贯入性

如图1，考虑两个物体B^I(I=A,B)互相接触的情况，物体所占据空间域为Ω^I∈R³。物体B^I的表面Γ^I由三部分组成：Γ_σ^I上面力已知；Γ_u^I上位移已知；Γ_c^I则是两个物体的接触面。

图1 接触体之间的法向间隙

接触物体在运动学方面需要满足不可贯入性要求，不可贯入性是指物体B^A和B^B的位形在变形和运动过程中不允许相互贯穿（侵入和覆盖），可用下式表达：

式中，xI(I=A,B)指的是物体BI上表面各点对应于变形后位形的坐标，即欧拉坐标：

式中，n^A是物体B^A表面的外法线方向单位矢量。假定接触边界描述了一个局部外凸的区域，我们可以将Γ^B上的任一点x^B与Γ^A上的某一点相关联。ẍ^A(x头上两点代表-，下同)是物体B^A表面上距离x^B最近的点，参见图1，二者距离用下式表示：

该距离可用于定义接触体B^A和B^B之间的法向间隙。如果ẍ^A已知，不可贯入条件可用以下不等式约束来表示

对于变形体与刚性表面接触的情况，上式仍然成立，此时x^A≡X^A，n^A≡N^A。

2、 法向接触力为压力

接触物体在动力学方面需要满足法向接触力为压力的条件。在不考虑接触面间的粘附的情况下，物体之间的法向接触力只可能是压力，不能为拉力。

如果gN=0，法向接触力p_N=p^A_N=p^B_N<0。如果gN>0，意味着两个物体间存在间隙，此时pN=0。即

上式就是无摩擦接触问题Hertz-Signorini-Moreau条件。在优化理论中，此种形式的约束条件称为Kuhn-Tucher条件或者Karsh-Kuhn-Tucher条件。

三、Coulomb 摩擦定律和摩擦机理

摩擦是由于两接触表面相互作用而引起的，其结果是产生运动阻力。当两个表面相对运动时，摩擦力将做负功，在接触表面上产生能量的损耗。在工程分析中，Coulomb模型因其简单和适用性而被广泛应用。摩擦的物理机制非常复杂，最终可追溯到原子尺度

1、Coulomb摩擦三定律

法国工程师Amontons于1699年提出了两条基本的摩擦定律。这两条定律已为实验所证实，能适用于大多数条件，但是也有一些显著的例外。

第一定律：摩擦力与两接触体表观接触面积无关。

实际上，任何表面从微观上几乎都是粗糙的，实际接触面积只占表观接触面积的很小一部分，摩擦力的大小仅与实际接触面积有关。

第二定律：摩擦力tT与法向载荷f_N成正比。

式中，μ为常数，即摩擦系数，必须指出，仅仅对于给定的一对接触滑动材料和一组给定的周围条件，摩擦系数才是常数。材料不同、周围条件（温度、湿度、真空度）不同，摩擦系数也不同。

此外，Coulomb于1785年提出了第三定律，即动摩擦力几乎与滑动速度无关。

目前有限元软件中最通用的切向摩擦本构关系是经典Coulomb摩擦模型，该模型体现了以上三个基本定律。也有研究者提出另外一些模型，能够考虑界面上微观的力学现象或者摩擦非局部特性。

2、摩擦机理

当两个表面相互压紧时，会在接触区的某些部分发生粘着，这是引起摩擦的表面作用的一种形式。如果两接触表面产生相对运动，表面的材料微凸体将发生变形和位移来适应相对运动，这种变形和位移将产生运动阻力，这是引起摩擦的表面作用的另一种形式。

当两个接触表面相对运动时，材料微凸体可能发生弹性变形、塑性变形或者断裂。塑性变形总是能带来能量的耗散，在大多数实际情况下，这种能量损耗占金属摩擦的大部分。当表面相互作用为粘着形式时则必然发生断裂，当互嵌的微凸体相对运动时也会引起断裂，然而，对于大多数金属，断裂引起的能量耗散小于塑性变形。

金属材料发生弹性变形所需要的能量绝大部分可以回复，因此弹性能量耗散与塑性能量耗散相比可以忽略不计。但是橡胶材料在发生弹性变形后，由于产生弹性滞后，显示出很大的不可逆的能量耗散，在某些情况下，这是摩擦的主因。

综上所述，引起摩擦的表面相互作用有两个来源：即粘着和材料位移。发生微观的弹性和塑性变形以及断裂，会引起能量的耗散。这些表面作用涉及到界面的原子间复杂的化学和电磁作用，因此一些研究者试图将承受切向载荷的接触界面考虑为一层独立的物质。

Bowden和Tabor提出了简单粘着磨擦理论，其出发点是：当金属表面压紧时，微凸体顶端相互接触。接触着的微凸体上压力很高，导致塑性变形，令接触面积增大到实际接触面积A恰好能法向支承载荷f_N为止。因此，对于屈服压力为p₀的理想的弹塑性材料，有

简单粘着摩擦理论认为，在接触表面紧密接触区会发生牢固的粘着。如为剪断节点所需的单位接触面积上的力，而t_T为摩擦力，则有

其中p_e是考虑到坚硬微凸体在较软的表面上“犁沟”所需的力而引入的附加项。此理论可以解释两条摩擦定律：摩擦力与表观接触面积无关；摩擦力与载荷成正比。

实验表明，对于高真空中的洁净金属表面，可能获得很大的摩擦力，表明实际接触表面比简单粘着理论所指出的要高得多。在简单粘着理论中，认为A决定于屈服压力p₀和法向载荷f_N。对于静态接触，这大致是正确的。但是在摩擦的情况下，还作用有切向力，屈服取决于正应力和剪应力的复合作用。所以Bowden和Tabor进一步考虑了复合应力对微凸体接点实际接触面积的影响，提出了修正的粘着摩擦理论：

_N/p₀为仅考虑法向载荷影响而得出的接触面积，而α(t_T/p₀)²表示由剪力或摩擦力产生的增量。

Green则提出微凸体塑性相互作用理论，后来由Edwards和Halling加以推广，他们考虑了相对滑动时微凸体上法向应力和切向应力随时间的变化，将摩擦系数定义为所有接触微凸体的瞬时剪力之和除以所有接触微凸体的瞬时法向应力之和。

微凸体相互作用理论与粘着理论并无原则性矛盾，事实上它们能够得出相同的摩擦系数表达式。微凸体相互作用理论更加完善，可进一步推广从而适用于具有实际微凸体高度分布并考虑加工硬化的表面摩擦。粘着理论的物理现实性稍差，但是它能够在各种条件下得到正确的摩擦系数值，且在分析上远比微凸体相互作用理论简单。