聶 赟，樂貴高，馬大為

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

0 引言

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

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

1 矢通量分裂法

1.1 控制方程

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

其中:

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

1.2 有限體積法和空間離散

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

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

對流通量Fc為

Gc也可寫成相似的形式。

為了離散流動控制方程:

采用應(yīng)用散度定理有

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

通過暫時略去某些下標(biāo)給出線性化通量向量:

1.3 空間數(shù)值方法

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

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

1≥Mn≥ －1 時:

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

通量向量線性化為

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

1.4 多步Runge-Kutta時間離散

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

對方程進行時間積分:

多步Runge-Kutta格式寫成:

式中 α1為權(quán)重系數(shù)為時間步。

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

1.5 斜率限制器

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

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

2 計算結(jié)果與分析

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

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

表1 流動參數(shù)Table 1 Parameter values of fluid flow

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

從圖2可見，超聲速外流在噴管出口處開始壓縮噴流，噴流近區(qū)馬赫盤消失。形成兩道間斷的斜激波。第一道斜激波后面有膨脹扇，第二道斜激波后面是射流激波，射流膨脹區(qū)域受到外流動的壓縮反射，反射激波與射流激波相交。計算結(jié)果與試驗紋影圖反映的流動特性很相似，能夠較準(zhǔn)確地反映發(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的密度等值線圖和馬赫數(shù)等值線圖。與工況2相比，變化的流動參數(shù)是超聲速來流馬赫數(shù)增大，射流的斜激波張角減小，激波反射點到噴管距離減小。

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

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

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

［4］李立沛，徐敏.大攻角側(cè)向噴流流場特性及噴流干擾效應(yīng)［J］.固體火箭技術(shù)，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.