王 帥，陳勁松

（北京航天發(fā)射技術(shù)研究所，北京，100076）

0 引 言

近年來隨著計算機技術(shù)的發(fā)展，數(shù)值模擬已成為燃氣流場研究的重要手段，在航天領(lǐng)域有著廣泛應(yīng)用。通過數(shù)值模擬方法能夠掌握燃氣流場的流動特性和場內(nèi)結(jié)構(gòu)的力熱環(huán)境，并且高精度的燃氣流場數(shù)值模擬技術(shù)是燃氣流場噪聲仿真的基礎(chǔ)，所以選取更適合燃氣流場模擬的計算格式，提高燃氣流場數(shù)值模擬的準確性和計算效率，具有重要的意義。

因計算流體力學(xué)（Computational Fluid Dynamics，CFD）商業(yè)軟件所提供的計算格式有限，要綜合比較多種計算格式只能通過編程實現(xiàn)，現(xiàn)選取計算格式中文獻引用較多的 Roe格式、AUSM+格式與 NND格式[1～4]編寫了燃氣流場計算程序，完成了二維前臺階激波反射算例、三維噴管燃氣內(nèi)流場算例和三維燃氣自由噴流流場算例的數(shù)值模擬。通過3種格式計算結(jié)果的對比分析，一方面能夠總結(jié)燃氣流場的流動特性，另一方面也能比較分析3種格式的特點和性能。

1 數(shù)值計算方法

程序采用有限體積法求解笛卡爾坐標系下三維Euler方程，其積分形式如下：

式中 q為守恒變量， q = [ ρ ,ρ u , ρ v, ρ w , ρ E ]T。f為對流通量，為單元網(wǎng)格體積；s為單元網(wǎng)格外表面；n為外表面微元法向單位向量；S為源項，在噴管流計算中源項為零。

根據(jù)有限體積法對上述方程在單元網(wǎng)格(I,J,K)上離散可得：

式中 下標“ I + 12”表示單元網(wǎng)格（I ,J,K） 上I方向正向面；“ I - 1 2”表示單元網(wǎng)格（I ,J,K） 上I方向負向面；F，G，H為各面上對流通量，數(shù)值上等于格心對流通量沿外表面法向分量與該面面積的乘積。

其中時間項采用顯式三階 Runge-Kutta法離散求解，而對流項分別采用3種迎風(fēng)格式處理。為便于比較，將3種格式以一種統(tǒng)一形式表示：

式中 第1項可看作中心差分項，第2項為耗散項。所有矢通量均取 I + 12面法向面通量。線性近似系數(shù)矩陣 A (q?)、單元界面馬赫數(shù) M aI+12等特征參量在各計算流體力學(xué)書籍中均有提及，這里不再展開[5]。

本文 Roe格式屬通量差分分裂（Flux Difference Splitting，F(xiàn)DS）格式，通過線性近似系數(shù)矩陣得到近似Riemann解，在具備高間斷分辨率的同時也存在非物理解出現(xiàn)的可能性，所以在計算中需要加入熵修正條件，但熵修正也會引入額外耗散。在噴管流場計算中，選用Harten型熵修正。

AUSM+格式將對流通量分裂為流動項Φ和壓力項P分別處理，從耗散項分析，其耗散項系數(shù)在差分算子?之外，類似于FDS格式，但其耗散項系數(shù)為標量形式，類似于矢通量分裂（Flux Vector Splitting，F(xiàn)VS）格式，所以 AUSM+可看作兩類格式的復(fù)合形式，AUSM+格式數(shù)值耗散小，間斷捕捉能力強，無需熵修正，也避免了大量的矩陣運算，具有較高的計算效率。

NND格式：

NND格式是由張涵信院士建立，具有TVD性質(zhì)的、空間二階精度的計算格式，在中國應(yīng)用廣泛。而在推廣到三維有限體積法的過程中，發(fā)現(xiàn)原始的通量型NND格式效果不太理想，所以采用了一種類似FDS的迎風(fēng)型NND格式[6]，該格式在計算中表現(xiàn)了更好的穩(wěn)定性。

Roe格式與AUSM+格式本身為空間一階精度，而NND格式為空間二階精度，所以在Roe與AUSM+格式中采用MUSCL重構(gòu)的方式將二者提高至二階精度，同時引入Van Leer限制器避免MUSCL重構(gòu)產(chǎn)生過大數(shù)值振蕩。

2 算法程序驗證

2.1 驗證算例計算模型

