爆轟胞格尺寸的統(tǒng)計分析

崔 皓，歸明月，張 振，張 輝

(南京理工大學 瞬態(tài)物理國家重點實驗室，南京 210094)

0 引 言

爆轟波在傳播時，其波陣面呈入射激波、反射激波和馬赫干組成的三波結構，該結構是非定常的，三波點的軌跡構成了爆轟的胞格結構。胞格的尺寸作為爆轟本征值，反映了可燃系統(tǒng)氣體動力學和化學反應動力學之間的非線性反饋作用，可以用來表征爆轟極限、爆轟臨界直徑和臨界能量等爆轟動力學參數(shù)[1-5]。然而，胞格尺寸的測量主觀性比較強。對于一個不規(guī)則的爆轟胞格，不同研究者測得的胞格尺寸可能相差較大。圖1是C3H8+5O2在直徑為52 mm的圓管內的實驗胞格，8個實驗人員對其進行了測量，得到的胞格尺寸分別為8、13.5、10、6.4、13、7.1、10、8.3 mm，最大值和最小值幾乎相差了一倍[6]。這種研究者的主觀因素帶來的誤差會明顯影響到可燃系統(tǒng)的爆轟動力學參數(shù)的確定。因此，需要對爆轟胞格的尺寸進行客觀的測量，以消除或減弱這種主觀因素引起的誤差，實現(xiàn)對爆轟過程的控制。

圖1 C3H8+5O2在直徑52 mm圓管的爆轟胞格結構Fig. 1 Detonation cell structure of C3H8+5O2 in a circular tube with a diameter of 52 mm

Shepherd[7]最早通過計算數(shù)字化胞格的能量譜，得到空間尺寸分布的頻譜，根據(jù)譜的峰值來確定相應的尺寸為胞格尺寸。隨后，Shepherd[8]采用該方法對不同燃料形成的爆轟的胞格尺寸進行了分析，比如氫氣、乙炔、乙烯等。這種方法的缺點是對噪點很敏感，在對煙跡板上的胞格進行數(shù)字化處理時，煙跡板上的積碳會導致背景灰度值波動較大，進而影響該方法對胞格尺寸的預測。Lee[9]采用自相關函數(shù)法對實驗胞格進行了分析，該方法的物理意義更加清晰，它可以提供x方向和y方向偏移量的相關性。作者對不同濃度氬氣稀釋的乙炔與氧氣的混合物形成的爆轟胞格進行了討論。對于不規(guī)則的胞格，該方法面臨噪音干擾的問題。為了降低這種干擾的影響，Hebral[10]使用Photoshop和Matlab對圖像進行降噪銳化處理，使跡線更加明顯，提高胞格尺寸的統(tǒng)計效果。近來，Shepherd[6]采用概率密度函數(shù)法和自相關函數(shù)法對數(shù)值爆轟胞格進行了初步的統(tǒng)計分析，其中三波點軌跡的捕捉采用正渦量和負渦量，這樣避免了三波點統(tǒng)計時的干擾。

以上研究表明，爆轟胞格尺寸的測量，尤其是不規(guī)則胞格，需要采用統(tǒng)計處理來擺脫人眼的主觀性。本文采用兩種統(tǒng)計方法：概率密度函數(shù)法和自相關函數(shù)法，對規(guī)則程度不同的數(shù)值胞格進行了統(tǒng)計分析，并對兩種統(tǒng)計方法進行了比較。由于爆轟的穩(wěn)定性對活化能比較敏感，本研究通過調整活化能的大小來形成規(guī)則程度不同的爆轟胞格。

1 物理模型和算法

1.1 控制方程和胞格的捕捉

對于氣相爆轟波而言，國內外普遍采用Euler方程進行數(shù)值研究，本文研究基于單步反應的Euler方程：

其中，

式中，u和v分別為x和y方向的速度分量，t、p、ρ和Z分別表示時間、混合物的壓力、混合物的密度和反應產(chǎn)物的質量分數(shù)。e為混合物單位質量的總能量，滿足以下關系：

