以超臨界二氧化碳為介質的徑向箔片軸承靜特性分析*

車國铚 楊啟超 滕 斌 魏志國 高志成

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

超臨界二氧化碳布雷頓發電循環是以超臨界狀態的CO2作為循環工質的熱力循環，二氧化碳的臨界壓力為7.38 MPa,臨界溫度為31.1 ℃[1]，在靠近臨界點的區域，S-CO2表現出密度大、黏度小的特點，這一物性特點使得占渦輪機輸出功比重較大的壓縮功減少，提高了S-CO2布雷頓循環的效率[2]。在同等條件(壓力8～20 MPa和溫度450～650 ℃)下，S-CO2布雷頓循環的效率在應用較多的蒸汽朗肯循環和氦氣布雷頓循環效率之上[3]。除此之外，S-CO2具有更高的介質密度，減小了系統設備的尺寸，使得系統結構緊湊[4]，在經濟性方面具有更好的潛力。因此超臨界二氧化碳布雷頓循環在太陽能、火電、核電等能源領域具有良好的發展前景[5]。美國桑迪亞(Sandia)國家實驗室[6]、韓國能源技術研究院(KIER)[7]、日本東京工業大學(Tokyo Institute of Technology，TIT)[8]等研究機構均率先對S-CO2布雷頓循環進行了理論及實驗研究，并在此基礎上開發了不同功率的循環系統。

壓縮機是S-CO2布雷頓循環系統中一個關鍵設備，S-CO2循環工質通過壓縮機提高至渦輪機的入口壓力，最終高溫高壓的S-CO2流體在渦輪機中做功。為了最大化系統效率，壓縮機轉子的轉速通常較高，可達到100 000 r/min[9],因此用于支撐高轉速轉子的軸承技術是超臨界二氧化碳布雷頓循環壓縮機的關鍵技術之一。美國Sandia國家實驗室曾選用滾珠軸承進行軸向力的承載，但軸承的運行壽命較短，其滾珠軸承壽命有20～2 000 h不等，后期采用氣體箔片軸承來進行長時間的測試[10]。FLEMING等[11]曾根據S-CO2循環的不同功率等級提出了S-CO2透平機械以及其中關鍵部件的推薦路線。眾研究機構設計實驗的超臨界二氧化碳循環功率多在100～300 kW范圍[12]，根據推薦路線，軸承推薦選擇氣體軸承。并且美國Sandia國家實驗室、TIT等在實驗測試樣機的軸承選擇上均選用以S-CO2直接潤滑氣體軸承，成功地實現穩定、長時間的轉子力的承載，驗證了氣體軸承的技術可靠性[13]。

以上眾研究機構在設計S-CO2壓縮機時，壓縮機的進口壓力設計在7.6～9 MPa范圍，進口溫度設計在305～330 K范圍。該進口溫度及壓力處在臨界點附近，S-CO2物性在臨界點附近呈現出劇烈的非線性變化。在利用S-CO2的特殊物性減少壓縮機耗功、提高循環效率的同時，也應該考慮物性的劇烈變化對軸承性能以及轉子穩定性的影響。傳統動壓徑向氣體軸承的研究主要以空氣為工質，發展至今，在數值模擬、理論分析及試驗研究方面都取得了較多的成果。2004年虞烈[14]通過引入彈性箔片靜、動變形，聯立求解氣體潤滑Reynolds方程，給出了以空氣為潤滑介質的箔片動壓軸承的完全氣彈潤滑解。2006年楊利花等[15]曾建立氣體軸承性能測試試驗臺，并對徑向箔片氣體軸承進行承載力和起飛轉速的實驗研究。2016年賈晨輝等[16]采用流體動力學軟件，以動壓氣體軸承為研究對象，研究氣體軸承瞬態非線性動力學行為，探討軸承-轉子系統的穩定性。2017年馮凱等人[17]建立了箔片徑向氣體軸承靜態特性和溫度特性的測試試驗臺，通過變載荷變轉速方式對起飛轉速、摩擦力矩和溫度進行測量。但目前對于S-CO2潤滑介質氣體軸承的研究還較少，近臨界區S-CO2密度、黏度的大幅變化，以及高壓下S-CO2內部流動的湍流狀態，這都將增加對預測、計算徑向軸承性能的難度。另外，基于空氣工質的氣體軸承研究及設計方法是否適用的問題也是研究學者關注的方向。

