高速重载弧齿锥齿轮点接触热弹流综合分析*

侯岩锟 刘振刚 胡敦珂 王秋菊 刘振侠

(1.西北工业大学动力与能源学院 陕西西安 710129；2.中国航空发动机集团航空发动机动力传输重点实验室 辽宁沈阳 110015)

弧齿锥齿轮有着承载力强、传动平稳等优点，被广泛应用于航空发动机中。在发动机工作状态转变的过程中锥齿轮的工作环境极其恶劣，拆解后时常出现的齿面点蚀胶合与齿轮的润滑状态密切相关。相关资料显示，齿轮超过50%的机械故障也是由于润滑失效引起的[1-2]。由于锥齿轮工作转速极高，润滑油供油温度较高，齿轮啮合的热效应不可忽略，因此对锥齿轮开展热弹流润滑研究是十分必要的。FRYZA等[3]通过薄膜比色干涉法测量了不同工况下点接触的油膜厚度，发现油膜厚度主要由卷吸速度与滑油黏度决定。ZHANG等[4]通过光学方法测量了球体滚动速度为0～51.15 m/s时的油膜厚度分布情况，发现速度较高时油膜中央不再平坦，而是呈现收敛楔形分布，且高速条件下的膜厚与Dowson-Higginson等温预测的最大偏差可达75%。ZHANG等[5]的实验与数值计算研究中也同样验证了这一现象。王延忠等[6]根据空间曲面共轭原理构建出齿面与刀具产形面之间的统一关系，发现等温条件下整个齿面上最小油膜厚度与最大压力出现在靠齿根方向的齿中点处。严宏志等[7]通过有限元方法计算得到了齿面压力分布，再基于Dowson-Higginson最小膜厚方程求解了弧齿锥齿轮的最小油膜厚度，发现弧齿锥齿轮的转速越高油膜厚度越大，并且齿轮啮入时的油膜厚度最小。XIAO等[8]将微弹流润滑点接触模型与黏着磨损模型应用在齿轮的运行和启停阶段，发现齿轮高转速对最小油膜厚度起着重要的作用，高转速相比载荷与粗糙度也会产生更大的温升。孙晓宇[9]发现卷吸速度沿接触椭圆短轴的分量对弧齿锥齿轮润滑膜厚起关键影响，由于螺旋锥齿轮接触齿面的椭圆率很大，限制了因卷吸速度沿接触椭圆长轴分量引起的润滑剂泄漏，所以卷吸速度沿接触椭圆长轴的分量对润滑膜厚的影响很小。韩兴等人[10]在针对椭圆点接触的研究中发现，随着速度参数的增加油膜显著增加，当该参数成倍增大时，二次峰不断向入口移动，且压力峰值高度逐渐增大，载荷参数对膜厚影响不大。MEZIANE等[11]通过实验验证了转速的增加会使最小膜厚增加，载荷的增加对最小膜厚影响不大但会使中心膜厚增加，同时外部温度的升高会使得中心膜厚与最小膜厚减小。

已有的有关弹流润滑研究大多集中在线速度30 m/s以下的齿轮以及工况变化并不剧烈的条件。由于航空发动机传动具有高速、高温、高载荷的特点，结合锥齿轮的传动特点，在其传动过程中必然会存在承载最高、工作环境最恶劣的位置，如锥齿轮分度圆位置啮入时的高速高载荷[12]。在工程实践中齿轮机构往往会受到安装结构、外部环境的限制，导致不能通过调整齿轮的部分结构参数与工况参数来改善润滑状态。但现有研究中大多定性研究不同参数(转速、载荷等)对齿轮油膜压力、膜厚、温升的影响规律，各输入参数对润滑特性的影响程度还属于未知，因此需要进行敏感性系数的计算，来获得各个输入参数对齿轮润滑效果的影响程度，在不需要改变齿轮结构的情况下，有针对性地改变某些参数来有效改善齿轮的润滑条件。本文作者以某型发动机的中央传动弧齿锥齿轮为研究对象，研究多工况多参数对齿轮啮合区域油膜压力、膜厚、温升分布的影响情况，并探讨不同参数对油膜特征参数(最大压力、最大温升、最小膜厚)的影响规律与影响程度，为航空发动机齿轮润滑设计与改进提供理论依据。

