一維與二維爆轟傳播的時空關聯(lián)特性數(shù)值研究1)

張文碩 楊鵬飛 姜宗林 劉云峰

(中國科學院力學研究所 高溫氣體動力學國家重點實驗室,北京 100190)

(中國科學院大學工程科學學院,北京 100049)

引言

爆轟是一種在預混可燃氣體中以超聲速傳播的極端燃燒現(xiàn)象,具有自持、自組織和增壓燃燒特征[1].在自然界中,爆轟可能發(fā)生于煤礦爆炸和化工廠事故,在Ia型超新星爆炸研究和航空航天新型動力的探索領域也得到越來越廣泛的關注[2-4].19 世紀末以來,人們從各個角度探索了氣體爆轟的起爆、傳播規(guī)律和精細結(jié)構(gòu),在理論和實驗認知方面都取得了重大進展.近年來,姜宗林[5]在前人研究的基礎上,結(jié)合自己的研究成果,提出氣體爆轟物理的統(tǒng)一框架理論,總體地概括了氣相爆轟波起爆和傳播的主要規(guī)律和關鍵物理參數(shù)的依存與競爭關系,為進一步分析爆轟現(xiàn)象建立了一個具有啟示性的理論基礎.

斜爆轟的研究起步較晚,由于應用需求推動其得到了越來越多的關注.斜爆轟波是一種駐定爆轟波(standing detonation wave),當高速可燃氣體經(jīng)過尖劈、鈍體等結(jié)構(gòu)時產(chǎn)生斜激波,波后溫度升高引發(fā)燃燒反應,燃燒帶與斜激波耦合構(gòu)成斜爆轟波.Li 等[6]的早期數(shù)值模擬研究揭示了斜爆轟波的基本起爆結(jié)構(gòu).Figueira da Silva 和Deshaies[7]通過數(shù)值模擬發(fā)現(xiàn)了兩種差異顯著的起爆結(jié)構(gòu),分別稱為光滑型和突變型.Teng 和Jiang[8]提出斜激波與斜爆轟波之間的角度差可以作為兩種起爆結(jié)構(gòu)的判據(jù).Choi 等[9]通過高分辨率的數(shù)值模擬發(fā)現(xiàn)斜爆轟波面也能夠演化出類似正爆轟的胞格結(jié)構(gòu),Gui 等[10]根據(jù)波面形狀特征將斜爆轟波面結(jié)構(gòu)分為3 個區(qū)域.Teng 等[11]的數(shù)值模擬結(jié)果表明高活化能、小楔角條件下更容易出現(xiàn)胞格結(jié)構(gòu).Fusina 等[12]模擬了斜爆轟波面在非均勻來流條件下的響應,說明斜爆轟波具有一定的抗干擾能力.Yang 等[13]的研究表明斜爆轟波在非定常來流條件下也能保持一定的穩(wěn)定性.

斜爆轟的形成和發(fā)展過程與正爆轟有一定的相似性,但其起爆結(jié)構(gòu)和胞格演化規(guī)律存在顯著區(qū)別,明確兩者的關聯(lián)和差異對深入理解和應用斜爆轟現(xiàn)象有重要意義.Ghorbanian 和Sterling[14]發(fā)現(xiàn)雙楔誘導斜爆轟的波系結(jié)構(gòu)與變速活塞的t-y圖有相似的特征,Daimon 和Matsuo[15]數(shù)值模擬比較了半無限長楔誘導斜爆轟與活塞驅(qū)動正爆轟的起爆結(jié)構(gòu),并通過活塞運動類比出多種楔面誘導的斜爆轟起爆結(jié)構(gòu).Higgins[16]指出,利用斜爆轟與活塞驅(qū)動正爆轟的相似性,可以直觀地理解沖壓加速器燃燒室內(nèi)的波系結(jié)構(gòu).上述這些研究結(jié)果表明,斜爆轟是駐定的,正爆轟是非定常的,考慮正爆轟的時間維度,兩種爆轟的流場結(jié)構(gòu)存在一定的時空關聯(lián)關系,即時空坐標下的非定常正爆轟與空間坐標下的定常斜爆轟有一定的相似性.這種時空關聯(lián)研究有助于理解斜爆轟與正爆轟在本質(zhì)上的聯(lián)系,深化對爆轟現(xiàn)象的認識.

