楊風波，馬大為，樂貴高，聶 赟

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

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

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

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

1 控制方程和離散方法

1.1 控制方程

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

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

1.2 空間離散

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

式中：Δx和Δy分別為軸向和徑向坐標方向上的空間步長。

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

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

定義：

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

將無粘通量Fi＋1／2，j分裂為對流通量項和壓力通量項進行處理，即

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

1.3.1 對流項通量

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

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

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

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

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

對流通量項可以表述為

1.3.2 壓力項通量

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

式中：

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

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

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

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

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

1.4 時間離散

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

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

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

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

圖1 計算區(qū)域

表1 流動參數(shù)

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

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

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

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

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

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

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

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

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

圖4 算例3的等值線

圖5 軸向壓力分布

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

圖6 噴管壁面壓力分布

3 結(jié)束語

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

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

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

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

［1］徐旭，張振鵬.用TVD方法模擬噴管內(nèi)的橫流流場［J］.推進技術(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］樂貴高，馬大為，臧國才.火箭底部流動的TVD數(shù)值模擬［J］.彈道學報，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］候中喜，易仕和，王承堯.超聲速開式空腔流動的數(shù)值模擬［J］.推進技術(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］.空氣動力學學報，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］鄭敏，張函信.無波動、無自由參數(shù)的耗散差分格式（NND）在噴流計算中的應(yīng)用［J］.空氣動力學學報，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］陳二云，馬大為，樂貴高，等.間斷有限元方法在彈尾超音速噴流計算中的應(yīng)用［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］沈孟育，張志斌，牛曉玲.滿足抑制波動原則的五階精度緊致格式［J］.清華大學學報，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］聶赟，樂貴高，馬大為.矢通量分裂法在噴管超聲速噴流計算中的應(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］.航空動力學報，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］.空氣動力學學報，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.