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

張建勛， 伍 林

(上海應(yīng)用技術(shù)大學(xué)機(jī)械工程學(xué)院，上海 201418)

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

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

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

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

1 模型介紹

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

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

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

2 模型求解

2.1 滑移修正控制方程

微尺度下氣膜間氣體分子的非連續(xù)型運(yùn)動(dòng)程度增強(qiáng)，階形板間存在上升和下降的氣膜壓力分布，應(yīng)用滑移系數(shù)固定的一階滑移邊界條件統(tǒng)一描述各個(gè)位置的氣體分子界面滑移狀態(tài)顯然是不符合現(xiàn)實(shí)物理規(guī)律的，新滑移邊界條件不僅滑移系數(shù)隨氣膜間當(dāng)?shù)乜伺瓟?shù)的變化而變化，而且結(jié)合分子自由程與氣膜厚度的關(guān)系將氣體分子的碰撞運(yùn)動(dòng)分為兩類進(jìn)行討論，可以更好地描述微尺度下的氣體流動(dòng)特性。

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

(1)

(2)

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

表1 滑移系數(shù)Table 1 Slip coefficients

(3)

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

(4)

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

2.2 限體積法網(wǎng)格劃分與方程的離散

因階形板模型在槽臺(tái)分界處存在氣膜厚度階躍變化的特性，若要對(duì)邊界壓力極大值進(jìn)行求解，需應(yīng)用有限體積法將計(jì)算參量轉(zhuǎn)移到有限體積邊界，解決氣膜厚度存在的奇異變化，具體的求解區(qū)域網(wǎng)格點(diǎn)劃分形式如圖2所示。

W為相鄰左側(cè)單元網(wǎng)格中心線；E為右側(cè)單元網(wǎng)格中心線；C為有限體積單元中心線；w、e分別為有限體積單元的左右邊界圖2 網(wǎng)格單元?jiǎng)澐諪ig.2 Grid cell division

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

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

(5)

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

當(dāng)Kn≤1時(shí)，有:

(6)

當(dāng)Kn>1時(shí)，有:

(7)

2.3 邊界條件與計(jì)算過程

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

圖3 計(jì)算流程圖Fig.3 Algorithm flow chart

(8)

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

3 計(jì)算結(jié)果分析

3.1 模型參數(shù)與程序驗(yàn)證

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

表2 模型計(jì)算參數(shù)Table 2 Model calculation parameters

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

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

3.2 氣膜壓力與克努森數(shù)

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

3.3 氣體流速分布

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

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

3.4 階形板承載力

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

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

4 結(jié)論

針對(duì)氣膜厚度階躍變化導(dǎo)致克努森數(shù)跨區(qū)間分布的階形板結(jié)構(gòu)，引入適用于任意克努森數(shù)的新滑邊界條件修正雷諾方程，數(shù)值模擬了不同環(huán)境克努森數(shù)下階形板間的氣體流動(dòng)特性，得到了氣膜壓力和克努森數(shù)分布，氣體流速分布和氣膜承載力等靜態(tài)特性的數(shù)值結(jié)果。得到如下結(jié)論。

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

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

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