前臺階激波反射流場是CFD研究中的經(jīng)典算例，有充分的文獻結(jié)果作為比對，能夠驗證計算程序的正確性，此算例采用無量綱化處理，計算模型見圖1。

圖1 前臺階流場計算模型Fig.1 Computational Model of Forward-facing Step Flow Field

計算網(wǎng)格為邊長0.01的正方形網(wǎng)格。計算包含3類邊界條件：左邊界處為速度入口，固定來流馬赫數(shù)為3。上下兩面與臺階前面均為剛性壁面，滿足滑移反射邊界條件。右邊界處為自由邊界，滿足自由輸出條件。初始流場取自由來流均勻流場。

2.2 計算結(jié)果驗證

文獻[6]給出了Roe格式和 AUSM+格式的馬赫數(shù)等值線圖，文獻[7]給出了NND格式的密度分布云圖。程序計算結(jié)果與文獻結(jié)果比對如圖2至圖4所示，其中橫、縱坐標表示流場各點的位置坐標。

圖2 Roe格式程序計算結(jié)果與文獻結(jié)果對比Fig.2 Comparision of Simulation Result and Literature Result Using Roe Scheme

圖3 AUSM+格式程序計算結(jié)果與文獻結(jié)果對比Fig.3 Comparision of Simulation Result and Literature Result Using AUSM+ Scheme

圖4 NND格式程序計算結(jié)果與文獻結(jié)果對比Fig.4 Comparision of Simulation Result and Literature ResultUsing NND Scheme

從程序計算結(jié)果來看，3種計算格式能夠清晰地捕捉到激波位置，且激波位置基本一致，表現(xiàn)出了良好的數(shù)值穩(wěn)定性和激波捕捉能力。與文獻結(jié)果對比來看，程序計算結(jié)果與文獻結(jié)果基本一致，尤其NND格式結(jié)果非常吻合。而文獻[6]所示Roe格式與AUSM+格式結(jié)果與程序計算結(jié)果相比具有較高的數(shù)值波動，這是由于文獻所用限制器為superbee限制器，由此可見程序所用Van Leer限制器具有更明顯的數(shù)值耗散作用，但同樣能夠清晰分辨整個流場結(jié)構(gòu)。

通過與文獻結(jié)果的比較，能夠確認算法程序的計算結(jié)果是準確可靠的。之后通過噴管內(nèi)流場和自由噴流的計算，詳細比較分析計算格式在燃氣流場計算中的特性。

3 噴管內(nèi)流場模擬

3.1 計算模型

噴管內(nèi)流場計算模型選取了某試驗用直錐噴管，收縮半角為45°，擴張半角為18°。根據(jù)噴管的對稱性，建立了雙對稱面的四分之一噴管網(wǎng)格模型，見圖5。

圖5 噴管計算網(wǎng)格模型Fig.5 Computational Grid of Nozzle Flow Field

計算過程需處理4種邊界條件：入口邊界為壓力入口，固定總壓為5 MPa，總溫為3000 K，速度方向垂直向內(nèi)；出口邊界分2種情況，當出口速度為超聲速時，邊界變量由內(nèi)點插值得到，當出口速度為亞聲速時，邊界取環(huán)境壓力，其余變量由黎曼邊界條件得到；壁面邊界為絕熱光滑反射壁面；對稱面邊界滿足標量法向梯度為零，垂直對稱面的速度分量為零。

3.2 噴管對稱面計算結(jié)果對比

3種格式計算所得靜壓分布云圖如圖6所示。

圖6 噴管對稱面靜壓分布云圖Fig.6 Contours of Static Pressure at the Symmetry Plane of the Nozzle Flow Field

續(xù)圖6

經(jīng)比較可見，3種結(jié)果靜壓分布規(guī)律基本一致。管內(nèi)靜壓沿軸線方向逐漸減小，且在喉部附近靜壓梯度最大。對比軸線與壁面靜壓，可見壁面靜壓在喉部變化更劇烈，且在喉部之后存在明顯的低壓區(qū)。3種格式計算結(jié)果在擴張段稍有差異，擴張段靜壓等值線存在一個凹陷，AUSM+格式與Roe格式捕捉的比較明顯，而NND格式捕捉的較平緩。

3種格式計算所得馬赫數(shù)分布云圖如圖7所示。

圖7 噴管對稱面馬赫數(shù)分布云圖Fig.7 Contours of Mach Number at the Symmetry Plane of the Nozzle Flow Field

續(xù)圖7

