高黏度油润滑弹塑性接触问题数值仿真*

张生光 刘 刚 胡文颖

(1.中国航空发动机研究院 北京 101304；2.中国黄金集团(上海)贸易有限公司 上海 200126)

随着高速机床、航空航天、铁路运输等工程技术的迅速发展，对高速、重载下的机械装备提出了更加苛刻的高精度、长寿命等要求。滚动轴承作为提供回转支撑的关键零部件，是机械系统中不可或缺的组成部分。在轴承制造和使用过程中，由于重载或振动等工况，滚动体与滚道接触时可能发生塑性变形，从而造成主轴精度下降，丧失精度保持性。而轴承需要工作在润滑条件下，其发生塑性变形时是典型的弹塑性润滑工况，因此，有必要对弹塑性润滑问题开展研究。

目前，针对弹流润滑问题，已经开展大量研究，并建立了有效的仿真计算模型，取得了丰富的成果[1-5]。在弹流润滑模型基础上，牛荣军和黄平[6]考虑材料塑性变形影响，研究了线接触粗糙表面弹塑性润滑问题；采用余弦函数表征粗糙表面，分析了粗糙度对弹塑性润滑的影响。随后，也有一些研究者对弹塑性润滑问题进行了研究，他们多采用基于各向同性强化准则下的J2流动理论求解塑性应变，应力-应变本构关系多采用线性强化准则或者Swift强化准则。REN等[7]结合统一雷诺方程和半解析法求解了点接触弹塑性润滑问题，通过对比弹塑性润滑与弹流润滑结果，发现考虑塑性变形后最大压力明显降低，并且压力范围明显增大，同时润滑油膜厚度也有所降低。REN等[8]还研究了单个微凸峰和凹坑通过接触区时对弹塑性润滑的影响，结果表明存在单个微凸峰和凹坑时，表层下应力集中较大，即使载荷较小时也较容易发生塑性变形。HE等[9-10]发展了REN等[7]的模型，研究了不同粗糙表面点接触情况以及有限长与无限长线接触情况，同时还考虑了塑性应变引起的材料硬化对弹塑性润滑的影响。在此基础上，HE等[11]研究了钢球单次和多次滚过表面后，塑性变形及其累计对接触行为的影响，给出了等效塑性应变的分布情况，为理解重载轴承滚道变形机制提供了有力的理论支撑。

上述研究成功建立了弹塑性润滑问题的数学模型并提出了合适的求解方法，这对于低黏度润滑油的求解是可行的。但如果采用高黏度润滑油，上述求解方法尚存在一定问题。由于塑性变形计算采用分步加载的方式完成，如果采用高黏度润滑油并分步加载，计算结果的接触区压力曲线会存在较高的第二压力峰，进而导致较大塑性应变的产生，与实际载荷产生的塑性应变存在较大偏差。因此，为解决高黏度润滑油工况下的弹塑性润滑仿真问题，本文作者采用半解析法建立弹塑性润滑接触模型，塑性变形采用特征应变方法进行求解，并且通过逐步增加润滑油初始黏度的方式来进行仿真计算，最终实现高黏度润滑油工况下的弹塑性润滑有效仿真分析。

1 弹塑性润滑建模与求解

润滑条件下，接触压力通过雷诺方程求解：

(1)

式中：x、y表示坐标方向；h为两表面间隙；p为接触压力；u为卷吸速度；t为时间；η和ρ分别为润滑油黏度与密度。

求得接触压力后，可进一步求得次表层应力分布。对于弹塑性接触，表层下应力由接触压力产生的施加应力和塑性应变产生的残余应力两部分组成，即：

σ=σ0+σp

(2)

式中：σp为塑性应变εp引起的残余应力，可通过特征应变的方法进行计算[12-13]；σ0为接触压力产生的施加应力：

(3)

文中研究，弹塑性模型本构方程遵循小变形假设，应力-应变本构关系采用图1所示的线性强化准则；而对于理想弹塑体，切向模量ET=0。基于各向同性强化准则下的J2流动理论，采用径向返回算法[15]求解塑性应变εp。

图1 应力-应变本构关系

对于点接触，两接触体间隙方程为

