三角形網(wǎng)格CE／SE方法在帶錐形尾噴管兩相脈沖爆轟發(fā)動機(jī)流場計算中的應(yīng)用

王研艷，翁春生

（南京理工大學(xué) 瞬態(tài)物理國家重點實驗室，江蘇 南京 210094）

0 引 言

脈沖爆轟發(fā)動機(jī)（Pulse Detonation Engine，簡稱PDE）是一種利用脈沖爆轟波來產(chǎn)生周期性沖量的新概念非定常推進(jìn)系統(tǒng)，具有熱效率高、結(jié)構(gòu)簡單、單位燃料消耗率低、工作范圍廣等優(yōu)點。但如何充分利用爆轟產(chǎn)物的高溫高壓能量來提高脈沖爆轟發(fā)動機(jī)推力仍為得到解決，方法之一就是在發(fā)動機(jī)尾部安裝尾噴管。常規(guī)發(fā)動機(jī)尾噴管增加推進(jìn)效能的原理是設(shè)計一個最佳噴管構(gòu)型，將燃燒室內(nèi)排出的高溫高壓燃?xì)獾臒崮鼙M可能多地轉(zhuǎn)化為動能，以增加發(fā)動機(jī)的推力。而脈沖爆轟發(fā)動機(jī)尾噴管中的流動是周期性的、非定常的，存在復(fù)雜激波和膨脹波，甚至存在爆轟波，在同一周期的不同時刻噴管內(nèi)流場不盡相同，因此很難找到一個能適合所有時刻的噴管構(gòu)型，這就使得建立在傳統(tǒng)定常流動基礎(chǔ)上的噴管設(shè)計理論和概念在這里不再適用［1］。Nobuyuki T等人［2］利用 Strang算子分裂法和迎風(fēng)TVD格式數(shù)值模擬了多循環(huán)下帶收斂擴(kuò)張噴管的PDE的內(nèi)外流場，研究了以H2為燃料時噴管進(jìn)出口面積比和燃料填充率對性能的影響，得出在不同飛行馬赫數(shù)條件下存在不同的最佳噴管面積比的結(jié)論。秦亞欣等人［3］對帶不同結(jié)構(gòu)噴管的PDE進(jìn)行了二維軸對稱單循環(huán)數(shù)值研究，指出常規(guī)噴管中擴(kuò)張噴管可以產(chǎn)生最大的瞬時推力峰值和比沖。王志武等人［4］試驗研究了錐形噴管的進(jìn)出口面積比和長徑比對吸氣式脈沖爆轟發(fā)動機(jī)推力的影響，指出不同頻率下各類型噴管對發(fā)動機(jī)的平均推力影響不同，但總體來講收斂型尾噴管的增推效果最好。李建中等人［5］采用基元反應(yīng)模型和顯式迎風(fēng)TVD差分格式模擬了帶共用尾噴管多管脈沖爆轟發(fā)動機(jī)各爆轟管內(nèi)以及共用尾噴管內(nèi)的流場。Caldwell等人［6］實驗研究了不同填充率下引射器擴(kuò)張角度對脈沖爆轟發(fā)動機(jī)流場和推進(jìn)性能的影響。總體而言脈沖爆轟發(fā)動機(jī)尾噴管的最佳設(shè)計與PDE的頻率、飛行馬赫數(shù)、噴管類型、高度、燃料填充率等工作參數(shù)緊密相關(guān)［7－8］，所以至今還未得出非定常脈沖爆轟發(fā)動機(jī)尾噴管的設(shè)計準(zhǔn)則。

