基于Additive Runge-Kutta方法的激波聚焦起爆高精度數(shù)值模擬

王成， 宋清官

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室， 北京 100081)

?

基于Additive Runge-Kutta方法的激波聚焦起爆高精度數(shù)值模擬

王成， 宋清官

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室， 北京 100081)

基于詳細(xì)氫氧化學(xué)動力學(xué)模型，建立了描述氫氧爆轟的多組分反應(yīng)歐拉方程組. 針對建立的反應(yīng)歐拉方程組，數(shù)值方法上采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分，采用5階精度的加權(quán)本質(zhì)無振蕩(WENO)格式對空間對流項(xiàng)進(jìn)行離散，自主研發(fā)了大規(guī)模高精度計(jì)算程序. 該程序能夠處理化學(xué)反應(yīng)源項(xiàng)引起的剛性問題，且能節(jié)省計(jì)算時(shí)間和計(jì)算內(nèi)存. 對半球型、半橢球型、圓錐型3種結(jié)構(gòu)形式凹面腔內(nèi)的激波聚焦起爆過程進(jìn)行了數(shù)值模擬，數(shù)值模擬研究得到了不同結(jié)構(gòu)形式凹面腔內(nèi)的激波聚焦起爆過程.

詳細(xì)化學(xué)動力學(xué)模型； Additive Runge-Kutta方法； WENO格式；激波聚焦

氣相爆轟是一個超音速燃燒過程，依靠波陣面后反應(yīng)區(qū)能量的快速釋放來自維持，因此氣相爆轟是燃料氧氣混合物釋放燃料化學(xué)能極其充分的方式，因而具有較高的能量利用效率[1]. 脈沖爆轟發(fā)動機(jī)(pulse detonation engine， PDE)利用了氣相爆轟的特點(diǎn)，在發(fā)動機(jī)性能提升方面擁有巨大優(yōu)勢，正逐漸成為了國內(nèi)外新型發(fā)動機(jī)研究熱點(diǎn). 爆轟波的點(diǎn)火起爆是脈沖爆轟發(fā)動機(jī)工作的關(guān)鍵條件，激波聚焦能夠在可燃?xì)怏w中產(chǎn)生高溫高壓區(qū)，進(jìn)而可以點(diǎn)火起爆形成爆轟波，是脈沖爆轟發(fā)動機(jī)點(diǎn)火起爆研究的一個重要方向.

國外在20世紀(jì)90年代已經(jīng)開始了激波聚焦方面的研究，且發(fā)展較為成熟. Izumi等[2]對惰性激波在拋物形或楔形凹面腔的反射問題進(jìn)行了實(shí)驗(yàn)和數(shù)值模擬研究；Levin等[3]提出了基于激波聚焦起爆的兩級脈沖爆轟發(fā)動機(jī)概念，解決了脈沖爆轟發(fā)動機(jī)的點(diǎn)火起爆問題，并使得提高脈沖爆轟發(fā)動機(jī)的工作頻率成為可能. 國內(nèi)激波聚焦問題研究起步較晚. 董剛等[4]數(shù)值模擬研究了氣體反應(yīng)活性、入射激波強(qiáng)度及反射壁面形狀對激波聚焦起爆的影響；曾昊等[5]開展了兩級脈沖爆轟發(fā)動機(jī)的相關(guān)研究，數(shù)值模擬了凹面腔內(nèi)的激波聚焦起爆過程，分析了環(huán)形射流噴口位置、環(huán)形射流初始壓力、入射噴口寬度、凹面腔結(jié)構(gòu)對激波聚焦起爆過程的影響，并進(jìn)行了初步的原理性實(shí)驗(yàn)研究，研究了導(dǎo)流角、氣流出口面積、擴(kuò)張角尾噴管、凹面腔與射流入口間距離等對激波聚焦的影響.

為進(jìn)一步研究凹面腔內(nèi)的激波聚焦起爆過程，本文基于詳細(xì)氫氧化學(xué)動力學(xué)模型[6]，數(shù)值方法上采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分，空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式[7]，分別對半球型、半橢球型和圓錐型3種結(jié)構(gòu)的凹面腔的激波聚焦起爆過程進(jìn)行了數(shù)值模擬，對比研究了不同凹面腔結(jié)構(gòu)的激波聚焦規(guī)律.

1 控制方程

