滑移邊界下氣膜厚度階躍變化的階形板間氣體流動

張建勛， 伍 林

(上海應用技術大學機械工程學院，上海 201418)

微尺度下的氣體流動中由于存在氣流的非平衡效應，氣體介質的連續型分布假設不再適用，通過Boltzmann方程求解微尺度下的氣體流動[1-2]存在數學上的復雜性，應用邊界速度滑移修正 Navier-Stokes(N-S)方程進行求解相對快速簡捷且模型適應性較好。

基于Maxwell基本假設提出的滑移系數為固定值的一階滑移[3]邊界條件，適用于氣膜間克努森數處于滑流區間時的微尺度流動。楊琴等[4-5]通過矩方法對Boltzmann方程求解并推導出了滑移系數隨克努森數變化的擴展速度滑移邊界條件[6]，在近過渡區內與Boltzmann方程解吻合性較好，但是在克努森數接近甚至大于1時偏差仍然較大。Wu[7]從分子動力學理論出發，細化了微尺度下氣體分子與固體壁面間的碰撞運動并將其分為兩類進行研究，推導出了滑移系數隨氣膜間克努森數變化且適用于任意克努森數的新滑移邊界條件，在氣膜厚度連續變化的楔形模型行中與F-K模型解呈現良好的吻合性。

基于階形板的動壓原理[8]衍生出了一系列溝槽類氣體軸承，如廣泛應用于干氣密封等工業領域的平面螺旋槽氣體止推軸承。Zhang等[9]、Wang等[10]應用不同的數值方法對平面螺旋槽模型進行了求解。但溝槽類氣體軸承氣膜厚度存在階躍變化，槽區和臺區克努森數跨區域分布，當臺區克努森數進入過渡區時槽區可能還處于滑流區，基于無滑移和一階滑移邊界條件得到的數值結果將會與現實結果產生較大偏差。作為溝槽類氣體軸承的基本構成單元，微尺度下階形板間的氣體流動特性的深入研究對溝槽類氣體軸承的發展具有重要意義。

代入新滑移邊界條件得到修正后的雷諾方程，應用有限體積法結合Newton-Raphson方法迭代求解了不同參數下階形板模型中的氣膜壓力和克努森數分布曲線，沿氣膜厚度方向的流速分布曲線，氣膜承載力隨速度參量的變化趨勢等數值結果，對比了不同邊界條件下所得計算結果之間的偏差并分析其產生原因，為相關微型溝槽類氣體軸承的設計提供參考。

1 模型介紹

在流體潤滑中為提高動壓軸承的相關性能，轉子或軸承上可以開有各種形式的槽，以提高軸承在工作過程中產生的動壓效應。階形板模型所產生的氣膜壓力分布曲線可以看作各種溝槽類動壓氣體軸承氣膜壓力分布的一個基本單元。圖1為階形板間的氣體流動模型示意圖，上板為臺階形平板，下板為平行板。槽區高度為h1，寬度為b1，臺區高度為h2，寬度為b2，假定上板固定不動，下板以速度u0沿X軸正方向移動。

圖1 階形板間氣體流動示意圖Fig.1 Schematic diagram of gas flow between stepped plates

關于階形板間的流體流動產生的動壓效應可以用流量守恒原則解釋是其物理原理，進氣端與出氣端的流量應該相等，槽區和臺區的流量也應該相等。階形板間氣體流速在槽區的分布為斜邊內凹的直角三角形，臺區為斜邊外凸的直角三角形。相應的槽區應存在沿X軸正方向的正壓力梯度，臺區應存在沿X軸正方向的負壓力梯度，在槽臺交界處取得壓力極值。圖1中氣體流速分布示意圖不考慮滑移邊界條件，僅作為階形板間動壓原理分析闡述說明。

2 模型求解

2.1 滑移修正控制方程

微尺度下氣膜間氣體分子的非連續型運動程度增強，階形板間存在上升和下降的氣膜壓力分布，應用滑移系數固定的一階滑移邊界條件統一描述各個位置的氣體分子界面滑移狀態顯然是不符合現實物理規律的，新滑移邊界條件不僅滑移系數隨氣膜間當地克努森數的變化而變化，而且結合分子自由程與氣膜厚度的關系將氣體分子的碰撞運動分為兩類進行討論，可以更好地描述微尺度下的氣體流動特性。

