基于齿轮系统动力学的油膜刚度分析

菅光霄,王优强,*,刘晓玲,张 平,李云凯

(1.青岛理工大学 机械与汽车工程学院，山东 青岛 266520；2.华东理工大学 机械与动力工程学院，上海 200237)

0 引言

齿轮系统属于非线性振动系统，振动会引起载荷和齿轮表面速度的波动。外激励和内激励是引起齿轮系统振动的根源，外激励的主要形式为转速和输入功率的变化，内激励则主要包括误差激励和刚度激励等。刚度激励是内激励最主要的形式，是产生动载荷的主要原因[1]。对于高速重载的齿轮传动系统，惯性力和阻尼对其传动特性有显著影响，必须进行动力学分析。

目前，对于齿轮瞬态弹流的研究，大多以简化的静载荷模型来进行。黄兴保等[2]对摆动工况下直齿轮的瞬态弹流进行了研究。王文中等[3]利用数值模拟的方法，研究了渐开线斜齿轮的非稳态弹流润滑状况。时高伟等[4]研究了齿轮加速和制动过程中，变卷吸速度过程对齿轮热混合润滑的影响，得到了完全数值的热解。以上研究均是基于静载荷模型进行求解的，并未考虑齿轮传动过程中冲击或振动的作用。也有部分学者考虑了振动或冲击引起的动载荷对齿轮润滑的影响，如WANG等[5]研究了啮入冲击对齿轮瞬态弹流润滑的影响，发现了油膜压力和油膜厚度振荡的滞后效应。Cruz等[6]建立齿轮的摩擦动力学模型，在多轴齿轮传动系统的瞬态热混合润滑方面做了大量有益的工作，其研究结果表明油膜厚度及油膜压力的振荡加剧了齿轮系统的振动以及噪声的产生。XUE等[7]基于齿轮的时变啮合刚度，利用数值模拟的方法，得到了渐开线直齿轮的动态特性及动载荷分布。以上研究的不足之处均为在考虑动载荷作用的同时，忽略了冲击或振动所引起的轮齿表面速度波动的影响，且未能将齿轮系统动力学与油膜的动态特性耦合起来加以研究。部分学者如HUANG和ZHANG等[8-10]对弹流润滑的油膜刚度特性进行了研究，但未能与齿轮系统动力学相联系。邹玉静等[11]基于动力学理论与载荷分担理论，对渐开线直齿轮的动力学行为与弹流润滑特性进行了耦合研究，考虑了油膜刚度对齿轮系统振动的影响，但忽略了动力学模型所反映出的油膜动态特性，对油膜刚度模型的描述也比较单一。苑士华等[12]考虑了油膜润滑和齿轮动力学的耦合效应，研究了渐开线斜齿轮混合润滑状态下的油膜压力和油膜厚度分布以及油膜刚度与润滑油黏度之间的关系，其结果对于后续的齿轮振动与润滑的耦合研究具有重大的指导意义。

根据以上分析，在进行弹流润滑的研究时，将齿轮的弹流润滑模型与其动力学模型相结合，研究齿轮的动力学行为及齿面的润滑特性，并分析不同工况下啮合力等动态响应在啮合过程中的变化规律，以及在动载荷、速度波动作用下油膜厚度和油膜压力的变化。同时，为了更好地表征油膜的动态特性，基于弹流润滑数值解，建立润滑油膜的刚度模型，研究不同载荷和齿轮转速下油膜刚度的变化规律，用以验证动力学模型的求解结果。

1 齿轮动力学模型

1.1 动力学模型

图1为渐开线直齿圆柱齿轮的传动模型，反映了其一对齿在K点处啮合的情况。N1、N2为理论啮入、啮出点，P为齿轮节点，rbp、rbg分别为主、从动轮的基圆半径，它们的角速度分别为ωp和ωg(下标p和g分别代表主动轮和从动轮，下同)。A为啮入点，B为单齿啮入点，C为单齿啮出点，D为啮出点。

考虑齿轮时变啮合刚度的影响，建立齿轮动力学模型，该模型仅考虑直齿圆柱齿轮的扭转振动。如图2所示，其中y轴为沿啮合线方向。图中Ip、Ig；θp、θg；Tp、Tg分别表示主动轮和从动轮的转动惯量、扭转角位移和转矩。km和cm分别表示啮合刚度和啮合阻尼，et为静态传递误差。

忽略制造误差和几何误差的作用，只考虑刚度激励，故静态传递误差为零，只存在动态传递误差。

齿轮系统的动态传递误差为

δ(t)=rbpθp(t)-rbgθg(t)，

(1)

轮齿间动态啮合力为

(2)

式中，Fd为轮齿间动载荷。

由齿轮系统的动力学模型，建立其动力学微分方程为

(3)

采用Weber能量法[6]，求解齿轮的时变啮合刚度。利用四阶Runge-Kutta 法求解动力学微分方程，计算齿轮啮合过程中的动载荷和齿轮副振动速度，并将结果代入弹流润滑模型，求解齿轮的动态服役特性及其润滑工况。

