連軍偉,張茂森,金志磊,賀雷

(1.中國航天科技集團有限公司 低溫液體推進技術實驗室,北京 100076;2.北京航天動力研究所,北京 100076)

1 概述

人類探測太空的能力主要受空間推進技術的限制,針對高速度增量的航天任務,傳統(tǒng)化學推進技術的比沖較低,太陽能電池陣—蓄電池技術嚴重依賴太陽光照,深空生存能力較差[1]?？臻g核電源系統(tǒng)具有能量密度大,容易實現(xiàn)大功率(數(shù)千瓦至數(shù)兆瓦)供電,機動性和隱蔽性好,不依賴太陽光照等優(yōu)勢,是深空探測不可替代的空間電源[1]。

動態(tài)轉(zhuǎn)換的放射性同位素發(fā)電系統(tǒng)在空間核電源系統(tǒng)中轉(zhuǎn)換效率最高,其中研究較多且已應用的是閉式布雷頓循環(huán)的動力轉(zhuǎn)換系統(tǒng)(圖1),旋轉(zhuǎn)單元是該系統(tǒng)的核心組件,其結(jié)構和性能直接影響系統(tǒng)的優(yōu)劣。止推箔片軸承作為一種柔性自適應動壓軸承,在旋轉(zhuǎn)單元中承受軸向載荷,與滾動軸承相比,省略了復雜的冷卻潤滑系統(tǒng),使系統(tǒng)更加清潔和高效。

(a) 閉式布雷頓循環(huán)系統(tǒng)原理

1854年,法國科學家提出將氣體作為潤滑劑的可能,1886年法國科學家雷諾推導出潤滑膜內(nèi)壓力分布的雷諾方程[2],對流體潤滑的認識上升到理論高度。止推箔片軸承作為階梯型動壓軸承的衍生結(jié)構軸承,具有波箔剛度非線性的特點,其理論研究晚于徑向箔片軸承且相關研究較少。1983年,文獻[3]對止推箔片軸承的靜態(tài)特性進行了詳細的計算分析,計算模型考慮了波箔剛度的非線性分布,第1次使用柔度系數(shù)α分析箔片軸承靜特性,該方法一直沿用至今。2000年,文獻[4]采用有限差分法與有限元相結(jié)合的方法,將可壓縮流體的雷諾方程求解與彈性變形耦合起來進行止推箔片軸承分析。2017年,文獻[5]采用非等溫雷諾方程評估了轉(zhuǎn)子靜態(tài)和動態(tài)角度不對中對止推箔片軸承性能的影響,并深入分析了軸承的承載力、摩擦功、氣膜剛度和阻尼系數(shù)的變化情況。

20世紀80年代,文獻[6]應用圓柱薄殼線性理論推導出波箔變形方程,這是國內(nèi)最早關于箔片軸承中波箔剛度的計算模型。文獻[7]采用有限單元法對平箔式止推箔片軸承進行求解,并分析了箔片軸承靜態(tài)特性的影響因素。文獻[8]對止推箔片軸承進行了啟停試驗,分析了箔片材料和表面涂層對軸承啟停磨損的影響。文獻[9]采用有限差分法得出,推力盤傾斜狀態(tài)下軸承的軸向承載力和黏性摩擦力矩均隨傾角的增大而增大,且傾角越大靜態(tài)特性越大。文獻[10]基于牛頓迭代法建立偏角不對中時止推箔片軸承的潤滑模型,與線性替換法的比較說明牛頓迭代法具有同樣的適用性。文獻[11]基于CFX的有限元計算與MATLAB的有限差分求解動壓氣體軸承氣膜壓力,二者計算結(jié)果一致,隨著偏心率的增加CFX有限元的計算結(jié)果比有限差分法計算的略小。文獻[12]采用ANSYS和CFX軟件計算了超臨界二氧化碳止推箔片軸承動態(tài)特性,得出動態(tài)剛度和阻尼大小主要取決于氣膜厚度和箔片結(jié)構,且隨氣膜厚度的增大而減小。隨著更多類型箔片軸承的研究,箔片軸承的性能逐漸被熟知,但仍缺乏精準的理論分析模型和實用的設計規(guī)范[13]。