1 基本方程和数值解法

1.1 基本方程

任何润滑问题都涉及到黏性流体在狭小间隙内的流动，描述这种物理现象的基本方程为Reynolds方程，是Navier-Stokes(N-S)方程的一种特殊形式。而针对接触面的高载荷，涉及到Hertz接触理论，因此须考虑材料之间的弹性变形。高载荷又会导致润滑剂的黏度与密度发生变化，因此也需要采用相应的黏度/密度-温-压方程。而接触过程的热效应也需要与之相对应的能量方程来进行求解。在整个计算模型建立的过程中应用到的方程如下。

(1)膜厚方程

弧齿锥齿轮啮合可以简化为两椭球体接触的点接触模型，两椭球体的接触可以进一步转化为一椭球体与刚性平面的接触。如图1(a)所示，Rx1、Rx2、Ry1、Ry2分别为两接触体在x方向与y方向的曲率半径。在两方向的综合曲率半径为

(1)

如图1(b)所示，包含表面弹性变形的点接触油膜厚度方程为

(2)

式中：h0为接触区域中心油膜厚度(m)；E′为两接触固体综合弹性模量(Pa)。

图1 点接触等效模型

(2)Reynolds方程

假设润滑油为牛顿流体，遵守牛顿黏性定律，与黏性剪切力相比，油膜受到的惯性力和其他体积力可以忽略不计。由于油膜厚度很薄，可以认为油膜压力沿厚度方向保持不变，也可以忽略由表面曲率引起的速度方向的变化，即x方向的速度远大于y方向的速度。基于以上假设的广义二维Reynolds方程[13]为

(3)

式中：ρ为润滑剂密度(kg/m3)；h为油膜厚度(m)；η是油膜黏度(Pa·s)；p为油膜压(Pa)；us为x方向卷吸速度(m/s)，us=(u1+u2)/2，u1、u2分别为两固体表面的速度。

在求解Reynolds方程时，其边界条件及约束条件为

p(xin/out,y)=p(x,yin/out)=0

