方宜申 胡宗民 滕宏輝姜宗林

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

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

圓球誘發(fā)斜爆轟波的數值研究1)

方宜申 胡宗民 滕宏輝2)姜宗林

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

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

斜爆轟發(fā)動機是飛行器在高馬赫數飛行條件下的一種新型發(fā)動機，具有結構簡單、成本低和比沖高等優(yōu)點.但是斜爆轟發(fā)動機的來流馬赫數范圍廣，來流條件復雜，為實現斜爆轟波的迅速、可靠引發(fā)，采用鈍頭體來誘發(fā).利用Euler方程和氫氧基元反應模型，對超聲速氫氣/空氣混合氣體中圓球誘導的斜爆轟流場進行了數值研究.不同于楔面誘發(fā)的斜爆轟波，球體首先會在駐點附近誘發(fā)正激波/爆轟波，然后在稀疏波作用下發(fā)展為斜激波/爆轟波.模擬結果顯示，經過鈍頭體壓縮的預混氣體達到自燃溫度后，會出現兩種流場：當馬赫數較低時，由于稀疏波的影響，燃燒熄滅，鈍頭體下游不會出現燃燒情況；而當馬赫數較高時，燃燒陣面能傳到下游.分析表明，當鈍頭體的尺度較小時，駐點附近的能量不足以誘發(fā)爆轟波，只會形成明顯的燃燒帶與激波非耦合結構；當鈍頭體的尺度較大時，流場中不會出現燃燒帶與激波的非耦合現象，且這一特征與馬赫數無關.通過調整球體直徑，獲得了激波和燃燒帶部分耦合的燃燒流場結構，這一流場結構在楔面誘發(fā)的斜爆轟波中并不存在，說明稀疏波與爆轟波面的相互作用是決定圓球誘發(fā)斜爆轟波的關鍵.

斜爆轟，氫氧，起爆，熄爆

引言

航空航天技術的飛速發(fā)展，要求飛行器的飛行速度越來越快，傳統(tǒng)渦輪發(fā)動機難以滿足要求.爆轟發(fā)動機是現代發(fā)動機的一種選擇，主要有斜爆轟發(fā)動機、脈沖爆轟發(fā)動機[1]和旋轉爆轟發(fā)動機[2]三種.斜爆轟發(fā)動機具有結構簡單、質量小、成本低、比沖高和適應高速飛行的優(yōu)點，在飛行馬赫數高達10的情況下仍能為飛行器提供有效推力[3].

爆轟波是一種強激波誘發(fā)快速燃燒的耦合流場結構.爆轟波相對于預混氣體是超聲速傳播的.在實驗室坐標系下，如果預混氣以超聲速流動，理論上可能實現爆轟波的駐定.但實際上，駐定正爆轟波很難實現，而利用斜激波誘發(fā)的斜爆轟波已經成功實現了駐定[4]，并且相對正爆轟波具有較小的總壓損失.而關于斜爆轟燃燒技術在高超聲速飛行器動力裝置中的應用，近年來已經開展了多方面相關的基礎理論與工程應用的數值和實驗研究[5-9].

斜爆轟發(fā)動機不僅具有超燃發(fā)動機的優(yōu)點，且燃燒過程近似等容，相比于超燃發(fā)動機的等壓燃燒具有更高的燃燒效率.Li等[10]通過數值分析發(fā)現，楔面誘發(fā)的斜爆轟流場結構包含斜激波、反應區(qū)、爆燃波以及之后的斜爆轟波.Viguier等[11]在1997年用實驗證實了這一斜爆轟波結構的存在.Choi等[12]發(fā)現，斜爆轟波也有如正爆轟波一般的胞格結構，不過與正爆轟波不同的是：斜爆轟波三波點的軌跡不存在相互交叉.Teng等[13-15]對斜爆轟波胞格結構的變化規(guī)律給出了定量化的研究結論.斜爆轟發(fā)動機的來流馬赫數范圍廣，來流條件復雜.Zhang等[16-17]研究發(fā)現來流當量比不同會影響斜爆轟波的特征長度.為了實現斜爆轟波的迅速、可靠引發(fā)，可以采用鈍頭體進行誘發(fā).但是鈍頭體誘發(fā)斜爆轟波的過程比楔面誘發(fā)更加復雜，涉及到正激波、斜激波和稀疏波與燃燒反應的復雜相互作用，其中的波系結構和起爆機理還不清楚.