本文采用多組分反應(yīng)Euler方程組描述氫氧預(yù)混氣體的爆轟化學(xué)反應(yīng)過程

(1)

(2)

(3)

(4)

式中：p，ρ，e，T，u，v，h分別為壓力、密度、單位質(zhì)量的總能、溫度、x和y方向的速度以及單位質(zhì)量的焓；Y1，Y2，…，YN-1分別為基元反應(yīng)各組分的質(zhì)量分?jǐn)?shù)；R為普適氣體常數(shù)；ω1，ω2，…，ωN-1分別為各組分的質(zhì)量生成率.

為準(zhǔn)確地描述氫氧預(yù)混氣體的爆轟反應(yīng)過程，化學(xué)反應(yīng)模型采用9組分20步基元反應(yīng)模型[1]. 參與化學(xué)反應(yīng)的各組分分別為 H2，O2，OH，O，H，H2O2，HO2，H2O，此外氫氧預(yù)混氣體用惰性氣體N2稀釋，N2不參與化學(xué)反應(yīng). 化學(xué)反應(yīng)方程可統(tǒng)一寫為

(5)

kfm=AmTnmexp-EamR0T，

kpm=exp∑9i=1υ″mi-υ′mis0iRi-hiRiT，

kem=kpmexpp0R0T∑9i=1(υ″mi-υ′mi)，

kbm=kfm/(kem).(6)

(7)

式中：Wi為第i組分的摩爾質(zhì)量；γTB為第三體效應(yīng)指數(shù).

2 數(shù)值方法

在數(shù)值計(jì)算方法方面，基于WENO-LF[8]對物理通量進(jìn)行分裂，方程(1)離散為

(8)

由于化學(xué)反應(yīng)模型采用了9組分20步基元反應(yīng)模型，對流項(xiàng)與化學(xué)反應(yīng)源項(xiàng)之間，以及各個化學(xué)反應(yīng)源項(xiàng)之間時(shí)間尺度差別很大，導(dǎo)致歐拉方程組出現(xiàn)嚴(yán)重的剛性問題. 為了解決剛性問題，本文采用了顯隱式Additive Runge-Kutta方法對時(shí)間方向進(jìn)行積分[10]. 該方法能夠解決下述形式的常微分方程組：

(9)

式中：F(U)為N項(xiàng)的總和，對于多組分反應(yīng)歐拉方程組，N=2，即方程右端F(U)包含對流項(xiàng)Fns(U)和化學(xué)反應(yīng)源項(xiàng)Fs(U).

AdditiveRunge-Kutta方法的主要思想是將剛性項(xiàng)Fns(U)和非剛性項(xiàng)Fs(U)分開處理，對非剛性項(xiàng)Fs(U)運(yùn)用顯式Runge-Kutta(ERK)方法進(jìn)行顯式積分，對剛性項(xiàng)Fns(U)運(yùn)用對角隱式Runge-Kutta(ESDIRK)方法進(jìn)行隱式積分. s步AdditiveRunge-Kutta法的具體使用方法如下：

(10)

3 計(jì)算模型及結(jié)果分析

3.1 一維氣相爆轟

預(yù)混氣體H2/O2/N2的體積比為0.50/0.25/0.25，初始溫度T0=300 K，初始壓力pa=0.101 MPa. 在左側(cè)封閉端設(shè)置一段高溫高壓點(diǎn)火區(qū)，點(diǎn)火溫度Tin=3 000 K，點(diǎn)火壓力pa=0.5 MPa. 數(shù)值方法上對于空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式.

為研究網(wǎng)格尺寸對計(jì)算結(jié)果的影響，分別取網(wǎng)格尺寸為2.00，1.00，0.50，0.25，0.10，0.05 mm.

表1 不同網(wǎng)格尺寸下的爆轟參數(shù)

從表1中可以看出，平均爆速、CJ爆壓和壁面壓力在不同網(wǎng)格尺寸下基本趨于相同值，即網(wǎng)格尺寸對其影響不大，且各物理量值與理論值基本吻合；而Von Neumann峰壓力隨著網(wǎng)格尺寸減小而增大，且隨著網(wǎng)格尺寸的減少而收斂. 綜上所述，網(wǎng)格尺寸取為0.1 mm時(shí)滿足計(jì)算精度，既能準(zhǔn)確的計(jì)算各爆轟參數(shù)，又能減少網(wǎng)格數(shù)，節(jié)省計(jì)算時(shí)間.

