郗雪辰， 王新文， 王克劍， 趙學民

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

0 引言

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

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

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

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

1 模型介紹

1.1 控制方程

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

(1)

(2)

(3)

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

(4)

(5)

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

q=λQ

(6)

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

(7)

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

(8)

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

(9)

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

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

(10)

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

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

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

2 結果與討論

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

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

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

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

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

(11)

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

表1 活化能EI和指前因子A的對應關系

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結論

本文采用兩步誘導- 放熱反應模型，對不同參數(shù)條件下的一維爆轟波進行了數(shù)值模擬。計算結果顯示，增大EI和kR，或減小ER，會促進一維爆轟波跨過穩(wěn)定性邊界，向不穩(wěn)定狀態(tài)發(fā)展；對于給定的EI，穩(wěn)定性邊界對應的ER和kR呈正相關，即在穩(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)的關系，尋找無衰減單周期振蕩模態(tài)對應的臨界參數(shù)，利用兩步反應模型中的幾個關鍵參數(shù)，構造了一維爆轟波的穩(wěn)定性邊界。研究獲得的穩(wěn)定性邊界，可以作為利用兩步反應模型研究多維爆轟波穩(wěn)定性問題的參數(shù)選擇依據(jù)，屬于建立爆炸預測模型的基礎性工作。建立的爆炸預測模型，將來可以應用于爆炸現(xiàn)場分析及還原、爆炸物銷毀現(xiàn)場方案制定、風險評估、搜排爆防護裝備的研發(fā)改進，以及隔爆抑爆技術研究等領域。