1972年，Lehr[18]利用彈狀鈍頭體誘發(fā)爆轟波，試驗結果表明鈍頭體誘發(fā)的爆轟波與楔面誘發(fā)的爆轟波在流場結構上有極大的不同.Kaneshige等[19]在1996年的氫氣/氧氣爆轟實驗中發(fā)現，利用球形鈍頭體誘發(fā)爆轟波的過程出現了伴隨著燃燒帶的激波、爆燃轉爆轟、直接起爆等三種流場結構.Ju等[20]通過理論和數值研究對起爆給出了定量化規(guī)律.Maeda等[21-25]在乙炔/氧氣爆轟實驗中，得到了燃燒帶與激波非耦合、草帽狀斜爆轟波、直接起爆爆轟波等三種典型的球形鈍頭體誘發(fā)斜爆轟波的流場結構.在非楔面斜爆轟中，由于稀疏波的存在，對斜爆轟波結構帶來一定的影響[26].

本文在溫度300K,壓力50kPa的初始條件下，改變來流馬赫數和球體直徑，研究球形鈍頭體誘發(fā)爆轟波的過程以及不同流場結構之間的臨界條件.通過數值模擬給出流場結構隨馬赫數與球體直徑變化的分布圖.

1 數學與物理模型

本研究對圓球誘發(fā)斜爆轟波過程進行數值模擬.預混可燃氣體在超聲速流動下，在鈍頭體頂部形成脫體激波，駐點附近為高溫高壓區(qū).不同的來流馬赫數及鈍頭體大小決定了駐點附近的流場參數，駐點附近的流場參數又決定了下游的流場結構.當駐點附近得到的能量高過誘發(fā)爆轟波的臨界能量時，流場出現爆轟波；反之，則不能.

現有的數值研究表明，在超聲速流動下，黏性對流場影響很小，數值研究主要采用無黏流假設.因此控制方程簡化為二維軸對稱Euler方程

其中，ωi為化學反應中第i種組分的質量生成率，它由化學反應模型決定.在上述方程中ρi是第i種組分的密度，總密度和v分別代表x方向和r方向上的速度.比內能e的計算公式為

其中，Ri代表第i種組分的氣體常數，T為溫度.

數值模擬中采用自適應無結構的四邊形網格[27]和MUSCL-Hancock算法[28]，Rienmann問題的求解采用HLLC算子.時間項為由CFL數控制的顯式格式，CFL數為0.4.化學反應采用時間分裂算法.氫/空氣化學反應模型[29]含有11種組元(H2,O2,O,H,OH, HO2,H2O2,H2O,N2,N,NO)和23個基元化學反應，詳細的機理見Chemkin程序說明，化學反應帶來的剛性問題由DVODE軟件包解決[30].氫氣/空氣混合物的化學當量比φ=1.0，即H2:O2:N2=2:1:4.球體邊界為固壁邊界，下邊界為軸對稱邊界，左邊界為來流，其他邊界為自由邊界.

2 數值結果與討論

參考Lehr[18]在1972年的實驗，給定流場的初始溫度為300K，壓力為50kPa.實驗中的鈍頭體直徑為15mm，為研究鈍頭體尺寸對于爆轟的影響，選取球形直徑D=5～15mm.由于球形直徑的不同，初始網格尺寸也不同，最大網格不大于1.2mm，加密后最小網格不大于0.1mm，此時繼續(xù)加密網格流場結構不發(fā)生變化.圖1所示為Ma=4.0,D=5mm，加密層數為3和4，即最小網格為0.125mm和0.06mm時的網格加密情況.可以看到加密處基本不變.圖2為這兩種加密情況下，球頭前的壓力溫度分布曲線，圖中的數字表示加密層數，兩種加密網格的壓力和溫度曲線基本重合.

