乏油工况下圆柱滚子线接触弹流润滑分析

王靖岳田 聪武旭东张 硕

(沈阳理工大学 汽车与交通学院，沈阳 110159)

齿轮、柔性梁及滚动轴承等线接触高副机构广泛应用于机械、建筑和航空航天等领域，高副接触在传递动力的过程中起着十分重要的作用，润滑油过量或乏油都是影响传递动力的主要因素。经典流体动力润滑理论适用于面接触低副机构(如滑动轴承)，通常在研究高副接触问题时，需将接触表面的润滑油黏度和弹性变形等因素加入经典流体动力润滑理论中。 由于润滑过程的复杂性，对线接触弹流润滑问题的研究多采用数值分析方法，同时辅助以实验研究。 构造计算精度高、运行稳定、收敛迅速的算法是解决弹流润滑问题的重点，近几十年来比较有代表性的算法主要有牛顿法、直接迭代法、逆解法、多重网格法等[1]。

Zhang 等[2]针对销轴和衬套之间的弹流润滑问题，采用球盘试验台进行光学干涉实验，结果表明，在乏油工况下销轴母线变化对润滑油膜厚度的影响很大。 Chippa 等[3]采用多重网格法求解圆柱滚子与平板在负载冲力下的有限长线接触弹流润滑问题，计算结果表明，润滑油膜在进口区产生颈缩且压力出现波动。 Liu 等[4]在文献[5 －6]采用的多重网格基础上，提出了一种基于压力限制算子的松弛迭代算法，简化了求解线接触弹流润滑模型中载荷平衡方程的松弛过程。 Bujurke等[7]在研究高载荷线接触弹流润滑问题时，提出了雅克比自由牛顿－克里洛夫子空间方法，并用于求解离散化的瞬态雷诺方程和膜厚方程。 Duan等[8]为研究滚动轴承在纯滚动条件下的润滑状态，提出了一种超声波法检测不同区域的润滑油膜厚度，结果表明，该方法可对弹性流体力学润滑和混合润滑区域进行有效识别。 Yang 等[9]建立了润滑剂为非牛顿流体、考虑热效应的弹流润滑线接触数学模型，并运用多重网格法求解，对乏油工况下摩擦副的润滑特性进行研究，结果表明，当润滑油膜的入口弯月面远离接触面时，弯月面上游油层的有效厚度基本稳定。 Awati 等[10]运用多重网格法求解弹流润滑模型方程，得到了变转速下椭圆接触等温稳态弹流润滑的压力和膜厚，确定了润滑油的最小膜厚、中心膜厚等参数变化对弹流润滑特性的影响。 卢洪等[11]运用多重网格法获得了不同工况下等温稳态弹流润滑问题的数值解，模拟结果表明，采用的网格层数越多，计算得到的最小膜厚与中心膜厚越接近准确值。

齿轮、轴承等工作时处于润滑状态下，且需承受较高的载荷，易产生乏油工况，目前对高副接触中乏油润滑问题的研究较少，故本文重点分析乏油工况下圆柱滚子轴承的弹流润滑问题。 首先建立乏油工况下圆柱滚子的等温线接触弹流润滑数学模型，然后运用多重网格法进行数值求解，得到压力与膜厚分布，并进一步研究供油膜厚及黏度等因素对弹流润滑特性的影响。

1 乏油工况下弹流润滑模型

1.1 弹流润滑模型的基本方程

假设弹流润滑系统中温度处处相等，在黏压方程和密压方程中略去温度项。

1)Reynolds 方程

Reynolds 方程用来描述两个表面之间润滑薄膜的流动，是润滑流体力学和弹性流体力学理论的应用基础[12]。 等温条件下一维稳态Reynolds方程为

式中:x为坐标变量；p、h分别为润滑油压力与膜厚；ρ为润滑油密度；η为润滑油黏度；uR为卷吸速度，uR＝(u1＋u2)/2，其中u1、u2分别为两个圆柱滚子表面上任意一点的运动速度；θ为润滑油膜厚与接触总间隙长度之比。 进入弹流润滑区的润滑油分别以h1、h2的油膜厚度在两个圆柱滚子表面流动，引入等效供油膜厚hoil表示弹流润滑入口的供油条件，hoil＝(h1u1＋h2u2)/uR。

式(1)的边界条件为

式中xin和xout为求解域的入口和出口边界坐标。

2)膜厚方程

弹流润滑问题中膜厚由刚性体的中心膜厚、初始间隙量及弹流变形组成。 实际油膜厚度方程为[12]

式中:h0为滚子的中心膜厚；R为曲率半径；E′为两接触表面的综合弹性模量。E′计算式为

式中:E1、E2分别为滚子1、2 表面的弹性模量；v1、v2分别为滚子1、2 的泊松比。

3)黏压方程

黏压方程采用罗兰德斯公式[12]，在等温假设下原公式经化简得到

式中:η0为润滑油的环境黏度；A1＝lnη0＋9.67；A2＝5.1 ×10－9Pa－1；z＝α/(A1A2)，其中α为黏压系数。

