混合弹流润滑系统的建模与摩擦学特性研究

董国忠, 张辉, 李月, 董光能

(1.西安交通大学现代设计及转子轴承系统教育部重点实验室, 710049, 西安;2.株洲中车时代电气股份有限公司, 412001, 湖南株洲)

混合弹流润滑系统的建模与摩擦学特性研究

董国忠1,2, 张辉1, 李月1, 董光能1

(1.西安交通大学现代设计及转子轴承系统教育部重点实验室, 710049, 西安;2.株洲中车时代电气股份有限公司, 412001, 湖南株洲)

为了研究粗糙表面微凸体对混合弹流润滑区摩擦学特性的影响,建立了考虑粗糙表面微凸体的弹流润滑数学模型,可呈现粗糙表面的局部接触状态。通过生成虚拟粗糙表面,分别利用K-E弹塑性接触模型和平均流量雷诺方程计算微凸体接触压力与流体动压力,并且利用快速傅里叶变换技术计算基体的弹性变形量;通过绘制润滑系统的Stribeck曲线,研究了虚拟微凸体、名义载荷、综合粗糙度和微凸体曲率半径对弹流润滑摩擦学特性的影响。结果表明:微凸体的接触压力产生基体弹性变形,使膜厚增加,导致流体动压力减小,微凸体承载比和摩擦因数增大;名义载荷增加导致低速时摩擦因数变小,润滑状态在更低速条件下从边界润滑过渡到混合润滑;综合粗糙度的减小会使Stribeck曲线向左移动;微凸体曲率半径的增大仅使润滑状态加快从边界润滑过渡到混合润滑,然而对从混合区域过渡到弹流区域几乎没有影响。

摩擦学;混合弹流润滑;粗糙表面;建模

机械配副,如齿轮、滚动轴承、凸轮等,在低速、启停过程中或受到冲击载荷时往往处于混合弹流润滑状态,此时零件表面容易因摩擦磨损而失效,因此,研究压力、黏度、表面粗糙度等因素对混合弹流润滑特性的影响具有重要意义。在理论上,混合润滑同时具备边界润滑和流体动压润滑的特征[1],因此需要综合考虑微凸体接触和流体动压润滑作用,建立相应的混合润滑模型。

Greenwood等通过研究单个微凸体来分析表面粗糙度的影响,提出了GW接触模型,用来计算微凸体之间的接触压力[2]。然而,该模型仅仅适用于纯弹性变形,具有一定的局限性。Chang等提出的CEB模型[3]和Zhao等提出的微凸体接触模型[4]弥补了GW模型的缺陷。Kogut和Etsion通过有限元方法计算了半球-刚性平面接触下各弹塑性阶段本构方程的不同参数,提出了K-E弹塑性模型[5],比上述模型具有更高的精度。因此,本研究采用K-E弹塑性模型计算微凸体之间的接触压力。

另一方面,随着混合弹流润滑技术的发展,Zhu和Cheng首次利用统计学方法来考虑表面粗糙峰对点接触弹流的影响[6],在利用Greenwood-Tripp模型[2]计算表面粗糙峰弹性变形的同时,使用Patir的平均流量因子法[7]修正了雷诺方程。Venner在光滑点接触弹流模型中引入实际表面轮廓,研究了粗糙度的影响[8]。Masjedi等利用平均流量雷诺方程和Zhao提出的微凸体接触模型[4],建立了混合弹流润滑模型,并提出了描述压力和油膜厚度的方程[9]。

为了弥补平均流量模型不能呈现粗糙表面的局部接触状态的不足,并进一步分析系统混合润滑状态下的摩擦学特性,本文通过结合自回归模型、滤波技术和快速傅里叶变换法生成虚拟粗糙表面,采用平均流量雷诺方程计算了粗糙峰对流体油膜的影响,采用K-E弹塑性模型计算了粗糙峰的局部接触压力;通过快速傅里叶变换技术计算了基体的弹性变形量,利用Eyring模型计算了流体的剪切摩擦因数;最后,以该混合弹流润滑模型为基础,研究了虚拟粗糙峰、名义载荷、综合粗糙度和微凸体曲率半径对系统摩擦学特性的影响。

