基于格子Boltzmann 方法的動壓氣體軸承黏性熱耗散研究*

鹿翔宇 許 波 陳振乾

（東南大學能源與環境學院 江蘇南京 210096）

動壓氣體軸承因其轉速高、 壽命長、 低摩擦、 無污染等優點， 被廣泛應用于透平膨脹機、 渦輪壓氣機、 飛機發動機、 制冷機等高速旋轉渦輪機械中， 在航空航天、 能源與動力等領域發揮著十分重要的作用［1］。 在動壓氣體軸承的運轉過程中， 雖然轉子和箔片不會直接接觸， 但是軸承間隙內流體的高速剪切流動依然會產生大量的熱量， 容易引起溫度的分布不均和局部過熱， 進而燒蝕箔片。 此外， 軸承的零部件也會因溫度的升高而導致剛度和阻尼發生變化， 產生熱變形， 對軸承的性能造成一定的影響［2］。 因此， 研究動壓氣體軸承的氣動熱和熱特性規律， 以及相關的冷卻方法具有重要意義。

國內外學者對于軸承熱特性的研究主要集中在近20 年。 SALEHI 等［3-4］使用基于Couette 近似的簡化能量方程和可壓縮雷諾方程獲得了箔片氣體軸承中的氣膜溫度， 并采用有限差分法和迭代法對簡化表達式進行求解。 PENG 和KHONSARI［5-6］考慮空氣的可壓縮性和黏溫特性以及軸承表面的柔度， 建立了用于預測箔片氣體軸承三維溫度場的熱流體力學模型， 并研究了冷卻流量的影響， 發現膜的溫升隨轉速的增加呈指數增長。 LEE 等［7-8］提出了一個包含箔片結構和轉子的軸承傳熱模型， 計算了潤滑氣膜、 頂箔、 波箔、 軸承套和轉子的溫度， 并采用廣義雷諾方程和三維能量方程預測了氣膜溫度。 SIM 和KIM［9］通過軸的二維軸對稱熱傳導、 氣膜中的三維能量傳輸、 箔片熱阻和軸承殼體中的三維熱傳導模型， 建立了箔片氣體軸承的熱流體動力學模型， 該模型還結合了箔片接觸熱阻和入口流混合的分析模型， 以提高模型精度。 SAMANTA和KHONSARI［10］分析了箔片氣體軸承的熱彈性不穩定問題， 推導出一個封閉的解析解來預測在波箔中發生熱彈性不穩定的臨界速度， 并驗證了波數和膜厚對失穩臨界速度的顯著影響。 ZDZIEBKO 和MARTOWICZ［11］開發了箔片氣體軸承的有限元模型， 分析表明， 波箔與頂箔和軸承套接觸次數最多的區域， 也就是頂箔的熱量傳遞到軸承襯套和其他軸承的部件最集中的區域。

以上對于箔片氣體軸承熱特性的研究， 主要采用的是雷諾方程結合能量守恒方程的宏觀方法。 然而軸承間隙的尺寸一般非常小， 通常在幾十微米的級別［12］， 間隙內的流場已經處于滑移區， 甚至過渡區，會出現一些微觀效應， 例如速度滑移、 溫度躍變、 黏性熱耗散等， 此時采用一些微觀方法將更加合適。 分子動力學（Molecular Dynamics， MD） 模擬和直接模擬蒙特卡羅（Direct Simulation Monte-Carlo， DSMC）方法等基于粒子的方法在微流場模擬方面取得了一定的進展［13-14］。 然而， 對于大多數實際應用來說， 它們的計算量過于龐大。 介觀的格子玻爾茲曼方法（Lattice Boltzmann Method， LBM）［15］自從被提出以來， 已經被廣泛應用于微尺度流動的研究， 其計算量與宏觀方法相當， 比較適合用于文中的研究。 與經典雷諾潤滑方程相比， LBM 可以捕捉更加微觀的現象，邊界條件實施簡單。 并且當軸承間隙內氣流處于過渡區時， 此時不滿足連續性假設， 宏觀雷諾方程不再適用， LBM 將成為十分有潛力的研究方法。 然而鮮有學者將LBM 用于動壓氣體軸承的黏性熱耗散研究。

