周期性非均勻介質中氣相爆轟波演變模式研究1)

陳 達 寧建國 李 健

(北京理工大學爆炸科學與技術國家重點實驗室,北京 100081)

引言

氣相介質通常存在非均勻性,區別于經典的氣相爆轟理論,非均勻氣相介質的爆轟機理研究面臨許多挑戰性的基礎理論問題.研究非均勻氣相介質的爆轟過程具有極為重要的實際應用價值,可以為工業安全防護和爆轟推進設計提供重要的理論和技術支持.

在可燃氣體泄露后的流場中,非均勻性以各個熱力學狀態參數的空間分布梯度的形式體現,學者們研究了氣體組分濃度梯度[1-3]或者溫度梯度[4-5]對爆轟波波陣面和傳播極限的影響.當前對于非均勻氣相介質爆轟波的研究主要集中于熱力學狀態參數的空間分布存在梯度這一模型,進而研究梯度分布對爆轟波波陣面和傳播極限的影響.在可燃氣體泄露后的流場和爆轟發動機的工作過程中,存在著多方向的組分濃度梯度或者溫度梯度.國內外很多學者對垂直于傳播方向存在組分濃度梯度的氣體混合物的爆轟傳播過程,特別是波陣面的演化過程,進行了實驗[6-8]和數值模擬[9-11]研究.Ishii 和Kojima[6]研究了波前介質非均勻對爆轟波陣面形狀和爆轟不穩定性的影響,結果表明,波陣面會形成一偏轉角,并且該角度隨局部梯度的增加而增加.Ettner 等[12]研究了直管中不同濃度梯度對爆轟波陣面、爆轟穩定性和壓力分布等的影響,結果表明,壁面的約束使得波陣面出現彎曲或者馬赫反射,但不會出現規則反射.藺偉等[13]研究了濃度梯度對瓦斯爆炸影響的數值模擬,以及爆轟波形成的胞格大小和三波點隨梯度的變化規律,發現由于濃度梯度的存在,波陣面上會出現馬赫反射現象,同時研究了馬赫反射對整個傳播過程的影響.Han 等[14]通過高精度數值模擬和一步反應模型研究了縱向存在組分梯度情況下爆轟波的波陣面結構、流場細節和傳播極限問題.

同時很多學者也對平行于爆轟波傳播方向存在組分濃度梯度的情況進行了研究.Kuznetsov 等[15-16]便對組分活性梯度降低的爆轟傳播進行了類似的研究,并研究了驅動段長度的影響.Thomas 等[17]則發現小的組分濃度梯度有助于爆轟波的透射,而大的梯度則可能使激波陣面和反應區分離,導致爆轟波熄爆.王春等[18]對爆轟波穿越惰性氣團時的透射激波進行了參數分析.他們發現當爆轟波穿越惰性氣團再進入可燃混合氣體時,透射激波的參數主要受到3 個基本過程的影響.

上述研究中,熱力學參數的非均勻性是來自整體空間分布上的梯度.這種梯度分布模型雖然更接近實際的情況,但是并不利于建立非均勻介質爆轟波的穩態傳播機理.另外一種思路則是局部上添加擾動,爆轟波此時會以一種近似穩態的模式傳播,波陣面近似平面(雖然存在小尺度的三波結構).例如,Mi 等[19]通過在均勻介質中添加隨機分布的離散單元實現非均勻性,離散單元的存在使得爆轟波波陣面更加扭曲、離散,同時也促進了爆轟波的傳播,使其不容易熄爆失效.Mi 等[20-21]和Lam 等[22]則利用一個一維熱方程研究了具有無限劉易斯數和空間離散放熱源的系統中火焰傳播的動力學問題,該方程的源項為反應過程變量,該變量的擴散率為零,模型揭示了火焰傳播模式的譜系.Li 等[23]通過大振幅、二維密度正弦波動引入非均勻性,對薄炸藥層中爆轟波傳播極限附近的爆速進行了計算.結果表明,非均勻介質中的爆轟波具有與均勻介質中類似的胞格結構.爆轟能夠在一個非常薄的炸藥層傳播而不熄爆,并且比相應的均勻介質呈現出更大地傳播速度.Higgins[24]通過在化學反應速率模型的基礎上乘以一個與當地坐標有關的周期函數來實現人工熱點,作為均勻反應速率的替代方法,分析了在含能材料中離散效應對爆轟波傳播的影響.但是這些研究存在一個問題,即他們使用的反應模型為壓力敏感型方程,這使得胞格結構不易顯現出來,這有利于研究波陣面曲率的影響,但是忽略了胞格不穩定性和橫波結構的影響,這種假設更適用于凝聚相爆轟.