從整體看，三者所表現(xiàn)出的馬赫數(shù)分布規(guī)律是一致的，馬赫數(shù)沿軸線方向逐漸增大。3種格式在收縮段與喉部附近結(jié)果非常接近，而在擴張段存在細節(jié)上的差異，Roe格式所得擴張段等值線較平滑，無明顯振蕩，AUSM+格式等值線在壁面附近存在小幅震蕩，而NND格式等值線在壁面附近同樣產(chǎn)生了波動，但波動產(chǎn)生于擴張段后半部分，且越靠近出口越明顯。

由圖6和圖7可見，3種格式在收縮段和喉部結(jié)果基本一致，主要差異存在于流動較復(fù)雜的擴張段，且馬赫數(shù)差異較明顯。為更細致的比較3種結(jié)果的差異，讀取了喉部、擴張段中部和出口處沿噴管徑向的馬赫數(shù)分布，通過馬赫數(shù)沿徑向的變化詳細比較3種格式的區(qū)別，并分析造成差異的原因。

3.3 沿噴管半徑計算結(jié)果對比

3種計算格式在喉部、擴張段中部和出口處馬赫數(shù)沿噴管徑向的分布如圖8至圖10所示。

圖8 喉部馬赫數(shù)沿徑向分布Fig.8 The Radial Mach Number Distributions at the Nozzle Throat

圖9 擴張段中部馬赫數(shù)沿徑向分布Fig.9 The Radial Mach Number Distributions at the Middle of the Nozzle Divergent Section

圖10 出口馬赫數(shù)沿徑向分布Fig.10 The Radial Mach Number Distributions at the Nozzle Exit

首先分析3種格式計算結(jié)果的相同點：馬赫數(shù)沿徑向的變化規(guī)律基本相同。在喉部，馬赫數(shù)沿徑向向外逐漸增大，在壁面附近達到最大值；在擴張段中部，馬赫數(shù)沿徑向先增加后減小，在徑向中段到達最大值；在出口處，馬赫數(shù)沿徑向先幾乎不變，然后逐漸減小，在中段梯度最大，在壁面附近變化較緩。

下面著重比較3種格式之間的差異，Roe格式在3處截面所得曲線均較平滑，表現(xiàn)出了良好的數(shù)值穩(wěn)定性；AUSM+格式在3處截面馬赫數(shù)曲線壁面附近均存在數(shù)值振蕩；而NND格式在喉部、擴張段中部曲線較光滑，而在出口處壁面附近存在較大波動?？v向比較來看，兩種格式數(shù)值振蕩產(chǎn)生的原因是不同的。

AUSM+格式在壁面附近均產(chǎn)生振蕩，可見振蕩是由壁面條件引起的。從格式本身來看，AUSM+的耗散項系數(shù)為標量形式，主要與單元界面馬赫數(shù)相關(guān)，而根據(jù)壁面邊界條件，壁面處流通量為零，再根據(jù)質(zhì)量守恒，貼壁網(wǎng)格壁面的對面流通量也幾乎為零，所以這兩面的單元界面馬赫數(shù)也近似為零，即在垂直壁面方向耗散項系數(shù)近似為零，所以可認為壁面附近網(wǎng)格計算中缺少了垂直于壁面方向的數(shù)值耗散，從而導(dǎo)致壁面附近產(chǎn)生了數(shù)值振蕩。

NND格式只在出口附近的壁面處存在數(shù)值波動，且越靠近出口，波動越明顯，可見波動主要由出口條件引起。雖然本文采用的是迎風(fēng)型NND格式，但NND格式本身對正負流通量的處理仍分別采用迎風(fēng)差分和中心差分，所以其受下游數(shù)據(jù)的影響程度要高于另兩種迎風(fēng)格式。流動達穩(wěn)態(tài)時，出口處為超聲速，數(shù)值上取出口邊界為內(nèi)點的插值，而對于真實流動，出口參數(shù)并不是簡單的內(nèi)點插值，這一誤差對NND格式造成了更大的影響。

除去上述數(shù)值振蕩的差異，3種格式在徑向還表現(xiàn)出了變化幅度的差別，以擴張段中部處馬赫數(shù)曲線最明顯，AUSM+格式變化幅度最大，Roe格式次之，而NND格式變化幅度最小，且之前在對稱面靜壓云圖的比較中同樣發(fā)現(xiàn)NND格式有較明顯的抹平現(xiàn)象。既然3種格式同為空間二階精度，那么可以認為這種幅值的差異主要是由格式不同的數(shù)值耗散造成的。相同條件下，3種格式中 AUSM+格式數(shù)值耗散最小，而 NND格式的數(shù)值耗散最大。

