多葉動壓氣體滑動軸承靜態特性的有限差分算法*

陳 陽 張功學 吳 垚

(陜西科技大學機電工程學院 陜西西安 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始終高于瓦間承載方式。

圖7 軸承數和長徑比對靜態性能的影響(ε=0.6，mp=0.4)

4 結論

(1)長徑比和軸承數對多葉動壓氣體軸承的承載能力影響顯著，長徑比的增大可以明顯提高軸承的承載能力。這是由于長徑比改變了瓦塊承載面積大小。軸承數增加對應的是軸頸轉速增大，轉速提高增強了氣體動壓效應，從而提高了多葉軸承的承載能力。

(2)預負荷系數可為多葉動壓氣體軸承的軸瓦提供預緊壓力，較大的預負荷系數可以提高軸承的承載力，但也會使軸頸表面周向摩擦因數增大，同時減小了軸瓦間隙，使單位時間內流入瓦2的氣體質量流量減小。

(3)瓦塊分布方式對多葉動壓氣體滑動軸承的靜態性能影響顯著，瓦上承載方式的承載能力和質量流量要高于瓦間承載方式，瓦間承載方式的軸承在高速運行時會出現氣體從瓦間溝槽泄漏的現象，從而使承載力和質量流量下降，因此在設計軸承-轉子系統時應充分考慮瓦塊分布位置的影響。