4)密压方程

密压方程采用Dowson-Higginson 公式，在等温假设下原公式经化简得到[13]

式中ρ0为润滑油的环境密度。

5)载荷平衡方程

线接触的载荷平衡方程为

式中w为外载荷。

1.2 弹流润滑方程的无量纲化

1)无量纲Reynolds 方程

2)无量纲膜厚方程

式中H0＝h0R/b2。

3)无量纲黏压方程

4)无量纲密压方程

5)无量纲载荷平衡方程

2 数值求解方法

本文采用多重网格法求解离散后得到的非线性代数方程组。 多重网格法的基本思想是:针对同一问题，轮流在粗网格和细网格上进行迭代，从而消除高低频偏差分量，最大程度减少数值计算工作量。 多重网格法的求解是一个递归过程。

2.1 网格划分及方程离散

运用多重网格法划分网格。 本文设计5 层网格，由最下层到顶层的网格节点数逐渐增多，在最高层即最密的网格上，进行迭代计算时采用W 循环，求解压力时迭代收敛精度(相对误差)设为0.01%。 采用有限差分法对模型方程进行离散，各方程的离散形式如下。

1)无量纲Reynolds 方程的离散形式

式中:ΔX为X方向的节点间距；εi＋1/2＝(εi＋1＋εi)/2，εi－1/2＝(εi－1＋εi)/2，i＝1，2，…，n－1，n为等距节点总数。

2)无量纲膜厚方程的离散形式

式中Ki，j为一维弹性变形刚度系数，表示在j点处单位压力作用下使i点处产生的弹性变形。Ki，j计算式为

3)无量纲载荷平衡方程的离散形式

2.2 无量纲缺陷方程

应用多重网格法求解非线性问题时，只能在最密集的网格上，即网格层数k等于网格最高层数m时，对离散后的原始方程进行直接松弛，而在低层网格中不能直接松弛和迭代，但原始方程的缺陷方程可以松弛和迭代。 无量纲方程的缺陷方程如下。

1)无量纲Reynolds 方程的缺陷方程

设第k层(k≠1)网格上无量纲Reynolds 方程的缺陷方程为

式中:Fk为右端函数向量，在k＝m时Fk为零向量；Pk为k层网格上的压力；Lk为差分算子。 在第k－1 层网格上Reynolds 方程的缺陷方程为

式中:Fki表示第k层网格上节点i的函数值；ΔkX表示第k层网格的等距步长；P~ki为第k层网格上节点i的压力近似解；Ikk－1为完全加权限制算子。

2)无量纲载荷平衡方程的缺陷方程

k层网格中无量纲载荷平衡方程的缺陷方程表达式为

式中:nk表示第k层网格的节点数；表示第k层网格上节点j的压力；gk表示第k层网格上的数值，其计算见下述说明。

在最稠密的网格上，即k＝m时，令gk＝π/2；在最稠密网格以下的网格中，即k≠m时，gk由下一层(即k－1 层)网格上的压力反复松弛迭代获得的近似解来确定，获得的近似解与k层网格上的压力计算相互独立。 当节点j的压力在第k层网格上迭代计算的近似解被限制到第k－1 层网格上时可得到

3)无量纲膜厚方程的缺陷方程

k层网格中无量纲膜厚方程的缺陷方程表达式为

式中fik为k层网格计算得到的函数，当k＝m时，即在最高层网格时计算式为

2.3 压力的松弛迭代

采用高斯－赛德尔迭代方法，通过局部线性化得到压力的迭代格式为

式中ωp为低松弛因子，一般在0.2 ~1.0 之间，本文取为0.32。δPki计算式为

其中

式中:为第k层网格上节点i的压力迭代初值；D0，0、D1，0为压缩后的变形矩阵元素。

3 计算结果与分析

以相对转动的圆柱滚子为研究对象，滚子工况参数取值如表1 所示[14]。 数值计算得到线接触弹流润滑压力和膜厚分布，分析供油膜厚、黏度和供油条件变化等对润滑特性的影响。

表1 滚子工况参数Table 1 Working parameters of the roller

3.1 线接触弹流润滑压力与膜厚分布

取xin/b＝1，此时滚子接触为乏油润滑状态。图1 和图2 分别为乏油工况下线接触弹流润滑的压力和膜厚分布，图中横坐标x/b表示啮合位置。

图1 线接触弹流润滑膜厚分布Fig.1 Film thickness distribution for linear a contact elastohydrodynamic lubrication

图2 线接触弹流润滑压力分布Fig.2 Pressure distribution for linear contact elastohydrodynamic lubrication

由图1 可以看出，当润滑油进入接触区后，膜厚开始快速增大，但随着滚子从啮入到啮出位置的变化，膜厚增大到一定程度后开始减小，且伴随有颈缩现象，离开接触区后，润滑油膜的厚度趋于平稳且小于进入接触区前的膜厚。

