彭寶林，郭良斌

（武漢科技大學(xué)機(jī)械自動(dòng)化學(xué)院，湖北武漢，430081）

高供氣壓力下的靜壓平行圓盤止推氣體軸承的氣膜入口流速為音速，隨著通氣截面半徑及截面面積的增大，氣流的馬赫數(shù)不斷增大，在等熵流中，隨著氣流馬赫數(shù)的增大，靜壓降低，氣體止推軸承的承載效率降低[1]。研究表明，通過改變和重新設(shè)計(jì)流道結(jié)構(gòu)，使氣流在氣膜入口處的流速降至亞音速，或降低氣膜內(nèi)超音速流的馬赫數(shù)，可以提高圓盤止推氣體軸承的承載力。

本文采用絕熱理想氣體和在氣膜入口處有微小平行段的變氣膜厚度模型，根據(jù)無徹力定常無黏流旋度保持不變的原理，將圓盤止推氣體軸承內(nèi)的流場(chǎng)簡(jiǎn)化為在半徑和氣膜厚度方向上為變量的二維無旋場(chǎng)，推導(dǎo)出氣膜間隙內(nèi)二維無旋場(chǎng)控制方程組[2]及其特征線方程和相容性方程[3]，對(duì)特征線方程和相容性方程用差分法離散，用歐拉預(yù)估－校正法計(jì)算內(nèi)點(diǎn)[4]，用等間隙氣膜內(nèi)一維超音速流動(dòng)的解析解進(jìn)行驗(yàn)算。

1 擴(kuò)散對(duì)稱流動(dòng)控制方程組

變間隙圓盤止推氣體軸承結(jié)構(gòu)簡(jiǎn)圖如圖1所示。圖1中：r1為供氣孔半徑；r2為圓盤外半徑；P＊為供氣壓力（滯止壓力）；T＊為供氣總溫（滯止溫度）；pa為環(huán)境背壓（環(huán)境大氣壓）；h0為平行段間隙高度常數(shù)；rs為平行段最大半徑；Vr、Vz分別為氣流速度的徑向分量和軸向分量。不平行段壁面方程為y＝y(tǒng)（x）。理想氣體以音速?gòu)臍饽と肟谶M(jìn)入圓盤間隙，在平行段內(nèi)是一維擴(kuò)散對(duì)稱無旋流動(dòng)，到達(dá)不平行段后是二維無旋流動(dòng)。

圖1 變間隙圓盤氣體止推軸承結(jié)構(gòu)Fig.1 Schematic diagram of the aerostatic thrust bearing of variable clearance

定常運(yùn)動(dòng)下的連續(xù)方程為

定常運(yùn)動(dòng)下的歐拉運(yùn)動(dòng)微分方程為

定常絕熱運(yùn)動(dòng)下的能量方程（沿流線）為

熱完全氣體的狀態(tài)方程為

等熵流動(dòng)中聲速關(guān)系式為

從式（2）和式（5）得

由式（1）得

將式（6）兩邊點(diǎn)乘V，得

由無旋定義?×V＝0，有

所以

從式（7）和式（8）中約去密度項(xiàng)，得到理想氣體等熵定常無旋流動(dòng)的控制方程為

用柱坐標(biāo)系表達(dá)，r軸指向半徑方向，z軸垂直于氣膜厚度方向，則式（10）可寫成

在變間隙氣膜內(nèi)依然有環(huán)向速度Vθ＝0及對(duì)稱性?／?θ＝0，式（11）中的三維流動(dòng)可以化簡(jiǎn)為氣膜厚度和半徑方向上的二維流動(dòng)。

由無旋定義

式（12）可化簡(jiǎn)為

令y＝z，x＝r，則u＝Vr，v＝Vz，式（14）和式（13）變?yōu)?/p>

式（15）、式（16）便是變間隙圓盤止推氣體軸承中氣體流動(dòng)控制方程組。其中，音速a是速度u和v的函數(shù)，即

所以控制方程組由兩個(gè)一階擬線性偏微分方程組成，因變量為u和v。

2 方程組特征線方程和相容性方程

u和v的全微分方程為

將式（15）、式（16）、式（18）、式（19）聯(lián)立，對(duì)ux、uy、vx、vy求解，有

對(duì)于式（20）所示的方程組，導(dǎo)數(shù)ux不確定的條件是方程組系數(shù)為零，即

展開式（21）得

式（22）兩邊同時(shí)除以dx2，并令dy／dx＝λ，則式（22）可改寫為

解得

式中：下標(biāo)“±”分別表示第Ⅰ族和第Ⅱ族特征線；M為馬赫數(shù)，當(dāng)M＞1時(shí)，特征線是實(shí)數(shù)。