圖1 Ma=4.0，D=5mm，加密層數3(左)，4(右)網格分布Fig.1Ma=4.0，D=5mm,encryption layer 3(left)and 4(right)grid distribution

圖2 Ma=4.0，D=5mm，y=0時的壓力和溫度變化Fig.2Ma=4.0，D=5mm,the pressure and temperature at liney=0 variation diagram

2.1 算例與討論

算例1 球體直徑D=5mm

當球體直徑為5mm時，馬赫數的增大不會引發(fā)爆轟波.當馬赫數為4.0的時候，駐點附近溫度達到自燃溫度，氣體燃燒并向下游傳播.由于鈍頭體為球形，氣體先受到球體的壓縮作用，后受到膨脹作用.由于稀疏波的影響，燃燒熄滅，鈍頭體下游不會出現燃燒帶，如圖3.

圖3 Ma=4.0,D=5mm，壓力(上)、溫度(下)分布圖Fig.3Ma=4.0,D=5mm,pressure(upper)and temperature(lower) distribution

當馬赫數為5.0的時候，駐點附近的燃燒延續(xù)到下游流場，形成燃燒帶與激波非耦合的現象，如圖4.

圖4 Ma=5.0,D=5mm，壓力(上)、溫度(下)分布圖Fig.4Ma=5.0,D=5mm,pressure(upper)and temperature(lower) distribution

繼續(xù)增大來流速度直至馬赫數為15，流場結構都為燃燒帶與激波非耦合的情況，馬赫數較大時只會導致燃燒帶變寬，流場結構不會出現變化.

為界定燃燒熄滅與非耦合情況的臨界馬赫數，分析馬赫數為4.12和4.14時的流場結構.結果表明當馬赫數大于4.1時，出現燃燒帶與激波非耦合的情況，且最大溫度也會有較大的變化，如圖5所示.

可以看到，當馬赫數從4.1變化到4.12的時候，不僅流場結構出現變化，駐點溫度大幅升高，因此當精度為0.02時，臨界馬赫數等于4.1.

圖5 最高溫度隨馬赫數的變化Fig.5 Maximum temperature with di ff erent Mach numbers

算例2 球體直徑D=15mm

當球體直徑為 15mm時，流場結構與直徑為5mm時的不完全一樣.馬赫數為4.0的時候，流場結構與直徑為5mm時相似，但是在馬赫數為5.0的時候，直徑15mm鈍頭體可以直接引發(fā)爆轟波，如圖6所示.

圖6 D=15mm，壓力(上)、溫度(下)分布圖Fig.6D=15mm,pressure(upper)and temperature(lower) distribution

為研究在馬赫數4.0～5.0之間的流場情況，馬赫數變化步長設為0.1.結果發(fā)現，在馬赫數為4.1的時候，流場并不出現燃燒帶與激波非耦合的情況，而是在駐點位置引發(fā)爆轟波，并誘發(fā)下游流場形成爆轟波.由于C-J爆速大于此時的來流馬赫數，爆轟波無法駐定，向上游傳播.

對比直徑為5mm和15mm的結果可以看到，在爆轟燃燒中，在同樣的流場條件下，鈍頭體的尺度是誘發(fā)爆轟波的一個重要參數.

2.2 流場結構隨來流馬赫數和球體直徑變化的分布規(guī)律

來流馬赫數和球體鈍頭體直徑都是影響流場結構的關鍵參數.為了深入研究其對流場的影響，設置直徑變化步長為2.5mm，馬赫數變化步長為0.1.流場結構如圖7所示.

從圖7可以看到，馬赫數的變化對于流場是否燃燒的影響較大.當馬赫數從4.0變化到4.2之后，流場燃燒熄滅的情況都消失了.而球體直徑的大小對于誘發(fā)爆轟波的影響較大.當球體直徑為5mm時，鈍頭體無法誘發(fā)爆轟波；當球體直徑增大到7.5mm時，鈍頭體才能誘發(fā)爆轟波.