式中q為反應熱?；瘜W反應采用單步不可逆反應，反應速率滿足：

式中k、Ea和R分別表示指前因子、活化能和氣體常數(shù)，其中k的選取是為了保證特征長度為半個反應區(qū)長度，即L1/2=1，其值與活化能有關，計算公式為：

為了捕捉流場中爆轟波的精細結構，本文采用分裂格式，方程（1）中的對流項采用5階WENO格式(Weighted Essentially Non-Oscillatory)離散[11]，時間項采用3階TVD型的Runge-Kutta法[12]進行積分求解，化學反應源項采用五階Runge-Kutta法進行求解。以上所提到的計算格式已經(jīng)在相關的算例中經(jīng)過了驗證[13-15]。

爆轟胞格作為爆轟波陣面三波點的軌跡，其形成機理存在多種觀點，Quirk[16]和Sharpe[17]認為實驗中煙跡片上產(chǎn)生的胞格痕跡是由三波點處的最大渦量引起的，故本文在數(shù)值研究爆轟胞格時，采用記錄每個網(wǎng)格點上達到的最大渦量 ω=(?u/?y-?v/?x)來獲得數(shù)值爆轟胞格。同時利用渦量存在正渦和負渦的特性，代表三波點不同的傳播方向。當跟蹤網(wǎng)格點上的最大渦量（正渦量）時，可得到一組相互平行的向左傳播的跡線（面向下游）；當跟蹤網(wǎng)格點上的最小渦量（負渦量）時，可得到一組相互平行的向右傳播的跡線（面向下游）。將最大渦量和最小渦量疊加起來就是常規(guī)的魚鱗狀的胞格結構，這種處理的方法是為了便于文章后面采用的自相關函數(shù)對跡線周期的提取。

1.2 爆轟胞格尺寸的概率密度法

概率密度函數(shù)是用來描述隨機變量在給定的樣本空間中出現(xiàn)的可能性。概率密度函數(shù)范圍介于樣本空間中樣本值的最大值和最小值之間，其峰值所對應的變量值表示最可能發(fā)生的事件。

隨機變量x的累計分布函數(shù)F(x)，存在非負函數(shù)f(x)，使得對于任意實數(shù)有：

則稱x為連續(xù)隨機變量，其中f(x)稱為x的概率密度函數(shù)，簡稱概率密度，存在：

對于爆轟胞格，通過跟蹤三波點軌跡的渦量數(shù)據(jù)，計算任意相鄰兩條平行跡線間距。最終通過繪制相鄰兩條跡線間距的概率密度分布圖，其中最大概率對應的三波點的間距可以定義為胞格的尺寸。

1.3 爆轟胞格尺寸的自相關函數(shù)法

自相關函數(shù)主要應用在信號處理中，用來描述信號在一個時刻的取值與另一個時刻的取值的依賴關系。當信號中有周期性分量出現(xiàn)時，自相關函數(shù)分布的極大值能夠很好的體現(xiàn)這種周期長度。

對于爆轟波的二維流場中，以渦量ω（x,y）為待處理的信號，其自相關函數(shù)R（dx, dy）可以定義為：

其中，dx和dy分別表示x方向和y方向的信號的偏移。上述相關函數(shù)的定義表明，未偏移信號ω(x,y)與偏移信號ω(x+dx,y+dy)之間的關聯(lián)程度。

2 結果與討論

2.1 不同活化能的爆轟胞格結構