可見特征線法僅用于超音速流場(chǎng)的計(jì)算。用式（20）中的常數(shù)項(xiàng)替換系數(shù)行列式中的某一列，并令其等于零，可以求得沿特征線上因變量u、v之間所滿足的相容關(guān)系。為避免高階行列式運(yùn)算，討論L1與L2的線性組合，即

式中：k1、k2為任意函數(shù)。

將式（15）、式（16）代入式（25），方程改寫為

將式（26）中第1個(gè)和第2個(gè)括號(hào)內(nèi)的表達(dá)式分別寫成全微分du／dx和dv／dx形式，并令dy／dx＝λ，即

由于原偏微分方程簡(jiǎn)化為常微分方程的曲線是特征線，因而給出的λ即是特征線的斜率

沿式（29）所示的特征線，式（26）可簡(jiǎn)化為u和v的全微分表達(dá)式，即相容關(guān)系式

通過求k1、k2非平凡解來確定特征線的斜率λ，同時(shí)導(dǎo)出k1、k2的關(guān)系，故將式（30）中的k1、k2約去。由式（29）有

若k1和k2有不為零的解，則下式必為零，即

解得

式（33）確定了變氣膜厚度內(nèi)定常二維超音速流場(chǎng)中的兩族特征線，由式（31）和式（32）可得

或

將上式代入式（30），得到相容關(guān)系

或

容易證明式（37）與式（38）是等價(jià)的相容關(guān)系。將特征線方程改寫為簡(jiǎn)單形式，用速度矢的模及夾角V、β表示速度分量u、v，用馬赫角α表示馬赫數(shù)M，即

將式（39）代入式（34），解得

式（40）說明特征線上各點(diǎn)的切線與該點(diǎn)流速方向的夾角為馬赫角α，按照觀察者的目光順著流速方向規(guī)定，伸向左邊的為左伸特征線C＋，或第Ⅰ族特征線；伸向右邊的為右伸征線C－，或第Ⅱ族特征線。

3 特征線法計(jì)算流場(chǎng)概述

如圖2所示，設(shè)平面xy中有一條非特征線的起始曲線AB，其上流動(dòng)參數(shù)為已知，利用兩族特征線相交形成一網(wǎng)格，該網(wǎng)格擴(kuò)展到一定的區(qū)域，計(jì)算出該區(qū)域的流動(dòng)參數(shù)。例如起始曲線AB上點(diǎn)1發(fā)出的右伸特征線C－與由點(diǎn)2發(fā)出的左伸征線特征線C＋相交于點(diǎn)4，它的位置可由特征線方程式（34）求出，沿線段1－4和2－4各有一個(gè)聯(lián)系du和dv的相容關(guān)系，利用此關(guān)系可以求出交點(diǎn)4上u和v的值。對(duì)于起始曲線AB上各離散點(diǎn)（包括端點(diǎn)A和B），均重復(fù)上述步驟，當(dāng)AB上離散點(diǎn)距離足夠小時(shí)，結(jié)果產(chǎn)生了一條參數(shù)已知的連續(xù)曲線DE，如此繼續(xù)下去，直到圖2中整個(gè)ABC區(qū)域算完為止。圖2中，點(diǎn)C是從A發(fā)出的第Ⅱ族特征線和從點(diǎn)B發(fā)出的第Ⅰ族特征線的交點(diǎn)，而對(duì)于邊界點(diǎn)來說一般只有一條特征線通過，對(duì)此必須補(bǔ)充速度與壁面相切的條件，將內(nèi)點(diǎn)和邊界點(diǎn)聯(lián)合起來計(jì)算，這樣特征線可以伸展到整個(gè)流場(chǎng)，從而計(jì)算出整個(gè)流場(chǎng)的參數(shù)。

圖2 二維特征線網(wǎng)格及內(nèi)點(diǎn)處理過程Fig.2 Two－dimensional characteristic line grid and interior point process

將特征線方程和相容關(guān)系寫成差分方程形式，變間隙內(nèi)定常二維無旋流動(dòng)特征線方程和相容關(guān)系的差分方程如表1所示，平均系數(shù)法確定的有限差分系數(shù)計(jì)算方程如表2所示。

表1 變間隙氣膜內(nèi)超音速流動(dòng)的有限差分方程Table 1 Finite difference equations of variable clearance flow within the supersonic gas film

表2 平均系數(shù)法確定的有限差分系數(shù)計(jì)算方程Table 2 Calculation equations of variable clearance flow within the supersonic gas film

4 單元處理過程

圖2中點(diǎn)1和點(diǎn)2為不在一條特征線上的已知點(diǎn)，根據(jù)特征線差分方程式（47）、式（48）和相容關(guān)系的差分方程，采用預(yù)估－校正法可計(jì)算出點(diǎn)4（x4，y4）及流動(dòng)參數(shù)。