選取網(wǎng)格尺寸為0.1 mm，通過對比各組分的體積分?jǐn)?shù)，容易發(fā)現(xiàn)各組分變化的量級差別很大，H2，O2，H2O變化的量級為10-1，OH，O，H變化的量級為10-2，H2O2和HO2變化的量級分別為10-4和10-6且衰減很快，如圖1(g)～1(h)，因此不同組分在氫氧預(yù)混氣體爆轟波的化學(xué)反應(yīng)中起的作用并不相同，甚至差別很大，可以忽略對反應(yīng)影響較小的組分H2O2和HO2，對基元反應(yīng)模型做進(jìn)一步簡化，簡化的化學(xué)反應(yīng)模型應(yīng)用到下面的激波聚焦起爆研究中.

3.2 激波聚焦起爆

為了研究不同凹面腔結(jié)構(gòu)內(nèi)的激波聚焦起爆過程，本文選取了半球型、半橢球型和圓錐型3種結(jié)構(gòu)的凹面腔作為計(jì)算模型，如下：

工況1：直徑為50 mm，開口直徑為25 mm的半球凹面腔；

工況2：長徑為50 mm，短徑為25 mm的半橢球凹面腔；

工況3：高為25 mm，底面圓直徑為25 mm的圓錐凹面腔.

凹面腔和尾噴管為內(nèi)區(qū)域(實(shí)線)，噴管外為外區(qū)域(虛線)，如圖2所示. 對于內(nèi)區(qū)域，填充體積比0.23/0.16/0.61的H2/O2/N2預(yù)混氣體，等截面尾噴管軸線方向長12 mm，環(huán)形射流入口為壓力入口邊界，環(huán)形射流入射寬度2.8 mm，壓力pin=0.45 MPa，溫度Tin=450 K；對于外區(qū)域，填充空氣，壓力pa=0.101 MPa，溫度Ta=300 K. 凹面腔及尾噴管為剛性、無滑移、絕熱壁面. 數(shù)值方法上采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分，空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式. 從初始網(wǎng)格尺寸為0.1 mm. 選取凹面腔中軸線為壓力觀測線，選取凹面腔中軸線上的點(diǎn)M作為壓力觀測點(diǎn)，如圖2所示.

從圖3中可以看出，高溫高壓的環(huán)形射流由入口噴入半球形凹面腔內(nèi)，初始的溫度、壓力不足以起爆預(yù)混氣體，在凹面腔內(nèi)形成了兩道激波，如圖3(a)所示. 隨著激波的傳播，激波的強(qiáng)度逐漸減弱，如圖3(b)所示，但是當(dāng)上下激波相遇碰撞時(shí)，壓力增加，激波的強(qiáng)度加強(qiáng)，此時(shí)仍未能起爆預(yù)混氣體，如圖3(c) 所示.

激波繼續(xù)向前傳播，在半球型壁面處發(fā)生反射，由于壁面收斂的作用，壓力逐漸增加，形成了更強(qiáng)的激波. 激波繼續(xù)向前傳播，在半球凹面腔中心處發(fā)生匯聚，如圖3(d)所示，當(dāng)激波在球形凹面腔中心處發(fā)生匯聚時(shí)，壓力密度等參量急劇增加，達(dá)到了氫氧預(yù)混氣體的臨界起爆條件，發(fā)生起爆，如圖3(e)所示. 起爆之后，爆轟波在球形凹面腔內(nèi)迅速傳播，如圖3(f)～3(h)所示，實(shí)現(xiàn)了脈沖爆轟發(fā)動機(jī)的起爆過程.

圖3～圖6可以看出，對于半球凹面腔，在t=310.72 μs時(shí)，在中軸線上發(fā)生激波聚焦起爆，起爆壓力約為3.0 MPa；對于半橢球凹面腔，在t=330.05 μs時(shí)，在中軸線上發(fā)生激波聚焦起爆，起爆壓力約為2.5 MPa；對于圓錐凹面腔，在t=307.77 μs時(shí)，在中軸線附近上發(fā)生激波聚焦起爆，起爆壓力約為2.0 MPa，同時(shí)發(fā)現(xiàn)初始起爆點(diǎn)明顯位于中軸線偏上的位置，本文中給出的入流條件和物理模型均是關(guān)于中軸線對稱的，由此推測這種現(xiàn)象是圓錐凹面腔頂點(diǎn)的奇異性導(dǎo)致的數(shù)值誤差造成的.

