楊理,岳連捷,張新宇
1.北京理工大學 宇航學院,北京 100081 2.中國科學院 力學研究所 高溫氣體動力學國家重點實驗室,北京 100190 3.中國科學院大學 工程科學學院,北京 100049
爆轟燃燒是一種化學反應快、壓縮性強、傳播速度高的燃燒模式[1],在高超推進領(lǐng)域受到較高關(guān)注,其中,斜爆轟發(fā)動機作為一種結(jié)構(gòu)簡單、燃燒室更短、高馬赫數(shù)下具備高性能的爆轟推進方式[2-6],是過去幾十年空天飛行領(lǐng)域的研究熱點之一,圍繞斜爆轟波(Oblique Detonation Waves, ODW)的起爆、演化和駐定等問題已開展了大量的理論、實驗和數(shù)值計算研究,對ODW的基本流場結(jié)構(gòu)、起爆條件、演化模式和波陣面的不穩(wěn)定產(chǎn)生機制有較為充分的認識。
對于斜劈誘導的延遲斜爆轟波(Delayed ODW),其典型的流場結(jié)構(gòu)如圖1所示,由斜劈頂端發(fā)出的斜激波(Oblique Shock Wave, OSW),在壁面上經(jīng)一定長度誘導形成的燃燒波(Combustion Wave, CW)和ODW形成三波點(Triple-point)。三波點處OSW向ODW的轉(zhuǎn)換模式分為兩種:光滑(Smooth)和突變(Abrupt),學者對兩種轉(zhuǎn)換模式開展了大量的定性分析和定量預測[7-9]。類似于正爆轟波的胞格結(jié)構(gòu),ODW波陣面亦存在多維非穩(wěn)定性,Choi等[10]通過數(shù)值計算研究活化能對斜劈誘導ODW的影響,結(jié)果表明局部壓力擾動和壁面反射激波是促使類胞格結(jié)構(gòu)產(chǎn)生的因素。Verreault等[11]提出非穩(wěn)定ODW演化存在兩個階段:① 平面ODW的前導激波出現(xiàn)擾動造成化學反應區(qū)厚度的變化,進而誘導左行橫波的產(chǎn)生;② 前導激波的空間擾動演變?yōu)闀r間上的不穩(wěn)定性,產(chǎn)生傳向下游的擾動并與現(xiàn)有的左行橫波作用,最終促使右行橫波的產(chǎn)生。同時,正爆轟波的過驅(qū)因子概念亦可作為判斷ODW非穩(wěn)定性的判據(jù)。Teng等[12]基于數(shù)值計算發(fā)現(xiàn)不同活化能下斜爆轟波陣面存在兩種類胞格結(jié)構(gòu),并提出兩種類胞格結(jié)構(gòu)伴隨著不同的橫波模式。Liu等[13]對有限長度斜劈誘導的快速斜爆轟波(Prompt ODW)開展數(shù)值研究,結(jié)果表明ODW波陣面存在化學反應區(qū)和前導激波的完全耦合和部分耦合現(xiàn)象;特別地,下游的局部爆炸會使得耦合狀態(tài)由部分耦合轉(zhuǎn)換為完全耦合。Kasahara等[14]通過實驗觀察到高速射入可燃氣體的錐形彈丸誘導ODW波陣面可分為:① 強過驅(qū)斜爆轟波(Strong Overdriven ODW,SOODW); ② 弱過驅(qū)斜爆轟波(Weak Overdriven ODW, WOODW); ③ 準Chapman-Jouguet(CJ)斜爆轟波; ④ CJ斜爆轟波(CJODW)。
目前,另外一種理論亦提供關(guān)于波陣面?zhèn)鞑ヒ?guī)律的描述——爆轟激波動力學(Detonation Shock Dynamics,DSD)理論,其是一個基于法向速度演化規(guī)律的波陣面?zhèn)鞑ヮA測概念,且其預測結(jié)果與直接數(shù)值模擬結(jié)果具有很好的吻合度[15]。DSD理論假設爆炸系統(tǒng)的波陣面曲率遠小于化學反應區(qū)長度的倒數(shù),對于較小尺度的曲率,準穩(wěn)定爆轟波波陣面的法向速度Un可看作是曲率κ的函數(shù)。Yao和Stewart[16]基于理想狀態(tài)方程得到Un-κ關(guān)系類似于極曲線形狀,存在一個臨界曲率的偏轉(zhuǎn)點;進一步對基于非理想狀態(tài)方程的凝聚相炸藥開展相關(guān)理論分析[17],發(fā)現(xiàn)Un-κ關(guān)系為“Z”形曲線,在前期得到的極曲線形狀的Un-κ關(guān)系下支還存在一個凝聚相炸藥反應系統(tǒng)特有的低速爆轟區(qū)域。相關(guān)研究表明,Un-κ關(guān)系可對爆炸反應系統(tǒng)的爆轟波傳播特性提供新的認識[18-19]。
目前對ODW波陣面的非穩(wěn)定性已開展了較為充分的數(shù)值計算研究,發(fā)現(xiàn)了ODW的兩種類胞格結(jié)構(gòu),但對整個波陣面的幾何特性以及波陣面的Un-κ關(guān)系較少關(guān)注。本文在不同化學反應參數(shù)(釋熱量、釋熱速率和化學反應區(qū)參考長度)條件下,開展斜劈誘導斜爆轟波的高精度數(shù)值計算。分析沿波陣面的波角變化趨勢和Un-κ關(guān)系,以期對斜爆轟波陣面的類胞格結(jié)構(gòu)和其不穩(wěn)定性的演化機制加深理解。
圖1 延遲斜爆轟波的數(shù)值計算物理模型和流場結(jié)構(gòu)示意圖Fig.1 Schematic diagram of physical model for numerical computation and flow field structure of delayed oblique detonation wave
本文計算用求解器為自主編寫,可開展多塊計算域的高精度數(shù)值計算,控制方程是帶兩步化學反應的二維歐拉方程,未考慮黏性效應,其在貼體坐標系下形式為
(1)
式中:
(2)
直角坐標系(x,y)到貼體坐標系(ξ,η)的轉(zhuǎn)換關(guān)系為
(3)
式(1)~式(3)中:W為守恒解向量;F和G為對流項;S為化學反應源項。具體表達式為
(4)
其中:ρ、u、v、E、p分別為密度、x方向速度、y方向速度、單位質(zhì)量總內(nèi)能和壓力;λI和λR分別為誘導區(qū)和釋熱區(qū)的反應進度數(shù)(取值范圍:0~1)。
所采用的狀態(tài)方程為
(5)
式中:Q為化學反應釋放的熱量;γ為絕熱指數(shù)。
兩步化學反應模型中,誘導區(qū)和釋熱區(qū)的反應速率ωI和ωR分別為
(6)
式中:Ea為活化能;R為氣體常數(shù);T為溫度;υ為反應級數(shù),一般取作0.5;H(λI)為Heaviside函數(shù),控制誘導區(qū)和釋熱區(qū)的切換;KI和KR為相應的指前因子。該兩步化學反應模型的詳細介紹見參考文獻[20]。
為減少計算機截斷誤差對數(shù)值計算結(jié)果的污染,上述所有物理量參照來流參數(shù)進行無量綱化:
(7)
空間離散方法基于WENO格式[21]重構(gòu)左右值,本文采取耗散較小的五階WENO-Z格式[22],并采用具備較好間斷捕捉能力的Roe-HLLE通量分裂方法[23-24]計算對流項通量;對于離散后的常微分方程,為克服爆轟物理計算過程中存在的剛性問題,采用附加Runge-Kutta方法進行時間離散,所用到的butcher table為ARK4(3)6L[2]SA,時間推進格式可達到四階精度[25-26]。
本文計算所用模型示意圖如圖1所示,來流沿水平方向,斜劈角度固定為29°,為平衡計算量和保證計算域中ODW波陣面的長度,對個別算例的斜劈長度做出相應調(diào)整(分別為20、30、40 mm)。上方和左側(cè)邊界為超聲速來流邊界條件,下方斜劈為滑移固壁,右側(cè)邊界為零階外推出口邊界條件。表1給出不同化學反應參數(shù)(釋熱量、釋熱速率和化學反應區(qū)參考長度)條件下的6個算例,工況Q40.0L2.5K1、Q50.0L2.5K1和Q52.5L2.5K1的釋熱量由40.0增加到52.5,工況Q40.0L2.5K2和Q50.0L2.5K2相應的把釋熱速率增大為原來的兩倍,工況Q50.0L3.5K1增大化學反應區(qū)參考長度為3.5×10-4m,所有工況的來流馬赫數(shù)、壓力和溫度均相同,分別為:Ma∞=7,p∞=1.0×105Pa,T∞=300 K。所有算例采用均勻網(wǎng)格,網(wǎng)格尺寸一般為10 μm(相應的爆轟波化學反應區(qū)內(nèi)網(wǎng)格數(shù)為36);對于反應速率增大的算例,由于反應區(qū)尺度相應減小,采用網(wǎng)格大小為5 μm,確?;瘜W反應區(qū)內(nèi)網(wǎng)格數(shù)不小于36,以保證數(shù)值計算能很好地捕捉流場結(jié)構(gòu)。
特別地,考慮到燃燒過程中比熱和絕熱指數(shù)的變化,可燃物的氣體常數(shù)或比熱Ψ的計算表達式為
Ψ(λI,λR)=ΨuλI+Ψs(1-λI-λR)+ΨbλR
(8)
(9)
表1 計算工況和初始條件Table 1 Computation cases and initial conditions
式中:cp和cv分別為定壓比熱容和定容比熱容。
借助經(jīng)典CJ爆轟波理論分析的思路[1],不考慮化學反應動力學的影響和多維結(jié)構(gòu)作用,將化學反應釋熱簡化為能量添加,斜爆轟波波后狀態(tài)參數(shù)的理論值可通過Rankine-Hugoniot關(guān)系確定,其波角-偏折角(β-θ)關(guān)系為[27]
(10)
(11)
基于斜爆轟波分類理論[28],處于脫體角(最大偏折角)之上的半支稱為SOODW,而下半支則包含位于CJODW點右端的WOODW和左端的弱欠驅(qū)斜爆轟波(Weak Underdriven ODW, WUODW),如圖2所示。CJODW的波后法向馬赫數(shù)Ma2n=1.0,此時波角最小;由于過驅(qū)斜爆轟波(SOODW和WOODW)的Ma2n<1.0,使得下游的馬赫波與爆轟波能相交,因此過驅(qū)斜爆轟波能實現(xiàn)穩(wěn)定駐定(特別地,惰性斜激波的Ma2n?1.0),需要指出的是,SOODW的Ma2亦小于1.0;而WUODW的Ma2n>1.0,此時下游馬赫波與爆轟波不相交,導致WUODW無法自持傳播。
圖2 不同放熱量下的β-θ極曲線圖Fig.2 Polar diagram of β-θ for various amounts of heat release
圖3為表1中對應工況的數(shù)值計算結(jié)果,圖3 給出了整個流場的溫度云圖和三波點附近局部放大的數(shù)值紋影,并繪制Ma=1.0和λR=0.05的等值線。由圖3可發(fā)現(xiàn),隨著化學反應參數(shù)不同,斜劈誘導ODW在三波點處的波系相互作用類型和結(jié)構(gòu)存在明顯差異。圖3(a)和圖3(d)的流場結(jié)構(gòu)與圖1所示經(jīng)典流場結(jié)構(gòu)類似,但燃燒波陣面出現(xiàn)凸起折皺;圖3(b)、圖3(e)和圖3(f)中則在三波點處出現(xiàn)橫波(Transverse Wave, TW),其原因為:燃燒波后狀態(tài)與斜爆轟波后狀態(tài)無法匹配,使得圖3(a)和圖3(d)的波系結(jié)構(gòu)無法維持,促使橫波產(chǎn)生,類似于經(jīng)典激波/激波相互作用中的Type VI到Type V的轉(zhuǎn)變[29]。特別地,圖3(e)中橫波與燃燒波作用產(chǎn)生馬赫桿(Mach Stem, MS);需要指出的是,圖3(c)三波點處流場完全不同于圖1,此時一個反應馬赫桿(爆轟波)替代原有的燃燒波。
由式(10)可知Q=40.0,50.0,52.5時,脫體角所對應的波角βdetach分別為69.1°, 68.5°, 68.4°;CJODW的波角βCJODW分別為39.1°, 45.0°, 46.4°。本計算考慮到燃燒誘導過程的絕熱指數(shù)變化,即由斜劈頂端發(fā)出的斜激波前后的絕熱指數(shù)不同,從而計算得到的波角βOSW約為37.4°。圖4 給出上述3類波角(βOSW,βCJODW,βdetach)和數(shù)值計算得到的波陣面沿程波角變化關(guān)系。從圖4中可發(fā)現(xiàn),當前數(shù)值計算結(jié)果在ODW起爆(三波點)前的惰性斜激波部分的波角與理論值吻合較好。本文重點關(guān)注ODW起爆點后的波陣面波角變化,可將其分為3個區(qū)域,對應流場結(jié)構(gòu)見圖1標示,分別為:① 區(qū)域I,ODW衰退段(Decaying ODW),此時ODW波角隨著波陣面發(fā)展逐漸減小,波強度減弱,此過程中前導激波與化學反應區(qū)仍然高度耦合(見圖3中數(shù)值紋影的波陣面和λR=0.05釋熱面),不存在解耦現(xiàn)象。② 區(qū)域II,ODW與反射激波(Reflected Shock Wave, RSW)相互作用段(ODW/RSW interaction),波角在ODW/RSW作用點處存在明顯躍升,之后逐漸減??;RSW是由從三波點發(fā)出的斜激波經(jīng)壁面反射、與剪切層作用后形成,最終與下游的ODW相互作用,使得衰減的ODW增強,由于不同化學反應參數(shù)條件下RSW強度不同,該波角躍升程度存在較大差別。③ 區(qū)域III,類胞格結(jié)構(gòu)(Cellular-like structure),此時波角沿波陣面出現(xiàn)較大幅度、有規(guī)律的振蕩;該階段的ODW呈現(xiàn)非穩(wěn)定性,誘導左行橫波或右行橫波產(chǎn)生,波陣面上出現(xiàn)多波點結(jié)構(gòu)。
圖3 不同化學反應參數(shù)下的溫度云圖和三波點處的數(shù)值紋影(紅色實線:Ma=1.0,紫色實線:λR=0.05)Fig.3 Temperature contours and numerical schlieren near triple-point for different chemical kinetic parameters (red solid line: Ma=1.0, purple solid line: λR=0.05)
當釋熱量從40.0增至50.0和52.5(圖4(a)~圖4(c)),波角振蕩峰值增大,且圖4(b)~圖4(c)中部分波陣面的波角大于βdetach,說明存在SOODW,表明隨著放熱量增加,過驅(qū)因子f=(Ma1·sinβ/MaCJ)2(MaCJ為CJ爆轟馬赫數(shù))由1.379減小到1.273和1.257,ODW的不穩(wěn)定性大大增強[11-12]。相應的規(guī)律亦可從圖4(d)~圖4(e)中得到。增大放熱速率(圖4(a)與圖4(d)、圖4(b)與圖4(e)),與增大放熱量規(guī)律類似,波角振蕩峰值增大,放熱速率增加亦促使SOODW的出現(xiàn)(圖4(a)與圖4(d)的區(qū)域I);圖4(d)~圖4(e)中沿程波角的振蕩更為劇烈,說明增加放熱速率使得波陣面上的類胞格結(jié)構(gòu)的空間尺寸減小。數(shù)值計算中化學反應區(qū)參考長度增加,對比圖4(b)和圖4(f),發(fā)現(xiàn)ODW的起爆位置明顯延后,但對波角變化趨勢無明顯影響。
圖4(c)~圖4(e)中波陣面存在多處波角小于βCJODW,即該波角值不存在于對應的ODW極曲線上,選取典型波谷位置P1~P9進行分析。圖5為P1~P9處的局部放大數(shù)值紋影。由圖5(a)可知,P1點位于區(qū)域II的ODW/RSW 相互作用點附近,P1點之前的類胞格結(jié)構(gòu)的橫波為單一的左行橫波,且能發(fā)現(xiàn)明顯的前導激波與化學反應釋熱面解耦現(xiàn)象,至P1處,解耦現(xiàn)象更為顯著,可認為此處的ODW退化為激波誘導燃燒(Shock Induced Combustion, SIC)或斜激波,但是由于RSW存在,使得這一解耦過程得到抑制,并使波陣面與化學釋熱面重新耦合。而對于P2~P9點處的過小波角,其形成原因則為另一種情況:由圖5(b)~圖5(d)可知,P2~P9點附近的類胞格結(jié)構(gòu)中存在左行橫波和右行橫波,其波陣面結(jié)構(gòu)呈現(xiàn)為解耦ODW-斜激波-高度耦合ODW(類似于正爆轟波多波點結(jié)構(gòu)中的馬赫桿)交替演化,其中,高度耦合ODW會沿著波陣面向兩側(cè)成長(如圖5(b)所示),并與附近的多波點碰撞,形成新的類胞格結(jié)構(gòu),P2~P9則出現(xiàn)在高度耦合ODW或解耦ODW與斜激波的連接處。
圖4 不同化學反應參數(shù)下波角沿波陣面的變化Fig.4 Variations of wave angles along wave front for different chemical kinetic parameters
圖5 P1~P9位置處的局部數(shù)值紋影(紫色實線:λR=0.05)Fig.5 Local numerical schlieren for P1-P9 locations (purple solid line: λR=0.05)
根據(jù)2.1節(jié)對ODW波陣面的劃分,圖6以3種形狀的符號給出波陣面法向速度Un與曲率κ的關(guān)系,κ值的正負對應于波陣面形狀的擴張和收斂。對于區(qū)域I (ODW衰退段,圖6中黑色方塊符號),在極小的κ范圍內(nèi)(0值附近),Un急劇減小,Un-κ關(guān)系呈現(xiàn)出準垂直直線變化趨勢;而圖6(b)~圖6(f)對應流場中存在SOODW,使得ODW受到下游的壓力擾動影響,造成圖6(b)~圖6(f)中ODW衰退段中κ值出現(xiàn)小幅度波動;由于此區(qū)域的前導激波與釋熱面高度耦合,可認為波陣面呈現(xiàn)此Un-κ特性主要是受波陣面幾何形狀影響。
關(guān)于區(qū)域III(類胞格結(jié)構(gòu)段,圖6中綠色三角符號),以圖6(a)為例,由于ODW波陣面上出現(xiàn)多個類胞格結(jié)構(gòu),由不完全解耦ODW-左行橫波-高度耦合ODW周期演化構(gòu)成,在Un-κ圖上出現(xiàn)多個高度自相似的“D”形曲線。圖6(a)中給出一個完整類胞格結(jié)構(gòu)的波陣面形狀,并標出波陣面上的相應點A1~A6,A1~A3間的Un-κ關(guān)系近似為極曲線,與Yao和Stewart[16]的理論分析類似,存在一個極大曲率點A2,位于左行橫波與ODW波陣面作用點附近,此處溫度壓力由于左行橫波的作用而急劇變化,使得前導激波與化學反應面高度耦合。可認為A1~A3段連接著不完全解耦ODW和強耦合ODW,為一過渡結(jié)構(gòu)。而對于A2~A4段的A3附近(強耦合ODW),κ存在較大跨度但Un值變化較小,Un-κ關(guān)系呈現(xiàn)為平滑水平線,與經(jīng)典DSD理論里理想爆轟的Un-κ關(guān)系類似。對于A4~A6段,特別是A5~A6段,Un-κ關(guān)系與區(qū)域I的趨勢類似,但Un和κ值變化范圍增大,Un-κ擬線性關(guān)系的傾斜程度更大,主要是由于前導激波與釋熱面解耦程度越來越大,使得ODW強度減弱相較于區(qū)域I更為顯著。
圖6 不同化學反應參數(shù)下的Un-κ變化關(guān)系Fig.6 Un-κ relations for different chemical kinetic parameters
隨著釋熱量的增大(圖6(a)~圖6(c)),由于SOODW的出現(xiàn),使得橫波作用點附近出現(xiàn)極小的亞聲速區(qū),局部ODW波陣面同時受到波后化學反應和下游流場擾動影響,Un-κ的典型“D”形關(guān)系不是十分規(guī)則,且多個類胞格結(jié)構(gòu)的自相似性減小,說明經(jīng)典的DSD理論引入的Un-κ關(guān)系對于強不穩(wěn)定爆轟可能不完全適用。當增大釋熱速率(圖6(a)和圖6(d),圖6(b)和圖6(e)),Un和κ值的范圍適當增大,特別是對于較小釋熱量(Q=40.0)的圖6(a)和圖6(d),其變化趨勢基本相似。反應區(qū)參考長度lref從2.5×10-4m增加到3.5×10-4m,Un和κ值的范圍變小,但對Un-κ關(guān)系影響較小。
對于區(qū)域II (ODW與RSW相互作用段),其為類胞格結(jié)構(gòu)之特例,即RSW扮演著橫波的作用,考慮到ODW和RSW激波強度變化、二者作用產(chǎn)生的ODW是否為SOODW、新產(chǎn)生的過驅(qū)ODW會經(jīng)歷類似于區(qū)域I的衰減以及波陣面在空間上的壓力擾動等諸多因素,該區(qū)域流場變化較為復雜,可看作區(qū)域I和區(qū)域III二者的耦合。在此區(qū)域內(nèi),圖6中存在κ值跨度不一的“D”形曲線,或較小(圖6(a)、圖6(d)和圖6(e)),或較大(圖6(b)和圖6(f)),該區(qū)域其他部分總體上則呈現(xiàn)一定程度的擬線性變化趨勢。需要說明的是,圖6(c) 中,由于ODW/RSW作用點出現(xiàn)在類胞格結(jié)構(gòu)下游,此處未作詳細討論。
初步研究結(jié)果表明,對于斜爆轟波波陣面的Un-κ存在多種變化關(guān)系,并且這種彎曲效應與類胞格結(jié)構(gòu)的演變密切相關(guān),隨著波陣面上橫波作用,使得解耦ODW演變?yōu)楦叨锐詈螼DW,對應著不同的Un-κ變化趨勢,并呈現(xiàn)周期性規(guī)律;后續(xù)研究中,可對經(jīng)典的DSD理論進行發(fā)展和改進,結(jié)合當前斜爆轟波波陣面的Un-κ規(guī)律,以描述波陣面的演化傳播。
本文通過求解帶兩步化學反應的二維歐拉方程,在不同化學反應參數(shù)(釋熱量,放熱速率和化學反應區(qū)參考長度)條件下,對斜劈誘導斜爆轟問題開展數(shù)值研究,得到的主要結(jié)論如下:
1) 沿波陣面的波角變化曲線可分為3個區(qū)域:區(qū)域I,波角平滑減小,為衰退的斜爆轟波;區(qū)域II,波角躍升后衰減,為斜爆轟波與反射激波相互作用;區(qū)域III,波角有規(guī)律振蕩,為類胞格結(jié)構(gòu)的周期演化。
2) 某些位置的波角會出現(xiàn)小于對應的CJ斜爆轟波波角,原因有兩個:① 此處的斜爆轟波退化為激波誘導燃燒或斜激波;② 在同時存在左行橫波和右行橫波的類胞格結(jié)構(gòu)附近,此處連接高度耦合斜爆轟波或解耦斜爆轟波與斜激波。
3) 法向速度Un和波陣面曲率κ的關(guān)系在區(qū)域I呈現(xiàn)出準垂直直線變化趨勢,由于強過驅(qū)斜爆轟波的存在,曲率值存在小幅擾動;區(qū)域III,Un-κ關(guān)系呈現(xiàn)出“D”形曲線,其伴隨著斜爆轟波的解耦、橫波作用使得斜爆轟波重新高度耦合和斜爆轟波的再次解耦的周期性演變,該“D”形曲線關(guān)系在弱過驅(qū)爆轟波陣面中存在較高的自相似性和更為規(guī)則,但對于強不穩(wěn)定爆轟則不是十分適用;區(qū)域II,可看作區(qū)域I和區(qū)域III二者的耦合作用。
4) 在本文研究框架下,化學反應參數(shù)中的釋熱量改變(Q從40.0增加到52.5)對波角和Un-κ關(guān)系影響極為顯著;放熱速率增大到原來的2倍,波角和Un-κ關(guān)系的變化趨勢基本不變,但在數(shù)量值上有少許變化;化學反應區(qū)參考長度增加到原來的1.4倍,斜爆轟起爆點位置明顯推后,但其波角和Un-κ關(guān)系無明顯變化。