王愛峰,趙 偉,姜宗林

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

斜爆轟是超聲速的可燃氣體通過一定角度的楔面經(jīng)過斜激波壓縮、燃燒誘導(dǎo)直至快速放熱的一種特殊的爆轟現(xiàn)象。由于斜爆轟波可以駐定,燃燒時間極短且擁有很高的燃燒效率,因此在超燃沖壓推進系統(tǒng)(supersonic combustion ramjet,scramjet)和沖壓加速器(ram accelerator)[1]等方面都有良好的應(yīng)用前景。

D.T.Pratt 等[2]主要應(yīng)用傳統(tǒng)的爆轟理論、簡化的物理模型,研究馬赫數(shù)、放熱量、楔面角對斜爆轟駐定的影響,結(jié)果表明:對于固定當(dāng)量比的可燃氣體(即固定放熱量),不同的馬赫數(shù),形成駐定斜爆轟需要楔面角滿足一定的范圍,這個范圍為CJ 線與極限角連線所夾的區(qū)域(如圖1 所示)。H.F.Lehr[3]在觀測高速彈丸在氫氧介質(zhì)中飛行時,發(fā)現(xiàn)彈頭產(chǎn)生斜激波誘導(dǎo)的燃燒現(xiàn)象,觀測得到不同的燃料當(dāng)量比、飛行馬赫數(shù)以及彈丸結(jié)構(gòu)等條件下出現(xiàn)的不同的斜爆轟波結(jié)構(gòu)。T.Fujiw ara 等[4]模擬了無粘條件下的斜爆轟結(jié)構(gòu),得到斜爆轟結(jié)構(gòu)中出現(xiàn)斜激波、爆燃波與斜爆轟波共存的復(fù)雜結(jié)構(gòu)(如圖2 所示),E.K.Dabora 等[5]則通過斜激波管實驗觀測到不同氣體條件下的斜爆轟波的結(jié)構(gòu)。由于楔面誘導(dǎo)的斜爆轟結(jié)構(gòu)中涉及爆燃波與爆轟波之間復(fù)雜的相互作用,隨著對斜爆轟現(xiàn)象認識的深化,近年來主要的研究集中在斜爆轟內(nèi)部流場的精細研究。L.F.F.da Silva 等[6]利用氫氧混合基元反應(yīng)模型,分析不同楔面角、不同馬赫數(shù)、不同氫氧當(dāng)量比條件下斜爆轟波的發(fā)展演化,并對誘導(dǎo)時間與誘導(dǎo)長度進行研究。C.Viguier 等[7]通過實驗,發(fā)現(xiàn)斜爆轟中的胞格結(jié)構(gòu)的橫波傳播特征(見圖3),J.Y.Choi 等[8]的數(shù)值研究發(fā)現(xiàn)了高活化能狀態(tài)下的特殊的胞格結(jié)構(gòu),同時發(fā)現(xiàn)橫波的單向傳播現(xiàn)象。

圖1 斜爆轟駐定匹配關(guān)系(Q～=1.96)Fig.1 Matching relationship of standing ODW

圖2 斜激波-斜爆轟波結(jié)構(gòu)示意圖Fig.2 Schematic diagram of OSW-ODW structure

為了探索斜爆轟形成機制并分析其中的胞格結(jié)構(gòu),本文中利用數(shù)值方法,分析在30°楔角條件下,臨界爆轟馬赫數(shù)附近6.8、7.0 與7.5 的3 種馬赫數(shù)狀態(tài)下的斜爆轟結(jié)構(gòu),觀察斜爆轟波陣面產(chǎn)生的橫波的傳播特征。

1 控制方程與數(shù)值方法

1.1 控制方程

忽略粘性、熱傳導(dǎo),可以采用二維歐拉方程耦合化學(xué)反應(yīng)過程描述斜爆轟波?；瘜W(xué)反應(yīng)模型采用了總包化學(xué)反應(yīng)模型,通過反應(yīng)進度常數(shù)λ反應(yīng)化學(xué)反應(yīng)的完全度。

控制方程為

圖3 斜爆轟波結(jié)構(gòu)以及橫波傳播的紋影圖[7]Fig.3 Schlieren photograph of ODW structure and propagation of transverse wave

式中:ρ為混合氣體密度,u 和v 分別為x 和y 方向上的速度,E 為質(zhì)量總能量,λ為反應(yīng)進程度。p 為氣體壓力,滿足以下關(guān)系

式中:γ為絕熱指數(shù),q 為化學(xué)反應(yīng)生成熱?？紤]單步不可逆反應(yīng)模型,化學(xué)反應(yīng)生成源項表達式為

