氣體箔片軸承靜態工作點求解

鄧志凱，程文杰，曹廣東，肖玲，李明

(西安科技大學 理學院，西安 710054)

氣體箔片軸承(Gas Foil Bearings，GFB)作為一種無油軸承，具有摩擦功耗低，工作溫度范圍寬，耐沖擊和裝配對中要求低等優點。GFB已經成功應用于極端環境(dmn值超過3.0×106mm·r/min和538 ℃的服役環境)的氣體透平機中， 并且在燃料電池、發電機、飛機推進器、氣體處理領域有廣泛的應用前景[1-3]。然而，GFB的設計需考慮多種工況，同時涉及多學科交叉的知識，其靜動態性能及性能極限的計算往往基于剛性表面的假定，該假定限制了流體動力潤滑理論在一定工況下的應用。

國內外諸多學者致力于建立更完善的模型來計算軸承的性能，以減少剛性表面軸承的假定對計算結果造成的誤差。文獻[4]在略去波箔之間的相互作用及庫倫摩擦的基礎上，采用柔度系數描述底層支承拱箔彈性變形對氣膜厚度的影響，為降低氣彈耦合計算的復雜性，對頂層箔片的變形未深入考慮。由于箔片結構在氣膜壓力作用下產生彎曲變形，軸承間隙隨之改變，氣膜厚度也相應增大，因此，氣膜與箔片結構之間存在強耦合作用[5]。文獻[6]采用二維Navier-Stokes方程求解流場效應，結合箔片結構的有限元法對流體-箔片-軸頸之間的相互作用進行了穩態、準瞬態和瞬態動力學分析，在分析中得到了耦合軸頸運動、流體、箔片結構的動態性能，并與剛性表面軸承的動態性能進行了比較。文獻[7]通過在有限元公式中耦合氣體流動和箔片的變形來分析軸承的靜態性能，研究發現氣膜出口區的膜厚變化是由于頂層箔片的脫離所引起。文獻[8-9]采用有限差分法以更簡單的方式研究了類似的流固耦合分析。諸多研究表明,對氣體箔片軸承的完整分析需要考慮庫倫摩擦的影響。文獻[10-11]的研究指出拱箔與頂層箔片之間的微滑動產生了摩擦阻尼，該研究中提出的半經驗模型分析了單個拱箔的摩擦阻尼及其與相鄰拱箔的相互作用。

以上對氣體箔片軸承的研究考慮了箔片結構與氣膜之間的耦合作用，使軸承靜動態性能更接近實際值。然而，關于軸承靜態工作點的確定卻鮮有文獻提及，已有文獻中給出的軸承性能多是基于給定偏心率(直接偏心法)來計算軸承的氣體流量、摩擦損耗和承載力。為確定任意給定載荷下氣體箔片軸承的靜態工作點，本文提出了一種“二分法搜索+不動點迭代”的求解策略，其中，二分法用于搜索偏心率，不動點迭代用于尋找偏位角。彈流耦合求解中雷諾方程采用基于質量流量守恒的差分法離散進行超松弛迭代求解，頂層箔片采用Kirchhoff薄板模型進行有限元求解，以整周式GFB及三瓦插入式GFB為算例進行求解及分析。

1 氣體箔片軸承模型及假定

1.1 氣體箔片軸承模型

本文重點研究彈流耦合過程對軸承靜態工作點的影響，由于頂層箔片變形會影響氣膜厚度分布，為方便計算，采用單一的厚箔片替代完整的箔片結構。從軸承承載機理上分析，箔片結構的總剛度矩陣K=Kt+Kb，其中Kt為頂層箔片剛度矩陣，Kb為波箔剛度矩陣，當采用厚箔片時，頂層箔片剛度矩陣可近似替代箔片結構的總剛度矩陣。GFB的結構形式如圖1所示，圖中O，O′分別為軸承中心與軸頸中心，Ω為轉子工作轉速，對于三瓦插入式GFB，定義α為每塊瓦的瓦張角，β為瓦位角(由軸承上方垂線到第1塊軸瓦進氣邊的位置角)，ξ為槽寬包角。

1.2 計算假定

在確定軸承靜態工作點的過程中，為減少耦合分析中每次迭代的時長，求解靜態氣膜壓力場及頂層箔片變形時基于以下假定：

1)氣體潤滑過程是恒溫的；

2)流場是定常的，軸承端部及開槽處的氣膜壓力為標準大氣壓；

3)不考慮頂層箔片的曲率效應。

2 靜態氣膜壓力的求解

2.1 求解理論