1.2 渐开线直齿轮的接触参数

渐开线直齿轮的传动模型如图1所示，啮合点K处的瞬时曲率半径为

(4)

式中，Rp、Rg分别为主、从动轮在啮合点的曲率半径。

则综合曲率半径为

(5)

齿轮传动过程中存在扭转振动，这必然会引起两轮齿表面的速度波动，两轮齿表面的切向速度为

(6)

式中，Up和Ug分别为主、从动轮表面的切向速度，两者的后一项表示由于动载荷的存在，齿轮表面速度的波动值。

两当量圆柱的卷吸速度为

U(t)=[Up(t)+Ug(t)]/2，

(7)

(8)

u为齿轮副振动速度，表示振动引起的卷吸速度的波动值。

2 弹流润滑模型

2.1 基本控制方程

2.2.1 Reynolds方程

选取牛顿流体，考虑时变效应的等温线接触弹流的雷诺方程[9]为

(9)

式中，p为油膜压力(Pa)；h为油膜厚度(m)。

2.2.2 载荷方程

记w为单位长度上的载荷(N/m)，则压力应该满足载荷方程：

(10)

式中，使用动力学模型时，w为齿轮动态啮合力,w=Fd。

其他润滑控制方程如黏压方程、密压方程、膜厚方程以及各方程的无量纲化见文献[8]。

2.2 油膜刚度模型

图3为全局法计算油膜刚度的模型，该方法是通过计算赫兹接触区内的油膜压力分布与膜厚分布的改变量来获得油膜刚度的。在赫兹接触区域内的每个节点处，可以将油膜假设成若干个并联的“小弹簧”，然后将这些“小弹簧”的刚度进行叠加，从而得到整个接触区域内的油膜刚度。

利用全局法计算油膜刚度的具体计算公式为

(11)

式中，n1、n2分别表示进口区压力产生、出口区压力消失的节点，Δfi表示在节点i处，计算载荷fi与下一节点较重的载荷fi+Δfi的差值。Δfi=BepiΔx，Be为轮齿齿宽，为网格节点间距，对应油膜厚度的改变量为Δhi=hi(fi)-hi(fi+Δfi)。该方法能够捕捉到赫兹接触区内各网格节点的压力和膜厚特征，具有一定的实用性。

相对于全局法，另一种计算油膜刚度的数值方法为平均膜厚法：

(12)

3 数值方法、结果与讨论

编程进行数值求解，采用Weber能量法，求解齿轮的时变啮合刚度。利用四阶Runge-Kutta法求解齿轮的动力学方程。考虑上一瞬时计算结果对下一瞬时计算结果的影响，啮合力和齿轮副振动速度的收敛条件是：下一瞬时的计算值与上一瞬时的计算值之差的绝对值不超过0.001。将动力学求解结果代入齿轮的弹流润滑模型进行分析。基于油膜刚度模型，分析油膜的刚度特性。

进行润滑求解时，将润滑控制基本方程量纲一化后进行离散。选取牛顿流体，油膜压力的求解采用多重网格法[14]，弹性变形和油膜厚度的求解采用多重网格积分法[14]。压力的计算域为Xin=-4.6，Xout=1.4，采用W循环，利用多重网格法进行计算，划分网格共6层，每层节点数分别为31、61、121、241、481、961。每一瞬时压力计算的迭代初值使用前一瞬时压力迭代的结果，每层网格上均使用Gauss-Seidel低松弛迭代，迭代收敛判据为每个瞬时压力和载荷的相对误差小于0.001。油膜刚度的收敛判据为每一瞬时油膜刚度的相对误差不超过0.000 1。

某航空用齿轮传动和润滑的基本参数见表1。

表1 齿轮传动和润滑的基本参数 Tab.1 Basic parameters for gear and lubrication

3.1 动力学特性分析

图4为齿轮传动过程中啮合刚度随时间的变化曲线。A为啮入点，B为单齿啮入点，C为单齿啮出点，D为啮出点，在单双齿啮合的交替处(B、C两点)，啮合刚度值均存在突变，从而使齿轮系统产生换齿冲击。

图5(a)反映了齿轮传动过程中动态传递误差随时间的变化规律。在单齿啮合区和双齿啮合区，动态传递误差的平衡值不同。传动误差动态响应的平衡位置即为静态传动误差，在齿轮刚度的突变处(B、C点)，由于存在刚度激励，产生了阻尼振荡效应，并且其振荡幅值随着时间的推移越来越小。

图5(b)所示为振动所引起的卷吸速度的波动值，即齿轮副振动速度的变化。从图中可以看出，与传递误差的情况相同，在单齿啮合区，齿轮副振动速度较大，卷吸速度波动变化大。这是因为在单齿啮合区，单对齿承载，所受载荷大，且啮合刚度小。

