动压轴承的莫顿效应建模和参数分析

仝晓萌，蔡茂林，石岩，许未晴，胡永

(1.北京航空航天大学 自动化科学与电气工程学院，北京 100191；2.沈阳鼓风机集团股份有限公司，沈阳 110869)

1 概述

动压油膜轴承因转速高、承载能力强等优点被广泛应用于航空航天、船舶重工、医药卫生、食品加工等领域。转子高速旋转时，油膜中产生的动压力将转子悬浮于轴承中，使得转子和轴承脱离接触，而润滑油的黏度剪切会导致油膜的温度升高，并间接对轴承和转子有加热效应，通常油膜越薄，加热效应越明显。在同步振动下，转子周向某处在旋转一周过程中的平均油膜厚度低于另外一侧，导致油膜对转子加热效应不均匀，使得转子周向产生温差。该温差和转子自身的热传导同时发生，并不相悖，且转子的热传导时间常数远大于转子的旋转周期，因此同步振动时转子的周向温度梯度无法因转子的热传导而消失。非均匀加热效应会造成转子出现周向温差和转子的热弯曲，进一步加剧了转子的振动，这种现象被称为莫顿效应，如图1所示。莫顿效应是一种由多物理场耦合引起的非线性转子失稳现象，传统的转子动力学理论无法准确预测其是否存在，文献[1]将其定义为一种未被充分理解的失稳现象。莫顿效应普遍存在于油膜轴承支承的悬臂转子系统中。随着设备转速和结构复杂度的提高，莫顿效应已日趋普遍，其引起的过量振动具有在时域上幅值不断增加，相位不断变化的特点，传统的机械平衡手段无法发挥作用，极易造成转子和轴承的意外损坏[2]。建立精准的轴承和转子模型，实现对莫顿效应的准确预测，对抑制或避免莫顿效应并提高设备的安全可靠性意义重大。

图1 莫顿效应示意图

当前国际上对莫顿效应的预测方法包括线性反馈法和非线性时域预测法。文献[3]首次提出了线性反馈法，建立了轴的热弯曲和初始弯曲间的传递函数，若开环传递函数的实部大于1，则认为转子的热弯曲不断增大，最终导致系统不稳定。文献[4-5]假设转子的振幅、周向温差和热弯曲导致的不平衡量均存在线性关系，若线性系数乘积的实部小于1，则认为系统稳定。线性反馈法采用短轴承理论或简化的能量方程求解油膜温度，无法对油膜、轴承、转子、环境间的不同热边界条件进行精确建模，也无法准确预测轴承和转子的温度场和热变形，容易导致莫顿效应的漏报和错报；同时，线性反馈法假设莫顿效应仅和旋转设备的稳定转速相关，无法预测轴承-转子温度场、形变场及其导致的莫顿效应随时间的变化，如文献[5]在研究一台离心式压缩机的莫顿效应时发现，在稳定转速为11.2×103r/min时设备出现振动失稳，但若将转速从10.5×103r/min迅速升高至13.6×103r/min，设备又恢复至稳定状态，未出现莫顿效应引起的振幅过大现象。文献[6-9]建立了非线性时域预测法，通过三维有限元模型求解流体(油膜、气体)轴承的热弹流模型，并通过建立精准的温度场和热变形的边界条件，求解莫顿效应引起的转子的周向温差、热弯曲和振动随时间的变化。文献[10]通过试验测量了可倾瓦轴承的莫顿效应，并将转子周向温差与预测算法比较，验证了预测的准确性。

莫顿效应对设备的运行参数非常敏感，工业界常通过改变轴承或转子的结构以达到抑制莫顿效应的目的：文献[5]通过更换联轴器减轻了约5 kg质量；文献[11]通过更换轴承将其长径比从0.39缩减至0.33；文献[12]将轴承间隙增加了33%，并将轴承长径比从0.476降至0.250；文献[13]通过研究不同转子悬臂质量对莫顿效应的影响，发现将悬臂质量降低40%可将莫顿效应的失稳转速从8.0×103r/min提升至13.6×103r/min。然而，受制于安装空间、成本、时间等因素，复杂的机械结构变动在很多场合下无法实施。

考虑到莫顿效应本质上是转子在轴承内受热不均导致的振动失稳现象，通过调整轴承的支点位置、间隙和供油温度等来抑制莫顿效应，可避免对转子或轴承的机械结构大幅调整，具有更大的工程参考价值。文中首先对莫顿效应进行精准建模，然后借助实际案例验证模型的有效性，最后对轴承支点位置、间隙和供油温度进行参数分析，探究不同参数对莫顿效应的影响机理。

2 莫顿效应建模

2.1 油膜压力场

通过有限元算法求解雷诺方程可以得到油膜的压力分布P，考虑到油膜黏度μ随温度变化，需联合能量方程求解温度场分布。油膜压力场的边界条件如图2所示，模型中轴承两端的压力设置为0，存在部分润滑油从两端泄漏。

图2 油膜压力场边界条件

本模型采用简化的油液混合模型，因此供给流量和泄漏流量并不完全相等。雷诺方程为