斜爆轟結(jié)構(gòu)具有多尺度特征[17],包含化學反應的控制方程求解存在剛性問題[18],斜爆轟波的起爆、胞格等局部精細結(jié)構(gòu)與流場中其他區(qū)域相比要求更高的時空分辨率,高超聲速邊界處理也有很高的計算方法要求.因此,對斜爆轟的二維流場開展大尺度、高精度的數(shù)值模擬需要較大的計算量和較復雜的計算方法.一維正爆轟的計算量小,而且有經(jīng)典理論和較豐富的實驗驗證.如果能夠借助正爆轟與斜爆轟的時空關聯(lián)特性,通過一維數(shù)值模擬獲得二維斜爆轟流場的主要規(guī)律,不僅具有工程應用價值,也能加深對斜爆轟起爆機理等基礎問題的理解.

然而,過去的研究[14-16]側(cè)重于利用斜爆轟與正爆轟的相似關系來理解斜爆轟現(xiàn)象,采用了理想的半無限長楔構(gòu)型和抽象的化學反應模型,主要開展定性的討論,缺乏定量研究,目前尚不清楚這種時空關聯(lián)的精確度,難以投入工程應用.為了深入探索正爆轟與斜爆轟現(xiàn)象的時空關聯(lián)特性,從定量上獲得兩種爆轟現(xiàn)象的關聯(lián)規(guī)律,本文在文獻[15]工作的基礎上,利用多組分Euler 方程和基元反應模型模擬了有限長楔誘導定常斜爆轟波和相同過驅(qū)動度條件下的活塞驅(qū)動非定常正爆轟波,通過兩者流場特征的對比研究,闡述兩種爆轟相似性和差異,并給出了時空關聯(lián)關系.

1 數(shù)值計算方法

一維爆轟理論是建立在定常流動的基礎上的,其預測結(jié)果和實驗符合良好,原因在于黏性對爆轟波流動結(jié)構(gòu)的影響較小[19].本文的重點在于時空關聯(lián)問題,因此采用多組分Euler 方程統(tǒng)一描述一維和二維爆轟波的起爆和傳播過程.該方程可表述如下

對于一維爆轟的數(shù)值模擬

式中ρi(i=1,2,···,n) 為第i個組分的密度,n為化學反應模型中包含的組分數(shù).混合氣體的總密度為,u和v分別為x方向和y方向的氣體速度,e為氣體單位體積的總能,與氣體的焓、壓力和動能有關,表達式為

式中hi為第i個組分單位質(zhì)量的焓值,本文將混合氣體的每一個組分處理為熱完全氣體,其定壓比熱Cpi為溫度的單值函數(shù),由多項式擬合得到

進一步積分得到各組分的焓值

式中的擬合常數(shù)aji(j=1,2,···,7)和b1i參考文獻[20].結(jié)合理想氣體狀態(tài)方程和道爾頓分壓定律,混合氣體的壓力p為

式中Ri為第i個組元的氣體常數(shù),T為混合氣體的溫度.式(1)中的化學反應源項Sc包括n個組分單位體積的質(zhì)量生成率,根據(jù)質(zhì)量作用定律

式中Mi為第i個組分的物質(zhì)的量,nr為化學模型包含的反應數(shù),,,αik分別為第k個反應中第i個組分在生成物一側(cè)的計量數(shù)、在反應物一側(cè)的計量數(shù)和第三體增強系數(shù),ci為第i個組分的物質(zhì)的量濃度,kfk,kbk分別為第k個反應的正、逆反應速率,一般通過Arrhenius 定律計算kfk

式中Ak,nk,Eak分別為第k個反應的頻率因子、溫度指數(shù)和活化能,Ru為普適氣體常數(shù).進一步結(jié)合反應平衡常數(shù)可計算出kbk

預混可燃氣體為化學當量的氫氣/空氣混合氣,反應過程采用高壓修正的Jachimowski 機理[21],包含H2,O2,H,O,OH,HO2,H2O2,H2O,N2共9 種組分和19 個化學反應.二維數(shù)值模擬采用非結(jié)構(gòu)四邊形網(wǎng)格,平均網(wǎng)格間距為50 μm,通過有限體積法和二階TVD 格式離散求解,界面通量通過HLLC(Harten-Lax-van Leer contact)近似黎曼求解器[22]計算.一維數(shù)值模擬采用結(jié)構(gòu)網(wǎng)格,網(wǎng)格間距為10 μm,采用頻散控制耗散(dispersion controlled dissipation,DCD)格式[23]離散求解,時間推進采用3 階Runge-Kutta法,時間步長由CFL 數(shù)控制,計算中CFL 數(shù)取0.1.