綜上可知,目前在止推箔片軸承的誤差研究方面,主要著重于軸承不對中誤差對其靜態(tài)性能的影響,而錐度誤差對軸承承載特性的影響以及不對中誤差下軸承動態(tài)特性變化的研究較少。在生產(chǎn)和裝配過程中,錐度和不對中誤差均難以避免,不對中誤差主要有平行不對中、角不對中和綜合不對中。由于旋轉(zhuǎn)單元轉(zhuǎn)子系統(tǒng)的徑向設計間隙約束及推力盤徑向尺寸大于止推軸承高度,分析過程不考慮平行不對中的影響,本文所提及的不對中特指角不對中或主軸傾斜。本文在柱坐標系下建立止推箔片軸承的雷諾方程,并采用牛頓迭代法和有限差分法對雷諾方程和氣膜厚度方程進行耦合求解,分析生產(chǎn)和裝配過程引起的徑向錐度誤差和不對中誤差對止推箔片軸承承載特性的影響。

2 理論分析

止推箔片軸承主要由頂箔、波箔和底板組成。波箔和頂箔組成多個扇形瓦,每個瓦的度數(shù)和傾斜面占比相同。止推箔片軸承結(jié)構和箔片參數(shù)示意圖如圖2所示,β為扇形瓦的瓦角,βp為扇形瓦節(jié)線位置角,σh為扇形瓦傾斜面高度,ω為推力盤角速度,h2為設計間隙。

(a) 止推箔片軸承 (b) 箔片參數(shù)

2.1 壓力控制雷諾方程推導

推導等溫條件下的雷諾方程時進行如下假設:1)忽略氣體慣性力;2)潤滑氣體介質(zhì)視為牛頓流體,符合牛頓黏性定律;3)氣膜間隙內(nèi)為層流;4)忽略氣膜法向流動;5)潤滑介質(zhì)視為理想氣體。

經(jīng)推導得到柱坐標下潤滑介質(zhì)為理想氣體的雷諾方程,即

(1)

式中:r為徑向坐標;h為氣膜厚度;P為氣膜壓力;θ為周向坐標;μ為氣體動力黏度。

(2)

式中:Pa為環(huán)境壓力;R2為軸承外半徑;Λ為軸承數(shù)。

2.2 氣膜厚度方程

每個扇形瓦的氣膜由傾斜面氣膜和水平面氣膜組成,根據(jù)圖2的結(jié)構參數(shù)得到單個扇形瓦的氣膜厚度方程,即

h=h2+h(θ)+w(r,θ) ,

(3)

式中:w(r,θ)為箔片變形后構建的變形區(qū)域氣膜厚度;h(θ)為扇形瓦傾斜面氣膜厚度。

定義b為斜面占比,b=βP/β,氣膜厚度的表達式為

(4)

在工程應用過程中存在加工誤差,推力盤不可避免地產(chǎn)生一定的徑向錐度βz,將直接影響氣膜厚度分布。錐度為負數(shù)時表示推力盤出現(xiàn)內(nèi)錐狀態(tài),為正數(shù)時出現(xiàn)外錐狀態(tài), 按內(nèi)錐度和外錐度這2種形式構建徑向錐度模型(圖3),Ri為錐度半徑,R1為軸承內(nèi)半徑,Δh為錐度引起的氣膜厚度變化。錐度βz的定義為

(a) 內(nèi)錐度

(5)

建立考慮推力盤徑向錐度的氣膜厚度方程,即

h0=h2+h(θ)+w(r,θ)+Δh。

(6)

盡管設計及生產(chǎn)過程要嚴格控制形位公差,但在裝配過程仍會存在不對中現(xiàn)象[14-15],而不對中誤差會導致止推箔片軸承氣膜間隙發(fā)生變化,從而引起軸承承載特性的改變。如圖4所示,建立不對中模型和空間坐標系。

(a) 坐標系及扇形瓦標記

假設轉(zhuǎn)子與止推箔片軸承的夾角(不對中角度)為φ,在空間坐標下可分解為繞x軸的θx和繞y軸的θy?？紤]轉(zhuǎn)子不對中的氣膜厚度方程為

h1=h2+h(θ)+w(r,θ)-rcosθtanφ。

(7)

2.3 波箔剛度等效模型

波箔單位面積剛度Kb的表達式為

(8)