式(1)和式(2)為根據簡化的N-S方程推導出的流速u和流量方程Qx：

(1)

(2)

式中:u為氣體流速；Qx為氣體流量；p為氣膜壓力；ρ為氣體密度；h為氣膜厚度；λ為分子自由程;η為氣體黏度系數；M、N為滑移系數，如表1所示。表1中α為分子切向動量協調系數，即在邊界層處發生漫反射的氣體分子所占的比例，克努森數Kn為分子自由程與氣膜厚度的比值，f=min[1/Kn,1]。

表1 滑移系數Table 1 Slip coefficients

(3)

根據式(3)定義可得到無量綱化的雷諾方程如式(4)所示。

(4)

式(4)中：P為氣膜壓力；h為氣膜厚度；x為模型特征長度；Λx為定義軸承數。

2.2 限體積法網格劃分與方程的離散

因階形板模型在槽臺分界處存在氣膜厚度階躍變化的特性，若要對邊界壓力極大值進行求解，需應用有限體積法將計算參量轉移到有限體積邊界，解決氣膜厚度存在的奇異變化，具體的求解區域網格點劃分形式如圖2所示。

W為相鄰左側單元網格中心線；E為右側單元網格中心線；C為有限體積單元中心線；w、e分別為有限體積單元的左右邊界圖2 網格單元劃分Fig.2 Grid cell division

根據階形板的幾何特性共分為槽區(G)，臺區(L)和槽臺交界處(G-L)三類有限體積單元。

根據流量守恒原理，求解域內為無源流場故有限體積單元內部流量變化為零，任意有限體積單元左側邊界氣體流入量Qw與右側邊界氣體流出量Qe相等，如式(5)所示。式(6)、式(7)分別為新滑移邊界條件下式(4)根據上述網格點劃分形式得到的控制方程離散格式，可以看出因為根據分子自由程與氣膜間隙之間的關系將氣體分子碰撞運動進行分類討論，進一步聯系氣膜間克努森數變化對滑移邊界的影響，過渡區內氣膜間克努森數大于1的方程離散格式參數相對復雜。

(5)

式(5)中：下角標w、e′表示該參量在網格單元中的位置。

當Kn≤1時，有:

(6)

當Kn>1時，有:

(7)

2.3 邊界條件與計算過程

階形板氣體流入和流出端的環境壓力值設為已知邊界條件，沿氣體流動方向依次分布的任意有限單元體可構成以Pw、Pc、Pe為未知量的非線性方程，結合兩端已知邊界條件在求解區域內可寫成非線性三對角矩陣方程組的形式，可用Newton-Raphson方法可對其進行數值迭代求解，如式(8)所示。圖3為算法流程圖。

圖3 計算流程圖Fig.3 Algorithm flow chart

(8)

式(8)中：n為迭代次數；i表示流場模型中任意網格單元中心。

3 計算結果分析

3.1 模型參數與程序驗證

微尺度下的氣體流動特性與高空稀薄氣體效應存在共性，微尺度效應是氣膜間隙的減小導致克努森數的增大，稀薄氣體效應是氣膜間分子平均自由程的增大導致的克努森數增大。在研究微尺度氣體流動特性時可以采用模擬增大環境分子平均自由程直接增大氣膜間克努森數的方式進行相似理論研究。直接應用環境克努森數的改變可描述不同模型參數下的微尺度氣體流動。階形板幾何模型計算參數取文獻[3]，具體數值如表2所示。

表2 模型計算參數Table 2 Model calculation parameters

與文獻[3]取相同計算參數，圖4(a)中一階滑移和無滑移邊界條件下的階形板氣膜壓力分布曲線為應用本文方法所得的數值結果，與文獻[3]中的給出的解析解數值基本一致，可作為本文方法可靠性的依據。

圖4 氣膜壓力與克努森數分布Fig.4 Distribution of gas film pressure and Knudsen number

3.2 氣膜壓力與克努森數