3.4 計算結(jié)果驗證

一維等熵噴管馬赫數(shù)曲線與3種格式噴管截面平均馬赫數(shù)曲線比較如圖11所示。

圖11 一維等熵噴管馬赫數(shù)曲線與截面平均馬赫數(shù)曲線Fig.11 The One-dimensional Isentropic Nozzle Mach Number Distribution and the Average Mach Number Distributions

由圖11可見，3種格式所得平均馬赫數(shù)曲線與理論曲線基本一致，數(shù)值誤差不超過馬赫數(shù)為0.045，說明計算結(jié)果與理論分析結(jié)果是相符的。同時可以發(fā)現(xiàn)主要偏差出現(xiàn)在收縮段，這是由于噴管收縮段斜率較大，氣流存在較大的徑向速度分量，與擬一維流動狀態(tài)存在較大差別。

另外，3種格式同樣可以互為校驗，經(jīng)過之前的比較分析，除去格式本身特性所造成的差異外，3種格式的計算結(jié)果表現(xiàn)出了高度的一致性，彼此差異很小，這同樣能夠驗證計算結(jié)果是可靠的。

4 噴管自由噴流流場模擬

4.1 計算模型

噴管自由噴流計算模型是在之前噴管模型的基礎(chǔ)上加入外流場得到的，軸向計算范圍由噴管出口前1 m到噴管后 10 m，徑向計算范圍為 5 m。網(wǎng)格模型如圖12所示。

圖12 自由噴流網(wǎng)格模型Fig.12 Computational Grid of Free Jet Flow Field

續(xù)圖12

實際計算中發(fā)現(xiàn)5 MPa、3000 K的入口條件要得到穩(wěn)定的外流場需要很長的計算時間，所以將入口總壓降為1 MPa，總溫降為2000 K。對自由噴流流場而言，只靠格式本身耗散難以得到穩(wěn)定流場解，所以計算中加入了黏性耗散項和標準 -kε湍流模型。

4.2 計算結(jié)果對比分析

圖13和圖14分別給出了燃氣自由噴流流場對稱面靜壓分布云圖和馬赫數(shù)分布云圖。圖13中數(shù)值表示等值線上靜壓值，單位為MPa，圖14中數(shù)值表示等值線上馬赫數(shù)值。

圖13 自由噴流對稱面靜壓分布等值線Fig.13 Contours of Static Pressure at the Symmetry Plane of the Free Jet Flow Field

自由噴流是經(jīng)典的燃氣流場算例，燃氣由噴管射出后，交替出現(xiàn)收縮擴張現(xiàn)象，呈葫蘆狀分布規(guī)律。由圖13和圖14可見，Roe格式與AUSM+格式計算結(jié)果充分反映了自由噴流的流動特性，且兩者表現(xiàn)出了高度一致性。主要區(qū)別在于AUSM+格式所得馬赫數(shù)云圖在燃氣流邊緣等值線存在波動，這同樣是其標量形式的數(shù)值耗散所造成的，結(jié)合之前噴管內(nèi)流場計算結(jié)果，可見AUSM+格式具有流動邊界處易產(chǎn)生數(shù)值波動的特點。

5 結(jié) 論

a）開發(fā)了基于Roe格式、AUSM+格式和NND格式的燃氣流場計算程序，3種格式在各算例中模擬結(jié)果基本一致，表現(xiàn)了良好的數(shù)值穩(wěn)定性和激波捕捉能力。前臺階激波反射流場計算結(jié)果與文獻結(jié)果相符，驗證了程序計算結(jié)果的可靠性。

b）在噴管內(nèi)流場算例中，由于AUSM+格式在壁面附近垂直壁面方向耗散項數(shù)值幾乎為零，導(dǎo)致計算結(jié)果在壁面附近產(chǎn)生數(shù)值振蕩，由于類似原因，在自由噴流算例中，其計算結(jié)果在燃氣流邊緣也產(chǎn)生了波動。

c）在噴管內(nèi)流場算例中，由于NND格式對正負流通量的處理分別采用迎風(fēng)差分和中心差分，其受下游數(shù)據(jù)的影響高于另兩種迎風(fēng)格式，導(dǎo)致計算結(jié)果在出口附近產(chǎn)生波動。

d）由于數(shù)值耗散不同，3種格式計算結(jié)果在噴管徑向數(shù)值變化幅度上存在差異，AUSM+格式變化幅度最大，NND格式變化幅度最小。