2017年國內溫建全[18]在連續性方程和動量方程的基礎上推導層流狀態下適用于S-CO2徑向軸承的雷諾方程，并考慮湍流速度脈動的影響，以修正系數的方式加入雷諾方程中。計算時S-CO2物性采用差值方法獲得，并用有限差分數值計算方法初步計算了徑向氣體軸承的靜態性能，但未分析工況參數及結構參數對徑向氣體軸承靜特性的影響。

2013年國外CONBOY[19]建立了一個彈性流體動力學模型，考慮潤滑層內的流體速度場、流體動壓和摩擦損失,計算了S-CO2潤滑的推力軸承中的壓力場，結果表明，計算域內所有離散點均為湍流區。除此之外，進行了不同轉速以及不同軸承外半徑下，空氣止推軸承與S-CO2止推軸承的承載能力對比。

此外，2020年中國科學院大學BI等[20]為研究S-CO2動壓徑向氣體軸承的動態特性，考慮了S-CO2湍流狀態以及非線性熱力學性質，在傳統擾動法的基礎上，引入擾動密度、擾動黏度、擾動雷諾數和擾動湍流系數，擴展了頻率擾動法，用數值方法計算了包含所有動力學變量的頻率相關剛度和阻尼系數。

本文作者考慮S-CO2的實際物性以及湍流狀態的影響，運用有限差分法數值求解雷諾方程，重點分析近臨界區S-CO2的物性突變對氣膜、承載力、摩擦力矩等的影響；另外分析工況參數以及軸承結構參數對徑向箔片氣體軸承靜特性的影響，為S-CO2布雷頓循環壓縮機中徑向氣體軸承的研究提供理論依據。

1 箔片動壓徑向氣體軸承結構

圖1為文中研究的波箔型動壓徑向氣體軸承結構示意圖，軸承座里包含以特殊材料制成的耐磨、耐高溫的彈性波箔片以及平箔片，平箔片和波箔片的一端固定在軸承座上，另一端可自由伸縮。軸頸在旋轉過程中，帶動潤滑氣體由間隙大端向間隙小端運動，形成動壓效應，潤滑氣膜因此產生承載力。軸承在工作時，彈性箔片之間的變形作用以及其之間的庫侖摩擦作用，可以吸收不平衡渦動能量[21]，即使在受到不穩定渦動時，氣體軸承也能保持較好的穩定性。

圖1中，O為軸承的中心；Ob為軸頸的中心；e為軸頸的中心相對于軸承中心之間的偏心距；R為軸頸的半徑；Ω為軸頸的角速度；hmax為最大楔形間隙處的氣膜厚度；hmin為最小楔形間隙處的氣膜厚度。

2 數學模型及求解

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

傳統的動壓氣體軸承多以空氣為工質進行研究，空氣的物性變化較小，但近臨界區的S-CO2物性變化劇烈，微小的溫度壓力變化將會引起密度與黏度的劇烈變化。除此之外，由于S-CO2密度大、黏度小，導致潤滑氣膜中會出現湍流狀態，在氣體潤滑模型中必須加以考慮。因此，一般形式的可壓縮氣體雷諾方程便不再適用于S-CO2動壓氣體軸承。文中采用文獻[18]中由連續性方程和動量守恒方程導出的變密度變黏度層流雷諾方程，在此基礎上，將湍流的影響以周向和軸向的湍流潤滑系數加入變密度變黏度層流雷諾方程中，其中kx為周向湍流潤滑系數，kz為軸向湍流潤滑系數，這是與雷諾數有關的非線性函數[22]，可得如下形式的雷諾方程：

(1)

(2)

(3)

式中:h為潤滑氣膜實際厚度(mm)；C為軸承間隙(mm)；pa為徑向軸承外部環境壓力(Pa);μa為軸承外部環境下的S-CO2動力黏度(Pa·s)；ρa為軸承外部環境下S-CO2密度(kg/m3)；L為軸承長度(mm)。

對式(1)量綱一化后得到：

(4)

2.2 氣膜厚度方程

箔片氣體軸承在工作時，不同于剛性氣體軸承，箔片的彈性變形會影響軸承的氣膜厚度大小，數值分析時需考慮箔片的彈性變形方程。文中所采用的彈性箔片結構如圖2所示，模型中不考慮平箔的剛度，忽略彈性箔片之間的相互摩擦作用以及波箔與軸承座的摩擦，則箔片柔度[23]為

(5)

式中：pa為軸承外部環境壓力(Pa);s為波箔單元的長度(mm);l為半波箔長度(mm)；Eb為波箔材料彈性模量(Pa)；νb為波箔片的泊松系數;tb為波箔片厚度(mm)。

