楊風(fēng)波，馬大為，樂(lè)貴高，聶 赟

（1.南京理工大學(xué) 機(jī)械工程學(xué)院，南京210094；2.北京機(jī)電工程總體設(shè)計(jì)部，北京100854）

超聲速射流問(wèn)題廣泛存在于現(xiàn)代航空航天和火箭導(dǎo)彈發(fā)射等科學(xué)技術(shù)領(lǐng)域，而超聲速伴隨使得已經(jīng)具有間斷波系結(jié)構(gòu)的超聲速射流問(wèn)題趨于復(fù)雜化。鑒于試驗(yàn)設(shè)備及試驗(yàn)費(fèi)用昂貴，不能提供流場(chǎng)波系結(jié)構(gòu)演化的具體過(guò)程，研究分辨率高、求解經(jīng)濟(jì)性好的數(shù)值格式，開展燃?xì)馍淞鞯臄?shù)值模擬工作，成為燃?xì)馍淞餮芯款I(lǐng)域的一種必要手段，也是試驗(yàn)研究的有效補(bǔ)充。

目前，在模擬超聲速流動(dòng)方面，穩(wěn)定性好、間斷分辨率高的格式主要有TVD、ENO、WENO、NND、有限元、緊致格式和矢通量分裂等［1－8］。文獻(xiàn)［2，6，8］分別采用TVD格式、間斷有限元格式和矢通量法對(duì)超聲速噴流進(jìn)行了數(shù)值模擬，捕捉的波系結(jié)構(gòu)也越來(lái)越清晰，但均沒有給出壓力等值線，對(duì)接觸間斷沒有給出具體解釋。另外，上述格式求解過(guò)程多數(shù)伴隨著特征矩陣的分裂和特征值的求解，過(guò)程繁瑣、求解經(jīng)濟(jì)性差。上世紀(jì)90年代，Liou和Stefen綜合矢通量分裂方法（FVS）［9］和通量差分分裂方法（FDS）［10］的優(yōu)勢(shì)，構(gòu)造了著名的混合通量差分AUSM［11］（Advection Upstream Splitting Method）格式。為了有效抑制“紅寶石”現(xiàn)象，消除間斷附近的數(shù)值偽振蕩和抹平現(xiàn)象，Liou后續(xù)推出了 AUSM＋［12］格式，Kim 提出了AUSMPW［13］和 AUSMPW＋［14］改 進(jìn)格 式。AUSM系列格式將數(shù)值通量分解為對(duì)流項(xiàng)通量和壓力項(xiàng)通量，對(duì)流通量包含了馬赫數(shù)、當(dāng)?shù)芈曀?，以及相?yīng)的流動(dòng)特征通量，壓力通量?jī)H僅含有壓力項(xiàng)；另外，該系列格式通量求解無(wú)需求解Jacobian矩陣，在計(jì)算效率方面具有很大優(yōu)勢(shì)。

以Liou提出的AUSM＋格式為基礎(chǔ)，針對(duì)超聲速伴隨射流特殊波系結(jié)構(gòu)及超聲速區(qū)和低速區(qū)同時(shí)存在的問(wèn)題，對(duì)該格式進(jìn)行了適宜改進(jìn)，構(gòu)造了空間和時(shí)間均具有二階精度的數(shù)值格式，并發(fā)展到二維軸對(duì)稱Euler方程組進(jìn)行數(shù)值求解。數(shù)值計(jì)算結(jié)果與試驗(yàn)紋影圖吻合良好；該格式能清晰捕捉入射激波、反射激波、λ激波，接觸間斷等現(xiàn)象，并且基本無(wú)“紅寶石”現(xiàn)象，表明該格式具有較強(qiáng)的激波和間斷捕捉能力。

1 控制方程和離散方法

1.1 控制方程

在忽略化學(xué)反應(yīng)、粘性和熱傳導(dǎo)的假設(shè)下，軸對(duì)稱流動(dòng)Euler方程組的守恒形式可以表述為

式中：t為時(shí)間變量；ρ為密度；u，v分別為x，y方向速度分量；p為壓強(qiáng)；R為氣體常數(shù)；T為氣體溫度；γ為比熱比；e為單位質(zhì)量?jī)?nèi)能；E為單位體積總能量；h為單位體積總焓。

1.2 空間離散