2 計算模型的驗證

選取兩個已知解析解的問題分別對一維和二維計算模型進行驗證.圖1 為一維爆轟波在光滑直管道中的傳播問題,計算得到不同時刻的壓力分布.31.4 μs 后,壓力波形幾乎不變,通過捕捉激波位置變化得到爆轟波速為1970.71 m/s,GASEQ 軟件計算出的Chapman-Jouguet (CJ)爆速為1979.33 m/s,兩者的相對誤差為?0.436%,說明一維計算模型比較可靠.

圖1 典型一維爆轟波的數(shù)值模擬結(jié)果Fig.1 Typical pressure profile of one-dimensional detonation in the numerical simulation

圖2 為半無限長楔誘導的斜爆轟問題,當楔面足夠長時,可以忽略起爆結(jié)構(gòu)的影響,在下游測量爆轟角.分別在10 個角度下進行了數(shù)值模擬,均在楔面上得到了駐定的斜爆轟波(見表1).通過化學反應平衡的斜激波求解方法[24]得到理論解,數(shù)值模擬和理論求解采用相同的化學反應模型,表1 中兩者的相對誤差幾乎都在1%以下,其中相對誤差

圖2 半無限長楔面數(shù)值模擬的爆轟角與爆轟極曲線對比Fig.2 Detonation angle of numerical simulation vs.detonation polar diagram

表1 半無限長楔面數(shù)值模擬的爆轟角與爆轟極曲線對比Table 1 Detonation angle of numerical simulation vs.that of detonation polar diagram under various wedge angles

從圖2 可以看出數(shù)值模擬得到的爆轟角βnum隨著楔面折角θ的變化規(guī)律與理論解符合,說明二維計算模型比較可靠.

3 計算結(jié)果分析

本文重點研究有限長楔誘導的定常斜爆轟波與非定常活塞驅(qū)動的正爆轟波的時空關聯(lián)特性,計算域及初始條件設置見圖3.二維數(shù)值模擬的計算域為12 cm×12 cm 的方形區(qū)域,楔角θ為25°,楔面長度2 cm,楔面末端出現(xiàn)折角,在折角下游,壁面與來流方向平行.

沿用Daimon 和Matsuo[15]的近似方法,建立壁面參考系,時空變換關系表述為:壁面法向流動對應一維流場,壁面切向流動對應一維時間推進.一維算例如圖3(b)所示,將活塞的非定常運動分為兩個階段.第一階段初始時刻,活塞以1450.73 m/s 勻速運動,該速度高于CJ 狀態(tài)的氣體速度1130.59 m/s,誘導產(chǎn)生過驅(qū)動正爆轟波.數(shù)值模擬中,正爆轟波的過驅(qū)動度與25°楔面產(chǎn)生的過驅(qū)動斜爆轟波相同,均為1.235,過驅(qū)動度f定義為

圖3 計算模型示意圖Fig.3 Schematic diagram of numerical model

式中Mn為來流在爆轟波面法向的速度對應的馬赫數(shù),MCJ為相同來流條件下CJ 狀態(tài)的激波馬赫數(shù)[25].來流為預混氫氣/空氣混合氣,來流壓力p0= 101.3 kPa,來流溫度T0=300 K,當量比為1,對應的MCJ=4.843.當正爆轟波面與活塞的距離和斜爆轟波面與楔面的距離相同時,第一階段結(jié)束,活塞停止.第二階段中,活塞保持不動,爆轟波繼續(xù)向前傳播.根據(jù)時空變換關系,上述活塞運動方式與二維有限長楔構(gòu)型的壁面變化規(guī)律一致.為闡述斜爆轟與正爆轟的時空關聯(lián)特性將分別從波系結(jié)構(gòu)、壁面參數(shù)、波形變化和時間尺度4 個方面對二維和一維數(shù)值模擬結(jié)果進行對比.