∇·(D1∇P)+∇D2·V+∂h/∂t=0，

(1)

V=[ωR,0]，

式中：h为油膜厚度；V为转子表面线速度；ω,R分别为转子的转速和半径。

2.2 油膜温度场

通过有限元算法求解能量方程可以得到油膜的三维温度场T，即

(2)

式中：ρ1,c1,k1分别为润滑油的密度、比热容和热传导率；u,w分别为油膜沿转子周向和轴向的速度场分布。

然后更新油膜黏度为

μ=μrefexp[α(Tref-T)]，

(3)

式中：μref，Tref分别为油膜的参考黏度和温度；α为黏温系数。

图3 油膜温度场的边界条件

图4 油液混合示意图

(4)

2.3 轴承和转子温度场

轴承和转子暴露于空气中的表面设置为热对流边界条件。通过有限元算法求解三维热传导方程可以得到轴承和转子的瞬态温度场T，即

(5)

式中：ρ,c,k分别为轴承或转子的密度、比热容和热传导率。

2.4 轴承和转子热变形

将轴承和转子的温度场作为输入，通过有限元算法可以求解轴承和转子的三维热变形。若将转轴全部划分为三维(3D)体单元，运算量过大，不适用于转子的瞬态分析。考虑到远离轴承的转子段温升和形变不明显，因此仅在靠近轴承的转子段使用3D体单元，远离轴承的转子段使用二维(2D)梁单元以提高运算效率，并增加若干虚拟梁单元用于连接转子3D温度段和2D梁单元段以传递弯矩，转子热变形模型如图5所示。通过该步骤得到的轴承和转子热变形会影响轴承油膜的厚度，且莫顿效应造成的转子热弯曲会影响转子的动力学状态。

图5 转子热变形模型

2.5 油膜轴承动力学模型

轴承的动力学模型如图6所示，轴瓦在油膜压力的作用下可绕支点转动δ和平动yG，忽略轴瓦在压力作用下的弹性形变，但允许轴瓦在温升下的热膨胀。

图6 轴承动力学模型

轴承的运动方程为

(6)

式中：mG,IG分别为轴瓦的质量和转动惯量；kpvt为支点刚度；ypvt为支点变形；xG,yG分别为质心G在x，y方向上的平动位移；pjy,fjy分别为油膜在轴瓦表面网格节点j处沿y方向的压力、摩擦力；Mpj,Mfj分别为pjy,fjy产生的压力矩、摩擦力矩；eG为支点A到质心G的距离。

转子的动力学方程为

[F]=[Fm+Fg+Fb+Fe]，

(7)

式中：Mr,Cr,Kr分别为转子的质量、阻尼和刚度矩阵；U为有限元节点处的位移；Ubow为莫顿效应造成的热弯曲；Fm,Fg,Fb,Fe分别为机械不平衡力、陀螺力、轴承力和其他外力。

2.6 莫顿效应预测算法

莫顿效应的预测算法如图7所示。在程序运行的起点需设置转子、轴承和油膜的初始温度，转子和轴承的初始位置；然后利用有限元方法求解轴承和转子的瞬态位置，即通过雷诺方程求解油膜压力场和速度场，通过能量方程求解油膜温度场和黏度，考虑到油液的黏度和温度是相互耦合的，因此需在时域上交替求解雷诺方程和能量方程并及时更新油膜温度场和黏度；若轴承和转子的运动状态收敛，则程序继续求解轴承和转子的三维温度场分布，否则程序返回至动力学求解器；随后程序根据轴承和转子的三维温度场分布计算轴承和转子的热变形，并更新油膜的厚度和转子轴向的热弯曲，作为下个循环的动力学求解器的输入。

图7 莫顿效应预测算法

3 实例验证

文献[5]在测试图8所示的压缩机时发现，当转速接近8 000 r/min时，莫顿效应造成转子振动失稳，在数分钟内，转子的振动幅度过大导致系统自动停机保护，该转子通过双轴承支承，悬臂端的可倾瓦轴承造成了莫顿效应。以图8所示压缩机转子-轴承模型为例，验证本文建立的莫顿效应模型的有效性，为简化转子动力学模型，远离悬臂端的轴承的刚度和阻尼被固定为常数，其温度效应被忽略，系统的主要参数见表1。

图8 压缩机转子-轴承模型

表1 压缩机转子-轴承模型的主要参数

轴瓦、油膜和转子在第139 s的等温线如图9所示，可以清晰看到转子的周向存在温度梯度，高温点和低温点相距约180°相位。转子在轴承处的振动如图10所示，可以看到转子的振幅和相位随时间不断变化呈螺旋状，因此莫顿效应无法通过传统手段进行机械配平。由图9和图10可知，此时系统振幅过大，表明7 000 r/min已接近系统的失稳转速，这和文献[5]观察的现象接近。

图9 第139 s时轴承中轴截面的等温线

图10 供油温度为50 ℃时转子在轴承处的振动曲线

4 参数分析

4.1 轴承的支点偏移系数

