非牛顿润滑的滚滑线接触表面摩擦特性

李 玲，张 东，王晶晶，云强强，杨亚兰，阮晓光+

(1.西安建筑科技大学 机电工程学院，陕西 西安 710055；2.西安陕鼓动力股份有限公司，陕西 西安 710611)

0 引言

滚滑并存的线接触粗糙表面广泛存在于机械设备中，如齿轮、滚滑复合导轨、圆柱滚子轴承等[1-3]。为减缓这些线接触粗糙表面的摩擦与磨损，通常在接触面间添加润滑油，这样不仅使表面的摩擦与磨损状况得到改善，还能降低机械系统的振动，从而提高设备的工作稳定性。此外，这些机械结构往往处于混合润滑的工作状态[4]。因此，研究混合润滑状态下滚滑线接触粗糙表面的摩擦润滑特性对机械结构表面的润滑设计与润滑状态预测尤为重要。

关于混合弹流润滑的研究，JOHNSON等[5]首次提出了法向外载荷由粗糙微凸体与润滑油膜共同承担的载荷分配思想，建立了混合润滑状态下的统计学模型。其中，微凸体的承载量由GREENWOOD等[6]建立的微观粗糙接触模型(GW模型)求解，油膜的动压力通过求解Reynolds方程获得；GELINCK等[7]基于该思想数值拟合出了GW模型下快速计算微凸体接触压力的公式，并用其预测了线接触润滑表面的Stribeck曲线。肖会芳等[8]基于该载荷分配机制建立了混合润滑状态下粗糙界面动摩擦特性的预测模型，研究了界面油膜厚度、摩擦系数和临界滑动速度随粗糙度参数、载荷以及润滑油属性参数的变化规律；此外，将润滑接触界面的法向刚度等效为粗糙固体接触刚度与油膜接触刚度的并联连接模型，提出了一种混合润滑粗糙界面法向接触刚度的计算模型[9]。李玲等[10]首先建立了粗糙表面固体的接触刚度模型，并利用固体刚度求解了液体油膜的等效厚度，从而建立了混合润滑界面的法向接触刚度模型，该模型解决了润滑界面油膜厚度测量难的问题。AKCHURIN等[11]通过测量粗糙接触表面的真实粗糙度，基于载荷分配思想建立了一种模拟混合润滑状态下滑动线接触表面摩擦系数的确定性模型。AKBARZADEH等[12]根据界面能量方程分别求解了表面温度与油膜温度，建立了考虑粗糙度的直齿轮热效应弹流润滑分析模型，模拟了载荷、速度、粗糙度以及温度对齿轮接触表面摩擦系数的影响规律。董国忠等[13]建立了考虑粗糙表面微凸体的弹流润滑计算模型，研究了微凸体曲率半径、载荷以及粗糙度对弹流润滑区摩擦特性的影响规律。刘明勇等[14]建立了一种有限长线接触混合润滑模型，分析了横向、纵向以及二维规则表面粗糙度的波长、幅值和工况变化对润滑摩擦特性的影响。

综合上述文献可知，目前关于混合润滑接触表面摩擦特性的研究，主要集中在微凸体的纯弹性变形。然而，混合润滑同时具备边界润滑和流体动力润滑的特征，微凸体会产生弹塑性和完全塑性的接触变形。此外，更为重要的一点是，随着现代机械装备向着高速、重载的方向发展，其工作条件越来越恶劣，使得实际工程中的润滑油往往呈现出显著的非牛顿特性[15-16]，这意味着润滑油的剪切稀化行为对接触表面摩擦特性的影响不可忽略。目前，关于非牛顿混合润滑线接触粗糙表面摩擦润滑特性的研究主要采用数值迭代法同时求解Reynolds方程、表面接触变形方程、载荷平衡方程以及描述油膜流变特性的方程，但数值法比较复杂、且难以收敛，导致所建模型不能快速预测实际工程中润滑界面的摩擦学特性，降低了理论仿真的实用性。