通過對比3種結(jié)構(gòu)的凹面腔的激波聚焦過程，發(fā)現(xiàn)在此初始條件下，3種結(jié)構(gòu)的凹面腔均實(shí)現(xiàn)了激波聚焦起爆過程，且聚焦中心基本均位于中軸線上，壁面收斂導(dǎo)致的激波增強(qiáng)效果不足以起爆預(yù)混氣體. 由于半橢球凹面腔的曲率要比半球凹面腔和圓錐凹面腔的曲率要大，因此半橢球凹面腔的激波聚焦效果最為明顯，聚焦中心位置最為靠前，如圖3～7所示. 此外，通過對比發(fā)現(xiàn)，3種結(jié)構(gòu)的凹面腔實(shí)現(xiàn)聚焦起爆時(shí)間不同，其中圓錐凹面腔最早實(shí)現(xiàn)聚焦起爆，半球凹面腔次之，半橢球凹面腔最晚實(shí)現(xiàn)聚焦起爆.

此外本文還研究了不同凹面腔結(jié)構(gòu)對脈沖爆轟發(fā)動機(jī)工作頻率的影響. 由分析可知，凹面腔內(nèi)預(yù)混氣體起爆之后產(chǎn)生大量高溫高壓氣體，氣體由尾噴管噴出形成推力，氣體噴出后在凹面腔內(nèi)形成了低溫低壓區(qū)，此時(shí)入口處的氣體總壓要高于凹面腔內(nèi)的氣體總壓，此時(shí)高溫高壓的射流再次由入口進(jìn)入凹面腔內(nèi)，再次實(shí)現(xiàn)聚焦爆轟過程，進(jìn)入下一個循環(huán)，因此凹面腔內(nèi)形成了周期性的脈沖爆轟波，實(shí)現(xiàn)了循環(huán)起爆. 從圖中容易看出，對于3種結(jié)構(gòu)的凹面腔，觀測點(diǎn)M的壓力均在400 μs后發(fā)生了周期性的變化，因此對于3種結(jié)構(gòu)的凹面腔，起爆之后均實(shí)現(xiàn)了循環(huán)起爆. 其中，半球凹面腔的工作頻率約為16.7 kHz，半橢球凹面腔的工作頻率約為15.0 kHz，圓錐凹面腔的工作頻率約為20.0 kHz，即圓錐凹面腔的工作頻率最高，半球凹面腔的工作頻率次之，半橢球凹面腔的工作頻率最低.

4 結(jié) 論

采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分，空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式，對脈沖爆轟發(fā)動機(jī)激波聚焦起爆過程進(jìn)行了數(shù)值模擬，重點(diǎn)研究了半球凹面腔、半橢球凹面腔和圓錐凹面腔對激波聚焦起爆過程的影響，并得到了以下結(jié)論:

① 3種結(jié)構(gòu)的凹面腔均實(shí)現(xiàn)了激波聚焦起爆過程，且聚焦中心基本均位于中軸線上，其中半橢球凹面腔的激波聚焦效果最為明顯，聚焦中心位置最為靠前，因此為獲得明顯的聚焦效果，應(yīng)使用半橢球凹面腔；

② 3種結(jié)構(gòu)的凹面腔實(shí)現(xiàn)聚焦起爆時(shí)間不同，其中圓錐凹面腔最早實(shí)現(xiàn)聚焦起爆，半球凹面腔次之，半橢球凹面腔最晚實(shí)現(xiàn)聚焦起爆，因此為較快實(shí)現(xiàn)聚焦起爆，應(yīng)使用圓錐凹面腔；

③ 3種結(jié)構(gòu)的凹面腔均實(shí)現(xiàn)了循環(huán)起爆，其中，圓錐凹面腔的工作頻率最高，半球凹面腔的工作頻率次之，半橢球的工作頻率最低，因此為達(dá)到較高的工作頻率，應(yīng)使用圓錐凹面腔.

[1] 張博，白春華.氣相爆轟動力學(xué)特征研究進(jìn)展[J].中國科學(xué)：物理學(xué) 力學(xué) 天文學(xué)，2014，44(7):665-681.