一定偏心率下軸頸靜平衡位置處的氣膜厚度由2部分組成，一部分為箔片變形前的幾何間隙，另一部分為箔片變形的變形間隙。量綱一氣膜厚度表達式為

(6)

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

采用有限差分數值方法，將軸承的計算域離散化為差分網格，用有限個網格節點代替連續的求解域，潤滑氣膜展開如圖3所示的m×n的離散網格。圖中Δθ、Δz分別為周向和軸向的網格間距，Δθ=2π/m，Δz=2/n，i、j代表周向和軸向的網格節點編號。

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

(7)

式中：i+1/2和i-1/2分別為i后半點處和前半點處；j+1/2和j-1/2為j后半點處和前半點處。

根據式(7)對式(4)進行離散化。

離散式(4)左邊第一項：

(8)

離散式(4)左邊第二項：

(9)

離散式(4)右邊項：

(10)

將式(8)、(9)、(10)代入量綱一化方程(4)整理得到如下矩陣方程組：

Pi,j=

(11)

式中：Ai,j、Bi,j、Ci,j、Di,j、Ei,j、Fi,j為系數矩陣，各系數矩陣計算式如公式(12)所示。

(12)

由于軸承的兩端(z=1,z=n+1)以及彈性箔片的固定端與軸承外部是相通的，且氣膜未展開之前，該部分(θ=1，θ=m+1)是重合的，所以這幾個區域的壓力均為徑向氣體軸承的工作環境壓力，故量綱一求解邊界條件為

(13)

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

(14)

通常數值計算的迭代過程中，常采用一定的相對收斂準則來判斷迭代是否達到精度要求，達到精度要求后，便可終止迭代。文中相對收斂準則為

(15)

2.4 計算對象及方法

2.4.1 計算對象

文中動壓徑向軸承的結構參數以及潤滑氣體參數見表1。

表1 動壓徑向氣體軸承相關參數

2.4.2 密度黏度處理方法

對于S-CO2的物性處理，文中采用美國國家標準與技術研究院的NIST Refprop 軟件獲取S-CO2流體熱物性，并通過軟件編程調用其物性。

2.4.3 計算流程

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

2.5 靜態性能計算

通過耦合變密度變黏度湍流雷諾方程與箔片變形方程，求解獲得潤滑氣膜的氣膜壓力分布，在此基礎上進行靜態特性的計算，其中包括承載力、偏位角、摩擦力矩等。

文中軸頸外載荷方向假定垂直向下，因此氣膜的合力方向垂直向上，與外載荷相互平衡。假定偏心方向上的力為Wn，垂直偏心方向的力為Wt，則其計算公式為

(16)

(17)

式中：負號方向代表力垂直于軸頸并指向里。

軸承的量綱一承載力為

(18)

有量綱承載力為

(19)

偏位角滿足如下關系

(20)

(21)

作用在軸頸表面的摩擦力矩其量綱一化形式為

(22)

則有量綱摩擦力矩為

(23)

3 數值結果及分析

3.1 程序驗證

根據圖4所示流程，用MATLAB 編寫數值計算程序。為驗證程序的正確性，計算剛性表面軸承不同工況下的承載力、偏位角，并與文獻[18]結果對比。對比結果如表2所示，其中L/D為軸承的長徑比；Λ為軸承數；ε為軸頸偏心率；Re為初始雷諾數。從表2可以看出文中計算結果與文獻[18]結果非常接近，最大誤差為2.9%，因此驗證了文中所用計算程序具有一定的準確性。

表2 文中計算結果與文獻[18]結果對比

3.2 近臨界區物性及軸承靜態特性計算結果

圖5和圖6所示分別為S-CO2的密度與黏度在不同壓力下隨溫度的變化規律。可以看出，S-CO2密度與黏度隨溫度的升高均呈現減小的趨勢，且越接近臨界點時，變化越顯著。當壓力為7.6 MPa，溫度從305 K升高1 K時，密度減小約38%，黏度減小約39%。而與壓力的關系則相反，當溫度一定時，S-CO2密度與黏度隨壓力的升高而增大。當溫度為306 K,壓力從7.6 MPa升高至8.0 MPa時，密度增大約73.8%，黏度增大約80%。氣體軸承的靜態特性與工作介質的黏度、壓力及溫度都密切相關，因此，研究物性變化對軸承靜態性能的影響至關重要。