另外,在爆轟發動機的工作過程中,學者針對于非均勻性對于爆轟發動機的影響也開展了相應的數值仿真[25]與實驗[26]研究.例如,分析了不同的非預混過程對組分濃度分布的影響規律[27],以及湍流效應在這一過程中的作用[28].此類研究通常不考慮點火以及其后的爆轟過程,也有很多研究涉及發動機非預混爆轟的整個過程,如點火過程[29]、波頭演化[30],進而研究非均勻性對發動機推力的影響[31].

研究非均勻氣相介質的爆轟過程具有極為重要的實際應用價值,可以為工業安全防護和爆轟推進系統設計提供重要的理論支持.在非均勻性存在的前提下,氣相介質的爆炸過程,包括直接起爆,穩態傳播和失效機制,都極其復雜,很多物理機制尚不明確.當前的很多研究中,由于實驗手段的限制,很難對周期性非均勻介質中氣相爆轟波傳播機制進行實驗研究,而數值模擬可以容易地實現這一目的,便于對物理機制的理解.基于此,本文使用反應歐拉方程和兩步反應模型對氣相爆轟波在非均勻介質中的傳播過程進行數值模擬,主要探究非均勻介質橫向均勻溫度擾動對氣相爆轟波的傳播機理進行系統的參數研究,特別是非均勻性將以何種形式影響爆轟波的直接起爆、傳播和失效以及波陣面多波精細結構,這種影響又是否與爆轟波自身的不穩定性和動力學參數存在某種定量的關系,將是本文研究的重點.相對于復雜的基元反應模型,兩步法可以容易地實現爆轟參數的解耦,方便系統地進行單一參數研究,而又不失一般性.

1 控制方程和數值方法

1.1 控制方程

在本文中控制方程使用耦合了兩步化學反應模型的歐拉方程.忽略了熱傳導、分子擴散以及黏性效應的影響,控制方程表示為

其中為U為守恒變量,F和G為通量,S為反應源項,分別寫為

對于理想氣體,總能和狀態方程可表示為

其中p為壓力,ρ為密度,E為單位體積的能量,u和v分別為x和y方向上的速度,Q為反應熱,yI與yR分別為兩步反應模型中表征誘導和放熱過程的進程參數.

1.2 兩步化學反應模型

化學反應采用兩步反應模型[32],即先進行誘導過程,誘導過程結束后再進行放熱反應過程.第一步誘導過程速率可表示為

其中EI為誘導活化能,Ts為穩態ZND (Zel’dovich-von Neumann-D?ring)爆轟波波陣面后的溫度,H(1-yI)為階躍函數,形式如下

誘導過程結束后,第二階段為放熱反應階段,反應速率表示為

其中ER為反應階段的活化能,指前因子分別為誘導階段和放熱階段的反應速率常數,可以控制誘導區與反應區的長度.相關的爆轟參數見表1,在本文的研究中,爆轟波的參考長度xref=0.001 m,通過選擇特定的指前因子,誘導區寬度 ΔI=0.001 m.通過改變指前因子的數值,可以得到不同反應區寬度的爆轟波.誘導區寬度固定,反應區寬度改變,爆轟波的不穩定性也因此改變.kR的值越大,反應區寬度越小,爆轟波也越不穩定,胞格結構也越不規則.馬赫數和反應熱的關系為

表1 爆轟參數Table 1 Detonation parameters

為了便于與文獻比較,這里定義新的活化能參數,即

一般來說,對于常見的碳氫燃料,活化能 εI的值在4～ 12 之間,并且滿足 εI?εR.

1.3 計算方法

為了分辨爆轟波陣面精細的多波結構,必須使用足夠高的網格分辨率.為了滿足這一要求,在本研究中使用了基于塊結構化自適應網格方法(SAMR)的AMROC 程序[33-34].AMROC 是一個基于有限體積法開發的開源框架,廣泛應用于燃燒、爆轟和空氣動力學的模擬研究,而且用戶可以在該框架下嵌入自定義的算法和反應模型.本研究在AMROC 中耦合了兩步化學反應模型,使其可以進行相關的爆轟研究.在整個計算區域,使用粗網格,在流場劇烈變化的區域使用細網格,使之具有更高的網格分辨率,以減少計算時間和內存需求.針對反應歐拉方法,使用兩階Strang 分裂方法解耦對流項和化學反應源項,然后采用基于Roe-HLL 求解器的顯式二階MUSCL-Hancock 方法求解對流項.反應源項采用半隱式的4 階通用Runge-Kutta 剛性常微分求解器(GRK4A)[35]進行數值積分.

