多叶动压气体滑动轴承静态特性的有限差分算法*

2023-11-08张功学

润滑与密封 2023年10期

陈阳张功学吴垚

(陕西科技大学机电工程学院陕西西安 710021)

随着旋转机械向高转速、高精度、高功率和低能耗的趋势发展，滑动轴承对高速转子的润滑性能已成为影响设备高效稳定工作的关键因素。相较于传统滚动轴承和油液润滑轴承，动压气体轴承因其结构简单、易于制造、对环境无污染及服役温度范围宽等优点[1-2]，广泛应用于空气压缩机、高速离心机和透平膨胀机等旋转机械中。多叶动压气体轴承采用多块相同几何参数的圆弧状轴瓦来支承转子，其轴向槽间隙可有效防止压力扰动沿周向扩散以提高轴承的稳定性，同时轴向槽的存在还起到收集灰尘和杂质粒子、维持瓦块/轴颈表面清洁的作用[3-4]。

目前，国内外学者针对动压气体轴承的性能开展了大量深入的研究，并取得了丰富的研究成果。赵三星等[5]发现轴向槽圆柱气体动压轴承的承载力和偏位角会随着开槽位置的改变呈周期性变化。LEE和KIM[6]通过耦合求解波箔型径向箔片轴承的气膜压力控制方程和能量方程，研究了轴承间隙、载荷大小和轴颈转速对转轴三维温度场分布的影响，结果表明转子的最高温度随轴承间隙的减小而升高，随载荷和转速的增大而急剧升高。田宇忠等[7]发现轴承的承载能力会随着轴向槽宽度与数目的增大而逐渐降低。RAHMATABADI等[8]通过有限元法求解微极性流体润滑Reynolds方程，发现承载力、偏位角和摩擦力随预负荷系数和微极性流体耦合数的增大而增大，端泄流量随预负荷系数的增大而减小。燕震雷和伍林[9]结合牛顿迭代法和有限差分法求解考虑气体稀薄效应的修正Reynolds方程，探讨了不同边界滑移条件对三瓦可倾瓦动压气体轴承承载力的影响，表明气体稀薄效应明显降低了轴承承载能力，且随着克努森数的增大，连续模型、一阶滑移模型和WU的新滑移模型计算的承载力依次降低。冯凯等人[10]建立了考虑接触面库仑摩擦的三瓣式气体箔片轴承气弹耦合模型，研究了转速和轴承预载对承载力的影响，发现承载力随轴承预载和转速的增大而显著提高。DAS 和ROY[11]通过数值求解非牛顿流体润滑的量纲一化Reynolds方程，对比了牛顿/非牛顿流体润滑剂润滑下二、三、四叶滑动轴承静态特性曲线，发现当非牛顿润滑剂的膨胀系数为1.3时三叶轴承拥有最大承载力，而牛顿流体润滑下的二叶轴承具有最大的姿态角和流量系数。LV等[12]建立了考虑轴颈倾斜、轴瓦弹性变形和瓦块数目的轴向槽轴承弹流润滑模型，运用等效支点位置探究了轴颈倾斜角和轴向槽数目对承载能力的影响，发现承载能力随偏斜角和沟槽数的增加而减小。LI等[13]结合线性摄动法和有限差分法求解静压气体轴承的非定常Reynolds方程，得到了不同供气孔径且不均匀分布的气膜厚度和气膜压力分布，指出承载力随孔口直径的减小而增大，随偏心率、转速和供气压力增大而增大，并发现在小偏心率时轴承具有更好的稳定性。化文灿等[14]基于有限差分法对2种波箔型动压气体轴承的定常Reynolds方程和气膜厚度方程联合求解，得到了承载力和摩擦力矩随转速与偏心率的变化规律。

由于轴颈/轴瓦间隙、气体可压缩性、气膜不连续性及瓦块分布位置联合作用的复杂性，针对多叶动压气体轴承润滑性能的研究文献相对较少。而研究多叶动压气体轴承的静态特性，对提高气体轴承-转子系统稳定性具有重要的指导意义[15-17]。