h(x,y)=h0+fg(x,y)+δ(x,y)+ve(x,y)+vp(x,y)

(4)

式中：h0为两接触体未考虑变形时的相互接近位移；fg为两接触体间未考虑变形时的几何间隙；δ为表面粗糙度；ve为表面弹性变形；vp表示塑性应变εp引起的塑性变形。

当材料为纯弹性体时，塑性变形vp=0，润滑模型由弹塑性流体润滑变为弹性流体动压润滑。

润滑油黏度按Barus黏-压关系计算：

η=η0exp(αp)

(5)

式中：η0表示环境黏度；α为黏-压系数。

密度通过Dowson-Higginson密-压关系计算：

ρ=ρ0[1+Cap/(1+Cbp)]

(6)

式中：ρ0为初始压力和初始温度下的密度；Ca和Cb为系数，Ca取6×10-10m2/N，Cb取1.7×10-9m2/N。

载荷平衡方程为

w=∬Ωp(x,y)dxdy

(7)

为计算方便，需要将各方程量纲一化处理。将计算域划分网格并用有限差分法离散量纲一化方程。数值求解时，弹性变形、塑性变形、表层下应力等采用FFT快速算法加速计算，压力求解采用Gauss-Seidal低松弛迭代逐行扫描法求解，具体差分格式及求解步骤可参见文献[16]。计算时，塑性变形和加载路径有关，应分步加载。由于塑性应变产生的塑性变形会影响接触压力，接触压力产生的应力又会影响塑性应变的求解，因此各参数相互耦合，需迭代求解直至收敛。

2 高黏度润滑油工况下计算方法分析

当选用的润滑油黏度较小时，如η0=0.03 Pa·s，载荷选取为100 N，采用半径20 mm的光滑刚性球与弹性模量206 GPa、泊松比0.3的半无限大光滑理想弹塑体接触，屈服极限取σs=600 MPa，弹塑性润滑接触压力计算结果如图2(a)所示。由于塑性变形的存在，最大接触压力明显小于最大赫兹压力pH，并且接触区域略大于赫兹接触区，计算结果比较符合实际情况。但是当润滑油黏度较大时，如η0=0.25 Pa·s，计算结果如图2(b)所示。在x>0区域压力存在波动，由于两表面均为光滑表面，不存在由于表面粗糙度引起压力波动的可能，因此，该压力波动的存在并不合理。而通过观察图2(c)所示的各分步加载下的压力分布可见，波动的位置恰好对应各分步载荷产生的二次压力峰位置。这是因为润滑油黏度较大，分步加载时，在各个分步载荷下，二次压力峰都比较大，会在其正下方产生较大塑性应变以及表面残余变形。随着分步加载载荷逐渐增大，接触面积增加，接触压力区域增大，原来二次压力峰位置的压力实际有所减小，但是由于上一步加载时该位置存在较大塑性应变，无法消除，因此会在该位置产生较大的表面残余变形，最终导致接触压力的波动。

图2 不同黏度下弹塑性润滑接触压力

为解决上述问题，文中采用了步进黏度法来计算弹塑性润滑问题。计算过程如下，分步加载时，令初始黏度η0=η0/n。根据黏度的大小，n可以取不同值，黏度较大时，n取值较大，黏度较小时n取值较小或取1。当载荷加至给定值后，在压力收敛后逐步增大初始黏度即令η0=η0×i/n，直至η0为所需黏度。求解过程中，为使方程简化，通常引入量纲一化参数，文中模型位移、压力以及油膜厚度的量纲一化参数为：X=x/a，Y=y/a，P=p/pH，H=h/a。其中a为赫兹接触半径，pH为最大赫兹接触压力。弹塑性润滑问题的具体计算流程如图3所示。采用步进黏度法计算的压力分布如图2(d)所示，接触区内压力分布光滑无波动，符合实际情况。

图3 弹塑性润滑计算流程

为验证模型的准确性，与文献[11]中的弹塑性润滑计算结果进行对比，对比结果如图4所示。通过对接触压力对比，验证了模型的准确性。

图4 模型验证

3 弹塑性润滑计算结果及分析

