影響氣相爆轟波穩(wěn)定性邊界的通用兩步反應(yīng)模型參數(shù)分析

郗雪辰， 王新文， 王克劍， 趙學(xué)民

(1.山西警察學(xué)院治安系， 山西太原 030401; 2.北京理工大學(xué)宇航學(xué)院， 北京 100081)

0 引言

對涉爆案件進行現(xiàn)場勘查，是公安機關(guān)刑事偵查工作的重要內(nèi)容之一，利用計算機技術(shù)復(fù)原爆炸現(xiàn)場，或?qū)ΜF(xiàn)場進行分析研判，可以為涉爆案件的偵破提供有效的幫助；從安全防范的角度來看，對可能發(fā)生的爆炸進行風(fēng)險評估，有助于相關(guān)的規(guī)范標準和應(yīng)急預(yù)案的制定。為了能夠?qū)崿F(xiàn)計算機技術(shù)在公共安全和公安涉爆領(lǐng)域的應(yīng)用，建立一套可靠的爆炸預(yù)測模型就具有重要的實踐意義。爆炸發(fā)生時，燃料及爆炸物劇烈的化學(xué)反應(yīng)，以燃燒波的形式由點火源開始向外傳播。通常根據(jù)傳播速度將這種自持傳播的燃燒波分為爆燃和爆轟。爆轟燃燒是化學(xué)反應(yīng)和氣體流動過程的耦合，這種極快速的能量釋放，表現(xiàn)出很強的破壞性。因此，對爆轟現(xiàn)象及其傳播規(guī)律進行研究，也是相關(guān)領(lǐng)域的重要課題之一。

早期的研究者提出CJ(Chapman-Jouguet)理論來定量預(yù)測爆轟速度，之后又建立了ZND(Zeldovich-von Neumann-Doring)模型描述爆轟波結(jié)構(gòu),模型假設(shè)爆轟波有一個穩(wěn)定的一維結(jié)構(gòu)，這個結(jié)構(gòu)由前導(dǎo)激波和其后的化學(xué)反應(yīng)區(qū)構(gòu)成[1]。自持傳播的一維ZND爆轟是不穩(wěn)定的，但是通常很難在實驗中觀察到這一現(xiàn)象[2]。隨著計算機技術(shù)的發(fā)展，可以以高精度和足夠的時空分辨率來對非線性帶反應(yīng)的歐拉方程進行積分運算，從而分析爆轟反應(yīng)區(qū)的精細結(jié)構(gòu)。對于一維的平面爆轟波，它顯示出與二維、三維數(shù)值模擬時相同的振蕩性質(zhì)[3]。這樣可以單獨研究一維縱向不穩(wěn)定性，然后推廣到更復(fù)雜的高維不穩(wěn)定性研究中。

一維爆轟波的不穩(wěn)定性僅表現(xiàn)為傳播方向的振蕩。常用的化學(xué)反應(yīng)模型主要有單步反應(yīng)模型和基元反應(yīng)模型。早期的數(shù)值模擬主要采用單步反應(yīng)模型[4-9]。單步反應(yīng)模型具有較高的計算效率，但是得到的反應(yīng)區(qū)是連續(xù)的。真實系統(tǒng)中的反應(yīng)區(qū)通常為一個明顯的誘導(dǎo)區(qū)和其后跟隨的快速放熱反應(yīng)區(qū)。Fickett[10]引入了一個簡單的代表一步鏈分支的反應(yīng)動力學(xué)模型研究不穩(wěn)定一維脈沖爆轟波，通過調(diào)整模型可獲得兩種可能的反應(yīng)區(qū)結(jié)構(gòu)。但是這種化學(xué)反應(yīng)模型很難確定控制不穩(wěn)定性的關(guān)鍵參數(shù)。而基元反應(yīng)模型雖然可以很好地模擬復(fù)雜的化學(xué)反應(yīng)過程，但是受限于現(xiàn)階段的計算能力和高昂的計算成本，并非研究此類問題的最佳選擇。為了兼顧計算效率和計算成本，Short[11]提出了一種三步反應(yīng)模型，對鏈分支反應(yīng)驅(qū)動的脈沖爆轟過程的控制參數(shù)進行了研究。該模型將化學(xué)反應(yīng)分為鏈初始反應(yīng)、鏈分支反應(yīng)和鏈終止反應(yīng)。鏈初始反應(yīng)的Arrhenius反應(yīng)需要很高的能量，生成低濃度的自由基。緊接著是低活化能的Arrhenius反應(yīng)，加速自由基的生成。最后是放熱的鏈終止反應(yīng)。此后，Short和Sharpe[12]將三步反應(yīng)模型簡化成了兩步反應(yīng)模型。該模型包含了一個由Arrhenius速率定律控制的熱中性鏈分支反應(yīng)誘導(dǎo)區(qū)。計算得到的解與三步模型的解定量上相同。Ng[13]在前人的基礎(chǔ)上，采用修正后的兩步鏈分支反應(yīng)模型，通過改變參數(shù)來覆蓋穩(wěn)定性邊界內(nèi)穩(wěn)定性參數(shù)χ的范圍，研究了一維爆轟的不穩(wěn)定性。