圖7 馬赫數和圓球直徑對起爆的影響Fig.7 The influenc of Mach number and ball diameter on initiation

對馬赫數為4.2，直徑為5～7.5mm進行進一步的數值研究發(fā)現，如果直徑變化步長為0.5mm,當直徑小于6.5mm時，流場結構為燃燒帶與激波非耦合；當直徑大于6.5mm時，預混氣體在經過一段時間后起爆；當直徑為6.5mm時，駐點位置的燃燒改變弓形激波在駐點附近的形狀，但是不足以克服稀疏波的影響，對激波的影響不能傳播到整個流場，這一流動是定常的.可以看到激波面由兩道斜激波構成，在y=20mm處存在一個明顯的波面拐點，如圖8所示.拐點上游是斜爆轟波，下游是斜激波，燃燒面與激波面發(fā)生了明顯的解耦.這一現象在楔面誘發(fā)的斜爆轟流場[13-16]中是不存在的.與楔面誘發(fā)斜爆轟波不同的是，氣體在流經圓球時，不僅受到激波壓縮作用，亦受到稀疏波的作用，圓球誘發(fā)的稀疏波是導致解耦的原因.因此可以得出結論，稀疏波與爆轟波面的相互作用是決定圓球誘發(fā)斜爆轟波的關鍵.

圖8 Ma=4.2，D=6.5mm，壓力(上)、溫度(下)分布圖Fig.8Ma=4.2,D=6.5mm,pressure(upper)and temperature(lower) distribution

3 結論

本文采用數值模擬對圓球誘發(fā)斜爆轟波進行了研究，重點分析了不同球體直徑對起爆的影響.當球體直徑為5mm時，流場結構不會出現爆轟現象.隨著馬赫數的增加，流場從燃燒熄滅變到燃燒帶與激波非耦合的結構，臨界馬赫數為4.1；當球體直徑為15mm時，流場中能夠觀測到爆轟現象.隨著馬赫數的增加，流場從燃燒熄滅變到誘發(fā)爆轟波，但是不會出現燃燒帶與激波非耦合的情況；當球體直徑在5～15mm之間，馬赫數高于4.1的時候，燃燒熄滅的情況不再出現，球體直徑對于爆轟波引發(fā)的影響明顯，且當直徑大于6.5mm的時候開始出現爆轟現象.

以前的學者對楔面誘發(fā)斜爆轟波進行了較多的研究，但是對球體誘發(fā)斜爆轟波還缺乏深入的研究.這種斜爆轟波的引發(fā)過程更加復雜，涉及到正激波/爆轟波、稀疏波與斜激波/爆轟波的相互作用.由于不同的直徑會誘發(fā)不同強度的稀疏波，它們與燃燒的相互作用成為能否成功起爆的關鍵.本文的研究證實大直徑圓球更容易誘發(fā)爆轟波，和以前的起爆實驗與理論[17]結果也是定性一致的.下一步需要在此基礎上，對臨界狀態(tài)下膨脹波與激波誘發(fā)燃燒的相互作用進行深入、系統(tǒng)的研究，從而建立適應工程應用的起爆模型.

1 董剛，范寶春，謝波.氫氣--空氣混合物中瞬態(tài)爆轟過程的二維數組模擬.高壓物理學報，2004，18(1)：40-46(Dong Gang,Fan Baochun,Xie bo.Two-dimensional simulation of transient detonation process for H2-O2-N2mixture.Chinese Journal of High Pressure Physics,2004，18(1)：40-46(in Chinese))

2 LinW,ZhouJ,LiuS,etal.Experimentalstudyonpropagationmode of H2/air continuously rotating detonation wave.International Journal of Hydrogen Energy,2015,40(4):1980-1993

3 Voland RT,Huebner LD,McClinton CR.X-43A hypersonic vehicle technology development.Acta Astronautica，2006,59(1-5)：181-191

