以高壓CO2為介質的推力箔片軸承靜特性分析

車國铚，楊啟超，張克龍，高志成

(1.青島科技大學 機電工程學院，山東 青島 266061；2.中國船舶重工業集團公司 第719研究所，武漢 430064；3.廣東智空動力科技有限公司，廣東 佛山 528216)

超臨界二氧化碳布雷頓循環作為傳統蒸汽動力循環的替代方案，使用超臨界區的CO2作為循環介質，比氦氣布雷頓循環以及蒸汽朗肯循環的效率更高[1]，而且超臨界二氧化碳(S-CO2)的密度更大，使系統設備尺寸減小，具有良好的經濟性，因此在核電、火電、太陽能發電等能源利用領域具有良好的發展前景[2]。

壓縮機作為超臨界二氧化碳布雷頓循環中的核心關鍵設備，其安全高效運行對布雷頓循環效率有至關重要的影響。對S-CO2壓縮機及其關鍵軸承零部件的設計是目前透平機領域研究的熱點[3]。

目前，對推力箔片軸承的研究多以空氣為潤滑介質，而對采用較高壓力的CO2為潤滑介質的研究較少，且性能還沒有太多明確的結果[4]。與空氣相比，高壓CO2為潤滑介質對預測和計算箔片軸承的性能提出了新挑戰，如湍流狀態、非線性熱力學性質等，這成為研究高壓CO2推力箔片軸承性能的動機，也成為近年氣體軸承研究的熱點[5]。

國外對S-CO2布雷頓循環壓縮機中氣體軸承的研究較早，主要包括承載力、溫度特性、摩擦功耗等，而國內近幾年才開展研究。文獻[6]率先搭建功率為250 kW的試驗臺研究壓縮機中采用CO2直接潤滑的氣體軸承，在壓縮回路中通過熱敏電阻測得推力軸承出口處氣體溫度為15～37 ℃，指出轉子腔的壓力應設計在2 MPa以下，理想情況約為1.4 MPa，以降低轉子和氣體軸承的風阻損失。文獻[7]在高達4.8 MPa的CO2中單獨運行氣體箔片軸承，研究其在高壓CO2中的功率損失，結果表明氣體箔片軸承的功率損失隨CO2壓力的增加而增大。文獻[8]考慮真實氣體效應和薄膜內流動湍流狀態，對混合推力軸承進行了三維熱流體動力學分析，計算了一定尺寸推力軸承的壓力分布、承載力、功率損失等。文獻[9]采用計算流體動力學的三維流體結構模型，研究了高壓CO2箔片推力軸承的彈性流體動力性能和運行條件對其性能的影響，比較了不同流動狀態下箔片推力軸承性能的變化。文獻[10]在連續性方程和動量方程的基礎上推導層流狀態下適用于S-CO2徑向軸承的雷諾方程，并考慮湍流速度脈動的影響，在雷諾方程中加入修正系數，采用差值法獲得S-CO2物性，并用有限差分數值法初步計算了推力氣體軸承的靜態性能，但未分析工況及結構參數對徑向氣體軸承靜特性的影響。

本文考慮CO2的實際物性及湍流狀態對潤滑氣膜的影響，運用有限差分法數值求解流體潤滑雷諾方程，分析結構及工況參數對推力箔片軸承性能的影響，并與空氣為工作介質的軸承性能進行對比。

1 推力箔片軸承結構

8個瓦塊的推力箔片軸承示意圖如圖1所示，其中R1為軸承內半徑，R2為軸承外半徑，β為扇形瓦張角，b為節距比，表示平箔楔形部分占整個瓦塊的比例，δh為楔形面高度，h2為平箔水平面與推力盤之間的間隙，也稱最小氣膜厚度，h1為楔形部分與水平面部分高度的和。瓦塊上的每個點都可以用徑向方向上的半徑長度r和圓周方向上的角度θ表示，lin和lout分別為進氣口和出氣口。

圖1 推力箔片軸承結構示意圖

