基于边界元法求解三维弹性摩擦接触问题

舒小敏， 李 坚

(中国航发湖南动力机械研究所，株洲 412002)

接触问题广泛存在于工程机械领域。在机械中各部件通常以接触的形式传递运动或载荷，而接触区常常应力集中，致使机械部件更早地出现疲劳断裂等问题，降低部件使用寿命，增加安全隐患。为降低这些影响及危害，接触问题得到大量研究。

早期主要通过理论解析法求解接触问题[1-3]，但解析法局限于简单几何体、简单加载及简单边界条件的接触问题。实际工程中，接触问题不仅模型复杂而且受力加载情况多变。工程中很少存在一个完全满足解析法的理想模型。因此，寻求数值算法求解接触问题，常用的接触数值求解算法包括有限元法[4-8]和边界元法[9-12]。

边界元法求解三维摩擦接触问题有30年的历史。Kong等[13]采用数学规划法求解，该方法将库仑摩擦锥近似为多边形的金字塔模型。Yamazaki等[14,15]采用罚函数法求解三维摩擦接触问题。Garrido等[16]采用直接法并利用增量技术求解。此后，Leahy等[17-19]进一步采用直接法并利用二次单元求解了弹性及弹塑性摩擦接触问题。国内也展开了相关算法研究，如数学规划法[20-22]、罚函数法[23]和直接法[24，25]。

利用边界元法求解三维摩擦接触问题，其中一个关键点在于如何确定滑移方向。即当出现相对滑移时，滑移方向如何确定。不同于二维摩擦接触问题，滑移方向在一条线上，只有两个切向滑移方向。三维中切向滑移方向位于一个平面内，难以确定。

边界元法求解摩擦接触问题常采用两种方式确定滑移方向。一种是利用切向面力得到滑移方向；
另一种是利用切向相对位移得到滑移方向。在第一次迭代时，文献[16,25]通常假设物体之间没有相对滑动，根据切向面力得到滑移方向，而后针对滑移点采用切向相对位移得到滑移方向。文献[16]采用该方法并不能保证收敛。同样在文献[17-19]中第一次迭代，假设物体之间没有滑动，根据切向面力得到滑移方向，但却没有说明在以后的迭代中如何得到相应的滑移方向。此外，文献[16-19]只考虑了由粘结接触点变成滑移接触点，却没有再对滑移接触点进行判断，也就是一旦接触点进入滑移就不再对其判断。针对这一点，文献[25]采用了摩擦耗能原理来判断是否仍是滑移点。

针对以上问题，本文采用直接法求解并利用滑移方向预测技术。该预测技术[26]在二维摩擦接触问题也得到采用，不过其采用无摩擦预测技术。即首先将摩擦接触问题假设为无摩擦接触问题，计算得到相应的切向相对位移，该切向相对位移就是以后的滑移方向。不过，假设为无摩擦接触可能导致预测的滑移方向失真。为此本文提出并采用无限摩擦预测技术，同时采用增量法及摩擦耗能原理识别非滑移点，避免一旦接触点进入滑移状态就不再对其进行判断的问题。最后，利用数值算例验证了本文算法的有效性和收敛性。

对于弹性问题，在不考虑体力的情况下，每个体的边界积分方程可表示为[27]

((P,Q)∈Γ)

(1)

式中uj和tj分别为位移和面力分量，Ui j和Ti j分别为位移和面力核函数或基本解，ci j(P)为关于边界布置系数矩阵，其中(i,j=1,2,3)。位移和面力基本解表达式为

(2)

(1-2ν)(njr,i-nir,j)

(3)

式中r为源点P到场点Q的距离，n为场点Q处的外法线，G和v分别为材料的剪切模量和泊松比。

在摩擦接触问题中，常采用库仑摩擦定律来描述界面的摩擦现象，如图1所示。

图1 三维摩擦库仑定律Fig.1 3D coulomb friction

其中库仑摩擦定律对接触区切向面力进行了约束，该约束可表示为

(4)

式中μ为摩擦系数，tτ为切向面力tτ1和tτ2合力，tn为法向面力，如图2所示。只要切向面力tτ足够小，小于摩擦极限，则出现粘结现象，如图1的点α。其相对切向位移为零，

(5)

图2 接触方向定义Fig.2 Contact problem definition

当摩擦力大于摩擦极限，则会出现滑移现象。此时，切向面力tτ的大小等于μ|tn|，其方向与切向相对运动方向相反，如图1的点β。切向摩擦力与切向相对运动方向相反，可表示为

(6)

此外，摩擦存在时，需满足摩擦耗能条件为

(7)