1 数值模型

1.1 虚拟粗糙表面的生成

从微观上看,真实表面是粗糙的,并且真实表面的微凸体分布是随机的,无法用一个确定的函数来表达。Whitehouse等的试验结果表明:许多工程表面的高度分布都满足高斯分布[10]。因此,基于微凸体的高斯分布,通过自回归模型(AR)[11]、滤波技术和快速傅里叶变换(FFT)[12]可以生成虚拟粗糙表面,如图1所示。

(a)三维图

(b)二维图(图(a)的中间对称截面)图1 虚拟粗糙表面的二维、三维图

表1列出了在不同密度的网格下生成的粗糙表面的统计参数特征值,可见仿真结果与理论值很接近,由此验证了生成的具有高斯分布的虚拟表面的有效性。

1.2 流体方程

为了求解弹流润滑中的油膜流体压力和膜厚,采用Patir的平均流量雷诺方程[7]

表1 粗糙表面统计特征值的仿真与理论值比较

(1)

式中:ph为平均油膜压力;h为各点的名义油膜厚度;hT为各点平均油膜厚度;η为润滑油黏度;φx、φy分别是沿x、y方向的压力流量因子,其计算方法可参考文献[13];u1、u2分别是2个摩擦表面的滑动速度。

黏-压方程[14]为

η=η0exp{(lnη0+9.67)[-1+

(1+5.1×10-9ph)Z]}

(2)

式中:η0为空载时的黏度;Z为黏-压指数。

密-压方程为

(3)

式中:ρ0为空载时润滑油的密度。

1.3 K-E弹塑性接触模型

采用K-E弹塑性接触模型[5]来计算微凸体的变形和接触压力。模型中将微凸体的接触分为弹性变形、弹-塑性变形和塑性变形3个阶段,并将2个粗糙表面的接触等效为一个具有半球形微凸体表面与一个平面的接触[15],如图2所示。

图2 粗糙表面与理想光滑表面的接触

对单个微凸体进行有限元模拟,结合经验公式表征微凸体在不同变形区域内的变形特性,得到微凸体的局部接触状况。接触压力和接触面积分别为

pc=Syc′(ω/ωcr)n′

(4)

ac=acrb′(ω/ωcr)m′

(5)

式中:Sy为屈服极限;ω为变形量;ωcr为纯弹性接触到弹塑性接触的临界变形量;acr为临界接触面积;b′、c′、m′、n′是由不同弹性变形阶段所决定的常量。

具体求解过程为:计算每个微凸体的变形量ω;通过ω/ωcr的取值[5],确定b′、c′、m′、n′的取值;分别通过式(4)和式(5)计算接触压力和接触面积。

1.4 膜厚方程

根据赫兹接触理论[16]和表面弹性变形可知,膜厚方程包含膜厚常量h00、基体的几何形状、表面微凸体高度和基体的弹性变形量4部分。采用快速傅里叶变换计算基体的弹性变形,以点接触承载的形式表示微凸体间的接触载荷,分别考虑线接触和点接触2种弹流润滑接触形式。每个接触点对基体产生的弹性变形量为U(xi,yi),膜厚方程和弹性变形计算式如下:

(6)

点接触时的膜厚

(7)

线接触时的弹性变形量

(8)

点接触时的弹性变形量

U(xi,yi)=

(9)

式中:δ(xi)和δ(xi,yi)分别为节点xi与(xi,yi)处附加的虚拟粗糙表面高度;Ci(xi)和Ci(xi,yi)分别为节点xi与(xi,yi)处的影响系数;P(xi)和P(xi,yi)分别为节点xi与(xi,yi)处的总压力,即接触压力与流体动压力之和,GPa;IFFT表示逆傅里叶变换。