4 林志勇，周進，張繼業(yè)等.預混超聲速氣流斜激波誘發(fā)脫體爆轟研究.航空動力學報，2009，24(1)：50-54(Lin Zhiyong,Zhou Jin, Zhang Jiye,et al.Investigation of detached detonation induced by oblique shock in premixed supersonic fl w.Journal of Aerospace Power,2009，24(1)：50-54(in Chinese))

5 Rudy W,Dziubanii K,Zbikowski M,et al.Experimental determination of critical conditions for hydrogen-air detonation propagation in partially confine geometry.International Journal of Hydrogen Energy,2016,1-8

6 Wang C,Dong XZ,Shu CW.Parallel adaptive mesh refinemen method based on WENO finit di ff erence scheme for the simulation ofmulti-dimensional detonation.Journalof Computational Physics, 2015,298:161-175

7 Zhang Z,Li Z,Dong G.Numerical studies of multi-cycle acetyleneair detonation induced by shock focusing.Procedia Engineering, 2015,99:327-331

8 王成，寧建國，雷娟.障礙物對氫氧預混氣體爆轟波傳播的影響//慶祝中國力學學會成立50周年暨中國力學學會學術大會，北京，2007(Wang Cheng,Ning Jianguo,Lei Juan.The influenc of obstaclesonthepropagationofH2-airdetonation//Celebrationofthe 50 Anniversary of the Founding of the Chinese Society of Mechanics and the Academic Conference of Chinese Society of Mechanics, Beijing,2007(in Chinese))

9 Wang C,Dong XZ,Shu CW.Parallel adaptive mesh refinemen method based on WENO finit di ff erence scheme for the simulation ofmulti-dimensional detonation.Journalof Computational Physics, 2015,298:161-175

10 Li C,Kailasanath K,Oran ES.Detonation structures behind oblique shocks.Phys Fluids,1994,6:1600-1611

11 Viguier C,Silva LFFD,Desbordes D,et al.Onset of oblique detonation waves:comparison between experimental and numerical results for hydrogen–air mixtures.Proc Combust Inst,1997,26:3023-3031

12 Choi JY,Kim DW,Jeung IS,et al.Cell-like structure of unstable oblique detonation wave from high-resolution numerical simulation.Proc Combust Inst,2007,31:2473-2480

13 Teng HH,Jiang ZL,Ng HD.Numerical study on unstable surfaces ofobliquedetonations.JournalofFluidMechanics,2014,744:111-128

14 Teng HH,Ng HD,Li K,et al.Evolution of cellular structures on oblique detonation surfaces.Combustion and Flame,2015,162(2): 470-477

15 Teng HH,Jiang ZL.On the transition pattern of the oblique detonation structure.Journal of Fluid Mechanics,2012,713:659-669

16 Zhang Y,Gong J,Wang T.Numerical study on initiation of oblique detonations in hydrogen–air mixtures with various equivalence ratios.Aerospace Science and Technology,2016,49:130-134

17 Wang T,Zhang YN,Teng HH,et al.Numerical study of oblique detonation wave initiation in a stoichiometric hydrogen-air mixture,Physics of Fluids,2015,27(9):096101

18 Lehr HF.Experiments on shock-Induced combustion.Astronautica Acta,1972,17:589-597

19 Kaneshige MJ,Shepherd JE.Oblique detonation stabilized on a hypervelocity projectile//26th Symp.(Int.)on Combustion,Pittsburgh, 1996.3015

20 Ju Y,Masuya G,Sasoh A.Numerical and theoretical studies on detonation initiation by a supersonic projectile.Symposium on Combustion,1998,27(2):2225-2231

21 Maeda,S,Inada,R,Kasahara,J,et al.Visualization of the nonsteady state oblique detonation wave phenomena around hypersonic spherical projectile.Proc Combust Inst,2011,33:2343-2349

22 Maeda S,Kasahara J,Matsuo A.Oblique detonation wave stability around a spherical projectile by a high time resolution optical observation.Combustion and Flame,2012,159(2):887-896

23 Maeda S,Sumiya S,Kasahara J,et al.Initiation and sustaining mechanisms of stabilized oblique detonation waves around projectiles.Proceedings of the Combustion Institute,2013,34(2):1973-1980