由边界积分方程(1)可知，在三维中每个节点可以配置3个边界积分方程(i,j=1,2,3)。对于非接触的节点，由于边界条件的存在，只存在3个未知量，3个边界积分方程，正好可解。然而对于可能接触区的节点，由于面力和位移均为未知量，共有6个变量，却仅有3个边界积分方程，系统方程不可解。需在每个接触节点上补充3个接触约束方程，保证系统方程可解。由于摩擦的存在，计算结果和加载历史相关，本文采用增量法及点对面(node -to -surface)[16,28,29]的方式施加接触约束。如图2所示，假设物体A上的点a和物体B上的点b接触，点b位于单元Q1Q2Q3Q4中(图3)。则可建立如下的接触约束补充方程。

图3 点对面施加接触约束Fig.3 Contact constraints by node -to -surface method

对于接触体A的粘结接触节点a，接触约束条件为法向间隙为零，切向力相对位移为零。

(8a)

(8b)

(8c)

对于接触体A的滑移接触节点a，接触约束条件为法向间隙为零，切向力满足库仑摩擦定律。

(9a)

(9b)

(9c)

对于接触体B的接触节点b，接触约束条件为法向和切向面力相等，

(10a)

(10b)

(10c)

对于非接触节点，不管是接触体A的节点a，还是接触体B的节点b,约束条件皆为所有方向上面力为零。

(11a)

(11b)

(11c)

(12a)

(12b)

(12c)

从式(8～12)可看出，任意接触状态的节点都可以补充三个约束方程，从而保证最后的系统方程可解。这种根据接触状态补充约束方程，进而直接求解的算法[16-19,24-26]，本文称为直接法。

采用滑移方向预测技术，可得到一个预测的滑移方向。在以后的迭代过程中若发生滑移，滑移方向一直取预测技术中得到的滑移方向。本文介绍两种预测技术。

5.1 无摩擦预测技术

该方法将摩擦接触问题假设为无摩擦接触问题，然后将所有载荷进行加载并计算求解。根据计算结果，可求出切向方向的相对位移，该相对位移就是滑移方向。该方法[26]在二维摩擦接触问题中应用，数值算例取得良好的结果。其中滑移方向或滑移角可表示为

(13)

不过，该方法存在一个问题，即在接触面为平面时可能会失效。因为无摩擦接触时，只有法向方向施加位移约束式(9a)，而切向约束式(9b，9c)为摩擦力为零约束，结果导致接触体位移约束不足，计算结果失真。如图4上面的接触体，只有法向位移(z方向)受到位移约束，切向没有位移约束(x和y方向)。即使最后方程可解，得到的位移也将失真，导致滑移方向预测失真。因此，该方法的应用受到一定的限制，必须要保证位移约束条件足够充分(限制刚体位移)才能使用。

图4 接触面为平面的接触问题Fig.4 Contact problem with a flat contact surface

5.2 无限摩擦预测技术

该方法将摩擦接触问题假设为无限摩擦接触问题，然后将所有载荷进行加载并计算求解。此时法向和切向都是位移约束条件,即式(8a～8c)，不会出现位移约束不足的情况。根据求解结果，可得到相应的切向合力方向，该方向的负方向就是相对滑移方向。其中滑移方向或滑移角方向可表示为

(14)

本文采用该方法进行滑移方向预测。

6.1 摩擦接触迭代收敛条件

(1) 对于可能接触区的非接触点，需满足非贯穿条件，即法向间隙gn≥0。否则非接触点变为接触点。

(2) 对于可能接触区的接触点，需满足压应力条件，即法向压力要求非正，即tn≤0。否则接触点变为非接触点。

只有上面四条同时满足，不出现违反的情况，摩擦接触求解才算收敛。对于第(4)条的收敛判断准则，文献[16-19]均没有采用。即认为一旦判断为滑移接触点，不管以后如何变化，一直为滑移接触点。事实上这种假设值得怀疑，在二维摩擦接触中，发现有部分滑移状态[26]，指的是在增量加载过程中，一开始为滑移点，随后又变为粘结点。因此，本文将耗能摩擦条件也加入到收敛性判断的准则中。

6.2 算法流程

三维摩擦接触算法流程如图5所示。在该算法中，根据压应力条件、非贯穿条件及接触点摩擦状态(粘结或滑移)是否变化，来判断当前增量加载是否收敛。由于增量步加载，当前步长不一定保证收敛，需判断是否需要细分。图5的迭代次数是否超过n，就是为了解决该问题。当迭代次数超过一定次数时，进行增量步细分。当所有增量载荷加载完，整个计算完成。本文在每个增量步加载前，采用无限摩擦预测技术得到滑移方向。

图5 三维摩擦接触算法流程Fig.5 Flowchart of 3D friction contact algorithm

为验证本文数值算法的有效性和收敛性，与有限元法进行比较。在有限元法中采用 surfure -to -surface接触及拉格朗日乘子法求解[8]。

7.1 立方块与基座的摩擦接触

如图6所示，弹性立方块和弹性基座接触。立方块冲头的尺寸为10×10×10，基座的尺寸为16×16×10。两者采用相同的材料参数为弹性模量E=200 GPa，泊松比v=0.3。基座的底面固定，冲头的上表面受到一个均布压力p=200 MPa。在边界元法中采用线性四边形单元进行离散，有限元法中采用线性六面体单元，并采用不同的摩擦系数μ=0.05，0.1和0.2。