采用改進(jìn)型的AUSM＋格式對(duì)控制方程（1）進(jìn)行離散，得到如下軸對(duì)稱Euler方程組的半離散有限差分方程：

式中：Δx和Δy分別為軸向和徑向坐標(biāo)方向上的空間步長(zhǎng)。

1.3 空間格式構(gòu)造方法

以單元交界面Fi＋1／2，j無(wú)粘通量 為 例，改進(jìn) 的AUSM＋格式數(shù)值通量的構(gòu)造過(guò)程如下。

定義：

式中：Fc，F(xiàn)p分別為對(duì)流通量項(xiàng)和壓力通量項(xiàng)；c為當(dāng)?shù)芈曀?；Q＝（ρρuρvρh）；Ma為馬赫數(shù)。

將無(wú)粘通量Fi＋1／2，j分裂為對(duì)流通量項(xiàng)和壓力通量項(xiàng)進(jìn)行處理，即

式中：定義Mai＋1／2＝ui＋1／2／ci＋1／2；L／R表示速度為正向，則取邊界L通量（左通量），速度為反向，則取邊界R通量（右通量）。

1.3.1 對(duì)流項(xiàng)通量

考慮到界面馬赫數(shù)Mai＋1／2，j受左右特征波影響的物理特性，將邊界馬赫數(shù)作如下處理，即：

式中：＝Ma＋（MaL），＝Ma－（MaR）。

為了改進(jìn)傳統(tǒng)AUSM格式的不足，更好地反應(yīng)特征波傳遞對(duì)流場(chǎng)的影響，消除激波后形成的“紅寶石”現(xiàn)象，給出如下形式的馬赫數(shù)分裂函數(shù)［15］：

在低速流動(dòng)區(qū)域，為加速收斂，且抑制數(shù)值振蕩，對(duì)邊界馬赫數(shù)進(jìn)行了修正［16］：

式中：δ0為小量，δ0＝δ1｜Ma∞｜，δ1∈［0.05，0.5］。

對(duì)流通量項(xiàng)可以表述為

1.3.2 壓力項(xiàng)通量

和對(duì)流通量類似，界面上引入壓力分裂函數(shù)：

式中：

1.3.3 高階格式構(gòu)造

為將AUSM＋格式推廣到高階格式，且防止在激波附近出現(xiàn)過(guò)沖或過(guò)膨脹，捕捉更為清晰的波系結(jié)構(gòu)，引進(jìn)Van Leer可微保單調(diào)限制器，對(duì)界面兩側(cè)原始變量進(jìn)行處理，以獲得重構(gòu)的無(wú)粘通量。

Van Leer可微保單調(diào)限制器［17］滿足：

式中：ΔW＋i＝Wi＋1－Wi，ΔW－i＝Wi－Wi－1，W為原始變量（ρuvp）T；ε是一個(gè)防溢出的小量。

則可構(gòu)造如下界面兩側(cè)原始變量的數(shù)值通量：

1.4 時(shí)間離散

為了與空間高精度匹配，時(shí)間離散采用二階精度Runge－Kutta法。

將式（2）右邊各項(xiàng)統(tǒng)一定義為Ri，j，則有：

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

圖1為發(fā)動(dòng)機(jī)噴管噴流計(jì)算區(qū)域，表示通過(guò)射流中心軸線一半計(jì)算區(qū)域。?。?∶1.5］×［0∶3.5］（徑向×縱向）的計(jì)算區(qū)域（將噴口外徑無(wú)量綱化為1，且噴口外徑長(zhǎng)度用de表示，如圖2～圖4所示），采用144×60正方形網(wǎng)格對(duì)該區(qū)域進(jìn)行離散。流動(dòng)參數(shù)見表1，r和R分別為噴管出口的內(nèi)、外徑，下標(biāo)j表示噴管出口的流動(dòng)參數(shù)，下標(biāo)∞表示超聲速外流的流動(dòng)參數(shù)。GH為入口邊界，DE為超聲速來(lái)流條件，EF和FG為鏡面反射邊界，CD為簡(jiǎn)單波邊界，CB為外推邊界，AB為軸對(duì)稱邊界，其它計(jì)算區(qū)域假定為超聲速自由來(lái)流的初始條件。

圖1 計(jì)算區(qū)域

表1 流動(dòng)參數(shù)

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

圖2 算例1的等值線與紋影圖