式中:Eb為波箔的彈性模量;t為波箔厚度;l為半波箔長度;s為波箔單位長度;ν為波箔的泊松比。

2.4 網(wǎng)格劃分及邊界條件

圖5 單個扇形瓦的網(wǎng)格劃分

2.5 流固耦合計算流程

采用有限差分法求解(2)式得到氣膜壓力,通過迭代由(9)式判斷是否收斂?；贛ATLAB平臺按圖6程序邏輯進行編譯計算。

圖6 止推箔片軸承計算程序邏輯圖

(9)

3 計算與分析

3.1 程序驗證

(10)

式中:W為承載力。

本文程序計算與文獻[3]的對比結(jié)果如圖7所示,最大誤差為2.53%,兩者吻合度較好,說明本文止推箔片軸承數(shù)學模型的準確性較高。

圖7 本文程序與文獻[3]計算結(jié)果對比

3.2 錐度誤差的影響

空間核電動力轉(zhuǎn)換系統(tǒng)的止推箔片軸承結(jié)構參數(shù)及其他參數(shù)見表1。

表1 止推箔片軸承結(jié)構參數(shù)及其他參數(shù)

錐度一定時,各扇形瓦氣膜厚度和氣膜壓力在同一半徑和周向角度上相同,為提高計算精度并減少計算時間, 選取單個扇形瓦進行分析。 在額定設計間隙下,內(nèi)錐(Ri=R1)和外錐(Ri=R2)狀態(tài)下量綱一的氣膜壓力和氣膜厚度分布如圖8和圖9所示:當內(nèi)錐度為-4×10-4時,最大氣膜壓力增大18.4%,對應軸向承載力提高47.8%;當外錐度為4×10-4時,最大氣膜壓力減小12.7%,對應軸向承載力降低33.8%。以圖9b中氣膜厚度A和B位置為基準,內(nèi)錐度為-4×10-4時B區(qū)域的氣膜厚度減小至10.6 μm,外錐度為4×10-4時A區(qū)域的氣膜厚度增大至29.4 μm。綜上可知:內(nèi)錐度減小了平均氣膜厚度,提高了承載力;外錐度增大了平均氣膜厚度,降低了承載力。

(a) -4×10-4

(a) -4×10-4

設計間隙不變時軸向承載力和摩擦功隨錐度的變化以及軸向載荷恒定(200 N)時最小氣膜厚度和量綱一的最大氣膜壓力隨錐度的變化如圖10所示:設計間隙不變時,軸向承載力和摩擦功隨著錐度的增加逐漸減小,由于錐度的增加增大了平均氣膜厚度,削弱了動壓效應進而減小了軸向承載力和摩擦功。在軸向載荷恒定時,隨著錐度絕對值的增加,最小氣膜厚度逐漸減小, 最小氣膜厚度達到0.67 μm,若繼續(xù)減小至零時程序計算很難收斂,由于模型未考慮表面粗糙度, 工程應用時可認為此時軸承已處于臨界碰磨狀態(tài)。說明恒定載荷工況時只有允許一定范圍的錐度誤差才能保證軸承正常運行。最大氣膜壓力與最小氣膜厚度正好相反,說明軸向載荷恒定時只要錐度絕對值增加就會增大最大氣膜壓力。此外,內(nèi)錐度和外錐度引起的最大氣膜壓力和最小氣膜厚度均關于錐度βz=0對稱。

圖10 止推箔片軸承靜態(tài)參數(shù)隨錐度的變化

3.3 不對中的影響

為簡化分析模型,掌握不對中對軸承性能的影響,只選取繞某一軸傾斜產(chǎn)生的不對中進行分析。對中和不對中時止推箔片軸承的氣膜壓力、氣膜厚度和箔片變形如圖11所示:由于不對中的產(chǎn)生,軸承各瓦的氣膜壓力和氣膜厚度分布不均勻,箔片產(chǎn)生的變形量也不同;2#瓦最大氣膜壓力(2.458)最大,高于對中狀態(tài)下各瓦的最大氣膜壓力(1.805),5#瓦的最大氣膜壓力最小(1.494);由于氣膜壓力作用,最大變形量為(1.205 μm)同樣出現(xiàn)在2#瓦上,是5#瓦最大變形量(0.412 μm)的2.925倍;由于箔片變形構建了新的氣膜承載區(qū)域,氣膜壓力越高時箔片變形量越大,體現(xiàn)了止推箔片軸承較好的自適應性。