式中:Ea 為化學(xué)反應(yīng)的活化能,k 為化學(xué)反應(yīng)速率系數(shù)。

1.2 數(shù)值方法

控制方程的對流項采用三階精度ENO 格式[9]離散,考慮迎風(fēng)效應(yīng),采用Steger-Warming 流通量分裂[10]。時間采用三階精度TVD Runge-Kutta 方法進行迭代,化學(xué)反應(yīng)項采用解耦的方式處理。

所選計算區(qū)域為圖4 中矩形區(qū)域,并按照楔面角度進行旋轉(zhuǎn)。計算域大小為10 cm×4 cm,θ為楔面角,β 為斜爆轟角。為了研究斜爆轟波的復(fù)雜的細部結(jié)構(gòu),采用2 000×800 的網(wǎng)格進行計算,同時在入口邊界增加3 個虛點,保證上游的入口超聲速流體不會在楔面上形成數(shù)值反射,楔面采用滑移邊界條件,出口邊界為無擾動外插邊界。

為了研究馬赫數(shù)對斜爆轟胞格結(jié)構(gòu)的影響,采用如下的物理參數(shù):初始壓力p0=101.325 kPa,初始溫度T0=300 K,楔面角θ=30°,馬赫數(shù)Ma=6.8,7.0,7.5,放熱量Q=965 kJ/kg,Ea=4.794 MJ/kg,k=7.5 ns-1,γ=1.29。所設(shè)定的馬赫數(shù)、放熱量與楔面角關(guān)系均滿足斜爆轟駐定條件[2]。

圖4 斜爆轟結(jié)構(gòu)示意圖Fig.4 Schematic diagram of ODW structure

2 數(shù)值驗證

2.1 斜爆轟角驗證

為了驗證不同狀態(tài)下斜爆轟角的變化,假定平面斜爆轟:(1)來流為預(yù)混未燃氣體,混合均勻,流動穩(wěn)定;(2)化學(xué)反應(yīng)層等效于放熱,放熱量為Q;(3)流動無粘絕熱。根據(jù)上面的假設(shè)條件,利用流體守恒方程,可以得到下面關(guān)于放熱量、楔角、馬赫數(shù)和斜爆轟角之間的關(guān)系

式中:(Ma)1 為斜爆轟波前馬赫數(shù),(Ma)1n 為波前法向馬赫數(shù),θ為楔面角,β 為斜爆轟角,為量綱一放熱量,cp為定壓比熱,T1為波前溫度,γ為絕熱指數(shù)。

楔面角θ=30°、馬赫數(shù)不同、放熱量不同時,斜爆轟角的理論解和數(shù)值結(jié)果及相對誤差見表1。由表可見,數(shù)值解比理論解略大,最大相差低于1.25%,說明數(shù)值解的結(jié)果是可信的。

2.2 斜爆轟波流場驗證

為了進一步驗證斜爆轟結(jié)構(gòu)中復(fù)雜的結(jié)構(gòu),將數(shù)值結(jié)果與實驗結(jié)果進行對比,圖5 是實驗紋影圖[5]與數(shù)值密度云圖的對比結(jié)果?？梢钥吹降湫偷男奔げ?、斜爆轟波以及斜激波、斜爆轟波與爆燃波交匯成的三波點。由于斜爆轟波與爆燃波波后的氣體密度不同形成一條向下游發(fā)展的接觸間斷面。因此所采用的數(shù)值計算能夠?qū)π北Z波結(jié)構(gòu)進行詳細的解析。

表1 斜爆轟角Table 1 ODW angles

圖5 實驗紋影圖[5]和數(shù)值密度云圖(p0=101.325 kPa,T 0=293 K)Fig.5 Schlieren image of the flow field and density field

3 數(shù)值結(jié)果與討論

3.1 不同馬赫數(shù)條件下的溫度流場

圖6 為放熱量Q=965 kJ/kg、馬赫數(shù)為6.8、7.0 和7.5 的溫度云圖。為了保證流場定常,每種狀態(tài)下的密度殘差都小于10-4。從圖中可以看出,馬赫數(shù)為7.5 時,上游的斜激波區(qū)域長度較短,斜激波與斜爆轟波過渡平滑,流場中的胞格結(jié)構(gòu)規(guī)則且胞格較小,波陣面呈現(xiàn)一條直線。而馬赫數(shù)為7.0、6.8時,斜激波區(qū)域長度相對較長,波后流場存在較大的擾動,下游的胞格結(jié)構(gòu)不再規(guī)則且胞格大小不一,斜爆轟波陣面總體趨向于直線,但局部出現(xiàn)彎曲現(xiàn)象。

