特征線坐標系下柱對稱超音速流膨脹波計算方法研究

郭良斌，趙瑞兵

(武漢科技大學機械自動化學院，湖北 武漢，430081)

柱對稱超音速流常見于盤狀等離子體發生器[1-3]和靜壓圓盤止推氣體軸承[4]中。氣體軸承間隙內出現超音速流動時，其流動規律已不能用可壓縮雷諾潤滑方程來描述[5]。在圓盤止推氣體軸承的超音速流中，當某一壁面有外折角時，轉折尖點處會產生傾斜的膨脹波系。由于氣體沿圓盤徑向流動時其馬赫數是變化的[4]，膨脹波上不同點處的來流馬赫數并不相同，因此軸承中的膨脹波呈曲線形。

目前對軸向超音速來流的研究，多是關于其對細長旋成體的繞流的分析[6-7]。而圓盤止推軸承中，來流方向為徑向來流，不同于細長體。對于圓盤止推氣體軸承中柱對稱超音速流膨脹波的計算方法，目前尚未見有針對性的研究成果。

左克羅等[7]給出了直角坐標系下定常二維無旋超音速流特征線法的嚴謹數學理論和計算方法，但其邊界條件的處理也較為繁瑣，缺乏物理直觀性。僅討論單壁面無膨脹波反射時，求解區域是物理平面上求解初值問題的三角形區域，按Courant 等[8]的理論，該問題屬于特征線坐標系下的非特征初值問題，采用有限差分法來求解4個特征方程可以得到問題的解，但邊界條件的處理較為復雜。為此，本研究借鑒Liepmann[9]引入兩個黎曼不變量的方法，對特征線坐標系下的相容性方程進行積分求解，以期使邊界條件的處理更加直觀和簡便。

1 物理模型的建立

單孔圓盤止推氣體軸承的結構如圖1所示。圖1中，rj為節流孔半徑；rs為圓盤半徑；平行段半徑為rO；圓盤內任意位置處半徑為rx；氣膜厚度為h；不平行段自O點轉折，轉折角為δ。

圖1 單孔圓盤氣體止推軸承結構示意圖Fig.1 Circular gas thrust bearing with single orifice

2 數學模型的建立

給定氣體軸承中二維柱對稱流動氣膜入口處的參數。該處環形面積是流道的最小截面，故在氣體流出節流孔進入圓盤時速度滿足臨界條件，該處馬赫數Ma1=1，平行盤內任意點處半徑rx與氣體入口(節流孔)半徑rj為已知，k為氣體的絕熱指數，則可求出平行圓盤內各點處的馬赫數Ma，其計算公式[4]為：

(1)

由式(1)即可求得O點發出第一條膨脹波前不同半徑處的Ma分布。

2.1 膨脹波位置的計算

膨脹波位置示意圖如圖2所示。圖2中，μ為O點發出第一條膨脹波的馬赫角，δ為壁面轉折角，dr和dy分別是軸承徑向和氣膜厚度方向上的微元段，A1為O點發出的第一條膨脹波上的一點，A為該膨脹波與上壁面的交點。

圖2 膨脹波位置示意圖Fig.2 Diagram of position of expansion wave

將由O點發出的第一條膨脹波劃分成無限多個微元段，每個微元段可看作直線馬赫波。不同rx處馬赫角不同，微元段OA1可近似為直線馬赫波，馬赫角μ可由O點的波前馬赫數求得，即μ=arcsin(1/Mao)。

由圖2可知，dy和dr之間的關系為dy/dr=tanμ。膨脹波不同位置處馬赫角μ的值不同，μ值可由所在點的Ma求得，從而得到dr和dy的關系，確定A1點的位置，再根據式(1)即可求出A1點的馬赫數。同理，由A1點沿膨脹波依次求解，可得出整條波OA上各點的具體位置及波前參數。

2.2 自然坐標系到特征線網格的轉換

為了方便使用特征線法，在此構建特征線網格，推導由自然坐標系s-n向特征線坐標系ξ-η(見圖3)的轉化過程。取任一函數f，沿η特征線特征參數ξ保持不變，沿η特征線函數變化Δf僅與η有關。圖3中，從O到O′的變化可近似寫為：

(a)特征線網格 (b)自然坐標系

