聶 赟，樂貴高，馬大為

(南京理工大學機械工程學院，南京 210094)

0 引言

飛行器以Ma＞1速度飛行時，其發(fā)動機噴管尾部噴出的燃氣流與超聲速氣流會產生強烈的相互干擾，形成間斷及稀疏波系組成的復雜干擾流場［1］，這種氣動干擾會直接作用于飛行器上，將影響其氣動特性。因此，針對超聲速射流問題開展相關研究是非常必要的。由于此類流動的風洞試驗的復雜性及昂貴的費用，促使人們努力尋求經濟性好、魯棒性強及具有高精度的數值計算方法。國內外近年來發(fā)展了許多能夠很好抑制偽振蕩的數值格式，如 TVD、ENO和 WENO等［2－4］，均有很好的分辨率。20世紀 80 年代，Steger等［5］提出了矢通量分裂法，該方法是二階精度、保單調的迎風格式，在空間離散時保持了迎風差分固有的耗散特性，能夠避免過度的耗散。

在Steger等的基礎上，本文將矢通量分裂法應用到軸對稱Euler方程的數值求解，并對發(fā)動機噴管超聲速射流問題進行了數值模擬。模擬結果與實驗紋影圖反映的流動特性基本吻合，與其他數值格式的計算結果相比較，發(fā)現(xiàn)該方法具有很好的間斷捕捉能力。

1 矢通量分裂法

1.1 控制方程

軸對稱無粘Euler方程組表示為

其中:

式中 Ω∈R2;T為時間變量;ρ為密度;u、v分別為x、y方向的速度分量;p為壓強，對于理想氣體滿足狀態(tài)方程p=(γ－1)pe;γ為比熱容比;E為單位體積的總能量。

1.2 有限體積法和空間離散

在數值模擬中，對方程組(1)矢量形式的空間離散形式采用有限體積法。將計算空間離散為具有有限體積的小單元，并將計算得到的所需變量值存儲在單元體的中心。

控制方程(1)采用矢量形式表示如下:

對流通量Fc為

Gc也可寫成相似的形式。

為了離散流動控制方程:

采用應用散度定理有

對于有離散體和區(qū)域組成的網格來說，方程(5)變化為

通過暫時略去某些下標給出線性化通量向量:

1.3 空間數值方法

矢通量分裂法基于Van-leer格式，矢通量分裂法對對流通量分裂為

定義表面法向馬赫數Mn=(Vn＇)/c，其中Vn＇為垂直于單元表面的速度。

1≥Mn≥ －1 時:

式中 nx、ny為表面法向矢量。

通量向量線性化為

然后確定QL和QR，L和R是單元表面的左邊和右邊的外推變量值。通過外推原始變量值，然后構造左邊和右邊的守恒變量值。

1.4 多步Runge-Kutta時間離散

離散方程的求解通過時間步進法實現(xiàn)，采用多步Runge-Kutta方法數值求解方程組。在加入限制器的條件下，當采用矢通量分裂法進行空間離散時，數值解的精度在空間上可以達到2階。因此，為了和空間上保持一致高階精度，時間方向采用多步的Runge-Kutta方法。首先將時間［0，T］剖分為

對方程進行時間積分:

多步Runge-Kutta格式寫成:

式中 α1為權重系數為時間步。

多步Runge-Kutta格式在m=1時為一階時間精度，m≥2時為二階時間精度。矢通量分裂法對時間步長有嚴格的穩(wěn)定性限制條件，用多步Runge-Kutta顯示時間步進格式進行時間離散時的時間步有CFL條件限制。

1.5 斜率限制器

在高階情況下，為了提高方法的穩(wěn)定性，在時間方向的迭代格式中引入局部斜率限制器，其思想是限制一階斜率［6］，表達式為

計算時M取值參考文獻［7］。

2 計算結果與分析

圖1為發(fā)動機噴管噴流計算區(qū)域，?。?∶5］×［0∶1.5］，噴管尺寸與流動參數見表1，其中 l為噴管長度，r和R分別為噴管出口的內外徑，下標j表示噴管出口的流動參數，下標∞為超聲速外流的流動參數，計算區(qū)域內的網格數為6 000個。

圖1 計算區(qū)域Fig.1 Calculation domain

表1 流動參數Table 1 Parameter values of fluid flow

圖2、圖3分別給出了工況1、2的密度等值線圖和馬赫數等值線圖及在相同流動條件下的試驗紋影圖。