与(3)阻(3)注(2)缕(3)树(6)鼓(2)暑(1)渚(1)俎(1)苦(4)醑(3)举(3)主(3)故(3)具(1)许(4)吐(1)睹(1)赋(2)股(1)露(3)度(3)素(2)絮(1)误(3)舞(1)谱(1)铸(1)羽(2)遽(1)据(1)渡(2)趣(3)付(1)富(2)雾(1)驻(1)暮(3)诉(1)户(1)悟(1)

1.5 载荷平衡方程及摩擦因数

通过在计算区域对流体动压力和接触压力进行积分,分别得到流体动压承载力和接触承载力,二者之和即是总承载力

(10)

由润滑剂的流变性能可知:当流体的剪切应力τ比较小,即τ<τ0时,τ随剪切率呈线性变化,表现为牛顿流体性质;当τ>τ0时,τ随剪切率呈非线性变化,表现为非牛顿流体性质。本文针对的是弹流计算,流体剪切应力较大,由Eyring模型可知

(11)

式中:τ0为Eyring剪切应力。所以,混合润滑的摩擦因数

(12)

式中:μC为边界摩擦因数;动压区面积AH=Anom-Ar,其中Anom为名义面积,Ar为实际接触面积。

1.6 计算程序

计算程序的基本步骤如下:

(1)设置初始参数,如名义载荷FT、黏度η0、粗糙度特征值(n、β、σs)等;

(2)根据赫兹接触理论,设置初始流体压力Ph;

(3)利用快速傅里叶变换法计算基体的弹性变形得到油膜厚度,计算黏压方程和密压方程,求得黏度η和密度ρ;

(4)根据K-E弹塑性模型计算微凸体的接触压力Pc,利用经验公式计算流量因子和接触因子;

(5)用平均流量雷诺方程反复迭代,直到输出的流体压力Ph收敛;

(6)根据载荷方程判断是否承载平衡;

(7)若满足载荷平衡条件,则计算摩擦因数μf并输出结果,否则重复第(3)～(6)步,直到承载平衡为止。

计算程序流程如图3所示。

图3 计算程序流程图

1.7 混合弹流润滑模型验证

试验采用UMT-2摩擦磨损试验机,选取球-盘旋转运动方式。

(1)试件材料及参数。一般情况下,摩擦配副要避免选用硬度相同的同种材料,因此试验中选用硬度不同的2种不同材料:上试件球的材料为GCr15,精度等级为IT10,处理方式包括退火、淬火、回火;下试件材料为45钢。具体参数见表2。

表2 试件的材料参数

Ra:表面粗糙度;HRC:洛氏硬度;E:弹性模量;ν:泊松比。

根据试件的材料参数,可计算塑性指数

根据G-W模型[17],当塑性指数Ψ小于1时,2个表面可基本上保持弹性接触,只会发生微小的塑性变形,所以磨损前、后表面粗糙度变化不大。

(2)工作载荷。工作载荷为4 N,可计算得最大赫兹接触应力为476 MPa,接触半径为63 μm。

(3)润滑油。由于完整的Stribeck曲线的横坐标范围比较大,所以分别使用32号低黏度机油和蓖麻油来验证模型处在边界润滑-混合润滑区域与混合润滑-弹流润滑区域的计算结果。32号机油的黏度为0.025 Pa·m,用来测试边界润滑-混合润滑区域;蓖麻油的黏度为0.6 Pa·m,用来测试混合润滑-弹流润滑区域。试验中采用蠕动泵供油,供油速率为2 mL/s。

(4)下试样盘旋转速度。选择不同的试验速度,使得系统在机油和蓖麻油润滑下分别处于边界润滑-混合润滑和混合润滑-弹流润滑状态。下试样盘转速选取:对于32号机油,当旋转半径为5 mm时,转速分别取12、30、60、120、240、和360 r/min;对于蓖麻油,当旋转半径为2 mm时,转速取12、24、36、60、120、180、360和540 r/min。

图4a是针对边界润滑-混合润滑区域,使用32号机油得到的摩擦因数随下试样盘速度的变化曲线。图中的离散点是试验测得的数据点,曲线是通过模型计算得到的仿真结果。从图中可以看出,试验测得的数据点在理论曲线附近上下波动,与Stribeck曲线在边界润滑-混合润滑区域的变化趋势具有相似性。