圖7和圖8所示分別為轉速48 000 r/min、偏心率0.6條件下的S-CO2潤滑氣膜的厚度分布與氣膜壓力分布。由圖7可以看出，氣膜厚度在軸承軸向上呈現出兩端小，中間大的分布；在圓周方向上呈現為先下降后增大的趨勢。軸向中間氣膜厚度大的原因是箔片未變形時此處壓力較高，彈性箔片發生變形，產生額外的間隙空間，氣膜厚度增大。

圖8中，量綱一最大氣膜壓力為1.01，約為7.6 MPa,量綱一最小氣膜壓力為0.987，約為7.42 MPa,氣膜壓力峰值出現在圓周方向150°左右，氣膜壓力在圓周角度呈現先增后減的趨勢。這是由于轉子帶動黏性氣體運轉的過程中，間隙由大變小，再由小變大。間隙大的地方氣膜厚度大，壓力較低。間隙小的地方氣膜厚度小，壓力較高。

氣膜的厚度和壓力在軸承楔形間隙中的分布與軸承的性能息息相關，且溫度、壓力引起的物性變化也最直接體現在潤滑氣膜的變化上。圖9和圖10所示分別為7.515 MPa下，氣膜厚度與氣膜壓力軸向中截面隨溫度變化的曲線。從圖9可以看出，隨著溫度的升高，氣膜壓力峰值逐漸降低。溫度為307 K時，量綱一氣膜壓力峰值約為1.015，而溫度為313 K時，量綱一氣膜壓力峰值約為1.013，壓力峰值降低約0.015 MPa。這是由于隨著溫度的升高，S-CO2黏性減小，動壓效應減弱，壓力有所降低。圖10中，最小氣膜厚度同樣隨著溫度升高而減小。溫度為307 K時，量綱一氣膜厚度約為0.75，而溫度為313K時，量綱一氣膜厚度約為0.65，降低了約0.1。這是由于黏性減小導致動壓效應減弱的緣故，壓力有所降低，箔片變形減小，從而產生的額外間隙減小。

圖11和圖12所示分別為轉速48 000 r/min、偏心率0.6條件下的近臨界區承載力、摩擦力矩隨溫度和壓力的變化關系。從圖11可以看出，在305～331 K的近臨界溫度范圍內，隨著溫度的升高，徑向軸承承載力先急劇降低，然后逐漸減小，最終趨于平緩。越靠近臨界點時，溫度的變化對承載力影響越顯著。當壓力為7.6 MPa時，溫度由305 K升至307 K，承載力減小了約22.8%;壓力為8.0 MPa時，溫度由305 K升至307 K,承載力減小了約19.3%。同樣地，壓力的變化對承載力的影響也非常顯著，承載力隨著壓力的提高而增大。如圖11所示，溫度為307 K時，壓力從7.6 MPa升高至8.0 MPa，承載力大幅提高約32%；溫度為309 K時，壓力從8.0 MPa升高至8.4 MPa，承載力提高約25%。相較于傳統滾珠軸承，動壓氣體軸承本身承載力較小，若改變工作條件，徑向軸承承載力降低會不利于壓縮機轉子的穩定運行。因此在設計S-CO2壓縮機時，需設計合適的進出口溫度及壓力。

由圖12可以看出，壓力一定時，隨著溫度升高，徑向軸承摩擦力矩先急劇降低，然后逐漸減小，最終趨于平緩。靠近臨界點的區域摩擦力矩有大幅減小現象，這是由于黏度的大幅減小，潤滑氣膜與軸頸之間的黏性力減小導致的。壓力為7.6 MPa時，當溫度由305 K升高至307 K時，摩擦力矩減小約24%；壓力為8.0 MPa時，當溫度由305 K升高至309 K時，摩擦力矩減小約33.3%。同樣地，壓力的變化對摩擦力矩的影響也非常顯著，摩擦力矩隨著壓力的提高而增大。307 K時，壓力從7.6 MPa升高至8.0 MPa，摩擦力矩增大約28.8%。309 K時，壓力從8.0 MPa升高至8.4 MPa時，摩擦力矩增大約32.5%。摩擦力矩的大幅提高一方面使得功率損耗增大,另一方面會使得軸頸高速運轉過程中所產生的熱量增多，以致于導致過高的溫升使物性發生劇烈變化。

3.3 偏心率對徑向軸承靜態性能的影響