本文采用的物理模型如圖2所示，無量綱熱力學參數(shù)分別為：初溫T0= 1，初壓p0= 1，比熱比γ= 1.2，反應熱Q= 50。以半個反應區(qū)長度L1/2為特征長度，采用的網(wǎng)格標準為64/L1/2（半個反應區(qū)里有64個計算網(wǎng)格），該網(wǎng)格精度足以描述爆轟的胞格特性[3]。計算區(qū)域縱向長度(y方向)為300L1/2，流向尺度(x方向)是爆轟波傳播方向的尺度。在爆轟胞格的數(shù)值模擬中，需要足夠長的時間和空間才能得到穩(wěn)定的流場結構。為了節(jié)省計算資源，本文采用移動計算窗口(x方向尺度為1 800L1/2)，即爆轟波接近計算區(qū)域右端邊界時，在計算區(qū)域最左端舍棄一部分網(wǎng)格，同時在右端增加一部分網(wǎng)格，這部分新增加的網(wǎng)格中的參數(shù)設置為新鮮的未燃氣體。需要注意的是，這部分網(wǎng)格數(shù)量的選取要適當，如太多，左邊界會影響爆轟波的結構，如太少，會導致舍棄網(wǎng)格和增加網(wǎng)格的操作太頻繁，從而降低計算效率。經(jīng)過本文的計算實驗表明，這部分網(wǎng)格數(shù)量為500L1/2。

圖2 計算示意圖Fig. 2 Schematic diagram of the computational domain

為了表征爆轟胞格的規(guī)則程度，活化能選取了四種不同的情況，分別為Ea= 15、20、25和27，其中數(shù)值模擬的爆轟動力學參數(shù)已經(jīng)在Sharp[6]的研究中驗證過。圖3為活化能為15時的不同時刻流場的變化圖，左圖為壓力云圖，右圖為反應物質量分數(shù)云圖。初始時刻，在管道左端的高溫高壓點火，如圖3（a）所示，形成的壓力波往周圍傳播，在管道上下壁面反射，形成反射激波，往管道右端傳播，與此同時，激波后面的未燃氣迅速燃燒，轉變成燃燒產(chǎn)物，如圖3（b）所示。隨后，激波在往前傳播的過程中，形成的橫波往管道上下壁面之間來回反射，此時，波陣面只存在一組三波點（分別向上下傳播的2個三波點），如圖3（c）所示。經(jīng)歷一段時間后，波陣面開始出現(xiàn)多個子波，即多組三波點，如圖3（d）所示。隨著波陣面繼續(xù)往前傳播，經(jīng)過上下壁面的多次反射，形成充分發(fā)展的爆轟波（具有完整的爆轟波波系的結構），如圖3（e）所示。此時的爆轟波并不是最終的爆轟波，其胞格尺寸也不同于最終的胞格尺寸，該尺寸與數(shù)值計算的初始條件有關。由于爆轟波是一個非定常結構，必須要經(jīng)過足夠長的時間才能得到穩(wěn)定的爆轟波，如圖3（f）所示。

圖3 不同時刻爆轟波傳播的流場圖Fig. 3 Flow fields of the detonation wave at different time instances

圖4是當爆轟穩(wěn)定傳播時，根據(jù)網(wǎng)格點上的最大渦量繪制的三波點的跡線，由于渦量存在正負，圖4（a、b、c）分別是根據(jù)正渦量最大值ω+、負渦量最大值ω-、正負渦量疊加(ω++|ω-|)/2繪制的三波點軌跡。渦量的符號代表三波點不同的傳播方向，正渦量代表一組相互平行的向左傳播的跡線（面向下游），負渦量代表一組相互平行的向右傳播的跡線，正渦量與負渦量的疊加得到魚鱗狀的胞格結構。

圖4 活化能Ea = 15時的三波點軌跡Fig. 4 Triple-point tracks for Ea = 15

圖5是活化能Ea為20、25、27時，根據(jù)網(wǎng)格點上的正負渦量疊加(ω++|ω-|)/2繪制的數(shù)值胞格。隨著活化能的增加，激波后的誘導區(qū)變長，反應區(qū)變短，能量釋放越快，此時，誘導區(qū)的溫度擾動對反應的影響越大，相應的爆轟越不穩(wěn)定，形成的胞格也越不規(guī)則，流場達到穩(wěn)定所需的時間也越長，數(shù)值胞格所選的區(qū)域也越靠后。為了排除初始流場的影響，數(shù)值胞格的統(tǒng)計處理需采用穩(wěn)定的流場信息。因此，活化能為20、25和27時，選取的胞格區(qū)域為7 200 ≤x≤ 9 000。由圖5，活化能為20時，胞格結構開始出現(xiàn)不規(guī)則，三波點的跡線出現(xiàn)相交的趨勢，活化能增加到25時，三波點的跡線出現(xiàn)合并，胞格尺寸的取值范圍明顯增加，當活化能增加到27時，胞格極不規(guī)則，其尺寸的取值范圍進一步擴大。