目前， 格子Boltzmann 方法應用在熱流動研究的主要有3 種模型： 多速度 （Multi-speed， MS） 模型［16-17］， 雙分布函數（Double-Distribution Function，DDF） 模型［18-23］和混合模型［24］。 MS 模型是等溫LB模型的直接延伸， 其利用密度分布函數的高階速度矩來描述能量方程。 DDF 模型利用了2 種不同的分布函數， 一種用于流場， 另一種用于溫度場。 在混合方法中， 流動模擬采用LB 法， 而溫度場的求解則采用傳統的數值方法， 如有限差分法。 1998 年， HE 等［18］基于DDF 模型引入內能分布函數， 第一次提出了包含黏性熱耗散和壓縮功的熱格子Boltzmann 模型， 隨后有學者在其基礎上進行了一定的改進［20，25］。 2007年， GUO 等［21］在單一速度分布函數的基礎上定義了代表總能量的分布函數， 基于DDF 模型提出了總能形式的含有黏性熱耗散和壓縮功的熱LBE 模型。 該方法在演化過程中不涉及復雜的導數項， 計算簡單，數值穩定性高， 文中將基于該模型研究動壓氣體軸承的熱特性。

LBM 雖然在研究熱流動方面已經發展了20 多年， 但是目前相關研究依然是針對一些比較簡單的物理模型， 像軸承間隙這種彎曲變截面的高速剪切微通道鮮有研究。 因此本文作者將LB 方法應用到該領域， 基于GUO 等［21］提出的總能形式的DDF 離散速度模型， 得到了可以應用在貼體網格中的有限差分格式的顯式LB 方法， 并引入Langmuir 滑移邊界， 實現對軸承間隙微通道的模擬， 并研究了不同溫差大小、 偏心率、 轉速以及溫度階躍對于黏性熱耗散分布規律的影響， 從更加微觀的層面探究了軸承間隙黏性熱耗散的機制。

1 數值模型

1.1 有限差分DDF-LBE 模型

文中使用的包含黏性熱耗散和壓縮功的總能形式雙分布函數離散速度模型［21］， 其基于Hermit 多項式展開， 并保留至二階項。 它由密度分布函數和總能分布函數構成， 其中總能包括內能和機械能， 2 個分布函數的演化方程為

式中：fi和hi分別是i方向的密度分布函數和總能分布函數；和是對應的平衡態分布函數；ci是粒子的離散速度；τhf＝1／τh-1／τf， 其中τf和τh分別為動能和內能松弛時間；Zi＝ci·u-u2／2，u是宏觀速度；Fi和qi是與外力有關的2 個項。

對于二維九速（D2Q9） 模型， 2 個平衡態分布函數為

式中：R為氣體常數；T0為參考溫度；ρ為宏觀密度，p0＝ρRT0； 權重系數w0＝4／9，w1＝w2＝w3＝w4＝1／9，w5＝w6＝w7＝w8＝1／36， 離散速度ci取值為

2 個外力項Fi和qi的表達式為

式中：a為質量力， 總能，cv為定容比熱容。

為穩定可靠地求解離散方程（1）、 （2）， 對其進行有限差分格式的離散， 在時間［tn，tn＋1］ 內進行求積， 并對演化方程中的碰撞項做如下處理：

式中：tn＋1＝tn＋Δt； 碰撞算子，，Si＝-τfZiΩf／τhf＋qi。

為了能夠顯式求解上式， 可以引入2 個新的分布函數，

式中：ωf＝Δt／τf，ωh＝Δt／τh。

將新的分布函數代入式（8）、 （9）， 就能得到可應用于貼體網格中的有限差分形式的玻爾茲曼方程：

動能和內能的松弛時間有如下關系：

式中：μ為動力黏度；p是壓力；κ為熱傳導系數；cp是定壓比熱容。

于是普朗特數為Pr＝μcp／κ＝τf／τh。

為了能夠將上述方程應用到曲線坐標系ξ中， 以方程（12） 為例， 對其空間梯度項進行如下離散［26］：

?fi／?ξβ的離散理論上可以采用中心差分格式或者二階迎風插值格式進行求解：

式中： Δξβ是ξβ方向的網格間距。

中心差分格式數值耗散少， 但不穩定， 而迎風插值格式雖相對穩定， 但是卻有較強的數值耗散。 為了解決這些問題， 可以將上述2 種方法組成一種混合的差分格式：

其中， 0≤?≤1 是調節中心差分格式和迎風插值格式權重的系數。 根據文獻的建議， 在后續研究中無特殊說明時取0.1。