由图2 可以看出，接触区压力分布与Hertz 压力[12]分布接近，随着滚子相互啮合，压力呈现出先增大后减小的变化趋势。 这是因为润滑油在进入接触区后，滚子相互啮合挤压润滑液体，故压力增大，随着滚子相互分离，润滑油流动空间增大，压力降低。 在入口与出口附近压力变化较平稳。

3.2 供油膜厚对润滑特性的影响

为反映两个相对滚动的圆柱滚子间的润滑状态，研究中心膜厚及有效供油膜厚与润滑油进入位置xin之间的关系，xin/b反映了润滑油的乏油程度。 当xin/b在0 ~20 之间变化时，中心膜厚hcen和有效供油膜厚hoil的变化规律如图3 所示。

图3 hoil和hcen随xin/b 的变化规律Fig.3 Variation pattern of hoil and hcen with xin/b

由图3 可以看出:供油膜厚hoil和中心膜厚hcen均随着xin/b的增加而增大；当xin/b在0 ~6范围内变化时，hoil和hcen增大较快，近似直线分布；当xin/b在6 ~20 区间内变化时，hoil和hcen增大变缓，并逐渐趋于稳定；当xin/b为20 时，hoil和hcen均达到稳定的最大值，说明此时已经达到了润滑油充分供油状态。 当xin/b在0 ~2 范围内时，滚子接触视为乏油工况。

图4 和图5 分别为不同供油膜厚(取为0.6、1.2、2.0、3.2 μm)时弹流润滑的压力和膜厚变化。

图4 不同供油膜厚下的压力变化Fig.4 Change in pressure at different oil supply film thicknesses

图5 不同供油膜厚下的润滑膜厚变化Fig.5 Change of film thickness at different oil supply film thicknesses

由图4 可见，随着供油膜厚增加，润滑油的压力峰值也随之增大；由图5可见，随着供油膜厚增加，润滑膜厚增加且伴有颈缩现象。 供油膜厚的变化对等温稳态下的弹流润滑特性有显著影响，增大供油膜厚，乏油程度相应得到改善。

3.3 黏度对润滑特性的影响

取xin/b＝2.0，此时处于乏油状态，润滑油黏度η分别取为0.08、2.0、10.0 Pa•s，计算得到相应的压力和膜厚变化，如图6 和图7 所示。 中心膜厚随黏度的变化如图8 所示。

图6 乏油条件下黏度对压力的影响Fig.6 Impact of viscosity on pressure under spent oil conditions

图7 乏油条件下黏度对膜厚的影响Fig.7 Impact of viscosity on the film thickness under spent oil conditions

图8 中心膜厚随黏度的变化Fig.8 Central film thickness variation

由图6 可见，随着黏度变化，乏油工况下的接触压力没有明显变化，不同黏度下的压力分布曲线几乎重合；由图7 可以看出，随着黏度增加，乏油工况下的润滑膜厚减小且伴有颈缩现象；由图8 可知，随着润滑油黏度增加，中心膜厚减小。 由于本文在建立模型时未考虑温度效应，计算结果与实际情况有所差别，总体来看，在乏油工况等温条件下，线接触弹流润滑的压力受黏度影响很小，膜厚受黏度影响较大。

3.4 不同供油条件对润滑特性的影响

xin/b分别取为1.10、1.50、8.0、12.0 时润滑压力和膜厚的分布如图9 所示。 图9(a) ~9(d)分别表示了严重乏油、由乏油到正常润滑的过渡状态及正常润滑等几种情况。

图9 不同xin/b 时的压力和膜厚分布曲线Fig.9 Pressure and film thickness distribution curves for different xin/b values

由图9(a)可见:在严重乏油条件下，润滑油膜的压力分布规律与Hertz 压力分布十分接近；膜厚在接触区变化不大，变化范围在0.16 ~0.2 μm，但在出口区则明显增大。

由图9(b)可以看出:由于供油量增加，接触区的膜厚较图9(a)明显增加，膜厚的峰值也增大，且接触区膜厚出现颈缩现象；油膜压力则出现第二峰值。

由图9(c)可见:油膜压力的第二峰值较图9(b)进一步增加，且第二压力峰的位置提前。图9(d)与图9(c)中压力和膜厚分布相差不大，润滑油膜厚在压力区约为0.45 μm，可见过量供油并不能使润滑油膜厚度增加。

4 结论

建立了圆柱滚子的线接触弹流润滑数学模型，采用多重网格法对离散方程进行求解，针对等温乏油工况下的弹流润滑特性进行了分析，得到如下结论。

1)随着供油膜厚增加，压力峰值升高，润滑油膜厚增加且伴有颈缩现象。

2)润滑油黏度对压力的影响很小，对膜厚影响较大，随着润滑油黏度增大，膜厚减小。

3)严重乏油时，膜厚未见颈缩现象，压力符合Hertz 压力分布，随着乏油状态的改善，压力分布出现第二峰值，润滑油膜的厚度增大并出现颈缩现象。 随着供油量增大到一定程度，压力和膜厚分布趋于稳定，此时过量供油不会增强润滑效果，反而造成能量损耗。