1.4 初始設定

在本研究中,爆轟波的比熱比 γ =1.44,初始CJ(Chapman-Jouguet)爆轟波的馬赫數MCJ=5.6,無量綱的誘導區和反應區活化能分別為 εI=4.8,εR=1.0,其他相關參數詳見表1.阿倫尼烏斯定律的指前因子kI,kR可以用于控制爆轟波的誘導區和反應區寬度.通過Ng 等[32]提出的不穩定因子χ=εIΔI/ΔR來表征溫度擾動的化學敏感程度.另外,在本次數值模擬中使用了多種不同穩定因子的爆轟波作為初始設定.

本研究針對非均勻介質中爆轟波發展分別進行了一維和二維數值模擬.一維計算域設定如圖1(a)所示,左側為一ZND 爆轟波,右側為溫度T1的熱區或溫度為T2的冷區.當待起爆區為熱區時,表達式為T1=T0+A,為冷區時,表達式為T2=T0-A,T0為基準溫度295 K,A為溫度幅值.二維周期性擾動初始設定則如圖1(b)所示,橫向長度為 1 600ΔI(1.6 m),縱向長度為 3 00ΔI(0.3 m).左側初始設定為具有不同穩定性的ZND 爆轟波,右側則在待起爆區沿縱向添加了人工溫度擾動,以便實現非均勻性,溫度擾動以矩形波的形式體現,其中溫度振幅A以均值295 K 為基準上下浮動,而波長為s則分別取值10 ΔI,20ΔI和40 ΔI.關于邊界條件,左右兩側為出流邊界,上下兩側則為滑移固壁邊界.

圖1 計算域設置Fig.1 Setup of the simulation domain

1.5 收斂性測試

如圖2 所示,首先對于一維爆轟波的計算進行網格收斂性測試.相同初始條件下,對于5 種不同的網格分辨率進行了對比,分別為4pts/ΔI,8pts/ΔI,16pts/ΔI,32pts/ΔI和64pts/ΔI.可以從圖2 中清晰地觀察到,無論是爆轟波傳播進冷區或熱區,在網格分辨率為32pts/ΔI時均收斂,因此本文一維算例網格分辨率均為32pts/ΔI.之后對于二維ZND 爆轟波向胞格爆轟波的轉變過程進行了測試,結果如圖3 所示,網格分辨率分別為16pts/ΔI,32pts/ΔI和64pts/ΔI.對比可得,隨著網格分辨率的增加,可以捕捉到更為清晰的爆轟波胞格結構,與分辨率為64pts/ΔI結果相比較,分辨率32pts/ΔI下的計算結果并未有太大不同,考慮到計算成本,本文針對于二維爆轟波的數值模擬,網格分辨率采用32pts/ΔI.

圖2 一維爆轟波網格分辨率測試Fig.2 Numerical resolution test for one-dimensional detonations

圖3 二維胞格爆轟波網格分辨率測試(無擾動,kR=2.0)Fig.3 Numerical resolution test for transition of ZND to cellular detonations (without disturbance,kR=2.0)

2 結果與討論

為了研究爆轟波在非均勻介質中的傳播問題,做了以下3 方面的工作.首先,通過對平面CJ 爆轟波進入“熱”或“冷”區轉變過程進行簡單的一維數值模擬研究,揭示了爆轟波在不同類型擾動下的傳播機理,特別是在 “冷”區內的失效極限和重新起爆過程進行了分析.其次,將研究拓展到二維,揭示了爆轟波陣面結構在存在間斷溫度擾動時的結構和演化過程,在這一部分中將溫度擾動限定為一個周期.最后,將擾動設定為多個周期,并且改變擾動的振幅,討論了爆轟波在溫度擾動下的傳播過程和波陣面結構.通過從簡單到復雜的研究思路,希望可以更清晰地展示爆轟波波陣面結構的變化,更容易解釋傳播機理.

2.1 一維爆轟波轉變研究

首先對于平面CJ 爆轟波進入“熱”或“冷”區進行了一維研究,討論了爆轟波的轉變過程和隨后的傳播過程.計算設定如圖1(a)所示,左側為一初始溫度為T0的CJ 爆轟波,右側基準溫度為T0,通過添加或減小一個溫度幅值變化A表征熱區和冷區.