為了讓LB 方法可以模擬微尺度流動現象， 需要處理2 個重要問題， 第一個就是解決松弛時間τf和流動Kn數（Kn＝λ／L，λ為分子平均自由程，L為特征長度） 之間的關系， 通過結合動力學理論和LBM的運動黏度表達式， 可得［27］：

式中：H為通道的高度。

由于軸承間隙內的流場壓力和溫度分布不均， 因此對上式進行如下修正：

式中：ρ0為參考密度；T0為參考溫度。

1.2 邊界條件

將LB 方法應用到熱微流動的另一個重要問題就是如何處理速度滑移和溫度階躍邊界條件。 文中采用了Langmuir 滑移模型［28］， 它反映了氣體分子和壁面之間的界面相互作用， 可以表示為

式中： 下標w 代表壁面處的局部值； 下標g 代表靠近壁面處的局部值；αv和αT是與氣體種類和壁面材料屬性有關的調節系數。

為了方便確定這2 個調節系數， 可以將Langmuir滑移模型與Maxwell 滑移模型進行聯立

式中：σv為動量調節系數；σT為熱調節系數；σv和σT在文中沒特殊說明的情況下均取1；γ是比熱容比。

于是

式中：N＝1／δx表示網格數量。

將Langmuir 模型的處理方法應用到非平衡外推邊界中， 最終可以得到適用于當前模型的速度滑移和溫度階躍的邊界格式：

式中：xw是邊界處的節點；xg是與邊界相鄰的節點。

1.3 求解步驟

數值計算程序的大致流程如圖1 所示。

圖1 計算流程Fig.1 Numerical calculation procedure

具體的實現步驟如下：

（1） 通過求解橢圓偏微分方程組生成貼體網格，并計算轉置速度。

（2） 根據物理問題， 初始化網格節點的密度、速度和溫度， 得到初始狀態的密度分布函數和總能分布函數及其平衡態。

（3） 根據Kn數確定τf， 進一步得到ωf， 同理根據Pr數確定τh， 進一步得到ωh。

（4） 密度分布函數的求解： 根據式（14） 計算密度和速度， 由結合式（10） 可得fni， 隨后根據式（19） 計算fni的有限差分， 最后根據式（12） 計算得到。

（5） 總能分布函數的求解： 根據式（14） 計算溫度， 由結合式 （11） 可得hni， 隨后根據式（19）計算hni的有限差分， 最后根據式（13） 計算得到。

（6） 判斷速度和溫度是否滿足收斂條件， 如果不滿足則根據溫度和密度調整Kn數， 并重復步驟（3） — （6）， 如果滿足則輸出結果并停止計算。

2 物理模型

簡化后的動壓氣體軸承的示意圖如圖2 所示。 軸頸和軸承套的半徑分別為r1和r2（r2＝r1＋C0， 其中C0為平均間隙）， 最大間隙為hmax， 最小間隙為hmin。軸頸以速度ω順時針旋轉，e為偏心距， 則偏心率為ε＝e／C0， 軸承的Kn數定義為Kn＝λ／hmin。 軸頸表面的溫度為T1， 軸承內表面的溫度為T2， 兩壁面的溫差ΔT（＝T1-T2） 可以采用無量綱數Eckert 數來表征，其定義為

圖2 動壓氣體徑向軸承簡化示意Fig.2 Schematic of gas-lubricated journal bearing

式中：U為軸頸表面的線速度。

軸承的具體參數和運行工況范圍見表1。

表1 軸承參數和運行工況Table 1 Bearing parameters and operating conditions

3 模型與可靠性驗證

以圓柱熱Couette 流為例， 對文中模型進行了驗證。 其幾何參數和貼體網格劃分如圖3 所示。 內外圓柱半徑比Rβ＝r1／r2＝3／5， 內圓柱以線速度U順時針旋轉， 溫度為T1， 外圓柱靜止， 溫度為T2。 考慮黏性熱耗散， 圓周方向和半徑方向的網格數分別為NX＝188，NY＝20。 圓柱熱Couette 流的量綱一溫度T*的解析解表達式［29］為

圖3 圓柱熱Couette 流示意Fig.3 Schematic of thermal cylindrical Couette flow

式中：r*＝r／r1， 為量綱一半徑。

文中在Pr＝0.7 的情況下， 分別模擬得到了Ec數為1、 5、 10 的3 種徑向溫度分布， 結果見圖4。 可以發現， 在相同網格分辨率下， 數值模擬結果與理論解吻合良好， 最大相對誤差不超過0.8%， 證明模型精度較高。