研究爆轟波的不穩(wěn)定性，有助于認識其傳播規(guī)律?？梢詾楸ǚ揽丶夹g(shù)中涉及的安全防護、隔爆抑爆等問題的研究提供理論依據(jù)。目前，基于兩步反應(yīng)模型的氣相爆轟波不穩(wěn)定性研究較少。本文針對一維爆轟波，采用無粘Euler方程和兩步誘導(dǎo)- 放熱反應(yīng)模型，固定參數(shù)放熱量Q和比熱比γ，分別調(diào)整誘導(dǎo)區(qū)活化能EI、反應(yīng)區(qū)活化能ER和反應(yīng)速率常數(shù)kR，計算無衰減單周期振蕩模態(tài)的臨界參數(shù)，通過尋找參數(shù)間的對應(yīng)關(guān)系，研究影響穩(wěn)定性邊界的兩步反應(yīng)模型參數(shù)。研究結(jié)果可為將來利用兩步反應(yīng)模型研究多維爆轟波不穩(wěn)定性的參數(shù)選擇提供參考。

1 模型介紹

1.1 控制方程

本文的計算忽略粘性項，控制方程采用無量綱形式的一維反應(yīng)Euler方程組：

(1)

(2)

(3)

假設(shè)氣體為具有固定比熱比的理想氣體，總能量和狀態(tài)方程為：

(4)

(5)

上述參數(shù)ρ、u、p、e、q、T分別表示流體的密度、速度、壓力、總能量、化學(xué)反應(yīng)放熱量和氣體溫度。化學(xué)反應(yīng)放熱量q的表達式為：

q=λQ

(6)

其中λ表示放熱反應(yīng)進程變量。上述所有變量都用波前未燃氣體狀態(tài)參數(shù)進行無量綱化：

(7)

其中下角標0表示波前的氣體參數(shù)。上述方程與化學(xué)動力學(xué)模型耦合，用來描述爆轟波結(jié)構(gòu)。采用Ng[13]在Short[12]研究基礎(chǔ)上修正的兩步鏈分支反應(yīng)模型，用兩個化學(xué)速率控制方程模擬化學(xué)反應(yīng)過程。第一步是熱中性誘導(dǎo)區(qū)或點火過程，反應(yīng)速率是對溫度敏感的Arrhenius形式：

(8)

其中，ξ表示誘導(dǎo)區(qū)反應(yīng)過程的變量；kI表示誘導(dǎo)區(qū)化學(xué)反應(yīng)速率常數(shù)，kI=-Uvn，Uvn是CJ爆轟波的激波前沿速度；EI表示誘導(dǎo)區(qū)活化能；Ts表示氣體經(jīng)過前導(dǎo)激波之后的溫度；H(1-ξ)是階躍函數(shù)：

(9)

第二步是緩慢放熱之后的快速熱量釋放過程，其反應(yīng)速率方程為：

=[1-H(1-ξ)]ρkR(1-λ)exp(-ER/T)

(10)

其中λ表示鏈重組反應(yīng)過程變量，kR表示放熱區(qū)反應(yīng)速率常數(shù)，ER表示放熱區(qū)活化能。上述公式中的EI和ER都采用RT0進行無量綱化，Ts采用T0進行無量綱化。

1.2 計算模型及參數(shù)