CE／SE方法最先由Chang提出，它將時間與空間統(tǒng)一處理，從守恒律積分方程出發(fā)，設(shè)立守恒元和（求解元，使求解格式局部和整體都保證物理意義上的守恒。由于CE／SE方法無需黎曼分解，具有計算格式簡單、精度高、捕捉間斷能力強(qiáng)等優(yōu)點，所以近年來CE／SE方法成為國內(nèi)外學(xué)者研究的焦點之一。Chang等人［9］在最初提出的CE／SE方法的基礎(chǔ)上進(jìn)行了改善，提高了它的計算精度和穩(wěn)定性。Gary C Cheng等人［10］在Chang的基礎(chǔ)上將CE／SE方法進(jìn)一步發(fā)展，針對粘性流采用長扁型網(wǎng)格發(fā)展了CE／SE的幾種計算格式，并進(jìn)行了算例驗證和效率分析。Wang等人［11］運用矩形網(wǎng)格CE／SE方法對兩相脈沖爆轟發(fā)動機(jī)內(nèi)流場進(jìn)行了數(shù)值研究。李昕等人［12］運用矩形網(wǎng)格CE／SE方法求解流體動力學(xué)N－S方程與麥克斯韋方程，實現(xiàn)了對電磁軌道炮中等離子體電樞的數(shù)值研究。劉建文等人［13］考慮到流場的不規(guī)則性運用非結(jié)構(gòu)網(wǎng)格CE／SE方法對單相多管脈沖爆轟發(fā)動機(jī)模型進(jìn)行了數(shù)值研究。

本文在Chang的基礎(chǔ)上推導(dǎo)了帶源項的結(jié)構(gòu)化三角形網(wǎng)格CE／SE方法計算格式，這方面的內(nèi)容在國內(nèi)還未見報道。而此結(jié)構(gòu)化三角形網(wǎng)格CE／SE方法與矩形網(wǎng)格CE／SE方法相比無需坐標(biāo)變換，網(wǎng)格分布均勻；與非結(jié)構(gòu)網(wǎng)格CE／SE方法相比推導(dǎo)過程簡化，計算格式簡單。本文應(yīng)用此CE／SE方法數(shù)值研究帶錐形擴(kuò)張尾噴管的兩相脈沖爆轟發(fā)動機(jī)的內(nèi)外流場。分析尾噴管內(nèi)不同填充狀態(tài)時PDE的內(nèi)外流場及推進(jìn)性能。

1 理論模型

計算采用軸對稱兩相燃燒轉(zhuǎn)爆轟控制方程［14－15］：

式中，α＝0為二維流動；α＝1為軸對稱流動。其中下標(biāo)g、l分別表示氣相和液相；φg、φl為體積分?jǐn)?shù)比，滿足φg＋φl＝1；ρ、u、v、p分別是密度、軸向速度、徑向速度、壓力。E是單位總能，E＝e＋（u2＋v2）。Id、Fd、Qd、Qc的含義和表達(dá)式見參考文獻(xiàn)［15］。

2 三角形網(wǎng)格CE／SE方法

2.1 三角形CE／SE方法中守恒元與求解元的確定

CE／SE方法［9－10，16］將整個空間－時間計算區(qū)域劃分為若干個求解元。在每個求解元內(nèi)，假設(shè)流場的變量是連續(xù)的，并可以用Taylor級數(shù)展開，穿過相鄰求解元的邊界，流場的變量可以不連續(xù)。而在每個網(wǎng)格點對應(yīng)的守恒元上，空間－時間的積分通量是守恒的。考慮如圖1（a）所表示的交替網(wǎng)格，每個三角形的中心點用實心圓和空心圓分別表示其交替網(wǎng)格點。令空間Ω表示E3上的所有網(wǎng)格點M（r，s，n）的集合。令空心圓點表示tn時刻，實心點表示tn＋1／2時刻。如圖1（b）所示，GMBFDHC表示tn時刻，G′M′B′F′D′H′C′表示tn＋1／2時刻，G″M″B″F″D″H″C″表示tn＋1時刻。那么點G′對應(yīng)的求解元S（G′）為GBB′G′，GDD′G′、GCC′G′和M′B′F′D′H′C′；對應(yīng)的守恒元為四棱柱GBFDG′B′F′D′，GDHCG′D′H′C′和GCMBG′C′M′B′。

2.2 計算格式

在方程（1）中，根據(jù)Chang對時間和空間統(tǒng)一處理 的思想［9－10，14，16］，令x1＝x，x2＝y(tǒng)，x3＝t作為三維歐氏空間的3個坐標(biāo)，且令S＝R－；h＝（F，G，U），利用高斯定理對其積分有：