具有不同穩定性系數(kR=2.0,kR=3.0,kR=4.0)的爆轟波的一維計算結果如圖4(a)～ 圖4(c)所示,最終的穩態爆轟波波速如圖4(d)所示.穩態傳播的CJ 爆轟波遇到新的溫度界面時,由于兩側的熱力學性質不同,會在界面處產生右行的透射波和左行的反射波.研究發現,平面CJ 爆轟波在遇到初始溫度較高的界面時,會產生一個壓力較低但速度較高的透射波向前傳播,如圖4(d)所示,與此同時,會產生一個反向傳播的反射波.而當平面CJ 爆轟波遇到初始溫度較低的界面時,會產生與熱界面完全相反的現象,即透射波的壓力增加,速度降低.

圖4 壓力歷史和速度曲線圖Fig.4 Pressure history curve

爆轟波完成轉變過程后,根據傳播界面溫度的不同,透射爆轟波會以不同的穩態速度和壓力傳播.另外,計算結果表明,熱界面溫度大大降低了爆轟波穩態傳播時的峰值壓力,略微加快了爆轟波傳播速度.與此相反的是,冷界面增大了爆轟波穩態傳播的峰值壓力,減慢了爆轟波的傳播速度.然而,隨著冷界面溫度的繼續下降,爆轟波會經歷一個熄爆再起爆的過程,即為臨界狀態,如圖4(c)中冷區溫度為280 K 時曲線所示.當溫度再次降低時,爆轟波會完全熄爆,如圖4(a)～ 圖4(c)中溫度為190 K 曲線所示.在一維計算中,冷界面的臨界溫度在240 K 左右,當低于此溫度時,爆轟波會熄爆.

對于不穩定爆轟波和穩定爆轟波,前者的化學反應率對溫度更為敏感,在遇到冷區擾動時,反應率下降更快,更容易衰減熄爆.圖4(d)中色塊表示透射過程中出現熄爆和重新起爆的范圍.可以看到,隨著冷區溫度的下降,不穩定爆轟波(kR=4.0)熄爆更早出現.但是對3 種爆轟波而言,完全熄爆的冷區溫度基本都在240 K 左右,沒有明顯的區別.這與實驗中的多維爆轟波的結論有區別,原因主要是一維爆轟波和多維爆轟波傳播機制的不同.

2.2 二維激波/爆轟波單周期轉變過程研究

在這一部分中,首先使用激波然后再使用爆轟波作為入射波傳入非均勻介質中,通過這樣的對比研究,可以更容易的對爆轟波的傳播機制進行研究.在二維情況下,非均勻區域由中心的熱區域和上下邊界處兩個冷區組成的階梯分布構成,溫度擾動設定為一個半周期.與本文的一維研究的類似,激波透射后發展到穩定態經歷了一個轉變過程.不同之處在于,由于界面兩側的初始溫度不同,導致流動速度不同,從而在界面處產生剪切層,因此在界面附近處可以觀察到Kevin-Helmholtz 不穩定性.由于本研究中使用了無擴散效應的反應歐拉方程,因此流動中顯示的漩渦并不是真正的物理黏性現象,而是由于數值黏性造成的結果.盡管如此,如Han 等[14]所分析的,目前所使用的方法在某種程度上依然能夠體現Navier-Stokes 方程的某些特性.

如圖5(a)～ 圖5(f)所示,隨著激波在橫向非均勻介質中的傳播,激波前沿演化成一個復雜的結構.在溫度較高的中心區域存在一輕微前凸的激波波陣面,而在上下溫度較低的冷區,可以發現由入射波、馬赫桿、反射激波和剪切層組成的三波結構.入射激波與中心彎曲的激波面相連,并被剪切層隔開.因此,存在兩個剪切層,即由溫度界面演化成的剪切層和三波結構造成的馬赫反射剪切層,而且這兩個剪切層在斜入射波后形成一個薄通道.激波面后的流場由三部分組成,第一部分為中心處前凸激波面后的高溫低壓膨脹區,第二部分為斜入射激波后兩剪切層之間的狹窄通道,最后一部分為邊界附近兩正入射激波后的高壓低溫區.

圖5(e)的局部放大圖中給出了波陣面的精細結構以及流線分布,設置流線時將坐標系建在了波陣面上.斜入射激波的存在使得流線向邊界處彎曲,導致中心區流場的膨脹以及邊界處流場的收縮.相鄰的兩對剪切層在激波面后一定距離處合二為一,在靠近上邊界的另一對剪切層沿與中心線方向相反的方向卷曲.隨著渦流尺寸的增大,兩剪切層在距離激波面后較遠距離處重合,并且在非均勻介質中形成了獨特的激波復合結構體.