弹塑性润滑模型采用图2相同的参数，假设半无限大体为理想弹塑体，接触载荷取150 N，对应的赫兹接触半径与最大赫兹接触压力分别为a=216 μm，pH=1 532 MPa；计算域为{-2

当润滑油黏度较低时，取η0=0.02 Pa·s，不同卷吸速度下沿中心线的接触压力、油膜厚度以及y=0平面的Mises应力分布如图5所示。其中图5 (a)—(c)示出了不同卷吸速度下的弹流润滑计算结果，可见卷吸速度较低时，油膜厚度较小，接触压力与赫兹接触相似；随着卷吸速度的增加，润滑油膜厚度逐渐增加，弹流特征逐渐明显；由于弹流润滑出口处二次压力峰的影响，对应位置的Mises应力略有增大。图5(d)—(f)示出了不同卷吸速度下的弹塑性润滑计算结果，与弹流润滑计算结果相比，压力主峰明显减小，二次压力峰变化较小；随着卷吸速度的增加，油膜厚度也逐渐增大，其数值与弹流润滑的结果接近；从Mises应力结果可见，表层下区域已发生部分屈服，由于材料设置为理想弹塑体，因此其最大Mises应力只能达到屈服极限600 MPa。

图5 η0=0.02 Pa·s时不同速度下沿中心线压力、膜厚以及y=0平面Mises应力分布

图6给出了润滑油黏度η0=0.25 Pa·s时不同卷吸速度下沿中心线的接触压力、油膜厚度以及y=0平面的Mises应力分布。润滑油黏度较高时，PEHL计算需采用步进黏度法进行处理。图6(a)—(c)示出了不同卷吸速度下的弹流润滑计算结果，与低黏度润滑油相比，油膜厚度明显增大，接触压力的二次压力峰也更加明显，表层下Mises应力分布与低黏度润滑油相似，数值上也相差不大。图6(d)—(f)示出了不同卷吸速度下的弹塑性润滑计算结果，可见采用步进黏度法后压力、膜厚都比较平滑，油膜厚度与低黏度结果相比具有较大提升；压力主峰相比于弹流润滑计算结果明显降低，二次压力峰相比于弹流润滑计算结果，不再是随着速度的增加一直增大，而是增大到一定程度后就逐渐平稳。

图6 η0=0.25 Pa·s时不同速度下沿中心线压力、膜厚以及y=0平面Mises应力分布

表1给出了η0=0.02 Pa·s和η0=0.25 Pa·s情况下中心膜厚与最小膜厚的结果。当黏度与速度均较低时，最小膜厚为0，发生接触；随着速度升高，中心膜厚与最小膜厚均增加，并且高黏度润滑油的油膜厚度比低黏度润滑油高出较多。但是弹流润滑的中心膜厚与弹塑性润滑中心膜厚十分接近，几乎没有变化；弹流润滑的最小膜厚比弹塑性润滑的稍高。

表1 弹流润滑(EHL)与弹塑性润滑(PEHL)中心膜厚与最小膜厚对比

不同黏度、卷吸速度下的表面残余变形及等效塑性应变计算结果如图7所示。图中表面残余变形为量纲一化量，表面残余变形量S与赫兹接触半径a的比值。塑性变形与塑性应变区域主要与载荷大小相关，因此对比图7(a)、(b)、(c)、(d)可见，在宏观上塑性应变区域变化较小。随着卷吸速度的增加，表面残余变形略有增加，并且由于高速条件下二次压力峰较明显，导致其下方对应位置会产生较大的应力集中，从而更容易发生屈服。对于较高的速度和黏度工况，塑性应变区域会向表面及接触区出口处移动。

图7 不同黏度、速度下表面残余变形及等效塑性应变

4 结论

(1)研究高黏度润滑油弹塑性接触问题数值仿真方法，利用所提出的步进黏度法计算高黏度润滑油弹塑性接触模型，消除了接触区压力曲线的波动。

(2)弹塑性润滑压力主峰相比于弹流润滑结果明显降低，二次压力峰相对变化较小。弹流润滑的中心膜厚与弹塑性润滑中心膜厚接近，最小膜厚稍大。

(3)对于较高的速度和黏度工况，塑性应变区域会向表面及接触区出口处移动。