本文作者通过数学变换将多叶动压轴承的气体润滑Reynolds方程转化为标准偏微分方程形式，利用有限差分法和超松弛迭代法进行数值求解，得到三叶气体轴承各瓦块的气膜膜厚及压力分布、承载力和摩擦因数等静态性能，并分析了偏心率、预负荷系数、轴承数、长径比和瓦块分布位置对多叶气体轴承静态性能的影响规律。

1 数学模型

图1所示为不同承载方式的三叶动压气体滑动轴承的结构示意图。三叶动压气体滑动轴承是由三块相同几何尺寸的轴瓦依次均匀地分布在支点圆上构成的，瓦块分布位置的改变可以影响主要承载区域，主要承载区域集中在瓦2上时为瓦上承载方式，位于瓦1与瓦2之间时为瓦间承载方式。图中ξ为瓦间角，α为瓦包角。

图1 三叶动压气体滑动轴承结构

为了计算轴承气膜的承载力等静态性能，必须计算其气膜压力分布。根据气体运动方程、连续性方程、气体状态方程和 Navier-Stokes方程，可得到可压缩气体润滑Reynolds方程[18]的一般形式为

(1)

式中：ρ为气体密度；h为气膜厚度；μ为气体动力黏度；U为轴颈表面速度，U=ωR；ω为轴颈转速；t为时间；x为轴承周向坐标；z为轴承轴向坐标。

对多叶动压气体轴承而言，其气膜厚度与轴颈中心相对于轴承中心的偏心距有关，在忽略轴瓦厚度的前提下，第i块轴瓦与轴颈表面的气膜厚度表达式为

hi=Cb-(Cb-Cp)cos(βi-φi)+ecos(φi-θ)

(2)

式中：Cb为轴承半径间隙；Cp为装配间隙；βi为第i块轴瓦的支点位置角；e为轴颈中心相对轴承中心的偏心距；θ为轴承偏位角。

对动压轴承进行分析计算时，常以量纲一化形式进行，可以将问题归纳为最紧凑形式，突出各有关因素的作用。则等温条件下气体润滑Reynolds方程的量纲一化形式[19]如下：

(3)

量纲一气膜厚度表达式为

H=1-mpcos(β-φ)+εcos(φ-θ)

(4)

式中：H为量纲一气膜厚度；ε表示轴颈中心相对支点圆中心的偏心率，ε=e/Cb；mp为轴瓦的预负荷系数，mp=1-Cp/Cb。

量纲一压力边界条件[20]为

(5)

式中：L表示轴承宽度；φ0为从角起线到轴瓦进气端的角度；φ1为从角起线到轴瓦出气端的角度。

2 数值求解

为了快速求解动压气体润滑Reynolds方程，需将其化为标准偏微分方程的形式进行求解。通过简单的数学变换[20]，令S=PH，则Π=S2=(PH)2，式(3)可变为

(6)

同理可得：

(7)

则：

(8)

通过移项，多叶动压气体轴承的静态润滑Reynolds方程(3)可整理为

(9)

利用有限差分法求解动压气体润滑轴承压力分布，就是在瓦块的气膜区域按照差分的步长Δφ、Δλ划分为许多网格。如图2所示，共划分有 (m+1)×(n+1)个网格点，每个节点的位置用(i，j)编号表示，节点(i，j)上的P值用Pi，j表示，得到一组近似表达气膜压力分布的离散压力数值，再利用相应的数值积分便可求得承载力等静态性能。

图2 多叶气体轴承瓦块的网格划分

采用二阶差分法和中心差分法对式(9)进行离散：

≈[(H/2)i+1/2，j·Πi+1，j-((H/2)i+1/2，j+

(H/2)i-1/2，j)·Πi，j+(H/2)i-1/2，j·Πi-1，j]/(Δφ)2