图6 网格离散模型Fig.6 Discrete model

从图7～图9的接触压力可以看出，两种计算方法得到压力的最大值都随着摩擦系数的增大而增大，最大相对误差为4.79%，如图9所示。此外图中压力趋势基本一致，表明接触压力计算结果的可靠性。

图7 摩擦系数μ=0.05时的接触压力Fig.7 Contact pressure when friction coefficient μ=0.05

图8 摩擦系数μ=0.1时的接触压力Fig.8 Contact pressure when friction coefficient μ=0.1

图9 摩擦系数μ=0.2时的接触压力Fig.9 Contact pressure when friction coefficient μ=0.2

随着摩擦系数从0.05增加到0.2，从图10～图12可以看出，粘结的区域逐渐变大(有限元法的红色部分，边界元法的红点)。

图10 摩擦系数μ=0.05时的滑移和粘结区域Fig.10 Slip and stick area when friction coefficient μ=0.05

图11 摩擦系数μ=0.1时的滑移和粘结区域Fig.11 Slip and stick area when friction coefficient μ=0.1

图12 摩擦系数μ=0.2时的滑移和粘结区域Fig.12 Slip and stick area when friction coefficient μ=0.2

同时两种方法得到粘结区的面积基本相同，表明粘结区和滑移区的大小判断有效。

此外图10～图12中，箭头方向表示冲头接触面中接触点的相对滑移方向。可以看出，不管摩擦系数如何变化，其具有高度的对称性。这是由于 图6 的计算模型高度对称，必然要求滑移方向具有高度的对称性。图10～图12边界元法滑移方向的对称性，表明了本文无限摩擦滑移方向预测技术的有效性。

7.2 圆弧形冲头与基座的摩擦接触

弹性圆弧形冲头和弹性基座模型如图13所示。冲头的几何尺寸为16×16×5，基座的几何尺寸为16×16×8。两接触体采用相同的材料参数，弹性模量E=200 GPa，泊松比v=0.3。基座的底面固定，冲头的上表面受到一个均布压力p=200 MPa。在边界元法中采用线性四边形单元进行离散，有限元法中采用线性六面体单元离散(图14)，并采用不同的摩擦系数μ=0.05，0.1和0.2。

图13 圆弧形冲头和基座几何模型Fig.13 Geometric model of arc-shaped punch and base

图14 网格离散模型Fig.14 Discrete model

从图15～图17的接触压力可以看出，两种计算方法得到压力最大值的相对误差不超过0.52%。此外，图中压力趋势基本一致，表明接触压力计算结果的可靠性。

随着摩擦系数从0.05增加到0.2，从图18～图20可以看出，粘结的区域逐渐变大(有限元法的红色部分，边界元法的红点)。同时两种方法得到粘结区的面积基本相同，表明粘结区和滑移区的大小判断有效。

图15 摩擦系数μ=0.05时的接触压力Fig.15 Contact pressure when friction coefficient μ=0.05

图16 摩擦系数μ=0.1时的接触压力Fig.16 Contact pressure when friction coefficient μ=0.1

此外图18～图20中，箭头方向表示冲头接触面中接触点的相对滑移方向。可以看出，不管摩擦系数如何变化，其都具有对称性。这是由于 图13 的计算模型对称，必然要求计算滑移方向具有对称性。图18～图20中边界元法滑移方向的对称性，表明本文无限摩擦滑移方向预测技术的有效性。

图17 摩擦系数μ=0.2时的接触压力Fig.17 Contact pressure when friction coefficient μ=0.2

图18 摩擦系数μ=0.05时的滑移和粘结区域Fig.18 Slip and stick area when friction coefficient μ=0.05

图19 摩擦系数μ=0.1时的滑移和粘结区域Fig.19 Slip and stick area when friction coefficient μ=0.1

图20 摩擦系数μ=0.2时的滑移和粘结区域Fig.20 Slip and stick area when friction coefficient μ=0.2

边界元法基于直接法求解三维弹性摩擦接触问题，其收敛性一直难以保证。针对这一问题，本文采用增量法求解，并利用无限摩擦滑移方向预测技术。该预测技术先将摩擦接触问题假设为无限摩擦接触问题，物体之间不发生滑动，计算得到相应的切向面力，该切向面力的反方向即为滑移方向。在以后的迭代过程中若发生滑移，滑移方向一直取预测技术中滑移方向。同时采用摩擦耗能条件对进入滑移状态的接触点进行判断，避免接触点一旦进入滑移状态就不再对其进行判断的问题。不同摩擦系数下的数值计算结果(接触压力、粘结区滑移区的大小及滑移方向)都表明本文算法的有效性和收敛性及滑移方向预测技术的有效性。