從圖2可見，超聲速外流在噴管出口處開始壓縮噴流，噴流近區(qū)馬赫盤消失。形成兩道間斷的斜激波。第一道斜激波后面有膨脹扇，第二道斜激波后面是射流激波，射流膨脹區(qū)域受到外流動的壓縮反射，反射激波與射流激波相交。計算結果與試驗紋影圖反映的流動特性很相似，能夠較準確地反映發(fā)動機噴管噴流的流動特性。

從圖3可看到，超聲速外流與噴管噴流作用后，形成一道斜激波與射流激波組成的曲線激波，同樣受到外流壓縮，射流膨脹區(qū)域反射后形成反射激波。與工況1相比，反射點到噴管的距離明顯增大，噴管出口處斜激波張角變大。在試驗紋影圖中得到了很好的驗證。

圖2 工況1的等值線和紋影圖Fig.2 Contours and experimental schlieren picture of case 1

圖3 工況2的等值線和紋影圖Fig.3 Contours and experimental schlieren picture of case 2

圖4是工況3的密度等值線圖和馬赫數等值線圖。與工況2相比，變化的流動參數是超聲速來流馬赫數增大，射流的斜激波張角減小，激波反射點到噴管距離減小。

圖4 工況3等值線分布Fig.4 Contours distribution of case 3

圖5為工況1、2的噴管壓強變化與實驗結果比較曲線，p/p∞表示噴管壁面壓強與自由來流壓強之比。噴管壁面壓強在出口處略微下降，數值結果與文獻［8］吻合較好。

圖5 壁面壓強變化曲線Fig.5 Pressure distribution along wall

圖6為工況2、3沿軸線的壓強分布，壓強變化趨勢比較吻合。差異的造成主要是由于2種方法對間斷點捕捉能力的不同。在激波發(fā)生反射后，反射點兩側物理量不連續(xù)。采用矢通量分裂法計算的壓強曲線在反射點后有較大的壓強梯度變化，而采用TVD格式在反射點計算的峰值被抹平，是由于該格式解的精度在該處降為一階。因此，采用矢通量分裂法模擬超聲速噴流是有效的，在反射點不會出現(xiàn)為振蕩或抹平現(xiàn)象。

圖6 軸線壓強變化曲線Fig.6 Pressure distribution along jet axis

3 結論

(1)將二維矢通量分裂法應用到軸對稱Euler方程的數值求解，并對發(fā)動機噴管超聲速射流問題進行數值模擬，可實現(xiàn)對噴管超聲速噴流問題的一個預測。

(2)用矢通量分裂迎風格式計算噴管超聲速噴流流場是一條可實現(xiàn)的途徑。利用矢通量分裂法對流動物理量和流場間斷的高精度、高分辨率計算，能夠促進噴管射流的應用和研究。

(3)利用本文的矢通量分裂法對流場間斷的計算是成功的。能夠有效的抑制間斷附近的數值振蕩，與其他格式的數值解相比，優(yōu)于TVD格式。

(4)以二維軸對稱噴管超聲速噴流計算結果為基礎，模擬結果與實驗紋影圖反映的流動特性基本吻合，與其他數值格式的結果相比，該方法能夠較準確的捕捉到激波在軸線反射點的位置，且具有較高的分辨率，表明該方法在超聲速噴流的數值計算中具有廣泛的應用前景，有必要對其進一步的研究。

［1］瞿章華，劉偉.高超聲速空氣動力學［M］.長沙:國防科技大學出版社，2001.

［2］Liu X D，Osher S，Chan T.Weighted essentially non-oscillatory schemes［J］.J.Comp.Phys，1994，115(1):200-212.

［3］喬陽，劉勇，徐敏.基于多塊對接網格的多側噴流瞬態(tài)干擾特性研究［J］.固體火箭技術，2007，30(2):98-101.

［4］李立沛，徐敏.大攻角側向噴流流場特性及噴流干擾效應［J］.固體火箭技術，2011，34(5):543-547.

［5］Steger J L，Warming R F.Flux vector splitting of the inviscid gas dynamic equation with application to finite difference methods［J］.J.Comp.Phys，1981，40:267-272.

［6］Cockburn B，Shu C W.The Runge-Kutta local projection P1-discontinuous Galerkin method for scalar conservation law［J］.M2AN，1991，337:337-361.

［7］Cockburn B，Shu C W.TVD Runge-Kutta discontinuous Galerkin method for conservation laws V:Multidimensional systems［J］.J.Comp.Phys，1998，144:199-244.

［8］Agrell J，White RA.An experimental investigation of supersonic axisymmetric flow over boattails containing a centered propulsive jet［C］.Sweden FFA Tech Note AU-913，1974.