以整周式GFB為例，無論是靜態氣膜壓力還是動態氣膜壓力的求解，氣膜厚度h都是氣膜壓力的決定性參數，其表達式為

h=C0+ecos(φ-θ)，

(1)

式中：C0為軸承間隙；e為偏心距；φ為軸承展開角；θ為偏位角。

在此基礎上，靜態氣膜壓力的求解可歸結為定常二維雷諾方程[12]的求解

(2)

式中：p為氣膜壓力；x,y,z分別為軸承水平方向、豎直方向、軸向方向的坐標分量；μ為氣體動力黏度；U為轉子表面沿x方向上的速度分量。橢圓型偏微分方程(2)式無法進行解析求解，大多數學者直接從雷諾方程出發，采用有限差分法構造每一項的差分格式后代入原方程迭代求解。本文基于質量流量守恒原理[2,12]建立包含所有節點壓力在內的差分方程，該方法在處理雷諾方程并不成立的特殊邊界節點時依然具有統一的表達形式。

將軸承沿周向展開并對周向和軸向進行網格劃分，如圖2所示，任意節點(i,j)區域所滿足的守恒關系為

圖2 任意節點處的質量流量守恒關系Fig.2 Conservation of mass flow at any node

Qa+Qb-Qc-Qd-Qe-Qf+Qg+Qh=0，

(3)

式中：Q為質量流量。

根據文獻[12]的推導，單位時間通過任一單位寬度周向和軸向截面的質量流量為

(4)

式中：P為量綱一的氣膜壓力；H為量綱一的氣膜厚度；Λ為軸承數；ψ，λ分別為軸承量綱一的周向長度和軸向長度；Patm為標準大氣壓；R為軸承半徑。

b～h點的推導與a點類似, 以a點為例，其質量流量、氣膜厚度、氣膜壓力以及偏導數的半步長差分格式為

(5)

式中：Δψ，Δλ分別為軸承周向、軸向微段長度；i，j分別為軸承周向、軸向網格節點；P0為氣膜壓力的迭代初始預設值；Pa0為節點a處的氣膜壓力初始值。

將(5)式代入(4)式，即可得到任意節點(i,j)處的氣膜壓力差分方程。由于軸承端部及軸瓦進出口壓力取環境壓力，(4)式的邊界條件為

P(ψ,0)=P(ψ,L)=1，

(6)

式中：L為軸承長度。

對已構造差分格式的所有節點，給出邊界條件后可采用超松弛迭代法進行差分方程組的求解。

氣膜力分量Fx，Fy及氣膜承載力W的計算公式為

(7)

以上關于整周式GFB的求解過程同樣適用于三瓦插入式GFB靜態氣膜壓力的計算，需要注意的是，后者在除軸承端部之外的開槽處也取環境壓力，其邊界條件在滿足(6)式的基礎上還要滿足

P(ψm,λ)=1;m=1,2,3，

(8)

此外，三瓦插入式GFB的W滿足疊加原則,即

(9)

2.2 GFB氣膜壓力求解算例

為驗證氣膜壓力計算程序的可靠性，以整周式剛性表面軸承為算例，首先給出本文計算與文獻[2]計算的氣膜壓力分布對比(圖3)，其中，軸承結構參數及其他性能參數見表1。仿真結果顯示，本文計算得到的量綱一的承載力為3.73，與文獻[2]的計算結果差值在5%以內。

圖3 整周式GFB靜態氣膜壓力分布對比Fig.3 Comparison of static gas film pressure distribution of single-pad GFB

表1 整周式GFB結構參數及其他性能參數Tab.1 Structure parameters of single-pad GFB and other performance parameters

為比較三瓦插入式GFB與整周式GFB的氣膜壓力分布特點，剛性表面假設下整周式GFB及三瓦插入式GFB的靜態氣膜壓力分布分別如圖4、圖5所示，其中，三瓦插入式GFB的瓦張角α=114°，瓦位角β=38°，槽寬包角ξ=18°，其余參數與表1相同。與整周式GFB相比，三瓦插入式GFB在軸向方向上的氣膜壓力分布更均勻。2種GFB在軸承中截面(λ=1.5)上的氣膜厚度和氣膜壓力分布分別如圖6、圖7所示。不難發現，三瓦插入式GFB中軸向槽的引入使氣膜厚度沿軸承周向分布不再是連續函數，同時，在軸承周向方向上出現3個氣膜壓力峰，與整周式GFB相比，三瓦插入式GFB氣膜壓力沿軸承周向分布較為均勻，但氣膜承載力有所下降。

圖4 整周式GFB靜態氣膜壓力分布Fig.4 Static gas film pressure distribution of single-pad GFB