通常，平箔和波箔的一端通過點焊固定在軸承座上。平箔由楔形部分和水平部分組成，與推力盤表面之間形成楔形間隙。轉子運行過程中氣體被吸入推力盤與推力軸承的間隙，沿間隙變小的周向方向運動，形成動壓效應并產生一定的承載力以平衡軸向力。與滾動軸承相比，氣體箔片軸承具有摩擦功耗小的優點，在轉子受到不穩定渦動時，通過具有柔性支承的彈性箔片之間的相互作用吸收部分渦動能量，使轉子運行穩定[11]。

2 數學模型及求解

2.1 變密度變黏度湍流雷諾方程

考慮密度、黏度變化以及高壓狀態下湍流狀態對潤滑氣膜的影響，修正雷諾方程為

(1)

為簡化計算，對(1)式量綱一化，令

(2)

式中：pa為軸承外部環境壓力；μa為軸承外部環境下CO2的動力黏度；ρa為軸承外部環境下CO2的密度。

(1)式變為

(3)

2.2 氣膜厚度方程

假設推力盤無傾斜時，每個扇形瓦塊結構相同，因此每個扇形瓦塊之間潤滑氣膜的厚度分布相同。氣體箔片軸承工作時，箔片的彈性變形會影響軸承的氣膜厚度，數值分析時需考慮箔片的彈性變形。本文所采用的彈性箔片結構如圖2所示，s為波箔單元長度，tb為波箔厚度,tp為平箔厚度,l為波箔波紋直徑。本文采用經典Heshmat模型，不考慮平箔的剛度，忽略彈性箔片之間的相互摩擦作用以及波箔與軸承座的摩擦，則箔片剛度kb為[12]

圖2 波箔和平箔結構

(4)

其量綱一化為

(5)

式中：Eb為波箔材料的彈性模量；νb為波箔材料的泊松比。

對于彈性箔片軸承，潤滑氣膜厚度為

h=h2+g(θ)+u(r,θ),

(6)

式中：g(θ)為間隙內的潤滑氣膜厚度；u(r，θ)為彈性變形量。

對h2，g(θ)，u(r，θ)分別進行量綱一化，令

(7)

(6)式變為

(8)

2.3 變密度變黏度湍流雷諾方程的數值求解

采用有限差分數值法將扇形瓦塊的計算域離散為差分網格，用有限個網格節點替代連續的求解域，潤滑氣膜展開為m×n的離散網格(圖3)，圖中Δθ，Δr分別為周向和徑向的網格間距，i，j分別為周向和徑向的網格節點編號。

圖3 計算域網格劃分示意圖

基于有限差分法，采用計算精度較高的中心差商，將控制方程中周向和徑向的壓力偏導數用相鄰節點函數值的中心差商近似表示，即

(9)

根據(9)式對(3)式進行離散化。

(3)式左邊第1項離散為

(10)

(3)式左邊第2項離散為

(11)

由于靜態計算無需考慮時間項，因此(3)式右邊第一項離散為

(12)

將(10)—(12)式代入(3)式得

pi,j=

(13)

(14)

扇形瓦塊的每個邊界都與外界環境接觸，因此每個邊界都為外界環境條件，故有

(15)

為加快迭代計算過程的收斂速度，采用松弛迭代法進行計算，迭代格式為

(16)

采用相對收斂準則判斷迭代是否達到精度要求，達到精度要求后終止迭代。本文相對收斂準則參照文獻[11]選取，即

(17)

2.4 計算對象及方法

2.4.1 計算對象

采用文獻[13]中空氣循環機推力箔片軸承尺寸，確定推力箔片軸承結構參數見表1。

表1 推力箔片軸承結構參數

2.4.2 密度和黏度處理方法

采用NIST Refprop 軟件獲取CO2流體熱物性，不同溫度及壓力下的密度、黏度都可通過軟件編程調用。

2.4.3 計算流程

整體迭代計算流程如圖4所示。

圖4 氣膜壓力計算流程圖

2.5 靜態性能計算

聯立(3)和(8)式，獲得楔形間隙內潤滑氣膜的壓力分布，在軸承計算區域對其積分得到承載力和摩擦力矩。

(18)

(19)

(20)

3 數值結果及分析

3.1 計算程序驗證