圖5(g)給出了最大壓力歷史云圖(對應于爆轟波模擬中的數值胞格圖),可以看出兩側臨邊界處壓力大,中心壓力低,但是中心低壓區(初始高溫區)面積擴張,兩側高壓區(初始低溫區)面積收縮.

圖6 比較了不同溫度擾動幅值下的激波波陣面的結構特征及其壓力歷史云圖,其中擾動幅度A分別為15 K,35 K 和105 K.可以發現在溫度擾動幅度較小時并沒有典型的三波點結構生成,而是在中心區域產生了一個微凸的激波以及上下邊界兩個微凹的激波.這是因為較低的溫度幅度使得中心與上下邊界流動狀態相差不大,這也延緩了剪切層中渦的演化.從波陣面上的最大壓力歷史中可以發現不同擾動幅度下,馬赫桿的高度和微凹激波的高度發生變化,溫度擾動幅值越大,馬赫桿越短,幅值越小,馬赫桿越長.同時發現較大幅度的溫度擾動下,波后湍流產生的壓縮波使得中心區域出現較大的壓力震蕩.

圖6 不同溫度擾動幅值下激波波陣面紋影圖和最大壓力歷史Fig.6 Schlieren photography of the shock front at the same instant for cases with different temperature perturbations and pressure history

下面使用ZND 爆轟波作為入射波傳入存在單周期溫度擾動的非均勻介質中.與激波相比,因為透射和傳播過程涉及到化學反應,使得問題更為復雜.圖7 給出了擾動幅值A為105 K 時的爆轟波前沿的演化過程,其中kR=2.0.從圖中可以觀測到,在上下邊界冷區內,爆轟波波陣面出現一定程度的解耦現象,即誘導區長度不斷變長.此外,在解耦區出現不斷試圖再次起爆的熱點,但均以失敗告終.即使在中心區域處前導激波與反應面也稍有解耦.在未有化學反應的冷區也發現了三波結構.由于冷區存在周期性的熄爆與再起爆過程以及流場的湍流特性,產生了未反應的區域.而此區域的燃燒產生了不同強度的壓縮波.可以發現這些壓縮波的相互作用可能在冷區產生局部的熱點,即使此時冷區初始溫度為190 K.聯想到一維數值模擬中冷區溫度為190 K時,透射爆轟波最終熄爆,可以說明一維和二維的情況存在不同,這也表明了第二個維度上的梯度對于此問題的影響.

圖7 二維爆轟波波陣面演化圖,溫度擾動振幅A=105 K,kR=2.0Fig.7 Evolution of detonation front structures during the twodimensional with A=105 K and kR=2.0

圖8 給出了兩個時刻沿中心線和邊界處壓力與溫度的分布曲線.與激波的例子相比,爆轟波波陣面由于化學反應的影響,從未到達與激波波陣面相同的穩態,而是一個不穩定的演化過程.沿中心線和邊界處壓力與溫度歷史曲線表明,冷區波陣面前沿的壓力高于中心波陣面前沿的壓力.圖8(a)中壓力曲線的第二個峰值清晰地顯示了重新起爆(熱點)的位置(t=0.25 ms).觀察溫度曲線圖,可以發現在x=0.25～ 0.35 m 之間溫度出現高頻率的震蕩,表明在渦流中存在大量未反應完全的氣團.這是未完全反應可燃氣體沿著上下壁面通道被渦卷入波后的中部區域造成的.

圖8 不同時刻下沿中心線和邊界線壓力與溫度分布和流線圖,溫度擾動振幅A=105 K,kR=2.0Fig.8 Pressure,temperature distribution and streamlines along the centering and boundary lines at two instants with A=105 K,kR=2.0