3.1 波系結(jié)構(gòu)

圖4(a)～ 圖4(d)為二維有限長楔誘導的斜爆轟數(shù)值結(jié)果,依次為楔面上的壓力、溫度流場和折角后的壓力、溫度流場,圖4 中坐標原點在楔面起點處,x軸與來流方向平行,y軸與來流方向垂直.在楔面上,為了方便與一維計算結(jié)果比較,對坐標軸進行了旋轉(zhuǎn),旋轉(zhuǎn)后x軸與楔面重合,y軸與楔面垂直.圖4(a)和圖4(b)采用了旋轉(zhuǎn)后的坐標系,來流經(jīng)過楔面被壓縮,壁面附近的壓縮波匯聚[26]為斜激波(oblique shock wave,OSW),波后溫度約為1260 K,誘發(fā)可燃氣體自點火.經(jīng)過一定的起爆區(qū)距離,可燃氣體發(fā)生劇烈的放熱反應,燃燒反應帶角度快速抬升并與斜激波相交,燃燒和激波面強耦合誘導出過驅(qū)動的斜爆轟(oblique detonation wave,ODW).同時三波點下游出現(xiàn)滑移線和反射波,反射波在壁面上再次反射產(chǎn)生膨脹波,膨脹波在滑移線上透射與過驅(qū)動斜爆轟波相交.在本文的計算條件下,斜爆轟波面較為穩(wěn)定,未發(fā)展出斜爆轟胞格結(jié)構(gòu)[10].在圖4(c)和圖4(d)中,折角處出現(xiàn)膨脹扇,氣流通過膨脹扇后與壁面平行.一段距離后,膨脹扇追上過驅(qū)動斜爆轟波使其激波角減小,最終在計算域下游形成近CJ 斜爆轟波,其激波角為37.63°,與CJ 斜爆轟波的激波角37.03°接近.氣流通過近CJ 斜爆轟波時熵增最小,總壓損失最小[27],此時燃燒室出口總壓恢復系數(shù)最高,有利于產(chǎn)生更大的推力,這是有限長楔構(gòu)型在斜爆轟發(fā)動機中應用的一大優(yōu)勢.

圖4 不同工況下壓力及溫度的數(shù)值模擬云圖Fig.4 Numerical contours of pressure and temperature under different conditions

圖4 不同工況下壓力及溫度的數(shù)值模擬云圖(續(xù))Fig.4 Numerical contours of pressure and temperature under different conditions (continue)

圖4(e)～ 圖4(h)為一維數(shù)值模擬結(jié)果,將時間作為一個維度與二維云圖比較.可以看出,一維數(shù)值模擬在時空坐標下的溫度、壓力流場具有二維流場的全部特征,包括激波、燃燒波、三波點、過驅(qū)動爆轟波、滑移線、反射波、膨脹扇等.這初步驗證了一維非定常活塞驅(qū)動正爆轟與二維有限長楔誘導定常斜爆轟的時空關聯(lián)特性.在一維流場中,活塞壓縮產(chǎn)生的高溫高壓區(qū)域從壁面向外擴張,活塞表面附近的氣體加熱時間最長,因此放熱反應發(fā)源于壁面向前傳播(interior detonation).從時空關聯(lián)的角度,可以判斷二維楔面上燃燒波從壁面引發(fā)的現(xiàn)象并不是邊界層黏性效應的影響,而是壁面壓縮效應的結(jié)果.

3.2 壁面參數(shù)

進一步對二維和一維數(shù)值模擬結(jié)果的時空關聯(lián)作定量研究,沿用Daimon 和Matsuo[15]的做法,用楔面末端氣流速度將楔面x坐標轉(zhuǎn)化為時間t,用靠近壁面的氣流出口速度將壁面x坐標轉(zhuǎn)化為時間t,從而與一維壁面參數(shù)進行對比,如圖5 所示.需要指出,本研究不考慮黏性,壁面上采用滑移邊界條件.考慮黏性影響時,楔面上發(fā)展出邊界層,邊界層內(nèi)的氣流在黏性的阻滯下速度較低,此時利用靠近壁面的氣流出口速度進行時空關聯(lián)并不合適.研究表明,黏性對斜爆轟流場的影響主要局限在邊界層內(nèi)部,對斜爆轟波系結(jié)構(gòu)和主流區(qū)的影響很小[28-29].因此,開展黏性數(shù)值模擬時,應選取邊界層外的氣流出口速度做時空轉(zhuǎn)換,此時本文研究的時空關聯(lián)特性仍是成立的.