馬赫數(shù)的差別不大,但流場的結(jié)構(gòu)存在較大的波動。進一步觀察三波點的結(jié)構(gòu)發(fā)現(xiàn),馬赫數(shù)為7.0、6.8 時斜激波區(qū)較長,且三波點附近的流場結(jié)構(gòu)較為復(fù)雜;而馬赫數(shù)為7.5 時,三波點附近的流場較為簡單,因此有必要對主流三波點附近的流場結(jié)構(gòu)進行細致的研究。

圖7 為馬赫數(shù)為7.0 時的三波點附近的流場圖,從圖中可以看出,三波點后會產(chǎn)生一道激波,這道激波傳播到楔面發(fā)生反射,并能透過滑移間斷面,與爆轟波陣面產(chǎn)生的橫波相互作用后下游流場受到擾動,使得胞格結(jié)構(gòu)呈現(xiàn)不規(guī)則的狀態(tài)。同時爆轟波陣面產(chǎn)生的橫波傳播到楔面后發(fā)生了反射,使得接觸間斷面下游出現(xiàn)一系列的渦結(jié)構(gòu)。

為了進一步分析,截取y=3 mm 的水平方向的壓力和化學(xué)反應(yīng)度曲線(見圖8)。圖中圓點部分為三波點后產(chǎn)生的一道激波,方點部分為激波在楔面上發(fā)生反射形成的反射激波。A 區(qū)壓力經(jīng)歷突變,最終化學(xué)反應(yīng)度為1,熱量完全釋放。B 區(qū)壓力下降,是爆燃區(qū)域。C 區(qū)為激波與反射激波包圍的混合區(qū)域。C 區(qū)之后由于爆轟前導(dǎo)反應(yīng)段面上產(chǎn)生的橫波在楔面上發(fā)生反射,形成振蕩的D 區(qū)。

3.2 波陣面結(jié)構(gòu)特征

前面的研究發(fā)現(xiàn),在馬赫數(shù)為7.0、6.8 時流場中出現(xiàn)不規(guī)則的胞格結(jié)構(gòu),而馬赫數(shù)為7.5 時胞格結(jié)構(gòu)卻非常規(guī)則。三波點后產(chǎn)生的激波能夠與波陣面上產(chǎn)生的橫波相互作用,并產(chǎn)生擾動影響著下游流場的結(jié)構(gòu),同時橫波傳播到楔面也能發(fā)生反射,使滑移間斷面呈現(xiàn)一系列的渦結(jié)構(gòu)。這些現(xiàn)象的出現(xiàn)都與橫波存在很大的關(guān)系,因此對波陣面上橫波傳播的研究是非常必要的。

分別對馬赫數(shù)為7.5 和6.8 等2 種情況下的爆轟波陣面的結(jié)構(gòu)進行細部分析,如圖9 ～10 所示。當(dāng)馬赫數(shù)為7.5 時,斜爆轟波陣面上產(chǎn)生的橫波均向同一個方向傳播(如圖11(a)所示),這與正爆轟的橫波的雙向傳播特性不同,但這種單向傳播的橫波與文獻[8]一致。當(dāng)馬赫數(shù)為6.8 時,由于上下游波陣面的結(jié)構(gòu)相異,因此分為上游A 區(qū)和下游B 區(qū)(見圖6(a)),研究表明,波后流場的上游橫波仍然是單向傳播(如圖11(a)所示),但是到了下游,隨著馬赫干的出現(xiàn),橫波的傳播不再單向,而是類似于正爆轟中雙向傳播的橫波結(jié)構(gòu)(如圖11(b)所示),這種橫波單向與雙向傳播的同時出現(xiàn)在以前沒有發(fā)現(xiàn)過。

圖6 相同放熱量、不同馬赫數(shù)條件下溫度云圖(Q=965 kJ/kg)Fig.6 Temperature contour at different Mach numbers

圖7 三波點附近的流場圖(Ma=7.0,Q=965 kJ/kg)Fig.7 Flow field near primary triple point area

圖8 沿楔面方向壓力與化學(xué)反應(yīng)度分布(y=3 mm,Ma=7.0,Q=965 kJ/kg)Fig.8 Pressure and degree of chemical reaction profile along w edge direction