不同環境克努森數Kn0下的階形板間無量綱氣膜壓力分布和克努森數分布如圖4所示。從圖4可以看出，氣膜壓力在槽臺交界處出現最大值，克努森數在槽臺交界處出現階躍變化，階形板槽區和臺區氣膜間克努森數跨區域分布。不考慮邊界速度滑移時的氣膜壓力數值明顯大于考慮邊界速度滑移時的計算數值，一階滑移邊界條件相對于新滑移邊界條件明顯過高地估計了階形板間的氣膜壓力分布數值，過低地估計了克努森數分布數值。在其他模型計算參數保持不變的情況下，隨著環境克努森數的增大特別是進入過渡區時，氣膜氣體可壓縮性下降導致氣膜壓力分布和克努森數分布數值在階形板臺區和槽區的線性化增強，邊界速度滑移效應增強且氣體分子碰撞頻率增大，氣膜中氣體分子運動矢能降低，氣膜壓力數值減小的同時氣膜間克努森數數值增大。

3.3 氣體流速分布

不同環境克努森數下階形板間槽區和臺區中點處沿高度方向的氣體流速分布如圖5所示，計算結果中的流速分布曲線外形與理論分析一致。從圖5可以看出，相較于無滑移邊界條件下的流速分布，考慮邊界速度滑移時的計算結果在靠近下板處的氣體流速滯后且小于下板運動速度，在靠近上板處的氣體流速超前且大于靜止的上板。隨著環境克努森數的增大，特別是進入過渡區后可以看出，邊界速度滑移效應的影響更加明顯，階形板槽區和臺區靠近上板和下板邊界處的氣體流速超前和滯后的程度相較于滑流區明顯增大。槽區靠近下板處的氣體流速受邊界速度滑移影響的程度要大于靠近上板處，臺區靠近上板處的氣體流速受邊界速度滑移的影響程度要大于下板處。隨環境克努森數的增大不同滑移邊界條件下流速分布的線性化程度增強，相同條件參數下新滑移邊界條件下的計算結果所呈現出的邊界速度滑移程度要大于一階滑移。

Ux為沿x軸方向的速度圖5 氣體流速分布Fig.5 Distribution of gas velocity

3.4 階形板承載力

不同環境克努森數下階形板無量綱承載力隨下板速度參量的變化曲線如圖6所示。從圖6可以看出，環境克努森數的增大導致階形板達到最大承載力的臨界速度不斷升高，相同環境克努森數下新滑移邊界條件的臨界速度大于一階滑移。階形板的最大承載力數值趨于穩定值，未達到臨界速度之前邊界速度滑移效應明顯降低了階形板承載力， 隨著環境克努森數的增大，氣膜間更多的分子動能減少表現為更低的承載力數值，一階滑移相較于新滑移邊界條件過高地估計了階形板的氣膜承載力數值。

圖6 階形板承載力Fig.6 Bearing capacity of step plate

4 結論

針對氣膜厚度階躍變化導致克努森數跨區間分布的階形板結構，引入適用于任意克努森數的新滑邊界條件修正雷諾方程，數值模擬了不同環境克努森數下階形板間的氣體流動特性，得到了氣膜壓力和克努森數分布，氣體流速分布和氣膜承載力等靜態特性的數值結果。得到如下結論。

(1)邊界速度滑移效應對階形板間氣體流動特性的影響程度隨著Kn0的增大而增大，微觀上體現在氣膜壓力數值的減小和氣膜間Kn的增大，靠近上板和下板邊界處的氣體流速超前與滯后于固體板邊界速度，宏觀上體現在氣膜承載力的減小與臨界速度的增大。

(2)Kn0分別為0.05和0.5時，臨界速度分別在6 000 mm/s和25 000 mm/s時得到接近的穩定承載力數值，增大下板移動速度可減小滑移邊界對階形板間氣膜承載力的影響程度。

(3)氣膜間Kn進入過渡區后一階滑移與新滑移之間產生明顯偏差，滑移系數固定的一階滑移邊界條件過低地估計了邊界速度滑移效應的影響程度，應用滑移系數隨Kn改變的新滑移邊界條件可以更加準確地描述階形板間氣體流動。