從圖2可以看出，在超聲速伴隨射流中，沒有出現(xiàn)燃?xì)馍淞髦谐Ｒ姷娜c(diǎn)等經(jīng)典波系結(jié)構(gòu)。超聲速來(lái)流在噴管拐角處強(qiáng)烈壓縮噴口射流，出現(xiàn)兩道斜激波、射流激波。從密度等值線可以看出，第二道斜激波和第一道入射激波中間存在兩道間斷，從圖2（b）中可看出，由于壓力等值線中無(wú)間斷產(chǎn)生，故該間斷為接觸間斷。文獻(xiàn)［2，6，8］均沒有給出壓力等值線，沒有對(duì)接觸間斷進(jìn)行判定。第一道斜激波后伴隨有膨脹扇產(chǎn)生，第二道斜激波下方的射流激波受壓縮遇到中心軸線發(fā)生反射，反射激波和接觸間斷相交，馬赫盤結(jié)構(gòu)消失，從以上分析可看出，超音速伴隨使得已經(jīng)具有間斷的超音速射流波系結(jié)構(gòu)更加復(fù)雜。另外，通過(guò)對(duì)比分析可以看出，流場(chǎng)結(jié)構(gòu)與相同流動(dòng)條件下三階DG－FEM格式［6］的計(jì)算結(jié)果是一致，優(yōu)于二階TVD格式［2］，數(shù)值模擬的波系結(jié)構(gòu)和流場(chǎng)特征和試驗(yàn)紋影圖［18］吻合良好，表明二階AUSM格式的數(shù)值模擬結(jié)果是可信的。

另外，為了說(shuō)明本文網(wǎng)格布局與數(shù)值求解結(jié)果的合理性，采用三種網(wǎng)格布局對(duì)工況1進(jìn)行求解與對(duì)比分析。計(jì)算結(jié)果如表2所示（軸線上的壓力用pz表示），從表2可以看出，網(wǎng)格過(guò)稀對(duì)計(jì)算結(jié)果的精度影響較大，而網(wǎng)格達(dá)到一定密度時(shí)，計(jì)算結(jié)果與網(wǎng)格基本無(wú)關(guān)，從表2可以看出，網(wǎng)格3在網(wǎng)格2的基礎(chǔ)上增加網(wǎng)格密度2／3，射流激波反射點(diǎn)的位置和反射點(diǎn)壓力基本一致。故本文采用144×60的網(wǎng)格密度（徑向×縱向）。

表2 工況1的三種網(wǎng)格布局對(duì)比

圖3給出了工況2的密度、壓力等值線圖和對(duì)應(yīng)的試驗(yàn)紋影圖。由于壓力比提高，射流中心區(qū)域明顯更大，馬赫盤位置向下游移動(dòng)，出現(xiàn)了由一道斜激波和一道入射激波共同組成的λ激波，與相同流動(dòng)條件下三階DG－FEM格式［6］的計(jì)算結(jié)果是一致的。和工況1類似，在斜激波和入射激波中間出現(xiàn)了一道弧線狀的接觸間斷。相對(duì)于工況1而言，工況2的反射激波張角更大，激波和接觸間斷均變成了一道，這些超聲速伴隨射流流動(dòng)特征與相同流動(dòng)條件下的試驗(yàn)紋影照片［18］反映的流場(chǎng)特征基本一致。從圖3也可以看出，密度和壓力等值線基本沒出現(xiàn)“紅寶石”現(xiàn)象，激波清晰且薄，說(shuō)明本文改進(jìn)合理，有較高分辨率。

圖3 算例2的等值線與紋影圖

從圖4可以看出，改進(jìn)的AUSM＋格式捕捉到的流場(chǎng)波系結(jié)構(gòu)很清晰，明顯優(yōu)于TVD格式［2］。和工況2相比，工況3中仍然能清晰看到入射激波和斜激波組成的λ激波，以及其間的接觸間斷。另外，對(duì)比工況2和工況3，可以看出，當(dāng)保持噴管出口馬赫數(shù)和出口壓力不變而增大出口溫度比時(shí)，出口斜激波張角幾乎保持不變，但馬赫盤的位置明顯更靠近噴口，射流激波在中心軸線上的正規(guī)反射點(diǎn)距噴口平面的距離明顯更小，這與文獻(xiàn)［2］的分析是一致的。