为此，本文分别采用KOGUT等[17]提出的微观接触模型(KE模型)描述干摩擦粗糙表面的接触，Carreau流变模型[18]表征润滑油的非牛顿特性，并基于法向外载荷由润滑油膜与微凸体共同承担的载荷分配思想，建立了一种可快速预测混合润滑状态下滚滑线接触表面的摩擦系数模型，分析了载荷、粗糙度、润滑油入口黏度以及微观粗糙接触模型对摩擦系数的影响规律。该模型充分考虑了微凸体的多种变形机制(包括弹性、弹塑性以及塑性变形)和润滑油的非牛顿特性，使得模型预测结果更符合实际工况。该模型具有较强的实用性，且求解过程简单、易于收敛、耗时短。

1 滚滑线接触表面润滑模型

1.1 问题描述

两干摩擦粗糙表面的接触可以简化为一个等效粗糙表面与一个刚性平面的接触[19]，故将润滑状态下两粗糙圆柱体线接触模型等效为如图1所示的模型，粗糙表面与刚性平面间填满润滑油，等效半径为R、等效弹性模量为E*的圆柱体在外载荷Fs的作用下与刚性平面接触，产生宽度为2r的接触区域，同时粗糙圆柱体以滚动速度ur和滑动速度us在刚性平面上做滚滑运动。

1.2 摩擦系数

基于JOHNSON等[5]提出的载荷分配机制，混合润滑状态下表面承受的法向外载荷Fs由油膜压力Fm与微凸体接触力Fa共同承担，

(1)

式中γ1和γ2分别为润滑油膜和微凸体的承载因子，且1/γ1+1/γ2=1。

同理，表面总摩擦力Ffs等于微凸体摩擦力Ffa与润滑油膜摩擦力Ffm之和，即

Ffs=Ffa+Ffm。

(2)

假设表面不同微凸体均具有相同的摩擦系数[20]，则微凸体的摩擦力可表示为：

(3)

式中fa为微凸体的平均摩擦系数。

文献[21-22]的研究表明，Carreau流变模型能够很好地描述实际工程中绝大多数润滑油的剪切稀化非牛顿流变特性，基于该模型获得的数值解更加贴合实际。其中，润滑油的广义黏度ηg与剪应变率γ0为幂函数的关系：

(4)

式中：j为幂指数;γ0=us/hc，hc为被入口区热效应因子校正后的非牛顿润滑中心膜厚，Gm为润滑油的有效剪切模量，η为润滑油在工作压力下的黏度，η0为润滑油入口黏度，r为接触半宽，α为待定常数，P为赫兹接触压力，Pmax为其最大值，x表示与接触区域中心的距离，L0为线接触长度。

通过在线接触区域A内对润滑油的剪切应力τ进行积分，便可获得润滑油的剪切摩擦力

(5)

则非牛顿润滑油膜的摩擦系数可表示为：

(6)

式中:A为线接触区域的面积，τ为润滑油的剪切应力。

对于牛顿润滑剂，其摩擦系数可表示为[23]：

(7)

式中hNct表示被热效应因子校正后的牛顿润滑中心膜厚。

(8)

于是基于载荷分配思想，非牛顿润滑接触表面的总摩擦系数可表示为:

(9)

1.3 弹流润滑模型

假设整个接触区域内的油膜厚度相同且等于中心油膜厚度，MOES[25]获得了等温条件下光滑表面线接触区的牛顿中心油膜厚度方程为：

(10)

式(10)仅适用于两光滑接触表面间中心油膜厚度的预测。为了考虑表面粗糙度的影响，基于载荷分配思想采用E*/γ1代替式中的E*，Fs/γ1代替Fs，获得了粗糙表面间中心油膜厚度为

(11)

油膜厚度的计算主要取决于接触区的入口条件[26]，为了引入入口区粘性剪切热效应对油膜厚度的影响，GUPTA等[27]提出了热效应因子Tt的概念，并将其定义为热油膜厚度与等温油膜厚度的比值，则有