圖5 三瓦插入式GFB靜態氣膜壓力分布Fig.5 Static gas film pressure distribution of three-pad insertion GFB

圖6 GFB氣膜厚度分布Fig.6 Gas film thickness distribution of GFB

圖7 GFB氣膜壓力分布Fig.7 Gas film pressure distribution of GFB

3 頂層箔片的彎曲變形

3.1 Kirchhoff板理論

在GFB的設計中，頂層箔片的結構形狀和宏觀變形是箔片軸承靜動態性能的主導性因素，就無預緊箔片軸承而言，頂層箔片所具有的間隙接近于名義氣隙的平均值[2]，因此，本文針對無預緊軸承進行研究。由于頂層箔片的厚度相比其長、寬尺寸小了幾個數量級，頂層箔片的建模可以歸結為橫向復雜載荷下的薄板彎曲。

采用有限差分法計算薄板的彎曲變形存在收斂條件較高，復雜載荷下難以收斂的問題，相反，采用有限單元法則極大地縮減了求解時間，這對耦合求解氣膜壓力場及箔片變形來說無疑加快了每一步的迭代過程。因此，本文采用有限元法自主編程求解頂層箔片的彎曲變形。

在已有關于彈性板的彎曲理論中，Kirchhoff板理論更適用于頂層箔片的力學建模，由于忽略了板的橫向剪切變形，板內所有的力學量都可用板中面的撓度ω(x,y)表示。頂層箔片簡化后的Kirchhoff薄板模型如圖8所示，該薄板模型中左右邊界分別采用固支和簡支。板中每個節點有3個廣義位移分量(撓度ω，繞x軸的轉角θx和繞y軸的轉角θy)。在此基礎上，由四節點矩形單元位移模式下的形函數得到薄板的單元剛度矩陣

圖8 Kirchhoff薄板模型Fig.8 Model of thin plate in Kirchhoff

(10)

式中：B為幾何矩陣；D為彈性矩陣；N為形函數，N=[N1N2N3N4]，其子矩陣N1～N4推導較為復雜，具體表達式可參考文獻[13]。

將各單元的剛度矩陣進行組裝，可建立薄板結構的整體剛度矩陣。若薄板單元受橫向氣膜力的作用，則等效結點力列陣為

(11)

式中：Fzi，Mθxi，Mθyi分別為作用在節點i上的節點力、繞x軸和繞y軸的節點彎矩。各節點的位移列陣為

(12)

式中：ωi為薄板節點的撓度；θxi，θyi分別為節點繞x軸和繞y軸的轉角。

3.2 箔片變形求解算例

圖9—圖11給出了軸承數Λ=1，長徑比L/D=1.5，偏心率ε=0.7，偏位角θ=38.08°下2種GFB的頂層箔片變形。其中，三瓦插入式瓦張角α=114°,瓦位角β=38°，箔片結構參數見表2。由圖9、圖10可知，盡管整周式GFB的箔片厚度是三瓦插入式GFB的2倍還多，但其徑向變形依然比后者大出一個數量級,原因為三瓦插入式GFB的多槽結構使氣膜沿周向分布相對均勻且總承載力有所降低,每塊瓦(頂層箔片)的長寬比更接近于1，若將圖8所示薄板模型簡化為二維歐拉梁模型，在橫截面不變的情況下，頂層箔片變形將與其板長呈三次非線性關系。因此，在箔片參數相同的情況下，相對整周式GFB，三瓦插入式GFB每塊瓦的周向長度更短，頂層箔片變形更小，其變形后的承載輪廓面更接近于剛性軸承。由于氣膜壓力分布的均勻化，三瓦插入式GFB更有利于轉子的旋轉穩定性。在三瓦插入式GFB中，當某塊箔片位于膜厚最小區域時，其徑向變形約是其余箔片的5倍。

圖9 整周式GFB氣膜壓力下頂層箔片的變形Fig.9 Deformation of top foil under gas film pressure for single-pad GFB

圖10 三瓦插入式GFB頂層箔片的徑向變形Fig.10 Radial deformation of top foil under gas film pressure for three-pad insertion GFB

圖11 三瓦插入式GFB軸向方向膜厚最小處頂層箔片的徑向變形Fig.11 Radial deformation of top foil at minimum gas film thickness in axial direction of three-pad insertion GFB

表2 頂層箔片的結構參數Tab.2 Structural parameters of top foil

由圖11可知，位于軸承端部的頂層箔片徑向變形較小，與文獻[2,14]得出的結果相似，這種變形特征減少了氣體在軸承端部的泄漏。