与双齿啮合区相比，单齿啮合区内动态传递误差和齿轮副振动速度的稳态值及振幅均较大，动态传递误差和齿轮副振动速度的这种变化规律可以反映出单齿啮合区内齿轮传动的传动平稳性差。

3.2 润滑特性分析

图6为给定工况下基于动力学模型与静载模型所计算的载荷谱。考虑振动时，齿轮的承载也呈现出衰减振荡的变化规律。在单双齿啮合的交替处，齿轮载荷谱产生突变，这是由齿轮刚度的突变引起的。与双齿啮合区不同，在单齿啮合区内，由于啮合刚度小，载荷的振幅较大。

图7(a)、7(b)分别反映了最小油膜厚度和油膜压力沿啮合线的变化规律。无论是基于何种载荷模型进行求解，最小油膜厚度和中心压力的分布基本是一致的。考虑振动时，油膜厚度和油膜压力不是一条平滑的曲线，而表现为一定的“振荡效应”。这种“振荡效应”通常会引起齿轮润滑供油量的“脉动变化”，对于齿轮润滑是极为不利的。在B、C点，油膜压力存在突变，A、B点油膜压力的振幅较大，因此，啮入点以及单双齿啮合的临界点是齿轮弹流润滑的危险点。

由以上分析，基于齿轮的动力学模型进行弹流润滑的求解时，油膜压力和油膜厚度表现出一定的动态特性，当改变载荷和齿轮转速时，油膜压力、厚度的振幅和振频应当不同。这种振幅和振频的变化可以反映出油膜刚度随载荷和齿轮转速的变化。

为进一步说明动力学模型所反映出的油膜的动态特性，采用不同的齿轮载荷和转速，使用已有的油膜刚度模型求解油膜刚度，同时与动力学模型的求解结果所反映出的油膜振荡特性进行对比。

3.3 油膜刚度的讨论

3.3.1 载荷对油膜刚度的影响

给定主动轮转速np=1 600 r/min,分别取3种载荷：F=100 kN/m、120 kN/m、150 kN/m，首先求解3种载荷下齿轮的时变啮合刚度，然后基于动力学模型求解3种工况下的油膜厚度。齿轮啮合刚度和中心膜厚的求解结果分别由图8～9。

从图8中可以看出，在保证不失效的载荷作用范围内，随着载荷的增大，齿轮的啮合刚度增大。齿轮动力学的求解是以刚度激励为基础的，刚度的变化影响动力学模型的求解结果，动力学的求解结果会进一步影响油膜厚度和油膜压力等润滑特性的变化。

从图9中可以看出，随着载荷的增大，中心膜厚不断减小；且其振幅也在逐渐减小；同时，振频逐渐增高，这表明转速一定时，油膜的刚度随着载荷的增大而增大。

为了验证上述结论，基于动力学模型与油膜刚度模型，求解不同载荷下的油膜刚度。齿轮节点处油膜刚度的求解结果如图10所示。kave与kglo分别表示基于平均膜厚法与全局法计算的结果。如图10所示，转速一定时，随着载荷的增大，油膜刚度不断增大。这表明随着载荷的增大，油膜难以被压缩，这与文献[15]中杨静等的试验结果也是一致的。

3.3.2 齿轮转速对油膜刚度的影响

给定载荷F=100 kN/m，分别取主动轮转速np=1 000 r/min、2 000 r/min、2 700 r/min，基于动力学模型，求解3种工况下的油膜厚度，结果如图11所示。

从图11中可以看出，油膜厚度随转速的增加而增大，但是其振幅却明显增大，这表明随着齿轮转速的增加，油膜刚度减小。

为了验证上述结论，基于动力学模型与油膜刚度模型，进行不同转速下油膜刚度的求解。kave与kglo分别表示基于平均膜厚法与全局法计算的结果。图12(a)描述了载荷一定而转速不同时，齿轮节点处油膜刚度随转速的变化规律；图12(b)描述了载荷一定而转速不同时，齿轮啮合线上各点油膜刚度的分布情况。从图中可以看出，油膜刚度随着转速的增加而减小。

4 结论

考虑齿轮时变啮合刚度的影响，基于静载荷模型与齿轮动力学模型，分析齿轮的弹流润滑特性；建立油膜的刚度模型，对动力学模型的求解结果加以验证，得到如下结论：

1) 对齿轮系统进行弹流润滑分析时，无论是采用简化的静载荷模型还是动力学模型，得到的齿轮载荷谱、油膜厚度以及油膜压力的分布大致相同。对于低速轻载的齿轮系统，可使用简化的静载荷模型分析，但对于高速重载的齿轮系统，必须进行动力学分析。

2) 基于动力学模型进行齿轮的弹流润滑分析时，油膜厚度和油膜压力表现出一定的“振荡效应”，其振幅和振频可以反映油膜的刚度特性，载荷和齿轮转速的改变会影响其振幅和振频。

3) 齿轮传动过程中，转速一定时，随着载荷的增大，油膜的刚度增大；载荷一定时，随着齿轮转速的提高，油膜刚度减小。