(12)

式中：hNct为热牛顿润滑中心膜厚，L为热载荷参数，L=(η0λ0ur)/Kl，λ0为润滑油的温度—黏度因子，Kl为润滑油的热传导系数，sr=us/ur，表示滑动速度与滚动速度的比值(即滑滚比)。

实际工程中润滑油的剪切稀化行为成为摩擦特性分析中不可忽略的因素，为获得非牛顿润滑的中心膜厚，BAIR[16]采用修正因子φ对牛顿润滑中心膜厚进行了修正，

(13)

式中参数Γ=(η0ur)/(hNctGm)。

1.4 粗糙表面接触模型

Kogut和Etsion基于有限元法，从物理角度出发提出一种更加准确地描述微凸体接触行为的模型(简称KE模型)，且由于其充分考虑了微凸体的弹性、混合弹塑性和完全塑性变形过程，以及接触行为表达形式的简单性而被广泛应用。因此，粗糙表面的接触特性用KE统计学接触模型来表征。这里假设整个接触区域内两表面间的距离是相等的，各微凸体间无相互作用，则粗糙表面的中心接触压力可表示为：

(14)

BEHESHTI等[28]基于KE统计模型拟合出了便于预测粗糙表面中心接触压力的方程：

(15)

同理，混合润滑状态下，载荷由润滑剂与微凸体共同承担，故在式(15)中引入微凸体承载因子γ2，获得修正后的中心接触压力方程为:

(16)

联立式(14)和式(16)，有

(17)

2 模型求解与验证

2.1 模型求解方法

本文采用迭代法联合求解非线性方程式(1)、式(13)和式(17)，获得未知参数γ1、γ2和hc，然后将γ2与hc代入式(9)中即可求得混合润滑下线接触表面的总摩擦系数，求解流程如图2所示，具体求解步骤如下：

步骤1给定载荷、速度、粗糙度以及黏度等已知参数。

步骤2合理假设微凸体承载因子γ2的初始迭代值。

步骤3将γ2的初始值代入方程1/γ1+1/γ2=1，即可求得润滑油膜的承载因子γ1。

步骤4将γ1代入式(11)，并联立式(12)和式(13)求解出被热效应因子校正后的非牛顿润滑中心膜厚hc。

步骤5将步骤4求得的hc和对应的γ2分别代入式(17)的左边与右边，并判断方程两边的绝对差值是否满足给定的误差要求。

步骤6判断是否满足误差要求。若满足，则迭代停止；否则继续执行步骤2～步骤6，直到满足收敛条件为止。

步骤7联立式(8)和式(9)求解非牛顿润滑接触表面的总摩擦系数。

2.2 模型验证

为验证本文模型，将模型仿真获得的Stribeck曲线分别与基于式(7)的牛顿润滑模型结果和文献[29]中的实验测试结果进行对比。实验中施加的法向载荷为667 N，线接触长度为25.4 mm，轴径的转速为2 rpm～500 rpm，模型其他仿真输入的参数与文献相同。

如图3所示为本文模型仿真结果分别与牛顿润滑模型结果以及实验测试结果的对比。图中横坐标Sommerfeld number=ηus/Fs，采用接触区域内的平均接触压力Pm=sqrt(πE*Fs/32RL0)代替式(4)中的P求解压力对润滑油黏度η的影响。从图中可以发现，所提模型的仿真结果与实验测试结果和牛顿模型仿真结果均具有良好的一致性，但本文的非牛顿润滑摩擦模型结果更接近实验数据，即更符合实际工况，这说明了本文模型的正确性与有效性。值得注意的是，当表面处于混合润滑状态时，由本文模型获得的摩擦系数大于牛顿润滑摩擦模型的预测结果，且随着滑动速度的增加，两者之间的差距越来越明显，这主要归因于实际工程中润滑油的剪切稀化行为，其在很大程度上减小了润滑油膜厚度，导致进入直接接触的微凸体数目增多，从而引起更高的摩擦系数。