4 靜態工作點求解

4.1 二分法搜索確定偏心率

對于任意給定的軸承參數及外載荷，偏心率對氣膜承載力的影響遠大于偏位角，因此，先采用二分法確定偏心率，在此之前，需預設初始偏位角θ0。

對于給定的偏心率搜索區間(ε1,ε2)，為判斷預設區間的合理性，需要分別計算(ε1,θ0)，(ε2,θ0)下的氣膜承載力與外載荷的差值W1，W2。若W1W2>0，則所求的偏心率不在搜索區間范圍內，反之，則令εx=(ε1+ε2)/2，開始在所給區間內由(εx,θ0)計算氣膜承載力與外載荷的差值Wx。為保證解的精度，需給出偏心率的允許誤差Δ，總迭代次數Nt要滿足

(13)

在搜索過程中，當W1Wx>0時，則令εx=ε1，反之，令εx=ε2，直至迭代次數n>Nt時迭代終止，此時，ε即為要找的偏心率。

4.2 不動點迭代尋找偏位角

對于任意給定的外載荷角αf，靜態偏位角與αf差值δ不能過大，需要滿足收斂條件

(14)

式中：Fx，Fy分別為(ε,θn)下的氣膜力在水平與豎直方向的分量；γ為偏位角的允許誤差，一般取γ<1.0×10-4。預設偏位角后，新的偏位角可按照如下迭代格式進行尋找

(15)

當δ<γ時，迭代終止。至此，偏心率、偏位角都已確定。

4.3 求解策略

在應用“二分法搜索+不動點迭代”求解軸承靜態工作點的過程中，氣膜壓力采用有限差分求解，薄板變形運用有限元法求解，其中網格大小相同，單元及結點也完全相同，于是，相同結點的氣膜力及薄板撓度的傳遞在程序上容易實現。此外，在迭代過程中，需要不斷用箔片變形更新氣膜厚度，具體求解步驟如下：

1)預設初始偏位角，采用二分法迭代尋找偏心率，給定偏心率的迭代區間；

2)由當前迭代步中的偏心率、偏位角計算靜態氣膜壓力場中的氣膜厚度h1，對壓力場積分獲得氣膜合力F1，將氣膜壓力場代入箔片變形場中，計算得到箔片變形ω1；

3)將第2步中得到的箔片變形ω1代入靜態氣膜壓力場中，修正氣膜厚度得到新膜厚h2，重新計算氣膜壓力場并積分獲得氣膜合力F2，此時將氣膜壓力場帶入箔片變形場得到箔片變形ω2；

4)重復第2步和第3步直到氣膜合力Fn收斂；

5)判斷氣膜合力Fn與外載荷是否滿足偏心率的收斂條件(滿足二分法根的允許誤差的最大迭代次數)，若滿足，進入偏位角的迭代計算，若不滿足，重復第2步到第5步；

6)在確定偏心率的基礎上，采用不動點迭代尋找偏位角，由當前迭代步中的偏位角和已確定的偏心率計算氣膜壓力場及頂層箔片變形，迭代過程與第2步至第4步一致；

7)當外載荷角αf與偏位角θn滿足偏位角的收斂條件(|αf-θn|≤10-4)時迭代終止，若不滿足，重復上一步的迭代。

4.4 GFB靜態工作點的求解算例

由于三瓦插入式GFB的求解流程與整周式GFB完全相同，因此，本文以具有代表性的整周式GFB為算例進行靜態工作點的求解。表3給出了整周式剛性表面GFB在給定量綱一外載荷Wf下的靜態偏心率和偏位角，其中Λ=0.6，L/D=1.0。

由表3可知：本文靜態工作點求解方法得到的靜態偏心率與文獻[2]的相對誤差最大為6.60%，隨著偏心率的增大，相對誤差顯著減小，且小偏心率的相對誤差并不會造成過大的軸承承載力誤差；偏位角與文獻[2]的相對誤差都小于3.80%。上述計算結果表明本文所采用方法的可靠性。需要說明的是，文獻[2]的承載力是給出轉子偏心率與偏位角后計算得到的，本文則是根據文獻計算的承載力來反向搜索軸承的靜態偏心率與偏位角。實際上，更關心某個給定質量的轉子運行穩定時的工作位置，按本文的算法，任意給定一個轉子及外載荷，都可以很快確定它的靜態工作點，這也是本文求解靜態工作點的目的。

表3 整周式剛性表面GFB靜態偏心率與偏位角Tab.3 Static eccentricity and attitude angle of single-pad rigid surface GFB

