矢通量分裂法在噴管超聲速噴流計(jì)算中的應(yīng)用①

聶 赟，樂貴高，馬大為

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

0 引言

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

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

1 矢通量分裂法

1.1 控制方程

軸對(duì)稱無粘Euler方程組表示為

其中:

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

1.2 有限體積法和空間離散

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

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

對(duì)流通量Fc為

Gc也可寫成相似的形式。

為了離散流動(dòng)控制方程:

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

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

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

1.3 空間數(shù)值方法

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

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

1≥Mn≥ －1 時(shí):

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

通量向量線性化為

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

1.4 多步Runge-Kutta時(shí)間離散

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

對(duì)方程進(jìn)行時(shí)間積分:

多步Runge-Kutta格式寫成:

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

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

1.5 斜率限制器

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

計(jì)算時(shí)M取值參考文獻(xiàn)［7］。

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

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

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

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

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

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

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

圖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相比，變化的流動(dòng)參數(shù)是超聲速來流馬赫數(shù)增大，射流的斜激波張角減小，激波反射點(diǎn)到噴管距離減小。

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

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

圖5 壁面壓強(qiáng)變化曲線Fig.5 Pressure distribution along wall

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

圖6 軸線壓強(qiáng)變化曲線Fig.6 Pressure distribution along jet axis

3 結(jié)論

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

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

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

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

［1］瞿章華，劉偉.高超聲速空氣動(dòng)力學(xué)［M］.長(zhǎng)沙:國(guó)防科技大學(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］喬陽，劉勇，徐敏.基于多塊對(duì)接網(wǎng)格的多側(cè)噴流瞬態(tài)干擾特性研究［J］.固體火箭技術(shù)，2007，30(2):98-101.

［4］李立沛，徐敏.大攻角側(cè)向噴流流場(chǎng)特性及噴流干擾效應(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.