3组试验的平均摩擦因数与仿真结果的相对偏差见表3,从中可以看出,在边界润滑-混合润滑区域仿真结果的相对偏差比较小,均值为2.3%,这说明混合润滑模型在边界润滑-混合润滑区域比较可靠,计算结果比较准确。

(a)机油润滑试验

(b)蓖麻油润滑试验图4 摩擦因数的仿真与试验结果对比

速度/mm·s-1摩擦因数试验仿真相对偏差/%0.500.1380.137880.091.250.1330.134731.32.500.1260.129592.85.000.1170.119742.310.00.1010.102991.915.00.0960.090515.7

图4b是针对混合润滑-弹流润滑区域,使用蓖麻油得到的摩擦因数随速度变化的曲线图。试验数据点在仿真曲线(理论Stribeck曲线)附近上下波动,可以观察到Stribeck曲线在混合润滑-弹流润滑区域的变化趋势。通过对比3组试验的平均摩擦因数与仿真结果的相对偏差,发现混合润滑模型在混合润滑-弹流润滑区域的计算结果也较为准确。

2 结果与分析

2.1 粗糙表面微凸体的影响

(a)平均法计算的压力分布(线接触)

(b)基于虚拟粗糙表面计算的压力分布(线接触)

(c)油膜厚度分布(线接触)

(d)统计平均法计算的压力分布(点接触y=0平面)

(e)油膜厚度分布(点接触)图5 线接触、点接触下的微凸体压力分布和油膜厚度分布

将粗糙表面用一些统计平均特征值来表征,如图5a所示。不考虑接触压力对基体弹性变形的影响,计算得到微凸体的承载比(微凸体承载力/总承载力)为19%。图5b所示为微凸体接触压力和油膜流体动压力的分布情况,相应的油膜厚度如图5c所示。图5b中的计算结果考虑了微凸体的接触压力对基体弹性变形的影响,微凸体的承载比为22.4%,高于未考虑基体变形的计算结果。这是因为当接触压力导致基体发生弹性变形时,油膜厚度会增加,流体压力会减小,即微凸体的承载比变大了。点接触情况下的压力分布(图5d)和油膜厚度分布(图5e)与线接触情况下的类似。

图6是本文混合弹流润滑模型与平均流量G-W模型所得Stribeck曲线的对比。二者具有相同的参数:σs=50 nm,nβσs=0.05,fC=0.13,n=1.45×1011,τ0=2.5 MPa,FT=500 N,速度范围为0.001～10 m/s。本文模型的微凸体承载比例更高,所以在混合润滑区域计算出的摩擦因数偏大;本文模型在弹流润滑区域采用的是Eyring模型,考虑了流体的流变剪切应力变化,所以得到的弹流摩擦因数更小,且变化更加平缓。

图6 本文模型与平均流量G-W模型的Stribeck曲线对比

2.2 压力的影响

引入Gelinck模型[18]中的承载比例系数γ1=FT/FC和γ2=FT/FH,以加快本文模型的计算收敛速度。选择名义载荷FT为0.5、2、8和16 kN,速度范围为0.001～10 m/s。其他参数为:fC=0.13,τ0=2.5 MPa,α=2.2×10-8Pa-1,η0=0.020 2 Pa·s,R=0.02 m,E′=116 GPa,B=0.01 m,σs=50 nm,β=0.01 mm,微凸体密度n=1011m-2。生成的Stribeck曲线见图7,图中λ=h/σs为膜厚比。

图7 载荷对Stribeck曲线和膜厚比λ的影响