圖4 不同Ec 數下的溫度分布Fig.4 Temperature distribution at different Ec number

在后續的模擬中， 為了排除網格數量對于軸承間隙黏性熱耗散計算結果的影響， 還需進行網格獨立性驗證。 文中對軸承半徑為8 mm， 間隙尺寸為25 μm，偏心率為0， 轉速為5×104r／min 的動壓氣體軸承間隙熱流動進行模擬，Ec數固定為10。 間隙流域圓周方向和半徑方向的網格數分別為10 053×5、 20 106×10、 30 159×15、 40 212×20、 50 265×25。 5 種網格分辨率下模擬得到的間隙流場中最高的量綱一氣膜溫度如圖5 所示。 可以看出， 隨著網格數量的增加， 黏性熱耗散產生的最高溫度也隨之升高， 但是升高的速率逐漸減小， 說明增加網格數對于提升模型精度的作用逐漸降低。 綜合考慮模型精度和計算資源的消耗， 文中最終選取的網格分辨率為30 159×15。

圖5 網格獨立性驗證Fig.5 Grid independence verification

4 結果與討論

4.1 不同Ec 數的影響

軸承在運轉過程中， 由于散熱條件的不同， 一般軸頸表面和軸承套內表面的溫度是存在差異的。 為了更好地探究黏性熱耗散的變化規律， 假設軸頸和軸承套兩側為定溫壁面， 且存在一定的溫差， 并用Ec數來表示溫差大小。 軸承轉速為5×104r／min， 偏心率保持0.6， 圖6 所示為不考慮兩壁面處溫度階躍時，不同Ec數下間隙內最高氣膜溫度處徑向量綱一溫度T*（＝（T-T2）／（T1-T2）） 的分布。 可以看出， 在Ec數較小時， 兩側溫差較大， 間隙溫度分布幾乎呈線性分布， 此時是以熱傳導為主， 黏性耗散的影響基本可以忽略。 而當Ec數較大時， 黏性熱耗散的作用明顯， 溫度呈拋物線分布， 流場中最高溫度出現在靠近軸頸側， 且流場中大部分的溫度都已經高于T1。 總體來看， 隨著Ec數的增加， 也就是溫差的減小， 黏性耗散作用增強， 氣膜溫度的最大值是逐漸遞增的。

圖6 不考慮溫度階躍時不同Ec 數下最高氣膜溫度處徑向量綱一溫度分布Fig.6 Radial dimensionless temperature distribution at the highest gas film temperature for different Ec number without considering temperature jump

圖7 所示為考慮了溫度階躍情況下的不同Ec數時的溫度分布。 當Ec＝1 時， 軸頸側的溫度小于T1，溫度階躍為正值； 而當Ec＝2 時， 軸頸側的溫度大于T1， 此時溫度階躍變成了負值。 隨著Ec數的增加，兩側的溫度階躍均呈增大的趨勢， 說明溫度階躍的影響將變得更加顯著。 因此， 軸頸和軸承套兩側溫差的增大， 將更有利于削弱黏性熱耗散的影響。

圖7 考慮溫度階躍時不同Ec 數下最高氣膜溫度處徑向量綱一溫度分布Fig.7 Radial dimensionless temperature distribution at the highest gas film temperature for different Ec number considering temperature jump

4.2 不同偏心率的影響

圖8 所示為轉速5×104r／min， 軸承平均間隙為25 μm，Ec數為10， 不考慮溫度階躍時不同偏心率下間隙內最高氣膜溫度處徑向量綱一溫度分布。 可以看出， 隨著偏心率的增加， 間隙中的氣膜最高溫度也逐漸升高。 這是因為， 當偏心率變大時， 軸承間隙尺寸在圓周方向的變化幅度增加， 意味著氣膜的擠壓作用增強， 黏性熱耗散的量也隨之增大。 此外， 出現最大氣膜溫度處的位置隨著偏心率的增加逐漸由靠近軸頸側向流場中間偏移， 表明間隙內的流動逐漸接近于平板Couette 流。

圖8 不考慮溫度階躍時不同偏心率下最高氣膜溫度處徑向量綱一溫度分布Fig.8 Radial dimensionless temperature distribution at the highest gas film temperature for different eccentricity ratios without considering temperature jump