而后研究了不同穩定性下的爆轟波演變過程,其中kR分別為2.0,3.0 和4.0,溫度幅值分別為15 K,35 K 和105 K,其紋影圖與數值胞格圖如圖9所示.從圖中可以清晰地觀測到均勻介質中典型的二維胞胳爆轟波結構,唯一的區別在于由于橫向溫度分布不均產生的剪切層.由于計算域的寬度僅為誘導區長度的30 倍左右,所以橫向只容許存在一個胞胳.另外,爆轟波的不穩定性僅影響了胞胳爆轟波的相位,并未顯著改變胞胳的大小.如圖9(a)所示,在溫度擾動幅值較小時(15 K),可以清楚的觀察到規則的胞格結構,溫度擾動并沒有對胞格形狀產生明顯的影響.如圖9(b)所示,在溫度擾動幅值增大到35 K 時,可以看出胞格的形狀發生了一定程度的改變,主要體現在長寬比的增加,這說明橫波的速度相對于前導波陣面的速度減小.而且在圖9(b)的數值煙膜中,爆轟波存在一定程度的解耦-失效,然后重新起爆過程,這體現了爆轟不穩定性的影響.如圖9(c)所示,在溫度擾動幅值增大到105 K 時,可以看出數值煙膜中不存在典型的胞格結構,這說明橫波發展受到了溫度擾動的抑制,波陣面上不斷變化的爆轟波三波結構基本消失,取而代之的是與激波例子中類似地穩態的波陣面結構,只有在中心處能夠看到非常弱的橫波產生的痕跡.同時也可以看到,由于邊界兩側存在熄爆和重新起爆不斷重復的過程,馬赫桿的長度也發生了周期性的變化.

圖9 不同穩定性下爆轟波波陣面紋影、溫度擾動幅度以及胞格結構Fig.9 Schlieren photography of the detonation front at the same instant for cases with different instability

2.3 非均勻介質中波陣面結構

當在二維管道中布置一維ZND 爆轟波作為初始條件時,均勻介質中平面爆轟波的形態在向前傳播的過程中能夠維持一段時間或者距離.但是數值擾動會誘發波陣面局部上小尺度的扭曲或者褶皺,并在爆轟不穩定性的作用下,波陣面的扭曲放大,進而在橫向上產生壓力、溫度和密度的梯度,導致弱橫波的產生.弱橫波放大受不穩定性因子的控制,最終會達到平衡態,即橫波的強度以及橫波的間距相對保持不變,而且該平衡態取決于可燃介質的熱力學狀態和化學動力學特性,這也說明橫波間距是一個特征尺寸,此時便形成了一個完全發展的自維持胞格爆轟波.圖10 給出了kR=2 時從ZND 爆轟波到胞格爆轟波的演化過程中的數值胞格、壓力歷史和傳播速度曲線.可以看出,波陣面上的擾動最早在200 ΔI出現,而后在700 ΔI處達到平衡.數值模擬結果發現,隨著kR的增大,這兩個數值均減小,顯然這兩個距離受到爆轟波不穩定性的影響.

圖11～ 圖13 給出了不同擾動波長s情況下ZND爆轟波演化發展過程的數值胞格.圖11 對應的壓力歷史曲線(取自中心線)如圖14 所示.這一部分的研究中使用了小振幅的溫度擾動,即A=15 K.從圖11 中可以看出,引入溫度擾動后,只經過了幾個誘導區的距離波陣面上便出現了扭曲變形,并形成非常弱且非常規則的橫波,這可以從圖中規則的小胞格中看出,也可以從圖14 中小振幅震蕩的規則壓力曲線中反應出來.圖11 中規則小胞格的尺寸幾乎等于擾動的波長10ΔI.但是,經過了一段傳播距離達到位置x=400ΔI時,橫向不穩定性開始放大增長,最終在x=700ΔI處達到平衡狀態,形成了充分發展的穩態胞格爆轟波.壓力歷史曲線的震蕩也開始不規則,并且震蕩增加,與無擾動時壓力的震蕩特性相一致(如圖10(b)所示).這說明,人工溫度擾動可以在一定的范圍內抑制胞格不穩定性的發展,但是并不能夠終止這一進程.這主要是因為人工擾動的波長(s=10ΔI)遠小于特征橫波間距(40ΔI～ 50ΔI),并且人工擾動的振幅也較小.

圖10 均勻介質中無擾動時ZND 爆轟波到胞格爆轟波的演化Fig.10 Transition from ZND to cellular detonation without artificial perturbations

圖11 存在人工溫度擾動時的數值胞格,A=15 K,s=10ΔIFig.11 Numerical smoked foils for cases with artificial perturbations of A=15 K,s=10ΔI

圖12 存在人工溫度擾動時的數值胞格,A=15 K,s=20ΔIFig.12 Numerical smoked foils for cases with artificial perturbations of A=15 K,s=20ΔI

圖13 存在人工溫度擾動時的數值胞格,A=15 K,s=40ΔIFig.13 Numerical smoked foils for cases with artificial perturbations of A=15 K,s=40ΔI