3 结果与讨论

研究发现，Carreau流变模型可准确描述矿物油的非牛顿特性。因此，本文选用矿物油来分析非牛顿润滑工况下接触表面的摩擦润滑特性，该润滑油的基本属性参数如表1[30]所示。基于表1的仿真参数，利用本文模型获得了不同润滑工况下线接触粗糙表面摩擦系数随着滑动速度的变化规律曲线。

表1 仿真初始参数

3.1 载荷对摩擦系数的影响

不同法向载荷下表面总摩擦系数随滑动速度的变化规律如图4所示。仿真结果表明，边界润滑区域内的摩擦系数最大，当滑动速度us=1×10-2m/s时，表面总摩擦系数等于平均微凸体摩擦系数0.13，这是因为在该润滑区域内，接触表面间存在着大量微凸体与微凸体间的直接接触，摩擦力主要来源于微凸体之间的摩擦；混合润滑状态下摩擦系数随着滑动速度的增大而不断减小，这是由于随着滑动速度的增大，表面间的油膜厚度不断增大，从而减少了微凸体间的直接接触，引起摩擦系数的降低；当滑动速度进一步增大到临界值时，此时润滑状态从混合润滑过渡为流体动力润滑，在流体动力润滑状态下摩擦系数随着滑动速度的增大而增大，但其递增速率很慢。从图中还可以发现，对于相同的滑动速度，法向载荷越大，摩擦系数越大，这可解释为较大的载荷使得进入接触的微凸体数目增多，从而导致微凸体的承载量与其摩擦力增大。此外，载荷对流体动力润滑状态下的摩擦系数有明显的影响，且对于较小的载荷，润滑状态更容易从混合润滑过渡为流体动力润滑，即载荷较小时，表面达到流体动力润滑时所需的滑动速度较小。

3.2 粗糙度对摩擦系数的影响

如图5所示为表面粗糙度对总摩擦系数的影响规律。由图可知，流体动力润滑区域内的摩擦系数基本不受表面粗糙度的影响，这是由于该润滑状态下摩擦系数主要由润滑油膜的剪切行为主导，表面粗糙度的影响可以忽略不计。当滑动速度固定不变时，粗糙度较大的表面间存在更多微凸体与微凸体间的直接接触，引起粗糙微凸体摩擦系数的增大，从而导致表面总摩擦系数的增大。值得注意的是对于粗糙度较小的表面(σ=0.15×10-6m)，摩擦系数随着滑动速度的增大而急剧下降，大概在滑动速度us=1.5 m/s时，总摩擦系数达到了最小值，此时表面进入了流体动力润滑区域，即表面越光滑，形成流体动力润滑状态的速度越快。

3.3 润滑油入口黏度对摩擦系数的影响

不同润滑油入口黏度对各个润滑状态下摩擦系数的影响规律如图6所示。由图可知，在边界和混合润滑状态下，对于相同的滑动速度，摩擦系数随着润滑油入口黏度的增大而减小，这是因为润滑油膜入口黏度越大，油膜厚度越大，导致进入直接接触的微凸体数目减少，最终引起表面总摩擦系数的减小。然而，在流体动力润滑区域内，摩擦力主要来源于润滑油膜的剪切行为，因而对于较大的润滑油黏度，其剪切应力或摩擦力较大，从而导致更大的摩擦系数。从图中还可以发现，较小的润滑油入口黏度延缓了线接触粗糙表面从混合润滑过渡为流体动力润滑的进程。

3.4 不同粗糙接触模型对摩擦系数的影响

在所有润滑工况参数固定不变的条件下(Fs=300 N，η0=0.036 Pa·s，σ=0.3 μm)，不同微观粗糙接触模型下，微凸体承担的外载荷不同，从而导致油膜承载量、油膜厚度以及摩擦系数的不同。为了分析不同粗糙接触模型下的摩擦系数，分别基于仅考虑微凸体纯弹性变形的GW模型[7]和考虑微凸体多种变形机制(弹性、弹塑性以及塑性变形)的KE模型计算微凸体承载量、油膜承载量及其厚度，最终获得了两种模型下表面总摩擦系数、微凸体摩擦系数以及润滑油膜的摩擦系数分别随滑动速度的变化关系曲线，如图7所示。