圖1 二維三角形網(wǎng)格CE／SE方法中時空幾何圖形Fig.1 Spatial scheme of two－dimensional CE／SE method using congruent triangular meshes

式中，S（V）是V在E3上任意時空區(qū)域的邊界，h＝（F，G，U）是時間－空間通量流密度矢量。

引入雅克比系數(shù)矩陣A、B：A＝因為F和G均是U的齊次函數(shù)，所以F＝AU，G＝BU。這樣h進(jìn)一步表示為：h＝（AU，BU，U）。

根據(jù)求解元和守恒元的定義，對于任意的（x，y，t）∈SE（r，s，n），U（x，y，t）、F（x，y，t）、G（x，y，t）可用其相應(yīng)的離散量U＊（x，y，t；r，s，n）、F＊（x，y，t；r，s，n）、G＊（x，y，t；r，s，n）近似代替。利用泰勒級數(shù)展開有：

所以方程（2）用離散通量表示如下：

為了方便計算引入一個新的坐標(biāo)系，如圖2。新坐標(biāo)系下取單位長度：Δη＝MH＝和Δξ，其中w＝

圖2 （x，y）和（ξ，η）坐標(biāo)系Fig.2 Cartesian coordinates and new coordinates

在各CE上對式（3）進(jìn)行求解，得到三角形網(wǎng)格CE／SE方法求解格式如下：

當(dāng)r＝4j＋1，s＝2k＋1，或當(dāng)r＝4j＋3，s＝2k（s＝1，2，3…；k＝1，2，3…）時：

當(dāng)r＝4j，s＝2k＋1，或當(dāng)r＝4j＋2，s＝2k（s＝1，2，3…；k＝1，2，3…）時：

上式中，I是單位矩陣。

和Uη＋求解如下［14］：

式中：α是一個可調(diào)的參數(shù)，通常取1或2。

式（6）中：

2.3 源項的處理

方程（1）中等號右側(cè)出現(xiàn)了化學(xué)反應(yīng)源項，由于化學(xué)反應(yīng)的特征時間遠(yuǎn)小于對流的特征時間，其源項是剛性的，因此需要采用龍格－庫塔法對方程（1）中的化學(xué)反應(yīng)源項進(jìn)行求解［14］。具體方法是（1）先考慮除化學(xué)反應(yīng)之外的其他源項，在CE／SE格式下求出（2）將作為初值，用龍格－庫塔法來求解方程＝R。龍格－庫塔法的時間步長為：Δt＝rk，N通常取5～20。

3 計算結(jié)果及分析

采用三角形網(wǎng)格CE／SE方法對帶擴(kuò)張尾噴管PDE的內(nèi)外流場進(jìn)行耦合計算，其軸對稱計算區(qū)域如圖3，AEFGHA為PDE爆轟管和尾噴管平面計算區(qū)域，尾噴管擴(kuò)張角度為30°。EDCBHGFE為PDE外流場計算區(qū)域。

初值條件：模擬的初始時刻PDE有兩種填充方案：方案1：爆轟管和錐形管作為整體，均填充化學(xué)當(dāng)量比的汽油空氣混合物；方案2：將錐形管僅作噴管使用，即爆轟管內(nèi)填充化學(xué)當(dāng)量比的汽油空氣混合物，但錐形管內(nèi)只填充空氣。同時外流場區(qū)域則布滿均勻的空氣。管內(nèi)外氣體速度均為0m／s，溫度為298K，壓力為0.1MPa，管內(nèi)液相汽油液滴半徑為50μm。

為了驗證運用三角形網(wǎng)格CE／SE方法模擬PDE流場的準(zhǔn)確性，先對不帶尾噴管的PDE內(nèi)外流場進(jìn)行數(shù)值計算，并與實驗結(jié)果進(jìn)行對比驗證。比較發(fā)現(xiàn)模擬得到的圖線變化趨勢和實驗結(jié)果吻合得很好（圖4）。數(shù)值計算的爆轟波波陣面很陡，只有2～4個網(wǎng)格，說明此三角形網(wǎng)格CE／SE方法能夠精確捕獲爆轟波等強(qiáng)間斷。

