高朋 侯磊 陈予恕
(哈尔滨工业大学航天学院,哈尔滨 150001)
随着转子系统向高速重载方向发展,对支承系统的动力学性能和热性能提出了更高的要求.中介轴承作为双转子系统中高、低压转子之间的重要支承和传动部件,其动力学特性[1]和热特性[2]相较于普通支承轴承更加复杂,甚至许多情况下是非线性的[3].然而,目前关于中介轴承非线性热行为的研究几乎处于空白状态.因此,对中介轴承非线性热行为机理进行深入地研究,阐明动力学参数以及热参数对其的影响是非常有必要的.
迄今为止,已经有许多学者研究了轴承的Hertz接触分数指数非线性以及径向游隙等非线性因素对转子系统动力学特性的影响.Yamamoto[4]研究了轴承径向游隙对轴承−转子系统共振特性的影响,发现共振峰的频率和振幅随径向游隙的增大而减小.Fukata 等[5]研究了球轴承在恒定径向载荷作用下的径向振动,发现轴承的Hertz 接触分数指数非线性以及径向游隙等非线性因素会引起超谐共振和亚谐共振等非线性共振.Holmes 等[3]在试验台上复现了某中型喷气发动机中双转子结构的跳跃现象、亚谐共振和组合共振等非线性动力学行为.Mevel 和Guyader[6]基于Fukata 的动力学模型研究了轻载球轴承通往混沌的两种不同路径,即次谐波路径和准周期路径.Tiwari 等[7-8]应用改进的谐波平衡法,从理论上模拟了轴承Hertz 接触和径向游隙对水平转子非线性动力学行为的影响,并进行了实验验证.Ghafari等[9]提出了考虑滚动体非线性刚度的集中质量−阻尼−弹簧模型,研究了径向游隙对轴承平衡点的影响.Bai 等[10]通过数值仿真和实验验证,建立了考虑轴承非线性因素六自由度动力学模型,研究了对称球轴承−转子系统的亚谐共振.Zhang 等[11-12]基于文献[7-8]中的分析方法HB-AFT,重点研究了轴承的Hertz 接触分数指数非线性以及径向游隙等非线性因素对转子系统的共振滞后特性的影响.然而以上研究工作中,都集中在轴承−转子系统的非线性动力学行为上,没有考虑轴承的热效应.
摩擦会阻碍物体运动并以摩擦热的形式造成能量损失.中介轴承运行中产生的摩擦热可使轴承温度升高,可采用摩擦力矩来度量.摩擦力矩的大小很大程度上取决于润滑方式,对于采用油润滑的圆柱滚子轴承,润滑油将占据轴承内部空间,从而阻碍了滚子的运动[13].摩擦力矩的大小与润滑剂的性能、填充量以及滚子运转速度都有关系,然而,这些因素通常是紧密耦合在一起的,很难区分开来[14].1945 年,Palmgren 等[15]通过对各种类型和尺寸的滚动轴承进行大量试验,得出了轴承摩擦力矩的经验公式,Palmgren 经验公式作为一种预测滚动轴承摩擦力矩的精确方法已被广泛认可和使用.根据热传递的基础理论,Harris 等[13]建立了滚动轴承稳态热传递能量守恒方程,理论预测了滚动轴承主要零部件滚珠、内圈和外圈的稳态温度.Winer 等[16]则更进一步建立了圆锥滚子轴承和轴承座间的热传递模型,同时给出了转轴、轴承以及轴承座之间的热阻计算公式.DeMul等[17-18]提出了一个五自由度轴承受力模型,同时采用矩阵法分析载荷与挠度之间的关系.Jorgensen 和Shin[19]提出了一个准三维传热模型来预测考虑热膨胀的主轴−轴承系统的稳态温度分布.基于Palmgren的经验公式,Stein 和Tu[14]提出了一个状态空间模型,用于预测由角接触球轴承的热膨胀而引起的预紧力,并分析了转速和初始预紧力的影响.Sun 等[20]开发了一种叶片损耗模拟方法,并建立了一个传热模型来估计叶片损耗以及轴承主要零部件的热膨胀数据.Takabi 和Khonsari[21]提出了油浴润滑的深沟球轴承非稳态热传递模型,从而用来研究轴承的瞬态温度.艾思源等[22]聚焦于脂润滑的双列圆锥滚子轴承的热行为,并提出了一种准静态模型,用以计算轴承载荷分布和运动参数.Than 和Huang[23]提出了一种将准静态模型与有限元相结合的方法,研究了高速主轴轴承在预载荷作用下的非线性热行为.然而,上述文献均未考虑系统动力学特性对轴承热行为的影响.
本文提出一种将传热学与非线性动力学相结合的理论方法,研究双转子系统的非线性动力学特性对中介轴承热行为的影响机制,并详细分析动力学参数和热参数的影响规律.根据双转子系统的非线性动力响应定义出中介轴承的动载荷,借助Palmgren 经验公式将其代入轴承热平衡方程,从而建立动载荷作用下的中介轴承非稳态热传递模型.通过数值模拟研究动载荷作用下的中介轴承热行为,发现轴承的运行温度不仅与转子系统的运行转速相关,也与转子系统的动力学特性息息相关,这也为轴承的设计提供了一种新的思路.
简化双盘四支承的双转子系统如图1 所示,每个转子有一个圆盘和一个转轴组成,圆盘固结在转轴上.中介轴承位于高、低压转子之间,其不同于通常用于支承的轴承,中介轴承外圈也会随着高压转子旋转.其中li(i=1,2,3,4,5) 为轴的长度,ki和ci(i=1,2,3) 表示弹性支撑的刚度和阻尼系数,ω1和ω2(rad/s) 是低压转子和高压转子转速.定义转速比为λ=ω2/ω1,由于高压转子转速大于低压转子,故λ>1 为双转子同向旋转,λ <−1 为双转子反向旋转.
图1 简化双盘四支承的双转子系统Fig.1 Simple dual-rotor system with two disks and four supports
简化双盘四支承的双转子系统动力学方程可采用Lagrange 第二方程推导,具体推导过程见文献[24-25],系统八自由度动力学方程如下
方程动力学参数物理意义及数值见文献[24-25].
对于航空发动机双转子系统,中介轴承一般采用径向圆柱滚子轴承.中介轴承具有间隙、Hertz 接触引起的分数指数以及变刚度[26]等非线性因素,会对中介轴承的弹性恢复力产生不可忽视的影响.中介轴承运动示意图如图2 所示.
第k个滚子的瞬时角位置可以表示为
其中保持架转速为ωc=(ω1ri+ω2ro)/(ri+ro),式中ri,ro分别表示中介轴承内、外圈半径.
根据几何关系,第k个滚子与滚道间的相对变形为
其中2δ0为中介轴承的径向游隙.
图2 中介轴承运动示意图Fig.2 Motion diagram of inter-shaft bearing
竖直方向和水平方向上的中介轴承恢复力为
其中H(·)为Heaviside 函数,Kb和Nb表示中介轴承的Hertz 接触刚度和滚子数目.
中介轴承采用德国FAG®公司生产的型号为NU1020 圆柱滚子轴承,其参数列于表1.
表1 NU1020 圆柱滚子轴承参数Table 1 Parameters of the radial cylindrical roller bearing NU1020
中介轴承动载荷定义为一个周期内竖直方向和水平方向上的中介轴承弹性恢复力的有效值[27],即方均根(root-mean-square,RMS),定义式如下
其中,T是周期,N表示周期内的离散化点数,和为竖直方向和水平方向弹性恢复力周期内的平均值.
中介轴承一般采用圆柱滚子轴承,总摩擦力矩M分为3 个部分:载荷摩擦力矩Ml、黏度摩擦力矩Mν和滚子端面−挡边摩擦力矩Mf.中介轴承不承受轴向载荷,因此可以忽略Mf,则总摩擦力矩M为
摩擦力矩单位统一采用N·mm.
本文引用中介轴承动载荷Fb来计算,而非以往文献采用的静载荷[23],则载荷摩擦力矩Ml为
式中,fl是一个与轴承类型和载荷相关的系数,对于带保持架的向心圆柱滚子轴承,fl=0.000 2~0.000 4,小值用于轻系列轴承,大值用于重系列轴承[13];Dm是中介轴承节圆直径.
黏度摩擦力矩Mν显然与润滑剂黏度ν 有关,采用运动黏度(cSt,即mm2/s)来表示.中介轴承内、外圈均随低、高压转子同步运转,因此采用内外圈转速差
替代原公式中的转速,单位为r/min.则黏度摩擦力矩Mν为
式中fν是一个与轴承类型和润滑方式相关的系数,带保持架的向心圆柱滚子轴承,脂润滑fν=0.6~1,油雾润滑fν=1.5~2.8,油浴润滑fν=2.2~4,小值用于轻系列轴承,大值用于重系列轴承[13].
中介轴承总摩擦热Q为
中介轴承载荷摩擦热Ql为
中介轴承黏度摩擦热Qν为
摩擦热单位统一采用W.
润滑剂黏度通常是温度的函数,由此可见,温度对于黏度影响很大.因此,在计算黏度摩擦热时,需要考虑润滑剂的黏温关系曲线.根据文献[23],润滑剂在30,40,50,60,70,80◦C 下的运动黏度分别为15,10,7.8,5.9,5,4 mm2/s.
Reynolds 黏温模型[28]如下
式中,γ 表示黏温系数,矿物油有γ=0.018~0.036◦C−1;T0表示初始温度;ν0表示初始温度下的初始运动黏度.
利用Reynolds 黏温模型拟合的实验数据[23],拟合曲线如图3 所示.拟合曲线的黏温系数为γ=0.026◦C−1,这表明Reynolds 黏温模型适用于该润滑剂.根据拟合出的黏温关系曲线,可以得到任意温度下润滑剂的运动黏度.此外,润滑剂运动黏度随温度变化范围很大,这进一步证明了考虑润滑剂黏温关系的必要性.
由于轴承钢GCr15 的Biot 数通常比较小(Bi<0.1),可以采用集总参数法[29]建立中介轴承热传递模型,忽略滚子、内圈和外圈的内部温差,从而大大简化了建模难度.中介轴承的热结点及热传递网络如图4 所示,图4(a) 给出了中介轴承传热过程中的6 个热结点,分别为中介轴承滚子、内圈、外圈、低压轴段、高压轴段和润滑剂.图4(b)给出了中介轴承热传递网络,Rri,Rro,Ri,Ro为热传导的接触热阻,RLr,RLi,RLo,RLP,RHP为对流传热的热阻.
图3 润滑剂黏温关系曲线Fig.3 Lubricant viscosity-temperature curve
图4 中介轴承的热结点及热传递网络Fig.4 Thermal nodes and heat transfer network of the inter-shaft bearing
润滑剂温度假设为TL=(Tr+Ti+To)/3,即中介轴承滚子、内圈和外圈温度的平均值[14].此外,总摩擦热假设为平均分配[14],即分配给滚子的摩擦热Qr=0.5Qt,分配给内圈的摩擦热Qi=0.25Qt,分配给外圈的摩擦热Qo=0.25Qt.对每个热结点列写热平衡方程,即可得到中介轴承热传递控制方程
式中,t为时间,csteel表示钢的比热容,mr,mi,mo,mLP,mHP分别是滚子、内圈、外圈、低压轴段和高压轴段的质量;内圈−低压轴段热阻Ri=ln(di/dL)/(2πksteelB),外圈−高压轴段热阻Ro=ln(dH/do)/(2πksteelB);滚子−内圈热阻Rri=滚子−外圈热阻Rro=式中nb为平均受力滚子数[30],Rone=1.13/(ksteel为单个滚子与内圈或外圈的热阻,其中A表示滚子与内圈或外圈接触面面积,Pe∗表示修正的Peclet 数[31].
而对流传热热阻的计算通常简化为Nusselt 数的计算,通过下式
式中,h为表面传热系数,k为流体热导率,L为特征长度,A为对流传热面积.
润滑剂−滚子的Nusselt 数[32]准则方程为
式中,Re=VL/ν 为Reynolds 数,Pr=ν/α 为Prandtl数,该准则方程适用于10−1<Re<105.
润滑剂−内圈或外圈的Nusselt 数[33]准则方程为
式中Ta=为Taylor 数,其中δio=(do−di)/2表示中介轴承内外圈的间距,r表示滚道半径.
高压轴段或低压轴段−环境Nusselt 数[34]准则方程为
由于中介轴承弹性恢复力存在不可忽视的非线性因素,导致双转子系统动力学方程也是非线性的.本文采用四阶Runge-Kutta 法求解系统动力学方程(1)~(8),得到了双转子系统的非线性动力学响应,根据中介轴承弹性恢复力表达式(9),中介轴承弹性恢复力同时也就得到了.图5 表示当转速ω1=715 rad/s时竖直方向水平方向上的弹性恢复力的时间历程曲线,双转子系统参数:转速比λ=1.2,低压转子偏心距e1=3 µm 以及高压转子偏心距e2=2 µm.竖直方向和水平方向的中介轴承弹性恢复力都是时变的,其中水平方向弹性恢复力Fy在0 附近周期性变化,而竖直方向弹性恢复力Fx在600 N 附近周期性变化,这是因为系统在竖直方向存在重力.中介轴承弹性恢复力在竖直方向和水平方向均是时变的,因此,当双转子系统以特定转速运行时,采用弹性恢复力的形式来描述中介轴承所承受的载荷显然是比较复杂的,这样并不利于中介轴承受力情况的描述和分析.
图5 中介轴承弹性恢复力Fig.5 Restoring forces of the inter-shaft bearing
图6 表示动载荷以及振幅随转速变化曲线.升速曲线表示转速由低到高的过程,而降速曲线表示转速由高到低的过程.与中介轴承弹性恢复力相比,当双转子系统以特定转速运行时,中介轴承动载荷变成了一个恒定的常数,因此,利用动载荷来描述中介轴承的受力状况显然比时变的弹性恢复力更加简单方便.
图6 动载荷以及振幅随转速变化曲线Fig.6 The dynamic load and the vibration amplitude against the rotation speed
由图6(b) 中双转子系统幅频曲线可以看出,在ω1=570~900 rad/s 转速区间内,系统存在两个共振区,分别由高压转子高低压转子的不平衡激励引起.升速过程中,振幅急剧增大,直至转速达到和振幅骤降,即发生了跳跃现象;降速过程中,当转速达到振幅骤增,即发生了跳跃现象.当转速处于ω1∈和ω1∈升速曲线和降速曲线没有重合在一起,即出现双稳态现象.从数学上,这意味着系统动力方程(1)~(8)在两个双稳态区间和内存在两个稳定解,系统动力方程收敛到哪一个稳定解取决于系统运动的初始状态.升速过程和降速过程中,系统在4 个跳跃点Adown,Aup,Bdown和Bup的运动初始状态不同,因此升速曲线和降速曲线在双稳态区间不重合.由图6(a)可以看出,中介轴承动载荷变化规律与系统振幅变化规律一致.动载荷在4 个跳跃点Adown,Aup,Bdown和Bup处也发生跳跃现象,而在两个双稳态区间和内也发生双稳态现象.总而言之,双转子系统动力学建模中,考虑了中介轴承径向游隙、Hertz 接触的分数指数关系以及变刚度等非线性因素,同时根据系统的非线性动力学影响定义了中介轴承的动载荷,因此,动载荷会发生跳跃和双稳态等非线性现象[35],这也说明中介轴承动载荷能够在一定程度上反映系统的非线性动力学特性.
中介轴承热传递控制方程(18)~(22) 本质上是一阶常微分方程组,本文采用四阶Runge-Kutta 法进行数值求解,计算精度默认为<10−6,初始温度等于环境温度T∞=20◦C,计算时间设置为10 000 s,计算时间应足够大,确保温度收敛于稳态温度.参数取值如下:转速比λ=1.2、低压转子偏心距e1=3 µm、高压转子偏心距e2=2µm.
当中介轴承处于稳态热传递状态时,其滚子、内圈以及外圈温度随低压转子转速变化曲线如图7(a)所示,图7(b)为中介轴承线性化时滚子、内圈以及外圈温度随低压转子转速变化曲线,中介轴承线性化方法及过程参见文献[25].由图7(a)可以看到,无论是升速曲线还是降速曲线,都有滚子温度Tr高于内圈温度Ti高于外圈温度To,即Tr>Ti>To.同时发现,升速曲线和降速曲线的Tr,Ti和To随转速的变化规律基本相同,因此,后文以Tr为例对中介轴承的非线性热行为进行详细分析.在共振区A和B内,升速曲线Tr急速升高,直到转速增至出现跳跃现象,即Tr骤降;降速曲线Tr缓慢降低,直到转速降至出现跳跃现象,即Tr骤增.当转速处于范围内时,升速曲线和降速曲线的Tr不重合,即出现了双稳态现象.对比图7(a)和图7(b),可以看出,当中介轴承线性化之后,其热行为的跳跃和双稳态现象均会消失不见!
图7 滚子、内圈和外圈的稳态温度随转速变化曲线(实线为升速曲线,虚线为降速曲线)Fig.7 Steady-state temperatures of rollers,inner race and outer race against rotation speed(solid lines represent run-up curves,dotted lines represent run-down curves)
为了便于后文进一步分析中介轴承的非线性热行为,现定义相关符号和参数:Adown,Aup,Bdown和Bup表示“跳跃点”;表示“跳跃点频率”;表示“跳跃幅值”;表示“双稳态区间”.
为了进一步揭示非线性热行为的内在机理,非常有必要进行摩擦热分析,包括总摩擦热分析、载荷摩擦热分析和黏度摩擦热分析.图8 为中介轴承总摩擦热、载荷摩擦热和黏度摩擦热随转速的变化曲线,其中实线表示升速曲线,虚线表示降速曲线.可以看到,升速过程中,总摩擦热Q、载荷摩擦热Ql和黏度摩擦热Qν在跳跃点Aup和Bup处均出现跳跃现象;升速过程中,三者在跳跃点Adown和Bdown处均出现跳跃现象.此外,在Aup和Bup处,Ql向下跳跃而Qν向上跳跃,但Ql跳跃幅值显著大于Qν跳跃幅值,因此Q与Ql跳跃方向相同;在Adown和Bdown处,Ql向上跳跃而Qν向下跳跃,但Ql跳跃幅值显著大于Qν跳跃幅值,因此Q与Ql跳跃方向相同.同时,三者在双稳态区间∆ωA和∆ωB内均出现双稳态现象.
图8 中介轴承摩擦热随转速变化曲线Fig.8 Friction heat of the inter-shaft bearing against rotation speed
总而言之,中介轴承温度出现非线性热行为,即跳跃和双稳态现象,直接原因是总摩擦热出现非线性行为,而根本原因是双转子系统非线性动力学特性所导致的中介轴承动载荷出现的非线性行为.之所以产生这一独特的现象是因为本文引入中介轴承动载荷代入Palmgren 经验公式来计算载荷摩擦热,如果仍旧采用静载荷则无法发现中介轴承的非线性热行为.
图9 不同转速比下滚子温度随转速变化曲线Fig.9 Variation for temperature of rollers with rotation speed under different rotation speed ratio
图9 为不同转速比下滚子温度随转速变化曲线对比,转速比分别取λ=1.1,λ=1.15,λ=1.2和λ=1.25.可见,转速比λ 不仅影响滚子温度Tr而且影响共振区A和B内中介轴承非线性热行为.随着λ 的增大,Tr显著增大;跳跃点频率基本不变,而则显著降低;跳跃幅值均明显增大;双稳态区间∆ωA和∆ωB长度也基本不变,但∆ωB略宽于∆ωA.
首先分析低压转子偏心距的影响,图10 为不同低压转子偏心距下滚子温度随转速变化曲线对比,低压转子偏心距分别取e1=2µm,e1=3µm,e1=4µm和e1=5 µm.可见,滚子温度Tr除共振区B外几乎重合,表明低压转子偏心距e1只影响共振区B.随着e1的增大,跳跃点频率基本不变,而则显著增大;跳跃幅值基本不变,而则明显增大;双稳态区间∆ωA长度基本不变,而∆ωB则变窄.
图10 不同低压转子偏心距下滚子温度随转速变化曲线Fig.10 Temperature and friction heat against rotation speed under different unbalances of lower pressure rotor
图11 不同高压转子偏心距下滚子温度随转速变化曲线Fig.11 Temperature and friction heat against rotation speed under different unbalances of higher pressure rotor
最后分析高压转子偏心距的影响,图11 为不同高压转子偏心距下滚子温度随转速变化曲线对比,高压转子偏心距分别取e2=2µm,e2=3µm,e2=4µm和e2=5 µm.可见,滚子温度Tr除共振区A外几乎重合,表明高压转子偏心距e2只影响共振区A.随着e2的增大,跳跃点频率基本不变,而则显著增大;跳跃幅值基本不变,而则明显增大;双稳态区间∆ωB长度基本不变,而∆ωA则变窄.
总而言之,低压转子偏心距只影响共振区B,而高压转子偏心距只影响共振区A,随着对应的偏心距的增大,对应的跳跃点频率和跳跃幅度均增大,而对应的双稳态区间则变窄.
图12 不同径向游隙下滚子温度随转速变化曲线Fig.12 Temperature and friction heat against rotation speed under different radial clearance
图12 为不同径向游隙下滚子温度随转速变化曲线对比,径向游隙分别取δ0=3 µm,δ0=5 µm,δ0=8 µm 和δ0=10 µm.可见,径向游隙δ 主要影响共振区A和B内中介轴承非线性热行为,随着δ 的增大,滚子温度Tr的最大值几乎不变;跳跃点频率均略微降低;跳跃幅值均略微增大;双稳态区间∆ωA和∆ωB均明显变宽.值得注意的是,当δ0=3µm 时,中介轴承跳跃现象和双稳态现象等非线性热行为“消失了”,这意味着径向游隙的存在是引起中介轴承温度出现非线性热行为的因素之一.
中介轴承径向游隙延迟了滚子与内圈和外圈的接触时间,这相当于以另一种方式降低了中介轴承的刚度,从而降低了跳跃点频率.此外,径向游隙本质上是一种非线性因素,适当减小径向间隙将显著抑制中介轴承非线性热行为,即减小跳跃幅值和缩窄双稳态区间,当径向间隙减小到一定值时,中介轴承非线性热行为甚至会消失.
图13 为不同中介轴承Hertz 接触刚度下滚子温度随转速变化曲线对比,径向游隙分别取Kb=8Kb0,Kb=10Kb0,Kb=15Kb0和Kb=20Kb0(Kb0=107N/m10/9).可见,中介轴承Hertz 接触刚度Kb主要影响共振区A和B内中介轴承非线性热行为,随着Kb的增大,滚子温度Tr的最大值几乎不变;跳跃点频率均显著增大;跳跃幅值均略微减小;双稳态区间∆ωA和∆ωB也略微变窄.
图13 不同中介轴承Hertz 接触刚度下滚子温度随转速变化曲线Fig.13 Temperature and friction heat against rotation speed under different Hertz contact stiffness of the inter-shaft bearing
图14 为不同滚子数目下滚子温度随转速变化曲线对比,径向游隙分别取Nb=12,Nb=16,Nb=20和Nb=24.可见,滚子数目Nb不仅影响滚子温度Tr而且影响共振区A和B内中介轴承非线性热行为.随着Nb的增大,Tr显著减小;跳跃点频率均显著增大;跳跃幅值均略微减小;双稳态区间∆ωA和∆ωB也略微变窄.
图14 不同滚子数目下滚子温度随转速变化曲线Fig.14 Temperature and friction heat against rotation speed under different roller numbers of the inter-shaft bearing
本文考虑中介轴承的径向游隙、分数指数非线性和参数激励等多种非线性因素,基于双转子系统非线性动力响应定义了中介轴承动载荷,同时考虑了润滑剂的黏温关系,根据Palmgren 经验公式建立了动载荷作用下的中介轴承热传递模型,研究了动载荷非线性行为对中介轴承温度的影响,并进行了详细的参数分析.主要结论如下:
(1)中介轴承的非线性动载荷将导致中介轴承温度出现非线性热行为,即表现出跳跃和双稳态现象.
(2)转速比和滚子数目对中介轴承温度和非线性热行为都有重要影响,降低转速比或增加滚子数目不仅能够有效降低温度,而且有利于抑制非线性热行为.
(3)偏心距、中介轴承径向游隙和Hertz 接触刚度均只影响中介轴承非线性热行为,降低转子偏心距、适当缩小径向游隙或增大Hertz 接触刚度均有利于抑制非线性热行为.
本文研究表明,双转子系统的非线性动力学特性会使中介轴承热行为表现出复杂非线性行为,未来的工作将集中在非线性热行为的实验验证上.