图7a所示为基于GW模型和KE模型计算的表面总摩擦系数随滑动速度的变化规律。由图可知，当表面处于边界润滑的状态时，基于KE模型的摩擦系数略小于GW模型结果，两者差别很小。当表面进入混合润滑区域后，该趋势随着滑动速度的增加而变得明显，基于KE模型获得的摩擦系数明显小于仅考虑微凸体纯弹性变形的GW模型预测结果，即考虑微凸体的弹塑性和塑性变形时，表面的总摩擦系数减小。这是因为在相同的润滑工况下，与纯弹性的GW模型相比，KE模型下微凸体承担的法向外载荷较小，而润滑油膜承担的外载荷则相应增大，导致微凸体与微凸体间直接接触的数目减少，进而引起表面总摩擦系数的降低。然而，当表面进入流体动力润滑区域后，此时润滑油膜将两粗糙表面完全隔开，法向外载荷由润滑油膜的动压力来平衡，表面间无微凸体的直接接触，摩擦力主要取决于油膜的剪切行为，故流体动力润滑区域内两种粗糙接触模型下的摩擦系数相同。

图7b与图7c分别为微凸体摩擦系数与润滑油膜摩擦系数随滑动速度的变化规律。由图7b可知，两种模型下微凸体摩擦系数的差别与图7a一致，但对于润滑油膜的摩擦系数，由于KE模型考虑了微凸体的弹塑性与塑性变形，使得润滑油膜的承载量及其厚度增大，从而导致更大的油膜摩擦系数。从图7c中还可发现，润滑油膜的摩擦系数随着滑动速度的增大而增大，当表面进入流体动力润滑区域时，油膜摩擦系数不再受粗糙表面的影响，摩擦系数相等。

4 结束语

本文提出了一种非牛顿润滑的滚滑线接触表面的摩擦系数模型，模拟了不同润滑工况参数和不同粗糙微观接触模型对摩擦系数的影响规律，为机械结构表面的润滑设计、润滑状态预测以及润滑优化提供了理论指导。该模型充分考虑了润滑油的非牛顿特性和粗糙微凸体的弹塑性与塑性变形过程，使得研究成果更接近于实际工程，且本文模型易于求解，具有较强的实用性，得出如下主要结论：

(1)实际工程中润滑油膜的剪切稀化行为在很大程度上减小了油膜厚度，引起微凸体与微凸体间直接接触的数目增多，从而导致预测的表面总摩擦系数大于牛顿润滑模型的预测结果。

(2)在相同的滑动速度下，总摩擦系数随法向载荷和表面粗糙度的增大而增大，但随着润滑油入口黏度的增大而减小。载荷对流体动力润滑区域内的摩擦系数影响较大，而表面粗糙度和润滑油黏度对此状态下的摩擦系数影响较小，且该润滑状态下两粗糙表面被润滑油膜完全隔开，摩擦特性主要由润滑油膜的剪切行为所主导。

(3)在轻载、表面粗糙度较小、润滑油黏度较大的润滑工况下，接触表面更容易从混合润滑过渡为流体动力润滑状态。

(4)当表面处于边界和混合润滑状态时，基于考虑微凸体弹塑性与塑性变形的KE微观接触模型获得的表面总摩擦系数和微凸体摩擦系数均小于基于仅考虑微凸体纯弹性变形的GW模型预测结果。然而，对于润滑油膜的摩擦系数，基于KE模型的结果大于GW模型结果。

在目前的工作中，仅考虑了入口温度对油膜厚度的影响。为更加准确表征混合润滑界面的摩擦特性，未来应考虑界面相对运动引起的温度热效应，以及工作温度对混合润滑状态下线接触表面动摩擦特性的影响。