考虑摩擦生热的扭动微动热力耦合数值分析

谭升旭，刘 娟，沈火明，陈扬帆

(西南交通大学 力学与工程学院， 成都 610031)

0 引言

扭动微动是指在交变载荷[1]作用下，接触面之间发生小幅度往复扭转而引起的相对运动，它发生在人体髋关节、膝关节等部位，也发生在许多工程部件中，如核电设备中的部件[2-4]、高速列车转向架的中心板[5-7]等。微动损伤是诸多复杂因素共同作用造成，许多问题尚未深入研究以了解损伤的本质和规律。扭动微动已有一些研究，如蔡振兵等[8-10]系统地研究了金属材料和复合材料的扭动微动磨损，使用润滑、涂层和气氛控制等条件减轻微动磨损。半解析方法[11]和ABAQUS数值分析[12]已应用于扭动微动的数值模拟，但关于微动的研究仍集中在切向模式上，其他模式的研究很少，尤其是扭动微动。目前，扭动微动磨损的力学行为和磨损机制仍有待揭示。

研究表明，影响材料微动磨损的参数有50多个[13]。法向载荷、摩擦因数、相对位移幅值等对微动过程的影响已得到广泛关注[12,14-15]，温度已被认定为影响微动磨损行为和损伤机制的重要因素[16-17],微动磨损过程中的温度会受环境温度和微动接触面间产生的摩擦耗散能的影响。Attia[18-19]的研究表明，摩擦引起的温升会显著影响微动行为。此外，Wen等[20]的研究表明，由摩擦耗散引起的温升，会对微动过程中氧化物碎片的形成、三体层的发展变化，以及表面的材料性质造成影响。微动磨损速度的变化是高温促进的氧化物形成而造成的[21]。试验中很难测量接触面之间的温度，所以需要预测技术来促进微动过程中温度作用机理的研究。许多研究人员(例如Jaeger[22],Archard[23],Tian[24]等)通过引入Péclet数的分析模型，证明了表面温度与材料的热性能、接触面的相对移动速度有关。文献[25-27]建立了数值模型，分析了温升和材料变形之间的关系。

综上所述，温升在微动机理研究中具有重要意义。因此，本文基于数值方法和摩擦生热机理，建立预测扭动微动热力耦合行为数值分析模型，研究表面Von Mises 应力分布、表面温升和次表面温升，以及表面塑性变形与温升的关系，进而探索温升对扭动微动磨损的影响机制，有助于进一步揭示扭动微动磨损损伤机理。

1 模型和方法

1.1 数值模型

图1 数值模型结构示意图

采用主-从接触算法来满足非穿透接触约束，在迭代过程中，定义最大的穿透范围，且只允许主表面上的节点侵入从表面，因此不会出现导致结果不准确的过度穿透。

法向载荷直接施加会导致局部节点刚体位移，影响模型收敛，因此定义一个参考点，将其和上接触体的顶面绑定约束。在参考点上施加恒定的法向载荷和水平的正弦扭转循环角位移，使半球面与平面发生扭转，实现微动循环的加载[28]，平面的下底部则进行完全约束。载荷加载历程如图2所示，设置较短时间的加载和卸载分析步使接触状态平稳过渡，增强模型的收敛性。

图2 加载历程曲线

此外，基于所有克服摩擦所做的功都以热源形式消散的假设，热量将从接触区域流入接触体，则可以表示为式(1)的热功率密度，该项与局部接触压力、位置和时间的函数成正比。

Q(t)=Ff|v(t)|

(1)

式中：Ff为摩擦力；v(t)为接触节点的滑动速度。描述运动的位置函数s(t)可以表示为：

s(t)=rφ(t)=rθsin(ωt)

(2)

式中：ω=2πf，f为微动频率；θ为角位移幅度；r为距接触中心的距离。速度如下：

(3)

摩擦力可以表示为：

Ff=Pμ

(4)

故总热源功率可定义为：

Q(t)=Pμrθω|cos(ωt)|

(5)

1.2 模型验证

为验证数值模型的准确性，与Johnson的解析解[28]进行了比较。在接触载荷的作用下，接触区压力分布表示为：

(6)

式中：

(7)

(8)

式中：a0为接触面积的半径；p0为最大接触压力；E1和E2为两材料的弹性模量；μ1和μ2为两材料的泊松比；R1和R2为两试样接触表面的曲率半径；x为距离接触中心距离。对于切向接触问题，基于切向力与压力没有相互作用的假设，剪切力可以表示为：

(9)

式中：H( )为Heaviside的阶跃函数；μf为摩擦因数，其值为0.6；c为黏着半径。

摩擦副GCr15/2024合金的材料属参数如表1所示，所列分别是弹性模量、抗拉强度、屈服强度、泊松比、密度，以及热膨胀系数和比热容。

表1 摩擦副的力学和物理性能参数