Zhang Bo， Bai Chunhua. Research progress on the dynamic characteristics of gaseous detonation[J]. Scientia Sinica: Physica， Mechanica & Astronomica， 2014，44(7):665-681. (in Chinese)

[2] Izumi K， Aso S， Nishida M. Experimental and computational studies focusing processes of shock waves reflected from parabolic reflectors[J]. Shock Waves， 1994，3(3):213-222．

[3] Levin V A， Nechaev J N， Tarasov A I. A new approach to organizing operation cycles in pulse detonation engines[C]∥Proceedings of High-Speed Deflagration and Detonation: Fundamentals and Control. Moscov: [s.n.], 2001:31-44.

[4] 董剛，唐敖，葉經(jīng)方，等.激波聚焦誘導(dǎo)點(diǎn)火和爆轟的數(shù)值研究[J].爆炸與沖擊，2006，25(5):437-444.

Dong Gang， Tang Ao， Ye Jingfang， et al. Numerical studies on initiation and detonation induced by shock wave focusing[J]. Explosion and Shock Waves， 2006，25(5):437-444. (in Chinese)

[5] 曾昊，何立明，章雄偉，等.環(huán)形射流噴口位置對激波聚焦起爆的影響分析[J].推進(jìn)技術(shù)，2011，32(3):437-442.

Zeng Hao， He Liming， Zhang Xiongwei， et al. Investigation on the influence of jets spout location on detonation initiation via imploding annular shock waves[J]. Journal of Propulsion Technology， 2001，32(3):437-442. (in Chinese)

[6] 趙慧，李健，郝莉.不同曲率彎管對氣相爆轟波傳播特性的影響[J].北京理工大學(xué)學(xué)報(bào)，2014，34(增刊1):176-180．

Zhao Hui， Li Jian， Hao Li. The influence of different bend curvature on gaseous detonation wave propagation[J]. Transactions of Beijing Institute of Technology， 2014，34(suppl 1):176-180. (in Chinese)

[7] 藺偉，宋清官，王成，等.濃度梯度對瓦斯爆炸影響的數(shù)值模擬[J].北京理工大學(xué)學(xué)報(bào)，2015，35(4):336-340.

Lin Wei， Song Qingguan， Wang Cheng， et al. Numerical investigation about the effect of concentration gradient on methane explosion[J]. Transactions of Beijing Institute of Technology， 2015，35(4):336-340. (in Chinese)

[8] Jiang G S， Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics， 1996，126(1):202-228．

[9] Chi-Shu W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[J]. Lecture Notes in Mathematics， 1998，1697(1):325-432.

[10] Kennedy C A， Carpenter M H. Additive Runge-Kutta schemes for convection diffusion reaction equations[J]. Applied Numerical Mathematics， 2003，44(1-2):139-181.

(責(zé)任編輯：劉雨)

Numerical Application of Additive Runge-Kutta Methods on Detonation Initiation with Convergence of Shock Waves

WANG Cheng， SONG Qing-guan

(State Key Laboratory of Explosion Science and Technology， Beijing Institute of Technology， Beijing 100081， China)

Based on detailed hydrogen and oxygen chemical kinetics model, the multi-species reactive Euler equations were established to describe the hydrogen and oxygen detonation. Using third-order Additive Runge-Kutta methods in time discretization, using fifth-order weighted essentially non-oscillatory (WENO) scheme in spatial discretization, a high-solution large-scale program was developed. It was verified that the program could solve stiff source term caused by multispecies chemical reaction, which could also save computation time and computation memory. With this program, the process of detonation initiation with convergence of shock waves was numerically simulated for 3 types of concave cavity, hemispherical concave cavity, semi-ellipsoidal concave cavity and conical concave cavity. The simulation results show the characteristic of detonation initiation in different types of concave cavity.

detailed chemical kinetics model; Additive Runge-Kutta methods; weighted essentially non-oscillatory scheme; convergence of shock waves

2015-10-30

國家自然科學(xué)基金資助項(xiàng)目(11325209，11272056)

王成(1972—)，男，教授，博士生導(dǎo)師，E-mail：wangcheng@bit.edu.cn.

O 383

A

1001-0645(2016)02-0137-07

10.15918/j.tbit1001-0645.2016.02.006