(10)

≈[Πi+1/2，j·Hi+1，j-(Πi+1/2，j+Πi-1/2，j)·Hi，j+

Πi-1/2，j·Hi-1，j]/(Δφ)2

(11)

同理可得：

((H/2)i，j+1/2+(H/2)i，j-1/2)·Hi，j]/(Δλ)2+

[(H/2)i，j-1/2·Hi，j-1]/(Δλ)2

(12)

(13)

(14)

则量纲一化Reynolds方程可化为

[(H/2)i+1/2，j·Πi+1，j-((H/2)i+1/2，j+

(H/2)i-1/2，j)·Πi，j+(H/2)i-1/2，j·Πi-1，j]/(Δφ)2+[(H/2)i，j+1/2·Πi，j+1-((H/2)i，j+1/2+(H/2)i，j-1/2)·

Πi，j+(H/2)i，j-1/2·Πi，j-1]/(Δλ)2

=[Πi+1/2，j·Hi+1，j-(Πi+1/2，j+Πi-1/2，j)·Hi，j+Πi-1/2，j·Hi-1，j]/(Δφ)2+[Πi，j+1/2·Hi，j+1-(Πi，j+1/2+

Πi，j-1/2)·Hi，j+Πi，j-1/2·Hi，j-1]/(Δλ)2+

(15)

将上式化简整理可得：

(Hi+1/2，j+Hi，j-Hi+1，j)·Πi+1，j+(Hi-1/2，j+

Hi，j-Hi-1，j)·Πi-1，j+[(Δφ/Δλ)2·(Hi，j+1/2+Hi，j-

Hi，j+1)]·Πi，j+1+[(Δφ/Δλ)2(Hi，j-1/2+Hi，j-Hi，j-1)]·Πi，j-1-[Hi+1/2，j+Hi-1/2，j+Hi+1，j+Hi-1，j-

2Hi，j+(Δφ/Δλ)2·(Hi，j+1/2+Hi，j-1/2+Hi，j+1+

Hi，j-1-2Hi，j)]·Πi，j=Λ(Si+1，j-Si-1，j)Δφ

(16)

经整理成：

Ai，jΠi+1，j+Bi，jΠi-1，j+Ci，jΠi，j+1+Di，jΠi，j-1-Ei，jΠi，j=Fi，j

(17)

其中各系数为

(18)

因此，Πi，j可表示为

Πi，j=(Ai，jΠi+1，j+Bi，jΠi-1，j+Ci，jΠi，j+1+

Di，jΠi，j-1-Fi，j)/Ei，j

(19)

任意一点的气膜压力Pi，j可表示为

(20)

轴承在垂直和水平方向的量纲一气膜合力可按下式计算

(21)

式中：Fx和Fy分别为轴承在垂直和水平方向上的气膜合力；P0为静态气膜合力。

轴承的量纲一承载力和偏位角[16]为

(22)

由于气膜的剪切流动与压力流动，轴颈会与轴承内表面会发生摩擦，轴颈表面的周向摩擦因数的表达式为

(23)

式中：“-”表示摩擦力方向与轴承旋转方向相反。

由于气体润滑剂的可压缩性，任一区域流入流出的体积流量不一定守恒，所以用质量流量表征其流动性更为合理。单位时间内流进轴瓦周向截面的质量流量可表示为

(24)

基于上述模型及数值求解方法，采用有限差分法和超松弛迭代法编制程序，对图1所示的三叶动压气体滑动轴承进行静态性能计算，具体计算流程如图3所示。

图3 静态性能计算流程

3 结果与讨论

为验证上述理论模型的正确性及编制程序的可靠性，文中计算了当Λ=1，L/D=2，β=38.9°；Λ=0.5，L/D=2，β=38.9°和Λ=0.5，L/D=0.5，β=10°时的轴承承载力，并与文献[20]的结果进行了对比，如图4所示。两者结果吻合良好，误差在2%以内，说明了文中有限差分算法和程序的正确性。