如圖12 所示,當人工溫度擾動的波長增大至s=20 ΔI時,ZND 爆轟波的發展過程與s=10ΔI時相似,主要的差別是橫波不穩定性的放大發展被滯后,這說明,波長的增大能夠加強對胞格不穩定性的抑制.壓力歷史曲線隨著擾動波長的增加,震蕩的范圍增大,并且頻率與波長匹配.當人工溫度擾動的波長增大至s=40ΔI時,由于該波長與均勻介質中特征橫波間距接近,橫波不穩定性與人工溫度擾動立刻發生同步,結果如圖13 所示.只是在接近計算域終點附近時,胞格形態出現了一定程度的不規則性,此現象在具有更大不穩定性的例子中呈現的尤為明顯.但胞格的尺寸幾乎與人工溫度擾動的波長保持一致,只是不穩定性參數的增大加重了胞格尺寸或者胞格長度的散布,從壓力歷史曲線中也可以觀察到這一點.

從上述的分析中也可以發現,隨著不穩定性參數的增大,上述ZND 爆轟波的發展演化過程可被大大提前.因此總結來說,ZND 爆轟波在人工溫度擾動下的響應(發展演化過程)主要受制于兩個競爭性的因素:一是爆轟波內在的不穩定性(主要受熱力學和化學動力學特性控制),二是人工擾動的波長,前者是內因,后者是外因.但是,由于上述算例中,人工溫度擾動均設置了相同的、相對較小的振幅(15 K),即便存在溫度擾動,最終也能夠得到一個充分發展的胞格爆轟波,并且這些受擾動胞格爆轟波的主要特征與在均勻介質中的未受擾動胞格爆轟波基本一致.

由于人工擾動的存在,壓力震蕩的發散程度減小,但是頻率增大.人工溫度擾動可以看作是橫向周期性排列的熱點和冷點區域.熱點區域能夠通過增大反應速率而增強爆轟波的不穩定性,其中冷點通過減小反應速度而限制爆轟波的發展.但是冷點的負面效應要強于熱點的正面效應,這是由于阿倫尼烏斯反應率與溫度不是線性關系.因此人工溫度擾動的存在,一方面能夠限制波陣面的發展,另一方面也能夠限制爆轟波橫波的發展,使得波陣面更加的離散,壓力震蕩的波動幅度減小.

上述的研究中,重點關注較小溫度擾動幅度下ZND 爆轟波的發展演化過程.下面研究在較大溫度擾動時,波陣面的結構的變化.圖15 中給出了4 種不同溫度擾動時ZND 爆轟轉變為胞格爆轟波時的數值胞格結構,管道中線處的壓力歷史曲線見圖16.在x=1.2 m 處波陣面結構及波后的壓力分布如圖17所示.從數值胞格中可以看出,隨著溫度擾動幅度的增加,胞格的形態發生了很大的變化,當A=15 K時,溫度擾動的影響非常有限,胞格的形狀與無擾動時基本一致.圖17(a)中能夠看到擾動使得波陣面結構局部上發生細微的影響,波陣面存在小的間斷,但是并沒有對波陣面整體的結構產生大的影響,橫波依然很清晰,壓力曲線的震蕩也在正常的范圍內.但在A=35 K 和A=55 K 時,數值胞格的形態發生了較大的變化,三波的軌跡局部上發生了較大的扭曲,這說明橫波在穿越溫度擾動時的速度發生了變化,產生了震蕩,使得胞格形狀發生了一定的變化.從圖17(b)和圖17(c)中可以看出,溫度間斷使得波陣面更加的離散化,除了存在很長的自然橫波之外,還存在局部上溫度間斷造成的短的三波結構,使得波陣面表現出更多的湍流特性.在溫度振幅A=105 K的例子中,非常強的溫度擾動幾乎完全抑制了波陣面的演化,只能觀察到強橫波消失的湍流結構的波陣面,這也可以從壓力歷史曲線圖17(d)中看出,從波后橫向的壓力分布中看出橫波強度得到了大大地削弱.上述討論說明適度的溫度擾動能夠增加爆轟波自身的不穩定性,促進爆轟波的傳播,但是過大的溫度擾動則會通過抑制橫波的發展而限制爆轟波的傳播.除此之外,由于波陣面更多的湍流特性,產生較多小的未反應氣團.

圖15 存在不同振幅人工溫度擾動時(kR=4.0,s=20ΔI)的數值胞格Fig.15 Numerical smoked foils for cases with artificial perturbations of kR=4,s=20ΔI

圖16 存在不同振幅人工溫度擾動時(kR=4.0,s=20ΔI)的壓力歷史曲線圖Fig.16 Pressure history for cases with artificial perturbations of kR=4,s=20ΔI