圖13—15描述了T=314 K、p=7.515 MPa條件下，不同轉速時承載力、偏位角、摩擦力矩隨軸頸偏心率的變化。由圖13可看出，同一轉速下，徑向氣體軸承的承載力隨著偏心率增大而增大；當轉速提高時，承載力逐漸變大；在小偏心率時，各轉速下的軸承承載力差距較小，隨著偏心率的提高，差距逐漸變大。在轉子運轉過程中，當偏心率增大時，軸頸與彈性箔片之間的楔形間隙變小，在這樣的情況下，軸頸擠壓氣膜的作用就顯得愈加強烈，由此一來，楔形間隙中的氣膜壓力變高，徑向軸承的承載力得以提高。

從圖14可看出，隨著偏心率的增大，偏位角逐漸減小；同一偏心率條件下，當轉速增大時，偏位角呈減小的趨勢。由于偏位角能表征系統的穩定性，偏位角越小，穩定性相對越好。因此轉速的增大，在一定程度上有利于提高系統的穩定性。

從圖15可看出，隨著偏心率的增大，摩擦力矩呈現出加速增長的趨勢；且轉速越大，摩擦力矩越大，這是因為當轉子的轉速增大時，軸頸與S-CO2潤滑氣膜之間的摩擦加劇，摩擦功耗增大。另外從計算結果可以看出，相同偏心率變化時，轉速越大時摩擦力矩相應變化也越大。

3.4 軸承數對靜態性能的影響

軸承數是表征氣體軸承工作狀態的重要參數之一，由于超臨界二氧化碳布雷頓循環壓縮機的徑向軸承工作環境壓力高，相較于傳統的空氣動壓軸承，軸承數的數值較小。圖16—18描述的是T=314 K、p=7.515 MPa、長徑比為1條件下，不同偏心率時承載力、偏位角以及摩擦力矩隨軸承數的變化。從圖16可看出，隨著軸承數的增大，徑向氣體軸承的承載力呈現出減速增長的趨勢；在小軸承數條件下，承載力較小，且不同偏心率時的承載力差值較小，而隨著軸承數的增大，其差值不斷增大，增大一定程度后趨于不變。

如圖17所示，隨著軸承數的增大，徑向氣體軸承的偏位角呈現出減速下降的趨勢。由于偏位角能表征系統的穩定性，偏位角越小，穩定性相對越好。因此軸承數的增大，有利于保持系統的穩定。

由圖18可以看出，隨著軸承數的增大，徑向氣體軸承的摩擦力矩呈現出增長趨勢，其中小偏心率下近乎線性增長，隨著偏心率的增大，增長趨勢變為減速增長。偏心率為0.6時，軸承數從0.001增大到0.01時，摩擦力矩增大約2.75倍。

3.5 長徑比對靜態性能的影響

軸承的長徑比(L/D)是氣體軸承的重要結構參數，文中研究長徑比對徑向軸承靜態性能的影響。圖19和圖20所示分別為T=314 K、p=7.515 MPa、Λ=0.01條件下,不同偏心率時承載力和偏位角隨長徑比的變化情況。

如圖19所示，隨著長徑比的增大，承載力呈現線性增長趨勢；長徑比為0.5時，即使在大偏心率情況下，承載力也較小，當偏心率為0.6時，長徑比為2時的承載力為長徑比為1條件的2.4倍，擁有較大的承載力。但長徑比較大時軸承的計算應該考慮溫度場的影響。為保證徑向氣體軸承能夠對高速轉子進行穩定支撐，在設計S-CO2壓縮機徑向氣體軸承時應選擇合適的長徑比，以保證合適的承載力。

由圖20可以看出，隨著長徑比的增大，偏位角呈現減速減小的趨勢。在長徑比從0.5增大到1時，偏位角急劇下降，如偏心率為0.6時，偏位角減小了約17.6°；然而長徑比從1增大到2時，偏位角變化卻很小，如偏心率為0.6情況下偏位角僅減小3.5°。

4 結論

考慮S-CO2實際物性，數值求解箔片動壓軸承在不同溫度、壓力工作環境下的承載力和摩擦力矩等靜態特性。主要結論如下：

(1)近臨界區密度、黏度的劇烈變化會引起承載力和摩擦力矩的大幅變化，且越靠近臨界點，變化程度越大；同一壓力下，承載力和摩擦力矩隨著溫度的升高而減小，同一溫度下，承載力和摩擦力矩隨著壓力的升高而增大。

(2)隨著偏心率的增大，S-CO2軸承的承載力和摩擦力矩都有不同程度的增大，而偏位角減小；隨軸承數的增大，S-CO2軸承的承載力和摩擦力矩增大，而偏位角減小；增大徑向軸承的長徑比，S-CO2軸承的承載力增大，而偏位角減小。