(a) 不對中狀態(tài)(θx=0.01°)

止推箔片軸承和各瓦的承載力和摩擦功隨不對中角度的變化如圖12所示:軸承的承載力和摩擦功隨著不對中角度的增大而增加,但各瓦不同;當推力盤按圖4坐標系傾斜時,x軸上方的3個扇形瓦平均氣膜厚度減小,動壓效應增強,氣膜壓力升高,承載力和摩擦功隨不對中角度的增大而增加,2#瓦平均氣膜厚度最小,其承載力最大且增速最快,1#瓦次之,3#瓦增速最慢,x軸下方3個扇形瓦則相反;不對中造成各瓦的承載力和摩擦發(fā)熱嚴重不均衡,不對中角度越大時各瓦的差異越大。

(a) 軸承的承載力和摩擦功

由于不對中的存在,止推箔片軸承不僅產(chǎn)生沿軸向的位移剛度,還存在沿徑向的彎曲剛度。圖4中軸承產(chǎn)生繞x軸的彎矩Mx; 由于1#, 6#瓦與3#, 4#瓦承載力的差異及2#, 5#瓦關于y軸不對稱,產(chǎn)生繞y軸的彎矩My,則軸承在3個自由度上有9個剛度系數(shù)和9個阻尼系數(shù),組成的剛度系數(shù)矩陣K和阻尼系數(shù)矩陣C分別為

(11)

(12)

剛度系數(shù)和阻尼系數(shù)隨不對中角度的變化如圖13所示:當繞x軸傾斜時,軸向主剛度系數(shù)和主阻尼系數(shù)比軸向交叉項的大,隨著不對中角度的增大,軸向主剛度系數(shù)和主阻尼系數(shù)變化較明顯;在繞x軸方向,交叉項Kzθx和Czθx隨不對中角度的變化較明顯,且幅值遠高于其他項,主要原因為不對中角度增加,承載力快速增大,造成承載力在該方向的交叉彎曲剛度系數(shù)和阻尼系數(shù)增大;在繞y軸方向,交叉項Czθy隨不對中角度變化較明顯且為負數(shù),剛度系數(shù)主項Kθyθy的變化最大,當θx為0.016°時,主項Kθyθy增幅為28.5%,交叉項Kzθy增幅為22.1%,說明在繞y軸方向,隨著不對中角度的增加,由彎矩My和承載力引起的剛度變化占主導。上述現(xiàn)象說明不對中角度對止推箔片軸承的剛度系數(shù)和阻尼系數(shù)的影響較大,而剛度系數(shù)和阻尼系數(shù)直接影響轉(zhuǎn)子系統(tǒng)的穩(wěn)定性。

4 結(jié)論

基于牛頓迭代法和有限差分法,耦合求解考慮徑向錐度和不對中誤差影響的氣膜厚度方程及壓力控制的雷諾方程,對閉式布雷頓循環(huán)系統(tǒng)的止推箔片軸承進行承載特性影響因素分析,得到如下結(jié)論:

1)在設計間隙一定的情況下,內(nèi)錐度會減小平均氣膜厚度,進而增大軸承承載力,當內(nèi)錐度為-4×10-4時,承載力增大47.8%;外錐度會增大平均氣膜厚度,進而降低承載力,當外錐度為4×10-4時,承載力降幅達到33.8%。在軸向載荷為200 N時,最小氣膜厚度隨著錐度絕對值的增加而減小,當錐度大于2.25×10-3時,軸承進入臨界碰磨狀態(tài)。生產(chǎn)過程應嚴格控制推力盤的錐度,避免因加工誤差過大而造成軸承失效。

2)不對中會改變軸承各瓦的氣膜壓力和氣膜厚度分布,當不對中角度θx=0.01°時,最大氣膜壓力增加36.2%。不對中角度造成各瓦的承載力嚴重不均衡,對軸承剛度系數(shù)和阻尼系數(shù)的影響較明顯,直接影響轉(zhuǎn)子系統(tǒng)的穩(wěn)定性。隨著不對中角度的增大,軸承的摩擦功呈指數(shù)形式增加,導致軸承的發(fā)熱量急速上升,裝配過程應盡可能將不對中誤差降到最低。