根據圖4的計算流程，用MATLAB編寫數值計算程序。以CO2為潤滑介質的推力軸承研究較少，本文計算空氣作為潤滑介質的剛性表面與彈性表面的軸承性能，并與文獻[15]的數據進行對比，結果見表2。計算數據為量綱一的值，k為0時為剛性表面，h1/h2為氣膜間隙比。由表2可知最大誤差為7.14%，驗證了本文所用計算程序的準確性。

表2 本文與文獻[15]計算結果對比

3.2 推力軸承靜態性能計算結果

轉速n為40 000 r/min, 環境溫度T為300 K，pa為1.4 MPa時，CO2為潤滑介質的推力箔片軸承氣膜壓力和厚度分布分別如圖5和圖6所示。圖5中壓力沿扇形瓦氣流進口到出口方向先增大后減小，在傾斜面與水平面的交界處附近達到最大，CO2氣膜壓力峰值為1.5 MPa,各邊界處壓力為環境壓力；圖6中氣膜厚度在徑向方向上恒定，最終在氣流出口降至環境壓力。以圖5氣膜壓力的計算結果為基礎，計算得到承載力為344.88 N，摩擦功率損失為108.2 W。

圖5 高壓CO2推力箔片軸承氣膜壓力分布

圖6 高壓CO2推力泊片軸承氣膜厚度分布

T為300 K，pa為1.4 MPa，不同轉速下瓦塊徑向中截面處的氣膜壓力分布如圖7所示，隨轉速升高，氣膜壓力呈增大的趨勢。

圖7 高壓CO2時扇形瓦塊徑向中截面處氣膜壓力曲線

采用同一推力箔片軸承，分別計算高壓CO2與空氣為潤滑介質時扇形瓦塊徑向中截面處氣膜厚度曲線如圖8所示，在相同轉速下，高壓CO2為潤滑介質時的氣膜厚度遠大于空氣，說明高壓CO2使彈性箔片發生了較大的變形；在不同轉速下，高壓CO2的軸承氣膜厚度差異較大，而空氣時差異較小。

圖8 高壓CO2與空氣時扇形瓦塊徑向中截面處氣膜厚度曲線

空氣與高壓CO2作為潤滑介質時的承載力對比如圖9所示，空氣產生的承載力遠小于高壓CO2產生的承載力；隨轉速提高，空氣產生的承載力變化較小，而對應的高壓CO2變化較大；空氣潤滑介質軸承承載力僅為1.4 MPa下CO2潤滑介質軸承承載力的41%。

圖9 高壓CO2與空氣時軸承承載力對比

3.3 節距比對推力箔片軸承靜態性能的影響

桑迪亞國家實驗室的壓縮機設計最高轉速為75 000 r /min，故以此為參考，n為75 000 r/min，T為300 K，pa為1.4 MPa時，不同扇形瓦張角下推力箔片軸承瓦塊節距比與軸承承載力的關系如圖10所示，隨節距比增大，軸承承載力呈先增大后減小的趨勢。

圖10 不同瓦張角下節距比對軸承承載力的影響

定義達到軸承最大承載力時的節距比為最佳，由圖10可知，最佳節距比隨扇形瓦張角的增大而減小，軸承最大承載力隨扇形瓦張角的增大呈先增大后減小的趨勢，本文計算條件下當扇形瓦張角為45°時軸承承載力出現了最大值，約為429.5 N，反映不同節距比時，軸承最大承載力對應的最佳瓦張角不同。

推力箔片軸承的摩擦功耗與節距比的關系如圖11 所示，摩擦功耗隨節距比增大呈先增大后減小的趨勢，這是因為軸承的摩擦功耗與氣膜壓力直接相關，與圖10中軸承承載力的變化規律一致。當扇形瓦張角較小，即設計更多的瓦塊時，間隙內氣膜與推力盤摩擦導致的摩擦功耗越大。當推力箔片軸承的瓦張角為30°、節距比為0.5時，摩擦功耗接近600 W,當轉速繼續提高時，摩擦功耗繼續增大。

圖11 不同瓦張角下節距比對摩擦功耗的影響

由圖10和圖11可知，扇形瓦張角為45°、節距比為0.5時軸承承載力較大、摩擦功耗較小，因此本文后續數值分析計算采用此條件。

3.4 楔形高度對軸承靜態性能的影響