本文采用的是色散可控耗散格式(DCD)[14]捕捉激波間斷，這是一種基于修正方程的控制而構(gòu)造的精度高、無振蕩的差分格式，屬于總變差衰減格式(TVD)中的一種，對流項的差分格式具有空間二階精度，時間推進采用三階Runge-Kutta方法。計算域網(wǎng)格尺寸為0.05×0.05，網(wǎng)格數(shù)量為40 000×3。所用兩步反應(yīng)模型的主要計算參數(shù)包括：放熱量Q，比熱比γ，誘導(dǎo)區(qū)活化能EI，放熱區(qū)活化能ER，放熱反應(yīng)速率常數(shù)kR。其中Q=50，γ=1.2，另外3個參數(shù)為可調(diào)參數(shù)。利用不同的壓力振蕩模態(tài)分析爆轟波的不穩(wěn)定性，獲得兩步反應(yīng)各參數(shù)對應(yīng)的穩(wěn)定性邊界。

2 結(jié)果與討論

2.1 兩步反應(yīng)模型參數(shù)與穩(wěn)定性的關(guān)系

研究兩步反應(yīng)模型參數(shù)與穩(wěn)定性的關(guān)系，其結(jié)果可用于定義穩(wěn)定性邊界兩側(cè)的參數(shù)特性，有助于將來進行涉爆案件分析研判、防護裝備設(shè)計等研究時，快速選取合理的參數(shù)區(qū)間。在單步和兩步反應(yīng)模型中，EI與穩(wěn)定性的關(guān)系是相同的，且單步反應(yīng)只有一個控制參數(shù)，便于討論，因此采用單步化學(xué)反應(yīng)模型計算一維爆轟波的傳播。

2.1.1 誘導(dǎo)區(qū)活化能EI與穩(wěn)定性的關(guān)系

通過調(diào)整活化能EI獲得了幾種不同的振蕩模態(tài)，找到了無衰減單周期振蕩模態(tài)所對應(yīng)的臨界活化能，并將這種振蕩模態(tài)定義為穩(wěn)定性邊界；同時，獲得的臨界還可為之后的計算參數(shù)選擇提供參考。對于衰減的單周期振蕩，被認為是小于該邊界的，可判定為一維爆轟傳播穩(wěn)定。對于增長的單周期振蕩，以及出現(xiàn)倍周期分叉后的多周期振蕩模態(tài)，則認為是大于邊界的，判定為爆轟傳播不穩(wěn)定。

對于單步化學(xué)反應(yīng)Arrhenius速率模型來說，由活化能EI決定的化學(xué)反應(yīng)的溫度敏感性是最重要的穩(wěn)定性參數(shù)。單步反應(yīng)的Arrhenius方程為：

(11)

其中，A、EI、R和T分別表示指前因子、活化能、混合氣體的氣體常數(shù)和混合氣體溫度?；罨芎椭盖耙蜃拥膶?yīng)關(guān)系如表1所示：

表1 活化能EI和指前因子A的對應(yīng)關(guān)系

圖1 EI分別為25.00，25.25，25.26和27.00的振蕩模態(tài)

不同活化能的爆轟波振蕩模態(tài)如圖1所示。圖1(a)顯示的是EI=25.00時的振蕩模態(tài)，此時呈現(xiàn)為快速衰減的單周期性振蕩，而且從時間坐標看，后周期的振幅明顯小于其前周期的振幅，這種衰減的單周期波形說明爆轟的傳播是趨于穩(wěn)定的。以前的研究表明，隨著活化能EI增大，爆轟的傳播會向不穩(wěn)定狀態(tài)發(fā)展，呈現(xiàn)出不同的振蕩模態(tài)。為了尋找Q=50，γ=1.2條件下單步反應(yīng)的穩(wěn)定性邊界，通過緩慢增加活化能，引導(dǎo)爆轟向不穩(wěn)定狀態(tài)發(fā)展。發(fā)現(xiàn)當(dāng)EI=25.25時，結(jié)果表現(xiàn)為臨界穩(wěn)定的單周期振蕩模態(tài)，即達到了單步反應(yīng)的穩(wěn)定性邊界，如圖1(b)所示。當(dāng)EI=25.26時就成了無衰減振蕩，如圖1(c)所示，這表明爆轟波已經(jīng)開始不穩(wěn)定。繼續(xù)增加活化能，發(fā)現(xiàn)初始周期振幅幾乎不變，周期振幅會隨著時間越來越大，使整個波形表現(xiàn)為擴張的單周期振蕩。圖1(d)為EI=27.00時刻的圖形，此時的爆轟波傳播已經(jīng)出現(xiàn)了向雙周期振蕩轉(zhuǎn)變的現(xiàn)象，且振幅顯著增加。當(dāng)活化能增加到一定值后，單周期振蕩會經(jīng)過倍周期分叉，呈現(xiàn)出雙周期、四周期等多周期振蕩模態(tài)，最終變成混沌狀態(tài)。顯然，高活化能情況下的爆轟是不穩(wěn)定的，這是由于反應(yīng)速率對溫度的敏感性強，較小的溫度擾動就會導(dǎo)致反應(yīng)速率的大幅波動。