从图7中可以清晰地观察到润滑的3个状态:低速(0.001～0.01 m/s)时为边界润滑区域,由微凸体承载;中速(0.01～1 m/s)时为混合润滑区域,由微凸体和油膜共同承载;高速(1～10 m/s)时为弹流润滑区域,由油膜承载。随着名义载荷的增加,低速时摩擦因数减小,润滑状态更容易从边界润滑转换到混合润滑,使混合润滑区域扩大,这是因为随着载荷的增加,油膜厚度不断减小的缘故。这同Schipper做点接触弹流试验[19]时观察到的试验现象是一致的。从图7中还可以发现,在λ为3左右时开始进入弹流润滑区域,名义载荷的变化对从混合润滑区域过渡到弹流润滑区域的影响较小。

2.3 粗糙表面特征参数

图8显示了综合粗糙度σs和微凸体曲率半径β对Stribeck曲线的影响。由图8a可以看出,随着σs的减小,Stribeck曲线向左移动,从而使系统在摩擦速度较低时便可从边界润滑过渡到混合润滑,再从混合润滑过渡到弹流润滑。

(a)σs的影响(FT=2 kN,β=10 μm,nβσs=0.05)

(b)β的影响(FT=2 kN,σs=50 nm)图8 综合粗糙度和微凸体曲率半径对Stribeck曲线的影响

由图8b可以看出,随着β的增大,微凸体高度的变化趋势变得更平缓,这有利于油膜的形成,使油膜的承载比增大。低速时摩擦因数变小,润滑状态可在更低速条件下从边界润滑过渡到混合润滑,从而扩大了混合润滑区域的范围,但是β的增大对从混合润滑区域过渡到弹流润滑区域几乎没有影响。

3 结 论

本文建立了考虑粗糙表面微凸体接触的混合弹流润滑模型,通过K-E弹塑性接触模型和平均流量雷诺方程分别计算出了微凸体接触压力和流体压力。通过比较该模型的计算结果与球盘试验结果,验证了混合弹流润滑模型的有效性。绘制了弹流润滑系统的Stribeck曲线,研究了虚拟粗糙表面、名义载荷、综合粗糙度和微凸体曲率半径对弹流润滑摩擦学特性的影响,得出以下结论。

(1)微凸体的接触压力会使基体产生弹性变形、膜厚增加,导致流体压力减小,微凸体承载比和摩擦因数增大。

(2)名义载荷的增加导致低速摩擦时摩擦因数变小,润滑状态可在更低速条件下从边界润滑过渡到混合润滑。

(3)综合粗糙度σs的减小会使Stribeck曲线向左移动;β的增大可以促进润滑状态从边界润滑区域过渡到混合润滑区域,但对从混合润滑区域过渡到弹流润滑区域几乎没有影响。

[1] ZHANG H, HUA M, Dong G, et al. A mixed lubrication model for studying tribological behaviors of surface texturing [J]. Tribology International, 2015, 93(1): 583-592.

[2] GREENWOOD J A, TRIPP J H. The contact of two nominally flat rough surfaces [J]. Proceedings of the Institution of Mechanical Engineers, 1970, 185(48): 625-633.

[3] CHANG W, ETSION I, BOGY D B. An elastic-plastic model for the contact of rough surfaces [J]. ASME Journal of Tribology, 1987, 109(2): 257-263.

[4] ZHAO Y, MAIETTA D M, CHANG L. An asperity microcontact model incorporating the transition from elastic deformation to fully plastic flow [J]. ASME Journal of Tribology, 2000, 122(1): 86-93.

[5] KOGUT L, ETSION I. A finite element based elastic-plastic model for the contact of rough surfaces [J]. Tribology Transactions, 2003, 46(3): 383-390.

[6] ZHU D, CHENG H. Effect of surface roughness on the point contact EHL [J]. ASME Journal of Tribology, 1988, 110(1): 32-37.

[7] PATIR N, CHENG H. An average flow model for deter mining effects of three-dimensional roughness on partial hydrodynamic lubrication [J]. Journal of Lubrication Technology, 1978, 100: 12-17.

[8] VENNER C. Multilevel solution of the EHL line and point contact problems [D]. Enschede, Netherlands: Universiteit of Twente, 1991.

[9] MASJEDI M, KHONSARI MM. Film thickness and asperity load formulas for line-contact elastohydrodynamic lubrication with provision for surface roughness [J]. ASME Journal of Tribology, 2012, 134(1): 011503.

