燃?xì)饬鲌鏊惴ǔ绦驅(qū)崿F(xiàn)及對(duì)比

王 帥，陳勁松

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

0 引 言

近年來隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬已成為燃?xì)饬鲌鲅芯康闹匾侄?，在航天領(lǐng)域有著廣泛應(yīng)用。通過數(shù)值模擬方法能夠掌握燃?xì)饬鲌龅牧鲃?dòng)特性和場內(nèi)結(jié)構(gòu)的力熱環(huán)境，并且高精度的燃?xì)饬鲌鰯?shù)值模擬技術(shù)是燃?xì)饬鲌鲈肼暦抡娴幕A(chǔ)，所以選取更適合燃?xì)饬鲌瞿M的計(jì)算格式，提高燃?xì)饬鲌鰯?shù)值模擬的準(zhǔn)確性和計(jì)算效率，具有重要的意義。

因計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）商業(yè)軟件所提供的計(jì)算格式有限，要綜合比較多種計(jì)算格式只能通過編程實(shí)現(xiàn)，現(xiàn)選取計(jì)算格式中文獻(xiàn)引用較多的 Roe格式、AUSM+格式與 NND格式[1～4]編寫了燃?xì)饬鲌鲇?jì)算程序，完成了二維前臺(tái)階激波反射算例、三維噴管燃?xì)鈨?nèi)流場算例和三維燃?xì)庾杂蓢娏髁鲌鏊憷臄?shù)值模擬。通過3種格式計(jì)算結(jié)果的對(duì)比分析，一方面能夠總結(jié)燃?xì)饬鲌龅牧鲃?dòng)特性，另一方面也能比較分析3種格式的特點(diǎn)和性能。

1 數(shù)值計(jì)算方法

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

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

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

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

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

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

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

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

NND格式：

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

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

2 算法程序驗(yàn)證

2.1 驗(yàn)證算例計(jì)算模型

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

圖1 前臺(tái)階流場計(jì)算模型Fig.1 Computational Model of Forward-facing Step Flow Field

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

2.2 計(jì)算結(jié)果驗(yàn)證

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

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

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

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

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

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

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

3.1 計(jì)算模型

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

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

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

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

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

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

續(xù)圖6

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

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

圖7 噴管對(duì)稱面馬赫數(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é)果非常接近，而在擴(kuò)張段存在細(xì)節(jié)上的差異，Roe格式所得擴(kuò)張段等值線較平滑，無明顯振蕩，AUSM+格式等值線在壁面附近存在小幅震蕩，而NND格式等值線在壁面附近同樣產(chǎn)生了波動(dòng)，但波動(dòng)產(chǎn)生于擴(kuò)張段后半部分，且越靠近出口越明顯。

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

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

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

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

圖9 擴(kuò)張段中部馬赫數(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種格式計(jì)算結(jié)果的相同點(diǎn)：馬赫數(shù)沿徑向的變化規(guī)律基本相同。在喉部，馬赫數(shù)沿徑向向外逐漸增大，在壁面附近達(dá)到最大值；在擴(kuò)張段中部，馬赫數(shù)沿徑向先增加后減小，在徑向中段到達(dá)最大值；在出口處，馬赫數(shù)沿徑向先幾乎不變，然后逐漸減小，在中段梯度最大，在壁面附近變化較緩。

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

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

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

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

3.4 計(jì)算結(jié)果驗(yàn)證

一維等熵噴管馬赫數(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，說明計(jì)算結(jié)果與理論分析結(jié)果是相符的。同時(shí)可以發(fā)現(xiàn)主要偏差出現(xiàn)在收縮段，這是由于噴管收縮段斜率較大，氣流存在較大的徑向速度分量，與擬一維流動(dòng)狀態(tài)存在較大差別。

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

4 噴管自由噴流流場模擬

4.1 計(jì)算模型

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

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

續(xù)圖12

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

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

圖13和圖14分別給出了燃?xì)庾杂蓢娏髁鲌鰧?duì)稱面靜壓分布云圖和馬赫數(shù)分布云圖。圖13中數(shù)值表示等值線上靜壓值，單位為MPa，圖14中數(shù)值表示等值線上馬赫數(shù)值。

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

自由噴流是經(jīng)典的燃?xì)饬鲌鏊憷細(xì)庥蓢姽苌涑龊?，交替出現(xiàn)收縮擴(kuò)張現(xiàn)象，呈葫蘆狀分布規(guī)律。由圖13和圖14可見，Roe格式與AUSM+格式計(jì)算結(jié)果充分反映了自由噴流的流動(dòng)特性，且兩者表現(xiàn)出了高度一致性。主要區(qū)別在于AUSM+格式所得馬赫數(shù)云圖在燃?xì)饬鬟吘壍戎稻€存在波動(dòng)，這同樣是其標(biāo)量形式的數(shù)值耗散所造成的，結(jié)合之前噴管內(nèi)流場計(jì)算結(jié)果，可見AUSM+格式具有流動(dòng)邊界處易產(chǎn)生數(shù)值波動(dòng)的特點(diǎn)。

5 結(jié) 論

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

b）在噴管內(nèi)流場算例中，由于AUSM+格式在壁面附近垂直壁面方向耗散項(xiàng)數(shù)值幾乎為零，導(dǎo)致計(jì)算結(jié)果在壁面附近產(chǎn)生數(shù)值振蕩，由于類似原因，在自由噴流算例中，其計(jì)算結(jié)果在燃?xì)饬鬟吘壱伯a(chǎn)生了波動(dòng)。

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

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