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

(南京理工大學(xué) 瞬態(tài)物理國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210094)

0 引 言

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

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

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

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

1 物理模型和算法

1.1 控制方程和胞格的捕捉

對(duì)于氣相爆轟波而言，國(guó)內(nèi)外普遍采用Euler方程進(jìn)行數(shù)值研究，本文研究基于單步反應(yīng)的Euler方程：

其中，

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

式中q為反應(yīng)熱?；瘜W(xué)反應(yīng)采用單步不可逆反應(yīng)，反應(yīng)速率滿足：

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

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

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

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

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

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

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

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

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

自相關(guān)函數(shù)主要應(yīng)用在信號(hào)處理中，用來(lái)描述信號(hào)在一個(gè)時(shí)刻的取值與另一個(gè)時(shí)刻的取值的依賴關(guān)系。當(dāng)信號(hào)中有周期性分量出現(xiàn)時(shí)，自相關(guān)函數(shù)分布的極大值能夠很好的體現(xiàn)這種周期長(zhǎng)度。

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

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

2 結(jié)果與討論

2.1 不同活化能的爆轟胞格結(jié)構(gòu)

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

圖2 計(jì)算示意圖Fig. 2 Schematic diagram of the computational domain

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

圖3 不同時(shí)刻爆轟波傳播的流場(chǎng)圖Fig. 3 Flow fields of the detonation wave at different time instances

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

圖4 活化能Ea = 15時(shí)的三波點(diǎn)軌跡Fig. 4 Triple-point tracks for Ea = 15

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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