圖3 PDE模型計算區(qū)域示意圖Fig.3 Scheme of the PDE computational domain

圖4 爆轟管尾部壓力隨時間變化曲線Fig.4 Pressure contribution vs.time at the end of PDE

圖5是采用填充方案1時PDE的內(nèi)外流場壓力云圖，圖6為對應(yīng)的PDE中心軸線上壓力隨x變化曲線。聯(lián)系圖6分析圖5，由圖5（a）知PDE內(nèi)已形成穩(wěn)定的爆轟波，爆轟波后是從封閉端發(fā)出的Taylor膨脹波扇，以滿足封閉端速度為零的條件。在封閉端和Taylor膨脹波波尾之間是均勻的平臺壓力區(qū)。從圖6中曲線5、6可以看出形成穩(wěn)定爆轟時爆轟壓力為1.97MPa，均勻區(qū)平臺壓力為0.48MPa。圖5（b）中爆轟波傳播到錐形管內(nèi)，燃料燃燒所釋放的熱量維持爆轟波在管內(nèi)繼續(xù)傳播。錐形管的擴(kuò)張作用使得爆轟壓力略有下降，爆轟波波面出現(xiàn)向右彎曲現(xiàn)象。由圖6曲線7知，至t＝0.58ms時，中心軸線上爆轟壓力值降為1.85MPa。圖5（c）中爆轟波傳出錐形管，失去能量支持后退化為球形激波在外流場擴(kuò)散。錐形管出口處產(chǎn)生一系列膨脹波向上游傳播，管內(nèi)受到外流場的干擾壓力降低。圖5（d），外流場球形波向左傳播至0.61m，向右傳播到1.44m，錐形管出口外側(cè)出現(xiàn)明顯的渦環(huán)區(qū)域，錐形管外側(cè)背壓發(fā)生變化，發(fā)動機(jī)推進(jìn)性能受到影響。發(fā)動機(jī)管內(nèi)膨脹波傳播至0.36m，管內(nèi)壓力降低，平臺壓力區(qū)縮短，由圖6曲線16可知隨著排氣過程的進(jìn)行，t＝2.91ms時發(fā)動機(jī)頭部壓力降低到0.14MPa。

圖5 采用填充方案1時PDE內(nèi)外流場壓力分布云圖Fig.5 Pressure distribution in and out of PDE with fuel filled in the conical tube

圖6 采用填充方案1時PDE內(nèi)外流場軸線上壓力分布曲線Fig.6 Pressure distribution on the symmetry axis of PDE with fuel filled in conical tube

為了進(jìn)一步分析錐形管內(nèi)外的流動情況，圖7給出了采用填充方案1時PDE內(nèi)外流場的速度矢量圖。圖7（a）中速度沿錐形管壁面分布解釋了擴(kuò)張管內(nèi)爆轟壓力降低的原因。圖7（b）氣流流出擴(kuò)張管后在出口外側(cè)出現(xiàn)速度渦環(huán)，這是因為當(dāng)曲面爆轟波上下邊界傳出擴(kuò)張管后，管內(nèi)欠膨脹的高溫高壓氣體尾隨激波進(jìn)入大氣環(huán)境中，其斜壓作用使得管口外側(cè)出現(xiàn)渦環(huán)。

圖7 采用填充方案1時PDE內(nèi)外流場不同時刻速度矢量圖Fig.7 Velocity vector distribution in and out of PDE with fuel filled in conical tube