4.1 預(yù)估步

設(shè)u＋＝u2，v＋＝v2，x＋＝x2，u－＝u1，v－＝v1，x－＝x1。由差分方程式（47）、式（48）有

4.2 校正步

如有必要，可以使用迭代校正步，直至

式中：qi代表x4、y4、u4、v4等參數(shù)；εi為對(duì)應(yīng)于每一參數(shù)的允許誤差。

數(shù)值計(jì)算表明，連續(xù)使用校正法對(duì)最終結(jié)果影響很小，一般迭代3次后結(jié)果較為穩(wěn)定。

對(duì)于待解點(diǎn)在壁面的情形（圖3），點(diǎn)2為已知點(diǎn)，由點(diǎn)2沿流動(dòng)方向往下引出1條特征線，于壁面相交于點(diǎn)4，這時(shí)待解點(diǎn)4的位置和流動(dòng)參數(shù)由兩組條件決定：一是通過2－4線段的特征線方程和相容關(guān)系方程，二是壁面曲線方程及氣流速度與壁面相切的條件方程，即

圖3 壁面點(diǎn)處理Fig.3 Wall point processing

由式（48）和式（63）可以求出x4、y4，由式（51）和式（64）可以求出u4、v4。其預(yù)估步、校正步的步驟與內(nèi)點(diǎn)的情形完全相同。

5 計(jì)算與驗(yàn)證

利用平行段一維流動(dòng)的解析解來驗(yàn)證特征線方程和相容性方程的正確性，采用文獻(xiàn)[1]中的例子，取r1＝10 mm，r2＝40 mm，P＊＝1 MPa，T＊＝293 K，此時(shí)氣膜入口總壓大于氣膜出口截面上產(chǎn)生正激波的特征滯止壓強(qiáng)，氣膜入口馬赫數(shù)為1，氣膜內(nèi)全是超音速流動(dòng)，任意截面上馬赫數(shù)與通氣截面半徑r之間的關(guān)系式[1]為

在特征線法中，給出入口處的條件為u1＝u2＝ccr≈312.251 8 m／s，v1＝v2＝0，x1＝x2。已知點(diǎn)間距y1－y2＝0.1 mm。在間隙內(nèi)聲速關(guān)系[5]為

式中：k為比熱比；R為氣體常數(shù)。

氣膜厚度h＝0.2 mm，由已知點(diǎn)向前直接步進(jìn)160次，也就是網(wǎng)格數(shù)是160×3。在特征線法中不能直接使用ccr的精確值，否則點(diǎn)1和點(diǎn)2在同一條特征線上，同時(shí)已知點(diǎn)間距不能過大，否則，式（66）根號(hào)內(nèi)會(huì)出現(xiàn)負(fù)值。

圖4 解析法與特征線法計(jì)算結(jié)果Fig.4 Contrast of analytical solution and the numerical solution characteristic line

編程后，分別用數(shù)值法和解析法（圖4）對(duì)160個(gè)點(diǎn)上的馬赫數(shù)進(jìn)行計(jì)算。為了使圖上的點(diǎn)清晰可見，僅畫出了步進(jìn)次數(shù)為奇數(shù)的點(diǎn)。從圖4可看出，解析法和特征線法得到的點(diǎn)基本重合。由程序計(jì)算結(jié)果可知，誤差小于0.01%。因此，改變起始點(diǎn)間距，總能得到預(yù)期數(shù)目點(diǎn)上的解，使得出的流動(dòng)參數(shù)接近于連續(xù)函數(shù)。

6 結(jié)語

通過改變止推氣體軸承等間隙的流道結(jié)構(gòu)，推導(dǎo)出了二維無旋超音速擴(kuò)散對(duì)稱流動(dòng)的控制方程及其特征線方程和相容性方程，用等間隙氣膜內(nèi)一維超音速流動(dòng)的解析解驗(yàn)證了特征線方程和相容性方程的正確性。

[1] 郭良斌，彭寶林.理想氣體條件下平行圓盤止推氣體軸承承載力特性的初步研究[J].武漢科技大學(xué)學(xué)報(bào)，2011，34（1）：62－68.

[2] 童秉綱，孔祥言，鄧國(guó)華.氣體動(dòng)力學(xué)[M].北京：高等教育出版社，1990.

[3] 單鵬.多維氣體動(dòng)力學(xué)基礎(chǔ)（第二版）[M].北京：航空航天出版社，2008.

[4] M J Zucrow，J D Hoffman.氣體動(dòng)力學(xué)（上下冊(cè)）[M].王汝涌，譯.北京：國(guó)防工業(yè)出版社，1984.

[5] 潘錦珊.氣體動(dòng)力學(xué)基礎(chǔ)[M].北京：國(guó)防工業(yè)出版社，1989.