图3(a)和(b)分别表示压力和剪切应力的解析解与数值结果。可以看出，数值结果和解析结果之间有很好的一致性。

图3 解析解与数值结果曲线

2 结果和讨论

2.1 Von Mises 应力

普遍认为，第一屈服准则对材料磨损机制具有重要指导作用。Johnson[28]已经证明了剪切应力、Von Mises 应力分布以及最大应力的确切位置。在扭动微动的半解析法(SAM)研究中，Liu等[12]已经证明，在扭动微动磨损机制中，表面的首次屈服主要表现为黏着磨损，接触面材料损伤加剧表现为明显的塑性变形。

为揭示表面塑性变形和磨损，在法向载荷P=100 N作用下，分别取角位移幅值θ为0.15°，1°，8°，呈现出扭动微动的3种不同滑移状态，如图4所示。Von Mises 应力分布的数值结果如图5所示：在部分滑动状态下整个接触面的应力较小，没有发生屈服[图5(a)]，而对于混合滑移状态和完全滑移状态，高应力几乎遍布整个表面，且发生了明显的屈服[图5(b)和5(c)]；在混合滑动状态，高应力主要表现在发生磨损的环形滑动区域[30]；在完全滑移状态，整个接触区域的显著屈服意味着严重的表面磨损，这与Cai等[31]的实验结果一致。在混合滑移状态的微动循环期间，Von Mises 应力值并不总是大于屈服强度，这表明表面磨损并不是主要的损坏机制；事实上，磨损机制主要是与疲劳相关的微裂纹，因为屈服主要发生在次表面[12]。

图4 3种不同滑移状态微动曲线

图5 微动周期内 Von Mises 应力的分布曲线

2.2 温升分布分析

角位移幅值是影响摩擦生热的主要因素之一，图6表示混合滑移状态和完全滑移状态下表面和次表面温升的预测结果，表面温度的距离接触中心x=0.2 mm的位置，次表面是距离表面监测点h=0.2 mm深度的位置。显然，温升受位移幅值的影响较大，表面和次表面温升随着扭转角位移幅值的增加而显著增加，且完全滑移状态的温升大于混合滑移状态的温升；次表面温升较接触表面温升显著降低；温升在很短的微动时间内显著增加。图7表示高频率(f=200 Hz)和低频率(θ=100 Hz)下10 000 s的温升随时间的变化。结果表明，在这2种情况下，前几个周期的温升都急剧增加，而上升速度随时间的增加迅速下降[如图7(a)]。图7(b)、(c)和(d)分别显示了3个不同阶段的表面温升，可以看出3个阶段的温升增加速度差异。在高频率f=200 Hz的扭动作用下，最大增长速度出现在A阶段，B阶段的增速大于C阶段。经过数千秒后，温升几乎达到一个稳定值。

图6 表面(x=0.2 mm)和次表面(h=0.2 mm)温升随时间变化的预测曲线(P=100 N, f=100 Hz)

图7 较长时间的温升曲线 (x=0.2 mm,P=100 N,θ=5°)

2.3 塑性变形

2.1节的结果表明，表面屈服引起了大范围的表面塑性变形，表面磨损是完全滑移状态的主要损伤机制。且在2.2节的分析中表明，微动频率高时表面温升大，特别是在完全滑移状态下，表面温升明显。相关研究认为，接触中的温度会影响磨屑形成稳定层和从接触区排出的方式[20]。

为了研究温升对扭动微动表面磨损的影响，本文预测了表面塑性变形与温升之间的关系，在P=100 N,θ=5°的工况下，探究温度变化对材料塑性变形的影响。如图8所示，温升影响了表面塑性变形行为，在循环一周时，塑性变形变化尚微小，而随着微动周期的增加，扭动微动中的塑性变形随着温度的升高显著增大。此外，图9显示了在距离接触中心0.2 mm处监测到的扭动微动中的塑性变形随随温度升高的变化。随着温度的升高，表面塑性变形显著增加，但增加的速度缓慢下降。

图8 不同微动循环下塑性应变分布曲线(P=100 N，θ=5°)

图9 表面塑性变形随温升的变化曲线(x=0.2 mm,P=100 N,θ =5°)

3 结论

基于摩擦耗散的摩擦生热机理，建立了预测扭动微动的热力耦合数值分析模型。在模型验证基础上，研究了表面Von Mises 应力分布、表面温升和次表面温升，以及表面塑性变形与温升间的关系。结果表明：Von Mises 应力分析结果显示表面磨损是完全滑移状态的主要损伤机制，而与疲劳相关的微裂纹是混合滑移状态的主要磨损机制，该结论与以往研究中的实验现象一致；表面和次表面温升均对角位移幅值和微动频率敏感，但表面温升比次表面温升更显著，且温升变化速度先增大后急剧减小，在数千秒后便达到稳定状态；表面温升显著影响表面塑性变形，进而影响表面磨损，因此摩擦生热引起的表面温升明显影响扭动微动中的表面磨损行为。