馬赫數(shù)越大,加入的動能越大,通過壓縮獲得的溫度越高,所得斜激波部分的長度越短,因此在一定尺度范圍內(nèi),三波點后產(chǎn)生的激波對波后流場影響也很小,因此馬赫數(shù)為7.5 時波后流場結(jié)構(gòu)均勻,同時橫波只出現(xiàn)單向傳播的機制。但是馬赫數(shù)為7.0、6.8 時,橫波是由單向轉(zhuǎn)為雙向傳播的,兩種傳播狀態(tài)并存,如果上述的這種橫波的傳播方式是斜爆轟中一般的傳播規(guī)律,那么馬赫數(shù)為7.5 時必然也會出現(xiàn)。如果這個假定成立,那么影響這個方面的因素必然是長度尺度,因此增加長度尺度來考察馬赫數(shù)為7.5 時的波陣面形成的橫波的傳播情況。

長度由10 cm 增加到15 cm、其余參數(shù)不變時的流場如圖12 ～13 所示。由圖13 可以看出,在11 cm后就開始出現(xiàn)馬赫干,隨后橫波的雙向傳播越來越明顯。因此橫波由單向向雙向傳播的轉(zhuǎn)變不是偶然的,它與長度尺度相關(guān),是其中復(fù)雜波系逐漸削弱的一種體現(xiàn)。

圖9 局部溫度云圖(Ma=7.5,Q=965 kJ/kg)Fig.9 Detailed temperature contour

圖10 局部溫度云圖(Ma=6.8,Q=965 kJ/kg)Fig.10 Detailed temperature contour

圖11 波陣面簡圖Fig.11 Schematic of w ave front structure

圖12 溫度云圖(Ma=7.5,Q=965 kJ/kg,L=15 cm)Fig.12 Temperature contour

圖13 局部溫度云圖(Ma=7.5,Q=965 kJ/kg,L=15 cm)Fig.13 Detailed temperature contour

4 結(jié) 論

研究了30°楔角條件下,臨界爆轟馬赫數(shù)附近的3 種來流狀態(tài)時斜爆轟的胞格結(jié)構(gòu)以及橫波的傳播機制,得到如下主要結(jié)論:

(1)馬赫數(shù)為7.0、6.8 時,誘導(dǎo)長度較長,一定尺度范圍內(nèi),斜爆轟的胞格結(jié)構(gòu)呈現(xiàn)不規(guī)則的狀態(tài)。

(2)主流三波點后形成一道激波,并在楔面上發(fā)生反射,穿越滑移間斷面與斜爆轟波陣面產(chǎn)生的橫波相互作用,使得橫波傳播發(fā)生偏向,對下游流場產(chǎn)生擾動。

(3)斜爆轟波陣面產(chǎn)生的橫波在上游呈現(xiàn)單向傳播狀態(tài),而在下游過渡到雙向傳播。增加長度尺度,馬赫數(shù)為7.5 時,在11 cm 處也開始出現(xiàn)雙向傳播的橫波,說明上述橫波由單向向雙向傳播的轉(zhuǎn)變不是偶然的,傳播的狀態(tài)依賴于長度尺度。

[1] Sasoh A.Laser-propelled ram accelerator[J].Journal de Physique Ⅳ,2000,10:11-14.

[2] Pratt D T,H umphrey J W.Morphology of standing oblique detonation w aves[J].Propulsion and Pow er, 1991,7(5):837-845.

[3] Lehr H F.Experiments in shock induced combustion[J].Astronautica Acta,1972,17(4-5):589-597.

[4] Fujiw ara T, Matsuo A.A tw o-dimensional detonation supported by a blunt body or a w edge[R].AIAA Paper 88-0098,1988.

[5] Dabora E K,Desbordes D,Wagner H G.Oblique detonation at hypersonic velocities[C]∥Borisov A A,Kuhl A L,Leyer J C,et al.Dynamics of Detonations and Explosions:Detonations,Progress in Astronautics and Aeronautics.Washington:AIAA,1991:187-201.

[6] da Silva L F F, Deshaies B.Stabilization of an oblique detonation wave by a wedge:A parametric numerical study[J].Combustion and Flame,2000,121(1-2):152-166.

[7] Viguier C,Gourara A,Desbordes D.Onset of oblique detonation waves:Comparison between experimental and numerical results for hydrogen-air mixtures[C]∥Proceedings of Twenty-Seventh Symposium(International)on Combustion.Pittsburgh:The Combustion Institute,1998:3023-3031.

[8] Choi J Y, Kim D W,Jeung I S, et al.Cell-like structure of unstable oblique detonation w ave f rom high-resolution numerical simulation[J].Proceeding of the Combustion Institute,2007,31(2):2473-2480.

[9] Shu C W,Osher S.Efficient implementation of essentially non-oscillatory shock capturing schemes Ⅱ[J].Computational Physics,1989,83(1):32-78.

[10] Steger J L, Warming R F.Flux vector splitting of the inviscid gasdynamic equations with application to finitedifference methods[J].Computational Physics,1981,40(2):263-293.