圖5(a)為楔面壓力、溫度與活塞運動階段表面的壓力、溫度對比.可燃氣體在激波后自點火產(chǎn)生的爆轟波造成溫度、壓力躍升,產(chǎn)生局部高溫高壓.氣流通過壁面反射的膨脹波后,溫度、壓力有所下降,之后保持平穩(wěn).圖中兩算例,氣體自點火的時間有一定差別,壓力、溫度線型有0.4 μs 左右的相位差,其原因在3.4 節(jié)中進行分析.其他的壁面參數(shù)相差很小,可以認為兩者定量符合.圖5(b)為折角下游壁面壓力、溫度與活塞停止后表面壓力、溫度的對比.在膨脹扇的作用下,圖5(b)的入口溫度、壓力低于圖5(a)的出口溫度、壓力,膨脹扇下游溫度、壓力保持平穩(wěn).在二維數(shù)值模擬中,折角處除膨脹扇外還產(chǎn)生一道二次激波,在二次激波作用下,壁面溫度、壓力有所上升,而一維流場沒有受到二次激波的影響.

圖5 壁面壓力、溫度變化Fig.5 Wall pressure and temperature profiles

3.3 波形變化

圖6(a)和圖6(c)顯示了壁面與激波面距離相等的一系列位置處,二維和一維壓力波形的變化,可分辨出圖4(a)、圖4(c)、圖4(e)、圖4(g)中的波系結(jié)構(gòu)變化特征,包括激波波形、起爆波形、過驅(qū)動爆轟波形和近CJ 爆轟波形.各位置處二維與一維波形定性相同,呈現(xiàn)出相同的變化規(guī)律,壁面附近壓力平臺的數(shù)值很接近,但活塞驅(qū)動正爆轟的von Neumann壓力峰值比斜爆轟高,在近CJ 狀態(tài)下,二維von Neumann 壓力峰值約為23,而一維壓力峰值約為28,這可能是由于一維數(shù)值模擬空間分辨率更高,得到的壓力峰值更準確.同時,一維數(shù)值模擬結(jié)果中過驅(qū)動爆轟向CJ 爆轟的演化比二維慢,在y=0.06 m后,有限長楔上的壓力波形基本保持不變,而一維活塞流場在y=0.09 m處的壓力波形與y=0.075 m 處相比仍存在明顯差別,其原因在3.4 節(jié)中進行分析.

圖6(b)、圖6(d)顯示了與圖6(a)、圖6(c)對應的溫度波形變化,同樣可以分辨出圖4(b)、圖4(d)、圖4(f)、圖4(h) 中的結(jié)構(gòu)變化特征,包括激波波形、起爆波形、過驅(qū)動爆轟波形和近CJ 爆轟波形.由于二維和一維數(shù)值模擬的起爆距離有一定差別,而起爆的時間尺度很短,起爆階段中相同位置處的波形有一定差異.除此之外,二維與一維的溫度波形變化定量一致.

圖6 壓力、溫度波形變化Fig.6 Variation of pressure and temperature profile

3.4時間尺度

活塞運動中產(chǎn)生的激波、過驅(qū)動爆轟波和近CJ 爆轟波處于定常狀態(tài),而這些狀態(tài)之間的轉(zhuǎn)變是非定常過程,經(jīng)過時空變換后對應二維流場的結(jié)構(gòu)特征.圖7 顯示了二維數(shù)值模擬結(jié)果中斜激波角β隨坐標x的變化,從圖7 中可以看出波系結(jié)構(gòu)變化的各個階段.與斜爆轟波相比,斜激波與來流方向的夾角較小,約為32.20°,與理論值32.73°接近.斜激波與燃燒波相交形成三波點,三波點下游波后溫度迅速升高,為了匹配波面前后的壓力、溫度,斜激波變強,激波角增大,在三波點附近激波角達到最大值48.82°.三波點處的反射波在壁面上再次反射,產(chǎn)生膨脹波,膨脹波透過滑移線作用在激波面上,使激波角減小,形成較穩(wěn)定的過驅(qū)動斜爆轟波,本文計算條件下過驅(qū)動爆轟波的激波角理論值為42.28°,數(shù)值模擬結(jié)果與理論值符合.楔面末端折角處的膨脹扇與過驅(qū)動斜爆轟波作用,進一步使激波角降低,計算條件下CJ 狀態(tài)的激波角為37.03°,由圖7 可以看出在計算域下游激波角趨近于理論CJ 激波角,達到近CJ 狀態(tài).一維算例中激波馬赫數(shù)隨時間的變化也具有上述變化特征.