圖5給出了3個(gè)工況軸線上的壓力pz分布規(guī)律。圖中給出了改進(jìn)的AUSM＋格式和二階TVD格式計(jì)算對(duì)比值，從圖中可以看出2種方法的計(jì)算結(jié)果吻合較好，但在間斷附近存在差異。當(dāng)噴流在膨脹區(qū)受超聲速來(lái)流壓縮，在中心軸線上發(fā)生正規(guī)反射后，反射點(diǎn)后的馬赫數(shù)降低而壓力迅速升高，壓力呈現(xiàn)大梯度變化現(xiàn)象。在計(jì)算區(qū)域內(nèi)，工況1和工況3中，改進(jìn)的AUSM＋格式的壓力計(jì)算值在反射間斷點(diǎn)處迅速上升，表現(xiàn)為在間斷點(diǎn)處具有較強(qiáng)的激波捕捉能力，而TVD格式在反射點(diǎn)捕捉到的梯度有抹平現(xiàn)象，捕捉間斷能力稍差。

圖4 算例3的等值線

圖5 軸向壓力分布

文獻(xiàn)［18］給出了和本文工況1和工況2同初始條件和結(jié)構(gòu)參數(shù)下流場(chǎng)中噴管壁面的壓力和遠(yuǎn)場(chǎng)壓力的比值參數(shù)，為進(jìn)一步驗(yàn)證改進(jìn)的AUSM＋格式在超聲速伴隨射流中的計(jì)算精度，圖6給出了工況1和工況2情況下超聲速伴隨流場(chǎng)中噴管壁面的壓力pw分布曲線和在相同流動(dòng)條件下測(cè)得的試驗(yàn)結(jié)果。pw／p∞為壁面局部壓力和無(wú)窮遠(yuǎn)處?kù)o止大氣壓力之比，圖6中黑色間斷點(diǎn)為文獻(xiàn)［18］中的試驗(yàn)結(jié)果。從圖6（a）和圖6（b）中可以看出，數(shù)值計(jì)算結(jié)果與相同條件下的試驗(yàn)結(jié)果吻合很好。

圖6 噴管壁面壓力分布

3 結(jié)束語(yǔ)

將改進(jìn)的AUSM＋格式發(fā)展到軸對(duì)稱Euler方程組進(jìn)行數(shù)值求解，給出了無(wú)粘數(shù)值通量的詳細(xì)構(gòu)造，對(duì)3種超聲速伴隨射流進(jìn)行了數(shù)值模擬，結(jié)果顯示，噴流出口壓力越大，馬赫盤位置越靠近下游，出口斜激波張角和出口溫度基本無(wú)關(guān)。

改進(jìn)的AUSM＋格式能有效消除“紅寶石”現(xiàn)象，捕捉到的波系結(jié)構(gòu)清晰，較好地預(yù)測(cè)了接觸間斷等超聲速伴隨射流特有的波系結(jié)構(gòu)。

對(duì)比已有二階TVD格式，改進(jìn)的AUSM＋格式對(duì)激波間斷具有更強(qiáng)的捕捉能力，能較好抑制間斷附近振蕩、前沖、抹平等現(xiàn)象，能較好反應(yīng)間斷附近變量的大梯度變化現(xiàn)象。

軸對(duì)稱工況的波系結(jié)構(gòu)等值線圖與相同工況下的已有文獻(xiàn)結(jié)果、試驗(yàn)紋影照片吻合較好。

［1］徐旭，張振鵬.用TVD方法模擬噴管內(nèi)的橫流流場(chǎng)［J］.推進(jìn)技術(shù)，1997，18（15）：53－57.XU Xu，ZHANG Zhen－peng.The numerical simulation of cross－flow in nozzle by using a TVD scheme［J］.Journal of Propulsion Technology，1997，18（15）：53－57.（in Chinese）

［2］樂(lè)貴高，馬大為，臧國(guó)才.火箭底部流動(dòng)的TVD數(shù)值模擬［J］.彈道學(xué)報(bào)，1995，7（1）：35－40.LE Gui－gao，MA Da－wei，ZANG Guo－cai.Numerical simulation of flow in rocket base by TVD［J］.Journal of Ballistics，1995，7（1）：35－40.（in Chinese）