圖8是采用填充方案2時PDE內(nèi)外流場的壓力分布云圖，圖9為對應(yīng)的PDE中心軸線上的壓力隨x變化曲線。聯(lián)系兩者分析：圖8（a）中爆轟波傳播到噴管后退化為無化學(xué)反應(yīng)的主激波，受噴管擴(kuò)展壁面產(chǎn)生的膨脹波影響，在靠近壁面處發(fā)生彎曲，形成曲面激波。此時錐形管內(nèi)流場比方案一中更加復(fù)雜。從圖9的曲線7、8、9可以看出，隨著主激波在擴(kuò)張噴管中的傳播，激波越來越弱，波后氣體的壓力越來越小。至t＝0.59ms時，波后壓力值為1.58MPa。圖8（b）中主激波傳出噴管，以球形激波向外流場擴(kuò)散，噴管內(nèi)已存在的膨脹波系和外流場新產(chǎn)生的膨脹波共同作用，降低管內(nèi)壓力。對比圖6和圖9發(fā)現(xiàn)，方案二的排氣過程比方案一時快。具體體現(xiàn)在，t＝1.27ms時，方案一的膨脹波降低PDE內(nèi)部平臺壓力至x＝0.34m處；而方案二中PDE內(nèi)部x＝0.30m處的平臺壓力已經(jīng)開始下降。

圖8 采用填充方案2時PDE內(nèi)外流場不同時刻壓力分布云圖Fig.8 Pressure distribution in and out of PDE at different times，without fuel in conical tube

圖9 采用填充方案2時PDE內(nèi)外流場軸線上壓力分布曲線Fig.9 Pressure distribution on the symmetry axis of PDE at different times，without fuel in conical tube

圖10是爆轟波傳入錐形管后中心軸線上的馬赫數(shù)分布曲線，比較（a、b）圖發(fā)現(xiàn)當(dāng)錐形管內(nèi)不填充燃料時，爆轟波傳入錐形管內(nèi)沒有化學(xué)反應(yīng)的支持退化為激波，但是波峰位置處的馬赫數(shù)反而變大，分析其原因：在錐形管內(nèi)不填充燃料時波峰位置處壓力變小，速度值也變小，但前者的變化幅度遠(yuǎn)大于后者，導(dǎo)致當(dāng)?shù)芈曀僮冃。移錅p小量大于速度減小量。

圖10 爆轟波傳入錐形管后中心軸線上馬赫數(shù)分布曲線Fig.10 Mach number distribution on the symmetry axis of PDE at different moments

考慮到脈沖爆轟發(fā)動機(jī)非定常特性，論文中采用表面力積分來分析帶錐形尾噴管PDE的推進(jìn)性能［17］。圖11為不同填充方案時PDE瞬時推力F、沖量I、燃料比沖Ispf隨時間變化曲線。圖中曲線1對應(yīng)填充方案1；曲線2對應(yīng)填充方案2。分析三幅圖發(fā)現(xiàn)，當(dāng)錐形擴(kuò)張管內(nèi)不填充燃料時，它起到噴管的作用，使爆轟管內(nèi)排出的高溫高壓燃?xì)庠阱F形管內(nèi)得到膨脹，轉(zhuǎn)化為動能，從而增加了發(fā)動機(jī)的推力，比較圖11（c）圖中兩條曲線發(fā)現(xiàn)，至t＝4.8ms時此方案發(fā)動機(jī)的燃料比沖比另一方案的增加了15.5%，說明這種方案提高了燃料的利用率；但當(dāng)錐形擴(kuò)張管內(nèi)填充燃料時，它作為發(fā)動機(jī)爆轟的一部分，充分利用了脈沖爆轟發(fā)動機(jī)的空間，增加了發(fā)動機(jī)的瞬時推力和沖量，比較圖11（b）圖中兩條曲線發(fā)現(xiàn)，至t＝4.8ms時此方案發(fā)動機(jī)的沖量比另一方案增加了11.1%。顯然，當(dāng)我們考慮在每個周期中增加PDE的推進(jìn)性能時，需要利用所有可以利用的空間，此時可以考慮將其錐形尾噴管作為發(fā)動機(jī)的一部分，在其中填充或者部分填充燃料，從而增加發(fā)動機(jī)的總體推力。

圖11 錐形管不同方案時PDE的瞬時推力、沖量和燃料比沖比較Fig.11 The comparison of thrust，impulse and fuel－specific impulse of PDEs with different filling conditions in conical tube

4 結(jié) 論

本文推導(dǎo)了帶源項三角形網(wǎng)格CE／SE方法，數(shù)值求解帶錐形尾噴管脈沖爆轟發(fā)動機(jī)的內(nèi)外流場。計算結(jié)果表明：