圖3特征線網格和自然坐標系

Fig.3Characteristiclinegridandnaturalcoordinatesystem

(2)

沿流線s可知：

(3)

由式(2)、式(3)可得：

(4)

根據圖3右側自然坐標系的坐標幾何關系，可將式(4)變換為：

(5)

同理可得ξ方向的導數：

(6)

2.3 普朗特-邁耶流動中相容關系式的推導

笛卡爾坐標系下普朗特-邁耶(P-M)流動的氣體動力學方程為：

(7)

(8)

在自然坐標系s-n下，速度矢量或分量可用速度模和方向(ω，θ)來表示，自變量用流線坐標(s，n)來表示，氣體動力學運動方程[9]可表示為：

(9)

(10)

由馬赫角μ與馬赫數Ma的關系式cotμ=Ma2-1可得：

或

(11)

將式(9)通乘tanμ，式(10)通乘tanμcotμ，再將式(11)代入可得：

(12)

(13)

將式(12)、式(13)相加、相減分別得：

(14)

(15)

令f=ν-θ和f=ν+θ，并分別將式(5)、式(6)代入，則式(14)、式(15)可變為：

(16)

(17)

式(16)、式(17)是P-M流動的氣體動力學方程在特征線坐標系下的表達形式，它們給出一個簡單的結果：函數ν-θ和ν+θ分別在η特征線和ξ特征線上不變。Liepmann[9]將函數ν-θ和ν+θ定義為黎曼不變量。

2.4 柱對稱流動中相容關系式推導

自然坐標系下柱對稱流動氣體動力學方程為[9]:

(18)

(19)

比較式(9)和式(18)可以看出，二維平面流動與柱對稱流動控制方程的差別在于sinθ/r項。特征線坐標系下柱對稱流動控制方程的推導過程與前述P-M流動中的推導過程類似，在此不作贅述，直接由方程式(18)、式(19)寫出在特征線坐標系下柱對稱流動的相容關系：

(沿η)

(20)

(21)

3 特征線法求轉折點后的流動參數

特征點網格位置如圖4所示，轉折點為點1，并將點1設置為坐標原點。轉折點處半徑為r1(即rO)。氣流方向角θ定義為氣流方向與圓盤軸線的夾角。軸承徑向與軸線夾角為90°，即1點處氣流方向角θ1為90°。在圖4中，從1點繞A點順時針旋轉，角度θ減小，反之則θ增大。壁面轉折角δ定義為軸承壁面與徑向的夾角，由圖4可知δ角為負值，1′點繞O點順時針旋轉，δ值減小，反之增大。由于篇幅所限，本文只討論單壁面無膨脹波反射情況。

圖4 特征線網格位置示意圖Fig.4 Diagram of position of characteristic grid

3.1 壁面點參數的求解

在求得平行圓盤間隙各點參數后，由前述關于膨脹波位置的計算方法可求出從轉折點O處發出的第一條特征線上的參數，即圖4中點1～5的參數。接著求解經轉折點膨脹后特征線上的參數。首先求解壁面點參數。已知點2參數以及壁面偏轉角度δ(見圖5)，從點2發出的特征線交壁面于點1′，待解點1'的位置由幾何位置關系決定。

圖5 壁面點求解示意圖Fig.5 Diagram of solution of wall point

由邊界條件可知，點1′處的氣體流動方向與壁面平行，其參數由特征線坐標系下的相容關系得到。具體解法如下：

線ξ21′的斜率為：

(22)

其點斜式方程為：

(23)

整理得其直線方程為：

(24)

轉折壁面的斜率為：

k=tanδ

(25)

轉折壁面的點斜式方程為：

(26)

整理得其直線方程為：

y=xtanδ+yO-xOtanδ

(27)

聯立式(24)、式(27)便可求解出1′點的位置坐標。

由兩點間的距離公式可得，2點到1′點的距離為

(28)

再由邊界條件計算出沿壁面流動的角度：

θ1′=θ1+δ

(29)

從2點到點1′沿特征線ξ由相容關系式(21)可得：

(30)

將式(30)右端括號內的量近似為常數，且具有點2的已知值，積分得：

(31)

無論在自然坐標系還是特征線網格中，ν、μ均為馬赫數Ma的函數，由普朗特-邁耶方程可知其關系為：