可倾瓦轴承的支点偏移系数为0.5，即该轴瓦为对称结构，转子可沿正向或逆向旋转，如图9所示；当支点偏移系数超过0.6时，轴承的支点向轴瓦末端偏移，转子只允许单向旋转，如图6所示。为探究不同支点偏移系数对莫顿效应的影响，选取最广泛的支点偏移系数0.5和0.6进行对比分析。2种偏移系数下转子在轴承处的振幅如图11所示，可以发现当偏移系数为0.6时，转子的振幅降低了约80%；不同偏移系数下转子的均温和周向温差如图12所示，可以发现当偏移系数为0.6时，转子的均温和周向温差均降低：表明将轴承的支点偏移系数从0.5增加到0.6，不仅可以有效降低转子温度，还可成功抑制莫顿效应。分析原因：1)在多数情况下，增加支点偏移系数可提高轴承的承载能力和刚度，导致系统的临界转速升高，考虑到莫顿效应和临界转速对应的模态振型(尤其是悬臂处振幅较大的振型)紧密相关，提高相应的临界转速，也提高了莫顿效应的失稳转速[8]，使该失稳转速远离当前运行转速；2)增加支点偏移系数增加了轴承承载区油膜的平均厚度，有利于减少油膜对轴承的加热效应，因此转子的均温和周向温差均下降。

图12 不同支点偏移系数下的转子周向温差和均温

4.2 轴承间隙

为研究轴承间隙对莫顿效应的影响，现将表1中的轴承间隙增加20%，其余参数保持不变。考虑到莫顿效应是由转子周向温差导致的振动失稳现象，且转子温差造成的热弯曲和转子振动在时域上相互耦合，造成转子振动相位不断变化，振幅随时间波动，相应的油膜厚度也随时间改变。不同轴承间隙下转子在轴承处的振幅和最小油膜厚度如图13所示，增加轴承间隙20%可有效降低转子振幅约85%，且最小油膜厚度从不足0.2Cb增加到0.4Cb(Cb=74.9 μm，为轴承设计间隙，见表1)。不同轴承间隙下转子周向温差和均温如图14所示，增加轴承间隙20%既可降低转子的均温又可以减少转子的周向温差。由上述结果可知，增加轴承间隙，可降低润滑油对转子的加热效应，有利于抑制莫顿效应，但过度增加轴承间隙，将显著改变轴承刚度和阻尼，并影响轴承-转子系统的临界转速、振动等特性，需谨慎实施。

图13 不同轴承间隙下的转子振幅和油膜厚度

图14 不同轴承间隙下的转子周向温差和均温

4.3 供油温度

为探究轴承的供油温度对莫顿效应的影响，现将表1中的供油温度从50 ℃降低至30 ℃，转子在轴承处的振动曲线如图15所示，莫顿效应造成转子的振动呈螺旋式增长，表现为振幅不断增加，相位持续变化；对比图10可知，转子的最大振幅随供油温度降低而增大。不同供油温度下的转子周向温差和均温如图16所示，不难发现，降低供油温度时，转子的平均温度降低，周向温差反而增大。不同供油温度下的最小油膜厚度和平均油膜厚度如图17所示，将供油温度从50 ℃降至30 ℃，油膜的平均厚度增加了约0.1Cb，而最小油膜厚度未出现明显变化。图17中平均油膜厚度在仿真初始阶段迅速降低，这是由于转子和轴承的初始温度被设置在40 ℃，高于常温时，转子和轴承的热膨胀将减小轴承的实际工作间隙。

图15 供油温度30 ℃时转子在轴承处的振动曲线

图16 不同供油温度下的转子温差和均温

图17 不同供油温度下的平均和最小油膜厚度

以上结果表明，降低供油温度可以降低轴承和转子的运行温度，提高轴承的实际工作间隙；但由于降低油温会显著增加油液黏度(在该例中约86%)，反而加剧油膜对轴承的非对称加热效应，导致转子的周向温差不降反升，转子热弯曲也更强烈：因此降低轴承的供油温度在该实例中无法有效抑制莫顿效应。

5 结论

本文对可倾瓦动压油膜轴承进行了精准的热弹流建模，通过编写三维有限元软件求解转子-轴承系统的动力学、温度和热变形，并在时域上预测莫顿效应，通过分析轴承的支点偏移系数、轴承间隙和供油温度对莫顿效应的影响，得到如下结论：

1)建立的基于三维有限元模型的莫顿效应预测算法可准确预测转子热弯曲造成的振动失稳。

2)将可倾瓦轴承的支点偏移系数从0.5增加至0.6，可有效降低转子的平均温度和周向温差，进而抑制莫顿效应和振动失稳现象。

3)适当增加轴承间隙，可降低油膜对轴承的加热效应，有利于抑制莫顿效应，考虑到该间隙和轴承的刚度阻尼特性息息相关，需谨慎调整轴承间隙，以防止对系统的临界转速和振动特性产生不利影响。

4)降低轴承的供油温度，导致润滑油的黏度升高，加剧了油膜对转子的非对称性加热效应，使得转子的周向温差升高，加剧了莫顿效应。