一維與二維爆轟傳播的時空關聯特性數值研究1)

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

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

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

引言

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

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

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

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

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

1 數值計算方法

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

對于一維爆轟的數值模擬

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

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

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

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

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

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

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

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

2 計算模型的驗證

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

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

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

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

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

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

3 計算結果分析

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

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

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

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

3.1 波系結構

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

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

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

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

3.2 壁面參數

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

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

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

3.3 波形變化

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

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

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

3.4時間尺度

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

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

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

表2 爆轟波各發展階段的時間尺度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外,二維與一維數值模擬結果在非定常尺度上存在顯著差別.根據理論分析,只有波后氣體在壁面切向的運動速度不變的情況下,二維楔面誘導斜爆轟波后的氣體運動方程才與一維活塞驅動正爆轟的波后氣體運動方程一致.然而在上述非定常過程中,均存在波系相互作用,氣體速度在壁面法向的梯度較大,選擇一個特定的速度將二維空間坐標x轉換為時間坐標t時存在較大誤差.因而通過活塞數值模擬無法準確預測楔面算例的非定常時間,這種差異本質上是流動的二維效應造成的.

4 結論

經過時空變換,有限長楔誘導的定常斜爆轟與非定常活塞驅動的正爆轟具有相同的波系結構,雖然兩者的非定常時間尺度存在一定差異,但爆轟參數、波形結構演化規律在定性和定量上均有較好的一致性,體現出二維斜爆轟與一維非定常正爆轟現象存在著時空關聯特性.

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