圖7 有限長楔面上的斜激波角的變化Fig.7 Variation of oblique shock wave angle over finite wedge

表2 將一維數(shù)值模擬中非定常過程的弛豫時間與時空變換后的二維計算結(jié)果作了對比.其中t1為激波轉(zhuǎn)爆轟時間,在計算條件下,可燃混合氣的點火延遲時間為2.59 μs,楔面和活塞數(shù)值模擬中t1分別為2.87 μs 和2.58 μs,根據(jù)Teng 等[30]的分類,這種情況屬于動力學控制(kinetics-controlled)起爆,起爆的時間尺度與點火延遲時間相當.t2為起爆和產(chǎn)生穩(wěn)定過驅(qū)動爆轟波之間的時間,即壁面反射膨脹波的作用時間,二維數(shù)值模擬得到的t2約為一維的2 倍.流動通過楔面折角或活塞停止時,產(chǎn)生的膨脹扇以當?shù)芈曀賯鞑?由于過驅(qū)動爆轟波波后是亞聲速的,膨脹扇經(jīng)過t3后追上過驅(qū)動爆轟波使其弱化,二維數(shù)值模擬得到的t3約為一維的1/2.t4為過驅(qū)動爆轟波與近CJ 爆轟波之間的弛豫時間,二維數(shù)值模擬得到的t4約為一維的1/3.目前對近CJ 狀態(tài)還沒有定量的界定,本文中近CJ 狀態(tài)定義為氣體流過波面時總壓恢復σ達到CJ 狀態(tài)總壓恢復σCJ的95%以上時的狀態(tài)范圍,如圖8 所示.

表2 爆轟波各發(fā)展階段的時間尺度Table 2 Time duration of different stages in the development of detonation wave

圖8 近Chapman-Jouguet 斜爆轟波的定量定義Fig.8 Quantitative definition of near Chapman-Jouguet oblique detonation wave

除起爆時間t1外,二維與一維數(shù)值模擬結(jié)果在非定常尺度上存在顯著差別.根據(jù)理論分析,只有波后氣體在壁面切向的運動速度不變的情況下,二維楔面誘導斜爆轟波后的氣體運動方程才與一維活塞驅(qū)動正爆轟的波后氣體運動方程一致.然而在上述非定常過程中,均存在波系相互作用,氣體速度在壁面法向的梯度較大,選擇一個特定的速度將二維空間坐標x轉(zhuǎn)換為時間坐標t時存在較大誤差.因而通過活塞數(shù)值模擬無法準確預測楔面算例的非定常時間,這種差異本質(zhì)上是流動的二維效應造成的.

4 結(jié)論

經(jīng)過時空變換,有限長楔誘導的定常斜爆轟與非定?；钊?qū)動的正爆轟具有相同的波系結(jié)構(gòu),雖然兩者的非定常時間尺度存在一定差異,但爆轟參數(shù)、波形結(jié)構(gòu)演化規(guī)律在定性和定量上均有較好的一致性,體現(xiàn)出二維斜爆轟與一維非定常正爆轟現(xiàn)象存在著時空關聯(lián)特性.

通過一維活塞驅(qū)動正爆轟計算,經(jīng)過時空變換方法,能夠獲得相同過驅(qū)動度下的二維斜爆轟流場.由于一維活塞驅(qū)動正爆轟理論完備,計算結(jié)果精度高,在斜爆轟發(fā)動機燃燒室設計等工程應用中能夠有效節(jié)約計算時間、成本和復雜度.另外,兩者的對比研究,可以區(qū)分壁面效應、邊界層效應和黏性效應的影響,為進一步研究斜爆轟的波系結(jié)構(gòu)和復雜流動提供了一種新方法.