相同載荷下剛性表面軸承與柔性表面軸承的靜態偏心率和偏位角分別如圖12—圖14所示，其中軸承的結構參數見表4。

圖12 剛性表面與柔性表面下軸承偏心率對比Fig.12 Comparison of eccentricity between rigid surface bearing and flexible surface bearing

圖14 極坐標系下剛性、柔性表面軸承偏位角、偏心率的對比Fig.14 Comparison of attitude angle and eccentricity in polar coordinate system between rigid surface bearing and flexible surface bearing

表4 軸承的結構參數Tab.4 Structural parameters of bearing

由圖12可知，隨著外載荷的增加，柔性表面軸承相對剛性表面軸承的靜態偏心率不斷增大，在載荷最大處其偏心率差值約為0.05，在重載情況下，即使偏心率差值較小，也會造成較大的承載力誤差。實際上，若箔片表現出更軟的特性，其變形越大,剛性表面軸承的靜態偏心率與實際差值也將越大。比較圖12、圖13可知，無論是剛性表面軸承還是柔性表面軸承，相同載荷條件下，偏位角隨著偏心率的增加都會逐漸減小，不同的是，隨著偏心率的增加，柔性表面軸承的偏位角遞減的更快一些。

圖13 剛性表面與柔性表面下軸承偏位角對比Fig.13 Comparison of attitude angle between rigid surface bearing and flexible surface bearing

圖14給出了更直觀的靜態工作點“軌跡圖”，不難發現，柔性表面軸承的靜態工作點軌跡表現出了更明顯的下沉，該現象與箔片變形產生的偏轉有關。需要注意的是，(2)式等式右端的項表示由氣膜旋轉效應所貢獻的氣膜承載力，本文采用的定常流體潤滑假設忽略了氣膜擠壓效應下的氣膜承載力，若計入這種效應的影響，箔片將會產生更大的變形，此時，剛性表面軸承與柔性表面軸承的偏位角差值將會更大。因此，在GFB的動態性能計算中，更完善的模型應考慮雷諾方程中關于時間的微分項。

圖15—圖17分別給出了Wf=0.75,ε=0.62，θ=61.90°時軸承間隙最小處的量綱一的膜厚、量綱一的軸承承載力、箔片變形隨迭代次數的變化曲線。其中，首次迭代計算出的氣膜厚度為剛性表面軸承的膜厚h1；二次迭代時則引入首次計算的氣膜合力F1下的箔片變形ω1，由圖可知，最小膜厚經過箔片變形的修正會顯著變大為h2(h2>h1)，此時氣膜合力減小為F2(F2h3>h1)；重復以上過程，膜厚和承載力最終收斂于定值。實際上，上述流體-箔片的耦合分析中，也可理解為氣膜和箔片同時經歷了振動的瞬態和穩態過程。

圖15 軸承量綱一的氣膜厚度隨迭代次數的變化Fig.15 Variation of dimensionless gas film thickness for bearing with iterations

圖16 軸承量綱一的承載力隨迭代次數的變化Fig.16 Variation of dimensionless load capacity of bearing with iterations

圖17 頂層箔片量綱一的變形隨迭代次數的變化Fig.17 Variation of dimensionless deformation for top foil with iterations

另外，以上計算結果顯示，氣體箔片軸承采用剛性表面假設會帶來不小的誤差，尤其是在大偏心率的情況下計算得到的結果很可能失真。由于本文算例中所采用的箔片尺寸相對較厚，箔片變形很小，在工程實際中箔片的厚度只有零點幾毫米，往往會出現偏心率大于1的情況，因此，計入箔片的彎曲變形對氣體箔片軸承的靜、動態性能計算非常必要。

5 結論

通過“二分法搜索+不動點迭代”的求解策略尋找氣體箔片軸承的靜態工作點，得出以下結論：

1)采用本文設計的算法可確定任意給定外載荷下軸承的靜態工作點，剛性表面軸承的偏心率與偏位角的計算結果與文獻[2]的計算結果吻合較好(偏心率相對誤差不超過6.60%，偏位角相對誤差小于3.80%)。

2)無論是剛性表面軸承還是柔性表面軸承，相同載荷條件下，偏位角隨著偏心率的增加都會逐漸減小，不同的是，相對于剛性表面軸承，柔性表面軸承偏位角遞減的更快一些。

3)在計算氣體箔片軸承時，采用剛性表面假設會帶來不小的誤差，尤其是在大偏心率的情況下計算得到的結果很可能失真。考慮到工程中往往會出現偏心率大于1的情形，計入箔片的彎曲變形來計算軸承的靜動態性能更加合理。