圖9 展示了考慮溫度階躍時不同偏心率下最高氣膜溫度處徑向量綱一溫度分布。 隨著偏心率的增加，軸頸和軸承套兩側的溫度階躍均增大， 并且軸頸側的溫度階躍量要小于軸承套側的溫度階躍量， 原因是軸頸側的溫度梯度要小于軸承套側的溫度梯度。 由于溫度階躍的發生， 當偏心率為0.8 時， 軸承套側的溫度已十分接近軸頸側的溫度， 流場中的溫度基本都在T1之上， 因此在大偏心率情況下， 溫度階躍的影響不容忽略。

圖9 考慮溫度階躍時不同偏心率下最高氣膜溫度處徑向量綱一溫度分布Fig.9 Radial dimensionless temperature distribution at the highest gas film temperature for different eccentricity ratios considering temperature jump

為了更好地看出溫度階躍對于黏性熱耗散的影響， 圖10 給出了不同偏心率下氣膜溫度最大值的變化。 可以發現， 最高氣膜溫度的值隨著偏心率的增加而線性遞增， 并且考慮溫度階躍時， 溫度增加的速率更大， 再次印證了大偏心率情況下， 溫度階躍效應的影響顯著。

圖10 不同偏心率下的最高氣膜溫度變化Fig.10 Variation of maximum gas film temperature at different eccentricity ratios

4.3 不同轉速的影響

圖11 所示為偏心率ε＝0.6， 軸承平均間隙為25 μm，Ec數為10， 不考慮溫度階躍時不同轉速下間隙內最高氣膜溫度處徑向量綱一溫度分布。 可以看出，當轉速從3×104r／min 增加到1.1×105r／min 時， 間隙內氣膜溫度的最高值增加明顯。 這是因為， 隨著軸承轉速的提升， 間隙內流體的速度梯度增加， 流體剪切力增大， 黏性熱耗散的作用增強。 此外， 雖然間隙的幾何尺寸未發生變化， 但是隨著轉速的變化， 間隙內出現最大氣膜溫度的位置也在隨之改變， 可能是因為離心力的變化導致的流體流動和傳熱特性的改變。

圖11 不考慮溫度階躍時不同轉速下最高氣膜溫度處徑向量綱一溫度分布Fig.11 Radial dimensionless temperature distribution at the highest gas film temperature at different rotational speeds without considering temperature jump

圖12 所示為考慮溫度階躍時不同轉速下最高氣膜溫度處徑向量綱一溫度分布。 可以看出， 隨著轉速的增加， 軸頸和軸承套兩側的溫度階躍均明顯增大。 當轉速為7×104r／min 時， 軸承套側的溫度已十分接近T1， 而當轉速為9×104和1.1×105r／min 時，軸承套側的溫度將大于軸頸側的溫度T1， 此時間隙內整個流場的溫度都高于T1， 這是由于溫度階躍和黏性熱耗散共同作用的結果。

圖12 考慮溫度階躍時不同轉速下最高氣膜溫度處徑向量綱一溫度分布Fig.12 Radial dimensionless temperature distribution at the highest gas film temperature at different rotational speeds considering temperature jump

圖13 給出了不同轉速下間隙氣膜最高溫度的變化趨勢。 氣膜溫度最大值隨著轉速的增加而增加， 并且升高的速率逐漸增大。 當考慮溫度階躍時， 最高氣膜溫度的增加速度和不考慮溫度階躍時比較接近， 說明不考慮溫度階躍對于不同轉速下黏性熱耗散計算帶來的低估影響大致一樣。

圖13 不同轉速下的最高氣膜溫度變化Fig.13 Variation of maximum gas film temperature at different rotational speeds

5 結論

（1） 軸頸和軸承套兩側溫差較大時， 間隙溫度幾乎呈線性分布， 以熱傳導為主； 當軸頸和軸承套兩側溫度較小時， 間隙溫度呈拋物線分布， 此時黏性熱耗散作用明顯。

（2） 隨著偏心率的增加， 最大氣膜溫度線性遞增， 軸頸和軸承套兩側的溫度階躍也隨之增加， 且偏心率越大， 溫度階躍對于黏性熱耗散量計算的影響越大。

（3） 隨著轉子轉速的增加， 最大氣膜溫度升高，變化的速率逐漸增大， 軸頸和軸承套兩側的溫度階躍增大， 但溫度階躍對于不同轉速下黏性熱耗散量計算的影響程度相差不大。