2.1.2 放熱區(qū)活化能ER與穩(wěn)定性的關(guān)系

與單步化學(xué)反應(yīng)模型不同，兩步誘導(dǎo)- 放熱反應(yīng)模型的振蕩模態(tài)不僅和EI有關(guān)，還受到放熱區(qū)活化能ER和反應(yīng)速率常數(shù)kR的影響。兩步反應(yīng)中的誘導(dǎo)區(qū)活化能，采用單步反應(yīng)模型求解得到的無衰減單周期振蕩模態(tài)對應(yīng)的臨界活化能EI=25.25。

圖2 不同ER的振蕩模態(tài)

對于給定的反應(yīng)速率常數(shù)，通過調(diào)整活化能ER，得到的爆轟波面的壓力振蕩歷程如圖2所示。可以看到，爆轟波在最開始會短暫經(jīng)歷計算帶來的數(shù)值振蕩，之后則表現(xiàn)出良好的周期性。圖中穩(wěn)定振蕩曲線表示ER取1.0、2.0和3.0時穩(wěn)定性邊界對應(yīng)的波面振蕩模態(tài)，此時kR分別為1.32、1.46和1.61。在穩(wěn)定kR的情況下，增大或減小ER的取值，可以知道兩步反應(yīng)中放熱區(qū)活化能與爆轟波不穩(wěn)定性的對應(yīng)關(guān)系。由結(jié)果可知，在小ER的情況下，波面的壓力振蕩表現(xiàn)出振幅逐漸增加的單周期振蕩模態(tài)，且擁有更大的振幅和更短的波長，這表明一維爆轟波已經(jīng)跨過穩(wěn)定性邊界向不穩(wěn)定狀態(tài)模態(tài)發(fā)展。增大ER可促進爆轟波向穩(wěn)定發(fā)展，這一結(jié)果與Ng的理論相符。

2.1.3 反應(yīng)速率常數(shù)kR與穩(wěn)定性的關(guān)系

給定放熱區(qū)活化能ER，調(diào)整反應(yīng)速率常數(shù)kR，考察其對振蕩模態(tài)的影響，并尋找各參數(shù)條件下穩(wěn)定性邊界對應(yīng)的臨界反應(yīng)速率常數(shù)。圖3顯示了放熱區(qū)活化能分別為1.0、2.0和3.0時，幾種kR參數(shù)的振蕩模態(tài)，所取ER覆蓋了常見燃料的參數(shù)區(qū)間。由計算結(jié)果可以看出，增大kR對爆轟波的失穩(wěn)過程起到一定的促進作用，誘導(dǎo)振蕩模態(tài)從單周期向多周期過渡?；罨芊謩e為1.0、2.0和3.0時，穩(wěn)定性邊界對應(yīng)的kR分別是：1.32、1.46和1.61。從圖中可以看出，穩(wěn)定后的振幅未發(fā)生顯著變化，但是發(fā)現(xiàn)波面振蕩的壓力區(qū)間會隨著kR的增大緩慢提高。隨著ER的增加，一維爆轟趨于穩(wěn)定，因此穩(wěn)定性邊界對應(yīng)的kR也在增大。與活化能不同，反應(yīng)速率常數(shù)kR是兩步反應(yīng)為了便于模型化研究而加入的控制參數(shù)。因此該參數(shù)并不直接影響爆轟波的不穩(wěn)定性，只是在一定程度上控制振蕩收斂的速率。

圖3 不同kR的振蕩模態(tài)

2.2 兩步反應(yīng)參數(shù)對穩(wěn)定性邊界的影響

圖4 不同EI對應(yīng)的穩(wěn)定性邊界