[10] WHITEHOUSE D, ARCHARD J. The properties of random surfaces of significance in their contact [J]. Proceedings of the Royal Society: A, 1971, 316(1524): 97-121.

[11] CANOVA F, CICCARELLI M. Panel vector autoregressive models: a survey [J]. CEPR Discussion Papers, 2013, 31: 205-246.

[12] RAO K, KIM D, HWANG J. Fast Fourier transform: algorithms and applications [M]. Berlin, Germany: Springer, 2010.

[13] PATIR N, CHENG H. Application of average flow model to lubrication between rough sliding surfaces [J]. Journal of Lubrication Technology, 1979, 101: 220-229.

[14] HUANG P. Numerical method and program for elastic deformation and viscosity-pressure equation [M]. Singapore: John Wiley & Sons Singapore Pte Ltd., 2013: 161.

[15] BEHESHTI A, KHONSARI M M. Asperity micro-contact models as applied to the deformation of rough line contact [J]. Tribology International, 2012, 52: 61-74.

[16] DINTWA E, TIJSKENS E, RAMON H. On the accuracy of the Hertz model to describe the normal contact of soft elastic spheres [J]. Granular Matter, 2008, 10(3): 209-221.

[17] 温诗铸. 摩擦学原理 [M]. 3版. 北京: 清华大学出版社, 2012: 222-223.

[18] LU X, KHONSARI M, GELINCK E. The Stribeck curve: experimental results and theoretical prediction [J]. ASME Journal of Tribology, 2006, 128(4): 789-794.

[19] SCHIPPER D. Transitions in the lubrication of concentrated contacts [D]. Enschede, Netherlands: University of Twente, 1988.

AMixedElasto-HydrodynamicLubricationModelforStudyingTribologicalPropertiesofRoughSurfaces

DONG Guozhong1,2, ZHANG Hui1, LI Yue1, DONG Guangneng1

(1. Key Laboratory of Education Ministry for Modem Design and Rotor-Bearing System, Xi’an Jiaotong University, Xi’an 710049, China; 2. Zhuzhou CRRC Times Electric Co., Ltd., Zhuzhou, Hunan 412001, China)

To investigate the influence of asperities on tribological properties in mixed elasto-hydrodynamic lubrication (EHL), a numerical model which can present the local contact state of the rough surface is developed. This model can produce a virtual rough surface in Gaussion distribution, and with the use of average flow Reynolds equation and the K-E elasto-plastic model, the hydrodynamic pressure and the contact pressure can be predicted, respectively. Fast Fourier transform (FFT) method is used to compute the elastic deformation of substrate. The influences of virtual asperities, nominal load, roughness and asperity curvature radius on the tribological properties of EHL are also discussed through plotting corresponding Stribeck curves and film thickness shapes. The results show that the contact of asperities induces the deformation of substrate, which subsequently increases the film thickness and reduces the hydrodynamic pressure. Thus, the load rate of oil film and the friction coefficient are increased accordingly. With the increase of nominal bearing load, the friction coefficient is decreased at low speed, implying that the transition point from boundary lubrication regime to mixed lubrication regime occurs at lower speeds. Smaller roughness implies the movement of the Stribeck curves towards left. The increase of asperity curvature radius accelerates the transformation from boundary lubrication to mixed lubrication. However, the transition from mixed lubrication to hydrodynamic lubrication seems not to be obviously influenced by the asperity curvature radius.

tribology; mixed elasto-hydrodynamic lubrication; rough surface; modeling

2017-05-11。 作者简介: 董国忠(1991—),男,硕士生;董光能(通信作者),男,教授,博士生导师。 基金项目: 国家自然科学基金资助项目(51475358,51705400)。

时间: 2017-10-20

网络出版地址: http:∥kns.cnki.net/kcms/detail/61.1069.T.20171020.1622.016.html

10.7652/xjtuxb201801016

TH117.2

A

0253-987X(2018)01-0107-08

(编辑 葛赵青)