图4 承载力计算结果对比

3.1 气膜厚度与气膜压力分布

在ε=0.6，mp=0.4，Λ=3和L/D=1.5时，计算得到不同瓦块分布位置的量纲一气膜厚度和气膜压力分布，如图5所示。瓦块分布位置对气膜厚度和压力分布有着明显影响，沿着轴颈旋转方向，气膜厚度先逐渐增大后减小再增大，且在最小气膜厚度位置处出现最大气膜压力。瓦上承载方式的最小气膜厚度(为0.29)大于瓦间承载的最小气膜厚度(为0.25)，且最大气膜压力高出瓦间承载方式近30%。这是因为瓦上承载方式的主要承载区域位于主要承载瓦2上，而瓦间承载方式的主要承载区域位于瓦1和瓦2间，在外载荷作用下润滑气体经轴向槽的泄漏量大于瓦上承载。

3.2 偏心率和预负荷系数对静态性能的影响

图6给出了当Λ=0.6，L/D=1.5时，承载力W、偏位角θ、轴颈表面周向摩擦因数-fJ和单位时间内流入瓦2的质量流量Qmx_in2随偏心率和预负荷系数mp(mp=0.1、0.3和0.5)的变化曲线。由图6(a)(c)可知，随着偏心率和预负荷系数的增大，W和-fJ迅速增大；同时，在不同预负荷系数条件下，瓦上承载方式的承载能力始终优于瓦间承载方式。从图6(c)中可以观察到，当偏心率小于0.4时，瓦上承载和瓦间承载方式的摩擦因数相差不大；在预负荷系数为0.3和0.5且偏心率大于0.4时，瓦间承载方式的-fJ略微大于瓦上承载方式；当预负荷系数为0.7时，小偏心率下的瓦上和瓦间承载方式的摩擦因数几乎一样；但当偏心率为0.8时，出现瓦上承载方式的摩擦因数高于瓦间承载方式的情况，这是由于高预负荷系数和大偏心率的条件使瓦上承载方式的轴颈表面更贴近主要承载瓦。从图6(b)(d)可知，瓦上承载方式的偏位角θ呈现出随预负荷系数的增大而逐渐减小，随偏心率的增大先增大后逐渐减小的趋势；随着预负荷系数和偏心率的增大，Qmx_in2逐渐下降，这是因为轴瓦间隙随着预负荷系数的增大而变小，所以质量流量减小。瓦间承载方式的偏位角θ和Qmx_in2随偏心率和预负荷系数的增大而迅速减小。

图6 偏心率和预负荷系数对静态性能的影响(Λ=0.6，L/D=1.5)

3.3 轴承数和长径比对静态性能的影响

图7示出了当ε=0.6，mp=0.4，L/D分别为0.5、1和1.5时，轴承数对承载力、摩擦因数和质量流量的影响曲线。由图7(a)可知，承载力随Λ和L/D的增大而呈现增长速率放缓的变化趋势，这是由于轴承数的增加提高了轴颈转速，增强了气体动压效应，从而使轴承承载力变大；同时，类似承载力随偏心率和预负荷系数的变化曲线，瓦上承载方式的承载力始终高于瓦间承载方式。由图7(b)可知，-fJ随轴承数和长径比的增加而线性增加，瓦上承载方式的摩擦因数与瓦间承载方式相差不大。由图7(c)可知，瓦上承载方式的Qmx_in2随轴承数的增大先减小而后缓慢增加，且随长径比的增加而快速增大，这是因为承载瓦块2的面积随长径比的增大而变大。瓦间承载方式的Qmx_in2随轴承数的增加而缓慢增大，随长径比的增加而明显增大。由于瓦间承载方式较大的气体端泄量，瓦上承载方式的Qmx_in2始终高于瓦间承载方式。