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

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

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

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

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

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

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

1 推力箔片軸承結(jié)構(gòu)

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

圖1 推力箔片軸承結(jié)構(gòu)示意圖

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

2 數(shù)學(xué)模型及求解

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

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

(1)

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

(2)

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

(1)式變?yōu)?/p>

(3)

2.2 氣膜厚度方程

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

圖2 波箔和平箔結(jié)構(gòu)

(4)

其量綱一化為

(5)

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

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

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

(6)

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

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

(7)

(6)式變?yōu)?/p>

(8)

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

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

圖3 計算域網(wǎng)格劃分示意圖

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

(9)

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

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

(10)

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

(11)

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

(12)

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

pi,j=

(13)

(14)

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

(15)

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

(16)

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

(17)

2.4 計算對象及方法

2.4.1 計算對象

采用文獻[13]中空氣循環(huán)機推力箔片軸承尺寸，確定推力箔片軸承結(jié)構(gòu)參數(shù)見表1。

表1 推力箔片軸承結(jié)構(gòu)參數(shù)

2.4.2 密度和黏度處理方法

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

2.4.3 計算流程

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

圖4 氣膜壓力計算流程圖

2.5 靜態(tài)性能計算

聯(lián)立(3)和(8)式，獲得楔形間隙內(nèi)潤滑氣膜的壓力分布，在軸承計算區(qū)域?qū)ζ浞e分得到承載力和摩擦力矩。

(18)

(19)

(20)

3 數(shù)值結(jié)果及分析

3.1 計算程序驗證

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

表2 本文與文獻[15]計算結(jié)果對比

3.2 推力軸承靜態(tài)性能計算結(jié)果

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

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

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

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

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

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

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

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

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

3.3 節(jié)距比對推力箔片軸承靜態(tài)性能的影響

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

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

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

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

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

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

3.4 楔形高度對軸承靜態(tài)性能的影響

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

圖12 不同轉(zhuǎn)速下楔形高度對承載力的影響

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

圖13 不同轉(zhuǎn)速下楔形高度對摩擦功耗的影響

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

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

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

3.5 最小氣膜厚度對軸承靜態(tài)性能的影響

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

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

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

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

4 結(jié)論

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

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

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

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