圖5 不同活化能的數(shù)值爆轟胞格結構Fig. 5 Numerical detonation cell structures for different activation energies

2.2 爆轟胞格尺寸的統(tǒng)計分析

對于規(guī)則胞格（Ea= 15），采用兩種不同的統(tǒng)計方法對胞格尺寸進行處理。首先，采用概率密度函數(shù)法對規(guī)則胞格進行處理?；趫D4（a）繪制的三波點跡線，利用Sharpe提出的全局閾值法提取胞格跡線，即對渦量的數(shù)值設置一個閾值，本研究采用的閾值為27，小于該閾值的渦量舍去，僅保留高于該閾值的渦量，得到圖6所示的胞格跡線，此時的胞格較規(guī)則，三波點的跡線沒有出現(xiàn)相交的情形。對相鄰的胞格跡線間的距離進行概率統(tǒng)計，繪制概率曲線，如圖7所示，此時，最大概率值所對應的跡線的距離約為28，該值即為胞格的尺寸。

采用自相關函數(shù)法對規(guī)則胞格（Ea= 15）進行統(tǒng)計處理，圖8代表自相關函數(shù)在不同偏移方向的變化，其中x方向表示爆轟波傳播的方向，y方向表示爆轟波結構中橫波的傳播方向，θ方向是三波點跡線的方向。根據(jù)爆轟胞格特征尺度的定義，胞格橫向的尺寸λ為爆轟胞格的尺寸（圖9）。因此，自相關函數(shù)在y方向的變化體現(xiàn)了爆轟胞格的尺寸。根據(jù)圖8所示，自相關函數(shù)在三個方向均呈現(xiàn)周期性的變化，且其峰值逐漸降低，在y方向出現(xiàn)的第一個峰值，也是最大值，其對應的偏移量為爆轟胞格的尺寸，該值約為29，與概率密度函數(shù)法得到的胞格尺寸接近。

圖6 活化能Ea = 15時，根據(jù)全局閾值法提取的三波點的軌跡Fig. 6 Extracted triple-point tracks based on the global threshold method for Ea = 15

圖7 活化能Ea = 15時，三波點軌跡間距離的概率分布Fig. 7 Probability density function of the spacing between triple-point tracks for Ea = 15

圖8 活化能為15時，自相關函數(shù)在不同偏移方向的變化Fig. 8 Autocorrelation functions along different shift directions for Ea = 15

當活化能Ea= 20時，對數(shù)值胞格分別采用概率密度函數(shù)法和自相關函數(shù)法統(tǒng)計方法得到的變化曲線見圖10所示。圖10（a）是概率密度分布曲線，其最大峰值明顯小于活化能Ea= 15時概率密度分布的最大峰值，這是因為爆轟胞格的不規(guī)則程度增加，胞格尺寸的取值范圍增加，胞格的特征尺寸在所有胞格尺寸中占的比例有所下降，此時，胞格的特征尺寸約為60。圖10（b）是自相關函數(shù)在y方向偏移的變化曲線，隨著不規(guī)則程度的增加，自相關函數(shù)的周期性減弱，其峰值對應的爆轟胞格的特征尺寸為64。此時，兩種統(tǒng)計方法得到的胞格尺寸基本一致。

圖9 爆轟胞格的尺寸示意圖Fig. 9 Schematic diagram of the detonation cell size

圖10 活化能Ea = 20時爆轟胞格尺寸的統(tǒng)計分析Fig. 10 Statistical analysis of the detonation cell size for Ea = 20