(32)

求解方程組式(29)、式(31)可以得到點1′處的ν1′和θ1′。再將ν1′代入式(32)，用二分法可以求出1′點的馬赫數Ma1′，進而求得：

μ1′=arcsin(1/Ma1′)

(33)

3.2 內部點參數的求解

內部點的處理是指由不在同一條特征線上的兩個已知點求解出這兩條特征線相交點的參數。如圖6所示，點3和點1′不在同一特征線上，點3順流而下發射出的特征線ξ與點1′發射出的特征線η交于點2′，已知點3和點1′的參數，則點2′參數可聯立兩相容方程求解，其位置由特征線方向所決定。下面求其數值解。

圖6 內部點求解示意圖Fig.6 Diagram of steps of solution of the interior point

由特征線性質可知，特征線ξ32′的斜率為：

(34)

點斜式方程為：

(35)

直線方程為：

(36)

特征線η1′2′的斜率為：

(37)

點斜式方程為：

(38)

直線方程為：

(39)

聯立方程式(36)、式(39)可以給出待解點2′的位置坐標(x2′,y2′)。

求解出以上各點位置，再求沿特征線變化的距離：

(40)

(41)

在特征線網格下利用特征線法求解2′點參數，由點3到點2′的關系式為：

(42)

從點1′到點2′的關系式為：

(43)

(44)

在式(43)中，μ、θ、r值為點1′的參數，可將其積分為：

(v2′-θ2′)-(v1′-θ1′)=

(45)

聯立式(44)、式(45)可得：

(46)

(47)

從而可以求出點2′的參數ν2′和θ2′。由式(32)、式(33)求出μ2′。

4 算例分析

為驗證本計算方法的可行性，將此法應用于P-M流動，并與P-M流動理論的計算結果進行比較。

算例：一個馬赫數為1.4的均勻空氣(絕熱指數k=1.4)繞外凸壁膨脹，氣流逆時針方向轉折不同角度(1°～80°)時，計算膨脹波后的最終馬赫數。

分別用P-M流動理論和本文的特征線法對單壁面無膨脹波反射的二維流動進行求解，求解時用Fortran 90編程，用Matlab繪制成可視圖，如圖7所示。從圖7中可知，采用特征線法的計算結果與P-M流動理論值相吻合，表明本文計算方法是可行的。

圖7P-M流動理論和特征線法計算結果

Fig.7ResultsofcalculationbyPlantMayerflowtheoryandthecharacteristicmethod

5 結論

(2)將二維無旋超音速流轉化到特征線坐標系下，引入Liepmann定義的黎曼不變量后，方程形式簡潔，邊界條件的處理直觀易懂。

(3)將特征線坐標系下的特征線方法應用于圓盤氣體軸承中柱對稱徑向流單壁面膨脹波的計算，計算方法具有可行性。

[1] Harada H,Tsuno K.Study of a disk MHD generator for nonequilibrium plasma generator(NPG) system[J]. Energy Conversion and Management,1998, 39(5):493-503.

[2] Okamura T,Okuno Y. Performance of disk MHD generator with high magnetic flux density[J]. Energy Conversion and Management,2001,42(7):855-866.

[3] Ishikawa M,Matsuo T. Stability analysis of MHD disk generators and application to power systems with CO2recovery[J].Energy,1997, 22(2-3):239-247.

[4] 郭良斌，彭寶林. 理想氣體條件下平行圓盤止推氣體軸承承載力特性研究[J].武漢科技大學學報,2011,34(1):62-68.

[5] 王祖溫，孫昂.靜壓氣體軸承超聲速現象的研究與發展[J].機械工程學報，2006,42(1):6-10.

[6] 童秉綱，孔祥言，鄧國華.氣體動力學[M].北京:高等教育出版社,2012.

[7] M·J·左克羅, J·D·霍夫曼.氣體動力學[M] .王汝涌，吳宗真，吳宗善，等譯.北京:國防工業出版社,1984.

[8] Courant R,Friedrichs K O.超聲速流與沖擊波[M]. 李維新，徐華生，管楚洤，等譯.北京:科學出版社,1986.

[9] Liepmann H W,Roshko A.Elements of gasdynamics[M].California:California Institute of Technology，1957.