（1）該三角形網(wǎng)格CE／SE方法可以有效捕捉到直管和變截面管中的爆轟波、激波等強(qiáng)間斷，很好地刻畫了PDE內(nèi)外流場的細(xì)節(jié)。

（2）將錐形擴(kuò)張管作為脈沖爆轟發(fā)動機(jī)的噴管而不填充燃料時，擴(kuò)張管中壓力峰值大幅下降，但馬赫數(shù)峰值增加，發(fā)動機(jī)的排氣速度增加，此時PDE可以得到更大的燃料比沖，發(fā)動機(jī)對燃料的利用率提高。

（3）而將錐形擴(kuò)張管與爆轟管當(dāng)成整體全部填充燃料時，錐形管中流場更簡單，壓力峰值比爆轟管中略有下降，此時發(fā)動機(jī)充分利用所有空間，瞬時推力和沖量增大，總體推進(jìn)性能提高。

［1］BROPHY C M，DAUSEN D F，SMITH L R，et al.Fluidic nozzles for pulse detonation combustors［R］.AIAA Paper，2012－1035.

［2］NOBUYUKI T，YUICHIRO K，HAYASHI A K，et al.Numerical study and performance evaluation for pulse detonation engine with exhaust nozzle［R］.AIAA Paper，2009－5315.

［3］秦亞欣，于軍力，高歌.脈沖爆震發(fā)動機(jī)噴管性能數(shù)值分析［J］.航空動力學(xué)報，2010，25（2）：366－372.

［4］WANG Z W，YAN C J.Experimental investigation of nozzle effects on a two－phase valveless air－breathing pulse detonation engine［R］.AIAA Paper，2008－991.

［5］李建中，王家驊，王春，等.共用尾噴管多管脈沖爆震發(fā)動機(jī)數(shù)值模擬研究［J］.空氣動力學(xué)學(xué)報，2008，26（1）：96－100.

［6］CALDWELL N，GUTMARK E，HOKE J，et al.Investigation of fundamental processes leading to pulse detonation engine／ejector thrust augmentation ［R］.AIAA Paper，2008－116.

［7］曾昊，何立明，章雄偉，等.噴管收斂－擴(kuò)張角對爆震發(fā)動機(jī)性能影響分析［J］.推進(jìn)技術(shù)，2011，32（1）：97－102.

［8］KAILASANATH K.A review of research on pulse detonation engine nozzles［R］.AIAA Paper，2001－3932.

［9］CHANG S C.A new approach for constructing highly stable high order CESE schemes［R］.AIAA Paper，2010－543.

［10］CHENG G C，VENKATACHARI B S，CHANG C L，et al.Comparative study of different numerical approaches in space－time CE／SE framework for high－fidelity flow simulations［J］.ComputersandFluids，2011，45：47－54.

［11］WANG G，ZHANG D L，LIU K X，et al.An improved CE／SE scheme for numerical simulation of gaseous and two－phase detonations［J］.ComputersandFluids，2010，39：168－177.

［12］LI X，WENG C S.2－D viscous MHD simulation of plasma armatures by the CE／SE method［J］.Chinese ScienceBulletin，2009，54（10）：1641－1647.

［13］劉建文，趙書苗，鐘誠文，等.CE／SE方法在多管爆轟流場并行計算中的應(yīng)用［J］.爆炸與沖擊，2008，28（3）：229－235。

［14］翁春生，王浩.計算內(nèi)彈道學(xué)［M］.北京：國防工業(yè)出版社，2006：317－326.

［15］馬丹花，翁春生.爆震管內(nèi)擾流片對爆震波影響的數(shù)值分析［J］.推進(jìn)技術(shù)，2011，32（3）：425－430.

［16］CHANG S C，WANG X Y，CHOW C Y.The spacetime conservation element and solution element method：a new high－resolution and genuinely multidimensional paradigm for solving conservation laws［J］.JournalofComputationalPhysics，1999，156：89－136.

［17］VENKAT E T，ANTHONY J D，NOBUYUKI T，et al.Performance of a pulse detonation engine under subsonic and supersonic flight conditions［R］.AIAA Paper，2007－1245.