隨著活化能的進一步增加，當活化能Ea= 25、27時，三波點的跡線開始出現(xiàn)合并，爆轟胞格非常不規(guī)則。圖11為活化能為25和27時的概率密度分布曲線，兩者的分布基本一致，分布范圍和峰值出現(xiàn)位置基本相同。其中活化能為25時的胞格尺寸為62，活化能為27時胞格尺寸為67。

圖11 活化能Ea = 25、27時爆轟胞格尺寸的概率密度分布Fig. 11 Probability density functions of the detonation cell size for Ea = 25, 27

圖12是活化能為25和27時的自相關函數(shù)在y方向的變化曲線。圖中第一個峰值分別發(fā)生在胞格尺寸為43（Ea= 25）和63（Ea= 27）的位置。值得注意的是，與規(guī)則胞格不同，這兩種不規(guī)則程度較高的胞格的自相關函數(shù)的第一個峰值不再是自相關性最高的值。在第一個峰值后，會再次出現(xiàn)峰值更高的自相關函數(shù)值，其值分別為121（Ea= 25）和120（Ea=27），約為第一個峰值出現(xiàn)位置的兩倍。根據(jù)自相關函數(shù)法對胞格尺寸的取值，自相關性最大的位置是胞格尺寸，因此，自相關法得到的胞格尺寸應是121（Ea= 25）和120（Ea= 27），其值約為概率密度法得到的尺寸的兩倍。該現(xiàn)象與Shepherd[6]的研究工作是一致的。之所以出現(xiàn)這樣大的偏差，這和兩種統(tǒng)計方法的特性有關。

圖12 活化能Ea = 25、27時爆轟胞格尺寸的自相關函數(shù)Fig. 12 Autocorrelation functions of the detonation cell size for Ea = 25, 27

概率密度法統(tǒng)計跡線間距時，忽略了跡線上能量的數(shù)值，也就是將跡線的能量以相同的權重進行統(tǒng)計，因此其峰值所對應的胞格尺寸是表觀上看到的胞格的最大概率尺寸。而自相關法在統(tǒng)計渦量信號時，跡線上的渦量大小信息作為自相關函數(shù)的參數(shù)加入計算，渦量大的跡線權重大。因此得到的自相關峰值是考慮了跡線渦量的大小的峰值。

在不規(guī)則的爆轟胞格中，出現(xiàn)了橫波的合并，使得三波點的渦量發(fā)生急劇變化，這種現(xiàn)象也體現(xiàn)在自相關函數(shù)值的變化上。而概率密度法僅能反映此時胞格大小出現(xiàn)的概率，無法體現(xiàn)這種橫波合并引發(fā)的物理現(xiàn)象。因此，兩種統(tǒng)計方法預測的胞格尺寸出現(xiàn)較大的偏差，而在規(guī)則胞格中，沒有橫波合并這種現(xiàn)象，故兩種統(tǒng)計方法預測的胞格尺寸基本吻合。

3 結 論

本文通過數(shù)值研究了不同活化能（Ea= 15、20、25、27）時的爆轟胞格，隨著活化能的增加，胞格的不規(guī)則程度明顯變大。同時，采用概率密度函數(shù)法和自相關函數(shù)法這兩種統(tǒng)計方法對數(shù)值模擬得到的胞格尺寸進行統(tǒng)計處理。當活化能為15、20時，這兩種統(tǒng)計方法得到的胞格尺寸較吻合。當活化能增加到25和27時，數(shù)值胞格出現(xiàn)橫波合并現(xiàn)象，自相關函數(shù)法得到的胞格尺寸大約是概率密度函數(shù)得到的胞格尺寸的兩倍。這是由于概率密度法沒有考慮三波點跡線上的物理量分布，僅是一種表觀上觀測到的胞格尺寸出現(xiàn)的最大概率，而自相關函數(shù)法考慮了跡線上的渦量分布，進而能反映橫波合并這種物理現(xiàn)象，因此，自相關函數(shù)法比傳統(tǒng)的概率密度函數(shù)法更能體現(xiàn)爆轟胞格尺寸這種爆轟的本征值。