之前討論了兩步反應(yīng)參數(shù)與穩(wěn)定性邊界的關(guān)系，同時獲得了穩(wěn)定性邊界對應(yīng)的兩步反應(yīng)參數(shù)。對幾種EI條件下的放熱區(qū)參數(shù)進行了統(tǒng)計，得到的穩(wěn)定性邊界線如圖4所示。該邊界線表示了EI分別為15.25、20.25和25.25 3個參數(shù)條件下，反應(yīng)區(qū)參數(shù)在臨近無衰減單周期振蕩時的取值變化?？梢钥闯?，在穩(wěn)定性邊界上，活化能和反應(yīng)速率常數(shù)呈正相關(guān)，這種關(guān)系與EI的取值無關(guān)。在誘導(dǎo)區(qū)活化能EI取值在15.25～20.25的區(qū)間內(nèi)，穩(wěn)定性邊界對kR的變化更加敏感，也就是說在此區(qū)間內(nèi)kR對穩(wěn)定性的影響占主導(dǎo)。隨著EI增大，取值在20.25～25.25區(qū)間內(nèi)，穩(wěn)定性邊界對kR的敏感程度大大降低，充分說明低活化能情況下爆轟非常穩(wěn)定。低活化能對應(yīng)的穩(wěn)定性邊界擁有高反應(yīng)速率，反應(yīng)速率加快對爆轟波向不穩(wěn)定發(fā)展有促進作用。因此在穩(wěn)定性邊界上，爆轟波在傳播方向上的壓力振蕩非常劇烈，如圖5所示。

圖5 穩(wěn)定性邊界的振蕩模態(tài)

以單步反應(yīng)穩(wěn)定性邊界對應(yīng)的臨界活化能為固定參數(shù)，通過調(diào)整兩步反應(yīng)中的放熱區(qū)活化能和反應(yīng)速率常數(shù)，獲得了臨近無衰減單周期振蕩模態(tài)的參數(shù)對應(yīng)關(guān)系。在此基礎(chǔ)上，又分別計算了誘導(dǎo)區(qū)活化能15.25和20.25的情況，并獲得了這兩個活化能條件下的參數(shù)對應(yīng)關(guān)系，統(tǒng)計數(shù)據(jù)如表2所示。

表2 穩(wěn)定性邊界對應(yīng)的模型參數(shù)

圖6 一維爆轟波穩(wěn)定性邊界

為了得到更加直觀的穩(wěn)定性邊界，利用Matlab對上述數(shù)據(jù)進行插值計算，繪制了由兩步反應(yīng)關(guān)鍵參數(shù)構(gòu)造一維爆轟波穩(wěn)定性邊界面，如圖6所示。位于邊界面以上的參數(shù)組合，得到的爆轟波是不穩(wěn)定的；位于邊界面以下的參數(shù)組合，爆轟波是穩(wěn)定的。從圖6可以明顯看出，活化能增大，穩(wěn)定性區(qū)間就越小。

3 結(jié)論

本文采用兩步誘導(dǎo)- 放熱反應(yīng)模型，對不同參數(shù)條件下的一維爆轟波進行了數(shù)值模擬。計算結(jié)果顯示，增大EI和kR，或減小ER，會促進一維爆轟波跨過穩(wěn)定性邊界，向不穩(wěn)定狀態(tài)發(fā)展；對于給定的EI，穩(wěn)定性邊界對應(yīng)的ER和kR呈正相關(guān)，即在穩(wěn)定性邊界上，ER越大，kR也越大；在合理的參數(shù)范圍內(nèi)，EI越大，ER和kR滿足穩(wěn)定一維爆轟波的取值范圍越小；EI在15.25～20.25區(qū)間內(nèi)，穩(wěn)定性邊界對kR的變化比較敏感；EI在20.25～25.25區(qū)間內(nèi)，穩(wěn)定性邊界對kR的敏感程度大大降低。

通過分析各參數(shù)與爆轟波面壓力振蕩模態(tài)的關(guān)系，尋找無衰減單周期振蕩模態(tài)對應(yīng)的臨界參數(shù)，利用兩步反應(yīng)模型中的幾個關(guān)鍵參數(shù)，構(gòu)造了一維爆轟波的穩(wěn)定性邊界。研究獲得的穩(wěn)定性邊界，可以作為利用兩步反應(yīng)模型研究多維爆轟波穩(wěn)定性問題的參數(shù)選擇依據(jù)，屬于建立爆炸預(yù)測模型的基礎(chǔ)性工作。建立的爆炸預(yù)測模型，將來可以應(yīng)用于爆炸現(xiàn)場分析及還原、爆炸物銷毀現(xiàn)場方案制定、風(fēng)險評估、搜排爆防護裝備的研發(fā)改進，以及隔爆抑爆技術(shù)研究等領(lǐng)域。