圖17 人工溫度擾動s=20ΔI,kR=4.0 時,波陣面的紋影結構Fig.17 Schlieren photography of detonation front under artificial perturbations with different amplitude with s=20ΔI and kR =4.0

3 結論

使用反應歐拉方程和兩步反應模型對氣相爆轟波在非均勻介質中的發展演化過程進行了系統的數值模擬研究,重點分析了不同波長、幅度的溫度擾動對波陣面波系結構的影響,特別是對自維持傳播起到關鍵作用的橫波強度和間距的影響.為更深入地分析問題,本研究遵循從簡到繁的思路,先后進行了一維和二維數值研究,同時計算了激波算例,進行了對比研究.

一維研究表明,平面CJ 爆轟波在遇到初始溫度較高的界面時,會產生一壓力較低但速度較高的透射波向前傳播,同時會產生一個反向傳播的反射波.而當平面CJ 爆轟波遇到初始溫度較低的界面時,會產生與熱界面完全相反的現象,即透射波壓力增加,速度降低.當溫度繼續降低時,根據自身不穩定性的區別,爆轟波透射后會熄爆-重新起爆或完全熄爆.在本研究中,二維數值模擬使用平面ZND 爆轟波而不是充分發展的胞格爆轟波作為初始條件,這樣做的目的是避免在一開始就引入橫波不穩定性,使得問題過于復雜.當ZND 爆轟波作為入射波傳入橫向存在單周期溫度擾動的非均勻介質中時,由于在冷區和熱區爆轟波傳播的速度不同,溫度界面附近會產生剪切層,導致不同區波陣面的演化存在差異,冷區甚至存在熄爆和重新起爆的現象,使得整體上波陣面的結構與在均勻介質中有所不同.在溫度擾動幅值較小時(15 K),可以清楚地觀察到規則的胞格結構,且溫度擾動并沒有對胞格形狀產生明顯的影響.但當溫度擾動幅值增大到35 K 時,胞格的形狀發生了一定程度的改變,主要體現在胞格長寬比的增加.而當溫度擾動幅值增大到105 K 時,數值煙膜中不存在典型的胞格結構,這說明橫波發展受到了溫度擾動較強的影響,波陣面上不斷變化的爆轟波三波結構,只有在中心處能夠看到非常弱的橫波產生的痕跡.

對于不存在人為溫度擾動的均勻介質中,平面ZND 爆轟波可以在傳播一定距離后逐步演變為充分發展的胞格爆轟波,這是由于數值截斷誤差會誘發波陣面的局部扭曲,并可在內在不穩定性的控制下發展出橫波和波陣面復雜的三波結構.在ZND 爆轟波如此的發展模式下,若一開始便考慮溫度擾動,在早期階段,溫度擾動可通過抑制胞格不穩定性的放大來控制橫波間隔,并且維持一段時間的穩態發展.但隨后,胞格不穩定性再次開始放大,并最終轉變為完全發展的胞格爆轟波.這表明人為溫度擾動的存在通過抑制橫波的發展,延遲了ZND 爆轟波向胞格爆轟波的演化,并且內在不穩定性的增加可以減緩此延遲現象.另外,人工的溫度擾動可以在一定的范圍內抑制胞格不穩定性的發展,但是并不能夠終止這一過程.ZND 爆轟波在人工溫度擾動下的響應主要受制于兩個競爭性的因素:一是爆轟波內在的不穩定性,二是人工擾動的波長,前者是內因,后者是外因.但是,由于上述算例中,人工溫度擾動均設置了相對較小的振幅(15 K),即便存在溫度擾動,最終也能夠得到一個充分發展的胞格爆轟波,并且這些受擾動胞格爆轟波的主要特征與在均勻介質中的未受擾動胞格爆轟波基本一致.同時還發現,如果人工擾動的波長接近均質介質中的胞格爆轟波的橫波間隔,則這兩個因素發生同步,即立即形成具有極其規則胞格圖案的胞格爆轟波,僅在經過很長的距離后,相位變化才使胞格模式變得稍微不規則.

非均勻介質中溫度的不連續性使得波陣面波更加離散,表現出更多的湍流特性,除了自然的橫波外,還存在局部弱的三波結構.在溫度擾動幅度較小時,會增加胞格不穩定性,從而改變爆轟波陣面的傳播機制.相反,較大的人工擾動會抑制橫波的發展,此時僅能觀察到強橫波消失的湍流結構的波陣面,從而降低爆轟波的不穩定性.綜上,人工溫度擾動的特性和爆轟波內在不穩定性的競爭支配著爆轟前沿的演化過程和傳播機制.