固定最小氣膜厚度h2為10 μm，研究不同楔形高度對軸承承載力的影響(圖12)，楔形高度反應了楔形間隙傾斜角的大小，楔形高度逐漸增大，即斜面傾斜角逐漸增大時，軸承承載力先增大后減小；當其他條件固定時，轉速為70 000 r/min時達到軸承最大承載力；隨轉速提高，軸承最大承載力對應的楔形高度有增大的趨勢。

圖12 不同轉速下楔形高度對承載力的影響

軸承在不同楔形高度和轉速下的摩擦功耗變化規律如圖13所示，摩擦功耗隨著楔形高度的增大而逐漸增大，且高轉速、大楔形高度對應的摩擦功耗較大。

圖13 不同轉速下楔形高度對摩擦功耗的影響

當n為40 000 r/min，T為300 K，pa為1.4 MPa時，不同楔形高度對應的扇形瓦塊徑向中截面壓力分布如圖14所示，隨著楔形高度的增大，軸承水平區域壓力呈先增大后減小的趨勢，楔形區域氣膜壓力呈下降趨勢。楔形高度增大使楔形區域氣膜壓力降低，此區域軸承承載力減小，承載任務主要集中在水平區域，將加劇該區域的摩擦，不利于水平區域的工作。

圖14 不同楔形高度下扇形瓦塊徑向中截面壓力分布曲線

假設楔形高度以50 μm為初始參考值，楔形高度分別設計為100，150，200，250 μm時，楔形區域氣膜壓力達到與50 μm設計條件下相同值的瓦張角分別約為13°，17.5°，19°，21°，瓦張角有后移現象，如圖中紅色圓圈所示。分析其原因為推力箔片軸承的楔形高度對動壓效應具有重要作用，當楔形高度增大時，動壓效應的產生區域會從入口向右移動。

3.5 最小氣膜厚度對軸承靜態性能的影響

推力箔片軸承的楔形高度通常在設計時已固定,即軸承結構參數已確定，而推力箔片軸承與推力盤間隙的最小氣膜厚度h2在裝配以及實際工作過程中隨著軸向力的變化而變化。固定楔形高度為50 μm，不同轉速下最小氣膜厚度變化對軸承承載力和摩擦功耗的影響分別如圖15和圖16所示。

圖15 不同轉速下最小氣膜厚度對軸承承載力的影響

圖16 不同轉速下最小氣膜厚度對摩擦功耗的影響

由圖15可知，當最小氣膜厚度增大時，軸承承載力呈減速下降的趨勢，且各轉速下的軸承承載力差值逐漸縮小；由圖16可知，在不同轉速下軸承的摩擦功耗均隨最小氣膜厚度的增大而減小，高轉速的摩擦功耗明顯高于低轉速。軸承間隙的設計需在滿足軸承承載力的前提下盡量減小摩擦功耗。

4 結論

通過求解修正雷諾方程，獲得高壓CO2在推力箔片軸承中的壓力分布，求解軸承承載力和摩擦功耗，分析后得到以下結論：

1)在本文計算模型下，高壓CO2與空氣為潤滑介質的軸承相比，高壓會產生較大的箔片彈性變形；在相同轉速下，高壓CO2為潤滑介質時的氣膜厚度遠大于空氣的，且在不同轉速條件下高壓CO2的氣膜厚度差異較大，而空氣時差異較小；高壓CO2為潤滑介質的軸承承載力遠大于空氣的，空氣在不同轉速條件下的軸承承載力差異較小，而高壓CO2差異較大。

2)瓦張角一定，節距比增大時，高壓CO2為介質的軸承承載力呈先增大后減小的趨勢；隨著扇形瓦張角的增大，不同節距比的軸承最大承載力出現先增大后減小的趨勢，不同的節距比，最佳瓦張角也不同；本文研究的推力箔片軸承，扇形瓦張角為45°、節距比為0.5時有較大的軸承承載力和較低的摩擦功耗。

3)當最小氣膜厚度固定，存在一個最優楔形高度或最佳斜面傾斜角使軸承承載力最大，本文研究的推力箔片軸承最優楔形高度約為100 μm，此時傾斜角約為0.57°。