24 Maeda S,Sumiya S,Kasahara J,et al.Scale e ff ect of spherical projectiles for stabilization of oblique detonation waves.Shock Waves, 2015,25(2):141-150

25 Choi JY,Maeda S,Kasahara J,et al.Calculation of drag coefficients for hypersonic spherical projectiles initiating oblique detonation wave or shock-induced combustion//50th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.2012

26 Liu Y,Han X,Yao S,et al.A numerical investigation of the prompt oblique detonation wave sustained by a finite-lengt wedge.Shock Waves,2016:1-11

27 SunM,TakayamaK.Conservativesmoothingonanadaptivequadrilateral grid.Journal of Computational Physics,1999,150:143-180

28 Toro EF.Riemann Solvers and Numerical Methods for Fluid Dynamics.(Second ed).Berlin:Springer,1999

29 Kee RJ,Rupley FM,Meeks E,et al.Chemkin-III:a fortran chemical kinetics package for the analysis of gas-phase chemical and plasma kinetics.UC-405,SAND96-8216,Sandia National Laboratories, 1996

30 BrownPN,ByrneGD,HindmarshAC.VODE:avariable-coefficient ODE solver.SIAM J Sci Stat Comput,1989,10(5):1038-1051

NUMERICAL STUDY OF THE OBLIQUE DETONATION INITIATION INDUCED BY SPHERES1)

Fang Yishen Hu Zongmin Teng Honghui2)Jiang Zonglin

(State Key Laboratory of High-Temperature Gasdynamics,Institute of Mechanics，Chinese Academy of Sciences，Beijing100190，China)

(School of Engineering Science，University of Chinese Academy of Sciences，Beijing100049，China)

The oblique detonation wave engine is a new kind of engine which has a simple structure,low cost,and high specifi impulse.In order to ensure the initiation,blunt body is used to induce the oblique detonation wave.The oblique detonation wave fl w fiel induced by spheres in supersonic hydrogen/air mixture is numerically simulated,based on the Euler equations and a detailed hydrogen-oxygen chemical reaction model.Unlike the oblique detonation wave induced by a wedge,the reacting fl w around a sphere is much more complex.First,a normal shock wave/detonation wave is formed, then oblique shock wave/detonation wave is developed in the presence of a rarefaction wave.The numerical simulation results show that after the gases being compressed by the blunt body and reaching the auto-ignition temperature,two kinds of fl wfiled will appear.When Mach numbers are low,the combustion will be quenched and can not appear downstream of the blunt body due to the influenc of the rarefaction wave.When Mach numbers are high,combustion can spreadto the downstream region.When the scales of blunt body are small,energy around the stationary point is not enough to induce detonation initiation and an obvious decoupling of combustion and shock wave is formed.As the sphere becomes large enough,decoupling of combustion and shock wave will not appear in the fl w and this feature is indpendent of the Mach number.By adjusting the spheric diameter,the fl w structures with partial coupling of shock wave and combustion zone was obtained which does not exist in a wedgy-induced oblique detonation.The present investigations suggest that the interaction between rarefaction wave and detonation wavefront is the key issue for detonation initiation induced by a spheric body.

oblique detonation,oxyhydrogen,initiation,quench

O381

A

10.6052/0459-1879-16-143

2016–05–25收稿，2016–12–09錄用,2016–12–09網絡版發(fā)表.

1)國家自然科學基金資助項目(11372333,91641130,11532014).

2)滕宏輝，副研究員，主要研究方向：激波與爆轟物理，爆轟燃燒推進，多相燃燒與爆炸.E-mail：hhteng@imech.ac.cn

方宜申，胡宗民，滕宏輝，姜宗林.圓球誘發(fā)斜爆轟波的數值研究.力學學報，2017,49(2):268-273

Fang Yishen,Hu Zongmin,Teng Honghui,Jiang Zonglin.Numerical study of the oblique detonation initiation induced by spheres.Chinese Journal of Theoretical and Applied Mechanics,2017,49(2):268-273