［3］候中喜，易仕和，王承堯.超聲速開式空腔流動(dòng)的數(shù)值模擬［J］.推進(jìn)技術(shù)，2001，22（5）：400－403.HOU Zhong－xi，YI Shi－h(huán)e，WANG Cheng－yao.Numerical analysis of supersonic open cavity［J］.Journal of Propulsion Technology，2001，22（5）：400－403.（in Chinese）

［4］徐文燦，黃振宇.高精度ENO格式在射流數(shù)值模擬中的應(yīng)用［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2001，19（4）：401－406.XU Wen－can，HUANG Zhen－yu.Calculations of jets flowfield with high resolutions ENO［J］.Acta Aerodynamica Sinica，2001，19（4）：401－406.（in Chinese）

［5］鄭敏，張函信.無(wú)波動(dòng)、無(wú)自由參數(shù)的耗散差分格式（NND）在噴流計(jì)算中的應(yīng)用［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1989，7（3）：35－40.ZHENG Min，ZHANG Han－xin.Application of non－oscillatory and non－free－parameters disspative finit difference scheme to the calculation of free－jet flows［J］.Acta Aerodynamica Sinica，1989，7（3）：35－40.（in Chinese）

［6］陳二云，馬大為，樂(lè)貴高，等.間斷有限元方法在彈尾超音速噴流計(jì)算中的應(yīng)用［J］.計(jì)算物理，2008，25（6）：705－710.CHEN Er－yun，MA Da－wei，LE Gui－gao，et al.Discontinuous finite element method for supersonic flow of a missile propulsive jet［J］.Chinese Journal of Computational Physics，2008，25（6）：705－710.（in Chinese）

［7］沈孟育，張志斌，牛曉玲.滿足抑制波動(dòng)原則的五階精度緊致格式［J］.清華大學(xué)學(xué)報(bào)，2002，42（11）：79－82.SHEN Meng－yu，ZHANG Zhi－bin，NIU Xiao－ling.Fifth－order accurate compact scheme that suppresses oscillations［J］.J Tsinghua Univ，2002，42（11）：79－82.（in Chinese）

［8］聶赟，樂(lè)貴高，馬大為.矢通量分裂法在噴管超聲速噴流計(jì)算中的應(yīng)用［J］.固體火箭技術(shù)，2013，36（1）：56－60.NIE Yun，LE Gui－gao，MA Da－wei.Application of flux vector splitting method to the calculation of supersonic flow on nozzle jet［J］.Journal of Solid Rocket Technology，2013，36（1）：56－60.（in Chinese）

［9］LIOU M S，STEFFEN C J.A new flux spliting［J］.Journal of Computational Physics，1993，107（1）：23－39.

［10］VAN L B.Flux－vector splitting for the euler equation［J］.Lecture Notes in Phys，1982，170：507－512.

［11］LIOU Meng－sing.A new flux splitting scheme［J］.Journal of Computational Physics，1993，107：23－39.

［12］LIOU M S.Progress towards an improved CFD methods：AUSM＋，AIAA95－1701－CP［R］.1995.

［13］KIM K H，RHO O H.An important of AUSM scheme by introducing the pressure－based weight functions［J］.Computers＆ Fluids，1998，27（3）：38－80.

［14］KIM K H，RHO O H.Methods for the accurate computations of hypersonics flows，I.AUSMPW＋scheme［J］.Journal of Computational Physics，2001，174：38－80.

［15］劉晨，王江峰，伍貽兆.預(yù)處理AUSM＋格式在全速域反應(yīng)流模擬中的應(yīng)用［J］.航空動(dòng)力學(xué)報(bào)，2009，24（8）：1 831－1 836.LIU Chen，WANG Jiang－feng，WU Yi－zhao.Application of preconditioned AUSM＋scheme on numerical simulation of reacting flows at all speeds［J］.Journal of Aerospace Power，2009 24（8）：1 831－1 836.（in Chinese）

［16］梁德旺，王可.AUSM＋的改進(jìn)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2004，22（4）：404－409.LIANG De－wang，WANG Ke.Improvement of AUSM ＋scheme［J］.Acta Aerodynamica，2004，22（4）：404－409.（in Chinese）

［17］ANDERSON W K，THOMAS J L.VAN L B.Comparison of finite volume flux vector splitting for the Euler equations［J］.AIAA，1986，24：1 453－1 460.

［18］AGRELL J，WHITE R A.An experiment investigation of supersonic axisymmetric flow over boattails containing a centered propulsive jet，F(xiàn)FA－TN－AU－913［R］.1974.