p(x,y)≥0 (xin

(4)

(3)黏-温-压方程

润滑油的黏-压-温关系采用Roelands关系式[14]：

(5)

式中：η0为润滑油初始黏度(Pa·s)；T0为环境热力学温度(K)；T为油膜实际热力学温度(K)；黏-压系数z=0.68；黏-温系数s=-1.1。

(4)密-温-压方程

润滑油的密度-温度-压力关系采用Dowson-Higginson关系式[15]：

(6)

式中：ρ0为润滑油初始密度(kg/m3)；温度密度系数D=-0.000 65 K-1。

(5)载荷平衡方程

数值求解弹流问题是在给定载荷条件下进行的，所以求出的压力必须满足载荷平衡条件(方程)，即压力在整个接触区域的积分需要与给定的载荷相同，载荷平衡方程为

(7)

式中:x0、x1、y0、y1分别为计算域起始与终止坐标；p(x,y)为载荷分布函数；w为点接触外加载荷(N)。

(6)油膜能量方程

在不考虑油膜的惯性力、体积力与热辐射的影响，忽略沿x和y方向的热传导，油膜厚度的尺度远小于接触区域在x方向与y方向的尺度，可认为∂p/∂z=0，故油膜的能量方程为

(8)

式中：cp是润滑油比热容(J/(kg·K))；k是润滑油导热系数(W/(m·K))。

(7)运动方程

在已知接触区域润滑油膜的压力与膜厚分布的情况下，才能求解油膜在x、y方向的速度，从而对油膜的能量方程进行求解。对运动方程进行积分即可得到对应流速：

(9)

(8)连续方程

由于油膜能量方程(8)中的流速w会给分析造成一定的困难，因此使用连续性方程对其进行替换：

(10)

(9)固体区域能量方程

求解油膜能量方程还需要已知上下固体接触面热界面条件，固体区域热界面方程为

(11)

式中：角标1、2分别代表接触区域的上下表面；ρ、c、k、u、Ti,0分别代表齿轮材料的密度(kg/m3)、比热容(J/(kg·K))、导热率(W/(m·K))、沿x方向的切向速度(m/s)、初始温度(K)。

固体区域能量方程的热边界条件为

(12)

式中：d为固体计算深度。

1.2 方程量纲一化及数值求解过程

在开始数值计算前，需要对上述理论方程进行量纲一化与离散。量纲一化公式见表1。

表1 量纲一化公式

计算区域如图2所示，X为卷吸速度方向，Y为端泄速度方向，a为接触区域短半宽，b为接触区域长半宽，其中 -2.5

图2 计算区域的划分

根据润滑油的流动情况分析，入口区可能存在一个逆流区。对于低速情况，逆流区很小，可以忽略不计。但在高转速或重载的情况下，逆流区的范围将逐渐扩大，并且逆流区逐渐向中心方向移动，此时逆流区对温度计算的影响不可忽略。因此根据流动的方向选择不同的差分方向，克服逆流给温度分析造成的困难，给出如下差分格式[16]：

(13)

图3示出了点接触热弹流润滑数值计算的基本过程。首先会采用复合直接迭代法[17]对油膜压力与膜厚分布进行求解至收敛，再采用逐行扫描对油膜温度分布进行求解，在这一过程中不断迭代更新油膜的黏度与密度，直至最后压力与温度的误差均小于1×10-6。收敛判定准则为

1×10-6

(14)

文中选择这一精度等级，是因为多次验证后发现，高于该精度一个数量级的计算结果与该精度等级的计算结果偏差小于0.1%，而低于该精度等级的计算结果则会出现大于1%的偏差。

图3 点接触热弹流程序计算流程

1.3 程序数值验证

由于文中的后续计算均基于啮合区富油的假设，因此选择文献[18-19]进行了点接触热弹流程序的准确性验证。首先采用文献[18]中图8所展示的实验工况。文献中提出只有润滑剂液滴尺寸大于100 μm时才处于完全润滑的富油状态，因此选择该液滴直径下的实验数据，将与其相同的钢球结构参数、润滑剂物性、外加载荷、卷吸速度等参数代入热弹流计算程序得到的结果见图4。如图4(a)(b)所示，啮合区域的中心油膜厚度为图中的凹陷区域，该处的油膜分布较为平坦为一平面，因此选择如图4(c)所示的x=0与y=0处的油膜厚度进行对比验证。文献[18]中实验测量的啮合区域中心油膜厚度约为1.413 μm，而文中程序计算结果为1.523 μm，计算结果相对于实验误差约为7.8%。

图4 基于文献[18]工况计算结果

由于文献[18]中采用的润滑剂黏度为5.41 Pa·s，与文中研究的润滑油黏度相差较远，因此再采用与文献[19]中图3—11(c)(d)相同的实验工况进行计算，验证低黏度下的程序准确性。计算结果见表2。可以发现在2种黏度、2种卷吸速度的情况下程序计算结果与实验数据最大误差不超过12%。由于文献[18-19]均采用光干涉法进行测量，实验数据会受测量精度的影响，同时数值计算也受到数值格式等求解方法的影响，两者会存在一定的误差。综上所述，可认为文中采用的数值方法结果可靠。

表2 2种黏度、2种卷吸速度下的数值验证

2 数值计算结果

2.1 工况参数

以航空发动机附件机匣中的中央传动齿轮作为研究对象，对其在发动机达到设计点工况时进行弹流润滑效果分析。主动锥(z=51)的设计点转速为15 000 r/min，从动锥(z=38)设计点转速为20 132 r/min，设计点功率为1 000 kW，航空发动机整个工作过程中润滑油系统的供油温度控制在120～141 ℃。研究的齿轮几何参数见表3。z为一对啮合齿轮的齿数，m为齿轮模数，α为压力角，B为齿宽，β为节圆处螺旋角。140 ℃下齿轮材料物性参数和润滑油物性参数见表4和表5。

表3 锥齿轮几何参数

表4 140 ℃下齿轮材料物性参数

表5 4106润滑油140 ℃下物性参数

油膜厚度是设计摩擦学系统时应考虑的一个很重要的参数，而最小油膜厚度则对能否实现良好润滑有着重要影响。常用油膜比厚来表征润滑状态，油膜比厚计算式为

(15)

式中：λ为油膜比厚；hmin为最小油膜厚度；σ为两接触表面的综合粗糙度；σ1、σ2为两接触表面的粗糙度。

当油膜比厚λ<1时认为接触出现干摩擦；当油膜比厚1≤λ<3时，认为该接触处于混合润滑的状态；当油膜比厚λ≥3时，则认为当前处于完全润滑的状态。

根据式(15)可以计算得到研究对象的综合粗糙度约为0.565 7 μm。

系统的输入变量在一定范围发生变化时，会使系统的输出变量产生相应的波动，通过计算敏感性系数可以评判不同输入变量对于系统输出变量的影响程度。敏感性系数越大，说明该输入变量对系统输出结果的影响越大。由于多个输入参数共同决定着齿轮的油膜特征参数，因此需要对各个输入进行量纲一化处理，文中采用的量纲一化形式敏感性系数的计算公式为

k=1,2,…,n

(16)

式中：S为敏感性系数，S为正则为正相关，为负则为负相关;a为不同的输入参数(转速、功率、黏度、弹性模量);M为对应的油膜特征参数(压力、膜厚、温升);下标k代表各参数的序号;上标*代表参考值，文中指设计点工况(20 132 r/min，1 000 kW，140 ℃)，当敏感性系数的绝对值为1时，说明输入参数每变化10%，对应的输出参数也会变化10%，因此可以将“1”作为敏感性系数的参考值。

2.2 数值计算结果

以从动锥为研究对象，在对该锥齿轮的设计点工况附近进行了数值计算后，发现油膜压力、厚度、温升的三维分布基本类似。因此以设计点工况供油温度140 ℃、功率1 000 kW、转速20 132 r/min为例，计算得到的节圆位置的压力分布、膜厚分布、平均温升分布如图5所示。图5(a)中并未出现弹流润滑中常见的2个压力峰[20]，这是由于较高的转速会使得由油膜颈缩产生的压力二次峰向啮合入口处移动并与Hertz压力峰重合，并且中央油膜平坦的现象并不显著，仅出现了轻微的油膜凹陷。在齿轮啮合过程中，油膜压力与油膜温升的趋势相似，而油膜厚度的分布则与压力温升的趋势相反。为了便于分析，从计算结果中提取油膜特征分布二维轮廓曲线(y=0)和特征参数进行最大压力、最小膜厚、最大温升的比较。

图5 设计点压力、膜厚、温升三维分布

2.2.1 齿轮转速的影响

当固定锥齿轮功率为1 000 kW时，以60%、70%、80%、90%、100%的设计点转速(以No表示)进行啮合区域润滑特性的研究，其中60%设计点转速的工况为发动机的慢车状态。转速分别为设计点转速60%、80%、100%的情况下，油膜压力、油膜厚度、油膜温升的二维轮廓曲线如图6所示。各个工况下油膜最大压力(p)、最小油膜厚度(h)、最大温升(Δt)的变化如表6所示。如图6所示，功率固定时随着齿轮转速的上升，齿轮油膜压力、油膜厚度、油膜温升均有所上升，设计点工况相对慢车工况，转速增加了66.7%，最小油膜厚度增加了30.4%，最大压力增加了33.4%，最大温升增加了202.2%。固定齿轮功率时转速的上升会导致扭矩下降从而降低接触载荷，表6中压力却随着转速的上升而增加，这是由于高转速带来的动压以及油膜颈缩现象共同产生的二次压力峰值远大于接触压力产生的第一峰值。对比慢车工况(12 709 r/min)与设计点工况(20 132 r/min)，当转速下降接触载荷上升时，可以看到压力二次峰的下降以及其向啮合出口处移动的现象，使得由接触压力产生的第一峰逐渐显现。从图6与表6中可以看到，温升与压力变化趋势相近，随着齿轮转速的上升油膜温升也逐渐升高，较高的压力与较高的转速都会带来较大的生热，因此较高的温升也出现在压力二次峰处。

图6 转速对油膜压力、膜厚和温升分布的影响

表6 油膜特征参数随转速的变化

根据表6与式(16)求得各油膜特征参数在转速变化时对应的敏感性系数见表7。可以看到，最小膜厚的敏感性系数较为稳定，最大压力的敏感性系数在小幅下降，因此转速可作为调整最小膜厚与最大压力的一个重要参数，而转速对油膜温升的敏感性系数始终大于1，因此转速的变化会对温升产生较大的影响，是改变油膜温升的一个关键参数，在转速为设计点80%时影响达到最大。

表7 油膜特征参数对转速的敏感性系数

2.2.2 齿轮功率的影响

以从动锥为研究对象，当固定锥齿轮转速为20 132 r/min时，供油温度140 ℃，以60%、70%、80%、90%、100%的设计点功率(以Po表示)进行啮合区域润滑特性的研究。功率分别为设计点转速60%、 80%、100%的情况下，油膜压力、油膜厚度、油膜温升的二维轮廓曲线如图7所示。各个工况下油膜最大压力、最小油膜厚度、最大油膜温升的变化如表8所示。如图7所示，当转速固定，随着功率的上升，齿轮油膜压力、油膜厚度、油膜温升均有所下降。油膜压力的降低是由于接触载荷的增加会抑制油膜颈缩的现象，最大温升受到最大压力下降的影响也有所下降，而最小油膜厚度并未随着最大压力的减小而增加，则是由于最大压力的下降导致润滑油黏度的降低使得润滑油的承载能力下降，同时接触载荷的增加使得弹性变形量有所增加等多种因素共同导致。

图7 功率对油膜压力、膜厚和温升分布的影响

表8 油膜特征参数随功率的变化

根据表8与式(16)求得各油膜特征参数在功率变化时对应的敏感性系数见表9。可以看到，在设计点工况最小厚度、最大压力、最大温升的敏感性系数均较小，因此可以通过功率细微调节油膜特征参数。

表9 油膜特征参数对功率的敏感性系数

2.2.3 润滑油黏度的影响

以从动锥为研究对象，当固定20 132 r/min，1 000 kW的工况时，通过改变进口温度来改变润滑油黏度进行啮合区域润滑特性的研究。80、100、120、140、160 ℃润滑油进口温度对应的润滑油黏度分别为7.108、4.595、3.191、2.345、1.789 mPa·s。润滑油进口温度分别为80、120、160 ℃的情况下，油膜压力、油膜厚度、油膜温升的二维轮廓曲线如图8所示。

图8 润滑油温度对油膜压力、膜厚和温升分布的影响

油膜最大压力、最小油膜厚度、最大油膜温升的变化如表10所示。表10中的ttotal为当前工况的初始润滑油温度与温升之和。当转速功率保持不变，随着润滑油温度的上升，油膜压力、油膜厚度、油膜温升均有所下降。润滑油温度的下降导致润滑油的黏度上升，使得油膜的承载能力增强，导致油膜压力与厚度的增加，同时较高的黏度带来了较大的摩擦生热，使得油膜最大温升也有所增加。较低润滑油入口温度带来的较高黏度使得接触压力大幅上升使得齿轮的弹性变形更加明显，压力峰值并未对应最小油膜厚度出现的位置，反而对应油膜出现凹陷的位置，这是由于最小油膜厚度出现的位置是弹性变形产生的凹陷将齿轮局部向出口处挤压以及颈缩现象所共同决定的。同时入口润滑油温度越低，润滑油温度的下降速度与压力相比出现了明显的滞后。

表10 油膜特征参数随黏度的变化

根据表10与式(16)求得各油膜特征参数在润滑油黏度变化时对应的敏感性系数见表11。可以看到，在设计点工况最小厚度与最大压力对温度的敏感性系数均大于表7中对转速的敏感性系数，因此润滑油进口温度(黏度)也可作为调整最小膜厚与最大压力的关键参数。最大压力的敏感性系数先增后减，因此在120 ℃的供油温度处改变温度对最大压力的影响最佳。最大温升的敏感性系数很大且有着很大的波动，总温度的敏感性系数增加较快，温度过高带来的润滑油黏度降低会产生较大的生热，因此依靠润滑油温度(黏度)来调整高温润滑油的最高温度来防止结焦较为困难。

表11 油膜特征参数对润滑油黏度的敏感性系数

2.2.4 齿轮弹性模量的影响

以从动锥为研究对象，设计点工况(20 132 r/min，1 000 kW，140 ℃)的齿轮弹性模量以1.0Eo(Eo=2.08 GPa)表示，现选取0.9Eo、0.95Eo、1.0Eo、1.05Eo、1.1Eo研究齿轮材料弹性模量对啮合区域润滑特性的影响。其中0.9Eo、1.0Eo、1.1Eo情况下，油膜压力、油膜厚度、油膜温升的二维轮廓曲线如图9所示。油膜最大压力、最小油膜厚度、最大油膜温升随弹性模量的变化如表12所示。可以发现，随着弹性模量的增加，最大油膜压力、最大油膜温升均升高，最小油膜厚度则下降，这是由于弹性模量的上升减小了齿轮的弹性变形量，从而使得油膜承载的压力有所上升，导致了温升的增加与膜厚的降低。

图9 弹性模量对油膜压力、膜厚和温升分布的影响

表12 设计点油膜特征参数随弹性模量的变化

根据表12与式(16)求得各油膜特征参数在弹性模量变化时对应的敏感性系数见表13。可以看到，在设计点工况最大压力与最大温升的敏感性系数均大于1且数据较为稳定，因此弹性模量可作为调整最大压力与最高温升的关键参数。而最小膜厚的敏感性系数则较小且数据波动较小，可以通过调整弹性模量来微调最小油膜厚度。

表13 油膜特征参数对弹性模量的敏感性系数

在得到各油膜特征参数对输入参数的敏感性系数后，通过取不同工况下敏感性系数的均值来进行排序，以此获得各工况参数对油膜特征参数的影响程度，如表14所示。

表14 工况参数对油膜特征参数的影响程度

可见，不同的工况参数对油膜特征参数的影响程度排序并不相同，因此当存在干摩擦时，可以通过调整对应的工况参数来获得良好的润滑。而各工况参数中润滑油黏度、转速以及材料弹性模量对油膜的各个特征参数的影响较大，功率的影响则普遍较小。其中材料弹性模量、黏度、转速每变化10%都会导致润滑油温升改变10%以上。黏度、转速每改变10%都会导致油膜厚度出现不小于5%的变化。综上，对油膜最大压力的影响程度由大到小为弹性模量、黏度、转速、功率；对油膜最大温升与最小油膜厚度的影响程度由大到小为黏度、转速、弹性模量、功率。

3 结论

(1)随着主动轮输入转速的升高，油膜压力与厚度均有所上升，而油膜的最大温升则有着显著的提升。而在较高的转速下，载荷的变化对油膜压力、膜厚、温升的影响均较小。齿轮材料弹性模量的提升会使得油膜压力、温升有所上升，使得油膜厚度降低。随着润滑油进口温度的提高，油膜压力、油膜厚度与油膜温升均下降。

(2)高转速使得弹流润滑中Hertz压力峰与二次压力峰出现了合并的现象，经典的中央平坦油膜现象并不明显，油膜分布变得尖锐，当润滑油黏度增加时，才逐渐出现了中央油膜平坦的现象。

(3)工况参数中对油膜最大压力的影响程度由大到小为弹性模量、黏度、转速、功率；对油膜最大温升与最小油膜厚度的影响程度由大到小均为黏度、转速、弹性模量、功率。

(4)当润滑油进口温度为140 ℃时，设计点工况齿轮油膜比厚为0.93，慢车工况齿轮油膜比厚为0.71，均会出现乏油干摩擦的情况。当润滑油进口温度为120 ℃时，设计点工况油膜比厚为1.1，此时齿轮处于混合润滑的状态。综上，降低润滑油的工作温度、选择较软的齿轮材料或是增加齿轮转速均可以有效改善齿轮的恶劣工作环境。