鄭凱，蔣軍成，邢志祥，吳凡，吳潔，余明高

（1 常州大學(xué)環(huán)境與安全工程學(xué)院，江蘇常州213164； 2 重慶大學(xué)煤礦災(zāi)害動力學(xué)與控制國家重點實驗室，重慶400044）

引 言

甲烷作為天然氣和煤礦瓦斯的主要成分，其在生產(chǎn)儲運過程中一旦發(fā)生泄漏或積聚形成爆炸濃度范圍內(nèi)的可燃氣體時，遇火源可發(fā)展為爆燃、爆轟以及從爆燃向爆轟的轉(zhuǎn)變等爆炸事故[1-3]。一般而言，管道中甲烷-空氣爆炸傳播是涉及湍流流動的化學(xué)反應(yīng)過程。隨著高性能計算技術(shù)的發(fā)展，大渦模擬(large eddy simulation, LES)已經(jīng)成為一種預(yù)測管道內(nèi)甲烷-空氣爆炸湍流反應(yīng)的重要手段[4-6]。相對于大渦模擬，雷諾平均模擬(Reynolds average Navier-Stokes，RANS)由于對Navier-Stokes (N-S)方程進行了時間平均，計算準(zhǔn)確性較差，難以精確預(yù)測爆炸火焰結(jié)構(gòu)。然而由于計算能力的限制，當(dāng)前大渦模擬中計算網(wǎng)格尺寸要遠大于待求解火焰厚度。為了解決這一問題，學(xué)者們提出了三種亞網(wǎng)格模型：火焰表面密度模型(flame surface density，F(xiàn)SD)[7]，增厚火焰模型(thickening flame，TF)[8-9]和G 方程模型[10-11]。

在TF模型中，通過提高燃燒組分?jǐn)U散系數(shù)并同比例降低反應(yīng)速率，可以使層流火焰速度保持不變但火焰厚度增大，以保證火焰在大渦尺寸網(wǎng)格上求解。由于火焰厚度增加導(dǎo)致的亞網(wǎng)格火焰表面的丟失可用亞網(wǎng)格尺度褶皺因子來評估，其為火焰總表面與可求解的表面積的比值[12]。為了避免增厚過程改變組分的輸運參數(shù)，采用火焰探測函數(shù)對火焰前鋒進行動態(tài)增厚[13]。當(dāng)前基于增厚火焰模型的大渦模擬技術(shù)已經(jīng)取得長足進展，且在預(yù)測預(yù)混和非預(yù)混火焰的可行性方面開展了大量研究。Selle等[14]采用簡化的甲烷-空氣兩步燃燒機理和增厚火焰模型開展了復(fù)雜幾何結(jié)構(gòu)中反應(yīng)流的大渦模擬，表明增厚火焰模型能夠預(yù)測復(fù)雜燃燒裝置的反應(yīng)特性。Wang 等[15]基于動態(tài)亞網(wǎng)格模型參數(shù)計算褶皺因子，結(jié)合動態(tài)增厚火焰開展大渦模擬分析了湍流本生燈火焰。張科等[16]采用基于釋熱率構(gòu)造火焰探測函數(shù)的動態(tài)增厚火焰模型，對甲烷-空氣同軸射流燃燒器的非預(yù)混燃燒進行可壓縮大渦模擬。Guo 等[17]基于改進的動態(tài)火焰探測函數(shù)，開展了耦合多步反應(yīng)機理的增厚火焰大渦模擬，分析了增厚因子對火焰結(jié)構(gòu)的影響。

對于管道內(nèi)甲烷-空氣爆炸的大渦模擬，當(dāng)前應(yīng)用較多的為FSD 模型，研究結(jié)果表明其能夠準(zhǔn)確地捕捉爆炸預(yù)混火焰結(jié)構(gòu)，精確預(yù)測超壓變化和火焰?zhèn)鞑チ鲌鎏卣鱗18-22]。但關(guān)于增厚火焰模型大渦模擬的當(dāng)前報道相對較少，且相關(guān)研究主要側(cè)重于分析障礙物引起的湍流對爆炸火焰加速過程的影響。Yu 等[23]將增厚火焰模型應(yīng)用到可壓縮燃燒預(yù)測并利用其捕捉了爆燃向爆轟轉(zhuǎn)變過程。Emami 等[24]利用增厚火焰模型和27步反應(yīng)機理，通過求解二維過濾后的Navier-Stokes 方程，研究了障礙物管道中預(yù)混氫氣-空氣的爆燃向爆轟轉(zhuǎn)變過程。Volpiani等[25]利用動態(tài)火焰褶皺函數(shù)對爆燃火焰通過連續(xù)障礙物的傳播過程進行了大渦模擬，結(jié)果表明大渦模擬能夠預(yù)測爆炸超壓等特征參數(shù)。

本文擬采用單步不可逆甲烷-空氣簡化燃燒機理開展三維大渦模擬，研究對象為出現(xiàn)郁金香結(jié)構(gòu)[26]且基本以層流形式傳播的管道內(nèi)甲烷-空氣爆炸火焰。模擬中采用Légier[13](TF1)與Durand[27](TF2)提出的兩種火焰探測函數(shù)對火焰面進行動態(tài)增厚?；谟嬎愕玫降幕鹧嫒S等值面、二維切片分析火焰動態(tài)傳播特征，將火焰?zhèn)鞑ニ俣?、超壓曲線與前期實驗結(jié)果進行對比以驗證模型可靠性。基于有效增厚因子分析火焰探測函數(shù)對大渦模擬結(jié)果的影響。

1 數(shù)值模型與計算步驟

1.1 增厚火焰模型

對于以預(yù)混火焰方式傳播的管道內(nèi)甲烷-空氣爆炸，可以通過增加擴散速率與降低反應(yīng)速率的方式實現(xiàn)火焰增厚。增厚后的組分輸運方程如下：

式中，ρ為密度；u為速度矢量；Yk是組分k的質(zhì)量分?jǐn)?shù)；Dk為組分k的擴散速率；F為增厚因子；(Q)為組分k的反應(yīng)速率，可以用Arrhenius 公式進行描述，其中Q表示用于計算反應(yīng)速率的物理量。

考慮到前鋒褶皺對火焰表面積的影響，將火焰前鋒的擴散系數(shù)和反應(yīng)速率乘以效率因子Ξ。本文中采用Zimont 提出的湍流火焰封閉模型(turbulent flame speed closure model，TFC)對效率因子進行評估[28]，該模型可改寫成如下形式：

式中，Ut是湍流火焰速度；A是常數(shù)，對于甲烷-空氣建議值為0.52；u'是亞網(wǎng)格脈動速度；Ul是層流火焰速度；δ是層流火焰厚度；lt為湍流長度尺度，可用下式計算得到：

式中，Cs是Smagorinsky常數(shù)，可用Smagorinsky-Lilly渦粘模型計算得到[29]；Δ是網(wǎng)格尺寸。

u'的計算公式為：

為應(yīng)變率張量:

效率因子Ξ，為亞網(wǎng)格湍流速度在FΔ和Δ尺寸上的比值：

LES過濾后的組分輸運方程為：

其中,-表示過濾值；～代表平均值。

1.2 火焰探測函數(shù)

為了避免非反應(yīng)區(qū)域由于擴散系數(shù)的增加而導(dǎo)致的組分混合與傳熱問題，只在火焰前鋒內(nèi)進行動態(tài)增厚，火焰前鋒的動態(tài)增厚系數(shù)由火焰探測函數(shù)Ω進行控制。理論上，在反應(yīng)區(qū)域內(nèi)Ω的值為1，在區(qū)域外則為0。本文中采用兩種火焰探測函數(shù)(TF1和TF2)。

TF1為[13]:

TF2為[27]:

其中，c為反應(yīng)進程變量：

基于火焰探測函數(shù)TF1 和TF2，可計算得到有效增厚因子：

1.3 反應(yīng)機理

甲烷-空氣爆炸的預(yù)混火焰反應(yīng)過程用單步不可逆燃燒反應(yīng)機理進行描述[15]：

反應(yīng)速率可用Arrhenius公式進行描述：

式中，A是指前因子；Ea是活化能；ρ為密度；nk是組分k的反應(yīng)指數(shù)；Mk是組分k的摩爾質(zhì)量；R為理想氣體常數(shù)；T為溫度。反應(yīng)機理的速率常數(shù)如表1中所示。根據(jù)該反應(yīng)機理計算得到的層流火焰速度和絕熱火焰溫度分別為0.40 m/s和2328 K。

表1 反應(yīng)機理的速率常數(shù)Table 1 Rate constants for the one step mechanism

1.4 參數(shù)設(shè)置

模擬在尺寸為1000 mm×100 mm×100 mm 的三維管道中進行?？刂品匠滩捎糜邢摅w積法離散，壓力-速度場采用SIMPLE 算法耦合。采用二階上迎風(fēng)格式與中心差分格式離散對流項與擴散項。整個計算區(qū)域內(nèi)初始網(wǎng)格尺寸為2.5 mm。為了確保最大限度求解湍動能，采用自適應(yīng)網(wǎng)格對火焰前鋒附近區(qū)域網(wǎng)格進行動態(tài)加密。加密過程是基于溫度梯度來進行的，其中網(wǎng)格粗化與細化的閾值分別為30 K/m 與100 K/m。每個網(wǎng)格最多加密2 次，加密之后最小尺寸為0.625 mm。增厚因子Fmax的值為10。計算對象為當(dāng)量比為1的甲烷-空氣混合氣體。各組分的比熱和黏度用五次多項式和Sutherland 公式進行計算[30-31]。組分輸運系數(shù)與熱導(dǎo)率分別采用Chapman-Enskog 公式與動力學(xué)理論計算得到[17]。混合反應(yīng)物的熱物理特性由理想氣體混合法則計算得到。初始溫度和壓力分別設(shè)定為295 K 和0.1 MPa。管道壁面采用無反射固體邊界條件。壓力監(jiān)測點位于管道點火端中心處。在未燃氣體中設(shè)置一個半徑為5 mm 的半球形區(qū)域，假設(shè)該區(qū)域內(nèi)甲烷/空氣完全反應(yīng)，區(qū)域內(nèi)溫度設(shè)置為2100 K。Wen等[4]的研究表明將點火半徑設(shè)置為5 mm 時火焰鋒面位置的大渦模擬與實驗結(jié)果吻合度較高。管道中心軸線設(shè)為X軸(單位：m)，半球形點火區(qū)域球心坐標(biāo)為(0，0，0)，即火焰沿X軸正向傳播。時間步長設(shè)置為小于1×10-5s，并隨著火焰速度增加逐漸減小。每個時間步長至少需要60 次迭代。動量方程的殘差小于1×10-5，能量和物質(zhì)輸運方程的殘差分別小于1×10-6。

2 結(jié)果與討論

2.1 火焰?zhèn)鞑ソY(jié)構(gòu)與流場

對于本文大渦模擬所采用管道中傳播的甲烷-空氣爆炸預(yù)混火焰，前期實驗結(jié)果表明其傳播會經(jīng)歷經(jīng)典“郁金香”火焰形成的四個階段，即：半球形、指形、火焰接觸管道壁面到形成平板形狀、郁金香形狀[32]。圖1 分別給出了TF1 和TF2 時郁金香火焰形成過程中溫度二維切片，代表了火焰?zhèn)鞑デ颁h結(jié)構(gòu)演變過程。從圖1 中可看出，無論是對于TF1 還是TF2，大渦模擬均能夠完整再現(xiàn)郁金香火焰形成過程。點火完成后，在11 ms 時可觀測到火焰以半球形傳播。此時，TF1 和TF2 火焰?zhèn)鞑ソY(jié)構(gòu)和距離相差不大。此后，由于管道壁面的限制，半球形火焰會轉(zhuǎn)變成指形。在28 ms時，TF1 和TF2 中均可觀測到火焰以指形結(jié)構(gòu)傳播。但是相對于TF1 (0.11 m)，TF2 (0.13 m)中火焰前鋒傳播的要稍遠些，代表了更大的火焰?zhèn)鞑ニ俣?。隨著火焰的進一步傳播，在接觸到管道壁面后，火焰前鋒面積開始減少，在圖1 中39 ms 可看到TF1 和TF2 中靠近管道壁面處火焰前鋒逐漸消失，表明進入火焰接觸管道壁面到形成平板形狀階段。此時，TF2 仍然出現(xiàn)了相對較大的傳播距離?；鹧媲颁h面積的進一步減少最終導(dǎo)致平板火焰前鋒的出現(xiàn)。在實驗中，平板火焰出現(xiàn)時間為44 ms，而模擬結(jié)果分別為47 ms (TF1)和45 ms(TF2)，實驗和模擬結(jié)果基本吻合。此后，火焰前鋒出現(xiàn)逆轉(zhuǎn)，火焰靠近管道壁面處速度大于管道中心位置，火焰前鋒出現(xiàn)凹陷，最終形成經(jīng)典的“郁金香”火焰結(jié)構(gòu)。圖1 中86 ms 時火焰切片代表了TF1 和TF2 預(yù)測的郁金香火焰結(jié)構(gòu)?！坝艚鹣恪被鹧娉霈F(xiàn)后會保持這一結(jié)構(gòu)傳播直到整個爆炸過程結(jié)束。

圖1 郁金香火焰形成過程中不同時刻的溫度二維切片F(xiàn)ig.1 2D-slices of the temperature field at different instants during the tulip shaped flame formation

圖2 為TF1 中不同時刻大渦模擬預(yù)測的火焰前鋒三維1200 K 溫度等值面結(jié)構(gòu)(左側(cè))與流場演化特征(右側(cè))。對于TF1 與TF2，除傳播速度外，大渦模擬預(yù)測的火焰?zhèn)鞑ソY(jié)構(gòu)與流場基本一致。在流場結(jié)構(gòu)中，兩條曲線為火焰探測函數(shù)Ω=0.1 時的等值線，以劃分已燃區(qū)、爆炸反應(yīng)區(qū)和未燃區(qū)。在圖2中11 ms 時刻，可以看到火焰以光滑的半球形傳播，此時火焰處于自由膨脹階段。因此，在流場中靠近火焰前鋒流線由反應(yīng)區(qū)呈扇形指向未燃區(qū)。在靠近管道壁面處流線方向發(fā)生改變，最終在遠離火焰前鋒的未燃區(qū)所有流線均平行于管道軸線。在28 ms火焰由半球形轉(zhuǎn)變成指形后，受到管道壁面的限制，三維結(jié)構(gòu)中火焰封面在靠近管道壁面時被壓平。此時，已燃產(chǎn)物膨脹作用推動火焰前鋒高速傳播。對應(yīng)流場中流線起源于靠近管道壁面處的火焰前鋒，指向管道軸線方向并流向未燃氣體區(qū)域。火焰前鋒接觸到管道壁面后，在39 ms 可以看出火焰在接觸到管道壁面處的鋒面消失，進而前鋒總面積大大減小?；鹧媲颁h面積的減小意味著反應(yīng)速率的降低，導(dǎo)致該階段火焰尖端傳播速度出現(xiàn)下降。此時，火焰前鋒附近流線主要起源于火焰前鋒與管道壁面的接觸點。隨著火焰前鋒面積的逐步減少，平板火焰在47 ms 時刻出現(xiàn)。流場中在未燃氣體區(qū)域流線保持正向流動，在已燃產(chǎn)物區(qū)出現(xiàn)兩個渦旋，其旋轉(zhuǎn)方向均為由管道壁面指向中心軸線。渦旋的旋轉(zhuǎn)推動靠近壁面處氣體向前但拖拽中心軸線處氣體向產(chǎn)物流動，最終導(dǎo)致了“郁金香”火焰的出現(xiàn)。圖中86 ms 為完全形成的經(jīng)典“郁金香”火焰結(jié)構(gòu)。在“郁金香”火焰形成后，在三維圖中可看到火焰呈現(xiàn)四個錐形結(jié)構(gòu)前鋒。在流場中未燃氣體保持正向流動，而已燃氣體區(qū)域基本上保持逆向流動。

圖2 TF1中不同時刻的溫度三維等值面和二維切片流場結(jié)構(gòu)Fig.2 3D iso-surface by temperature of 1200 K and flow field with TF1 at different instants

2.2 火焰?zhèn)鞑ニ俣扰c超壓

對于管道中可燃氣體爆炸事故，火焰?zhèn)鞑ニ俣群统瑝菏呛饬勘ㄎｋU程度的重要參數(shù)。本文中大渦模擬可靠性可通過作者前期實驗與本文模擬得到的火焰?zhèn)鞑ニ俣扰c超壓曲線定量對比來進行驗證[32]。

圖3中給出了火焰?zhèn)鞑ニ俣入S時間變化曲線的實驗和大渦模擬對比結(jié)果。爆炸發(fā)生后，實驗與大渦模擬中都可觀測到火焰?zhèn)鞑ニ俣鹊闹笖?shù)形式增長。這一現(xiàn)象主要是由于火焰在半球形和指形階段前鋒面積的指數(shù)增長所引起[33]?；鹧娴募铀賯鞑ミ^程因前鋒接觸到管道壁面而結(jié)束。圖3 中，實驗觀測到的火焰?zhèn)鞑ニ俣确逯岛拖鄳?yīng)時間分別為17.5 m/s 和35 ms，TF1 和TF2 計算得到的對應(yīng)值分別為15.4 m/s、35 ms 和15.7 m/s、32 ms。顯然，TF1和TF2 預(yù)測的火焰?zhèn)鞑ニ俣确逯祷疽恢?，稍小于實驗結(jié)果，但實驗與大渦計算結(jié)果最大相差僅12%。而TF1 預(yù)測的到達火焰?zhèn)鞑ニ俣确逯禃r間與實驗結(jié)果相同，TF2 相應(yīng)值則稍大。在火焰接觸到管道壁面后，鋒面面積減少，火焰?zhèn)鞑ニ俣纫仓饾u降低，并在形成平板火焰時達到最小值。伴隨著郁金香火焰的形成，火焰再次出現(xiàn)加速。最終伴隨著火焰?zhèn)鞑ニ俣鹊闹饾u減小，火焰前鋒以郁金香形狀傳播直到結(jié)束?；鹧?zhèn)鞑ミ^程中，實驗和大渦模擬得到的火焰?zhèn)鞑ニ俣惹€基本吻合，但TF1 預(yù)測結(jié)果準(zhǔn)確度相對較高。

圖3 火焰?zhèn)鞑ニ俣入S時間的變化Fig.3 The time curves of flame speed

圖4對比了實驗與大渦模擬得到的超壓隨時間變化曲線。相對于火焰?zhèn)鞑ニ俣?，超壓的實驗與大渦模擬曲線在“郁金香”火焰階段存在較大差異。在圖4中，與火焰?zhèn)鞑ニ俣茸兓厔菹嘁恢碌氖?，爆炸超壓在火焰?zhèn)鞑コ跗谝渤手笖?shù)增長。這是因為火焰?zhèn)鞑ニ俣鹊闹笖?shù)增長表明燃燒反應(yīng)速率呈指數(shù)增長，導(dǎo)致產(chǎn)物生成速率和放熱速率也相應(yīng)提高，最終造成了超壓的指數(shù)上升。因此，實驗中觀測到的較大火焰速度導(dǎo)致超壓相對于大渦模擬結(jié)果要上升的更快，且兩者之間差值在火焰?zhèn)鞑ニ俣冗_到最大值時達到最大。在火焰鋒面接觸管道壁面后，超壓曲線斜率出現(xiàn)下降，該過程持續(xù)到平板火焰的形成時刻。此后，在“郁金香”火焰形成后，超壓基本上呈線性增長。TF1 和TF2 預(yù)測結(jié)果均與實驗基本吻合，但后者增長要稍微快于前者。

圖4 超壓隨時間的變化Fig.4 The time curves of pressure

2.3 有效增厚因子

在增厚火焰模型中，有效增厚因子是影響大渦計算結(jié)果的重要參數(shù)。圖5 給出了39 ms 時TF1 和TF2的Ω、F、Ξ以及Ξ×F分布二維切片，代表了火焰?zhèn)鞑ミ^程中兩種火焰探測函數(shù)對計算結(jié)果的影響效果。根據(jù)Guo 等[17]的觀點，一個理想的Ω值應(yīng)該在火焰前鋒內(nèi)為1，而在火焰前鋒外則為0。從圖5中可看出，TF1 和TF2 預(yù)測結(jié)果均能達到該標(biāo)準(zhǔn)。在靠近火焰前鋒中心位置，Ω的值基本接近1。隨著遠離火焰前鋒中心，Ω的值逐漸減小，并在已燃區(qū)和未燃區(qū)中減小到0。顯然，TF1捕捉到的反應(yīng)區(qū)域要明顯大于TF2。TF1 預(yù)測的Ω值大多為1，只有在靠近已燃與未燃區(qū)域邊界處少量值小于0.5，但總體上是從1突然下降到0的。而在TF2中，只有在預(yù)測區(qū)域中心位置處少量值為1，且從1 下降到0 的過程是漸進的。

在大渦模擬中，效率因子Ξ為在過濾網(wǎng)格上火焰總表面與已求解火焰表面積的比率，其值與火焰增厚導(dǎo)致的亞網(wǎng)格尺度火焰面積丟失有關(guān)[15]。根據(jù)流場分析結(jié)果可知，火焰在以指形結(jié)構(gòu)傳播過程中管道內(nèi)基本上保持層流流動。所以在圖5中整個管道大部分區(qū)域內(nèi)Ξ均為1。只有在火焰前鋒及靠近火焰前鋒附近的已燃區(qū)域內(nèi)Ξ的值大于1。在已燃區(qū)域內(nèi)，Ξ主要是受到火焰前鋒與管道壁面接觸后導(dǎo)致該區(qū)域流線方向改變的影響，在該區(qū)域內(nèi)Ξ的值較小(小于1.18)。而在火焰前鋒尖端，Ξ的分布則主要受到火焰拉伸效應(yīng)的影響，故其在該區(qū)域內(nèi)Ξ的值較大(接近于1.5)?？傮w上，TF1 和TF2 兩種工況中計算得到的Ξ分布基本保持一致，表明火焰捕捉函數(shù)對Ξ基本上沒有影響。

由于TF1和TF2中Ξ的分布基本一致，可知Ξ×F的分布只受到F的影響。而根據(jù)有效增厚因子計算公式可知，當(dāng)Fmax一定時，F(xiàn)是Ω的函數(shù)。因此，圖5中F與Ξ×F的分布規(guī)律均與Ω相似，且TF1和TF2預(yù)測的F最大值均小于Fmax= 10。對于TF1和TF2，受到Ξ分布的影響，Ξ×F的最大值均分布在火焰前鋒上。在反應(yīng)區(qū)域內(nèi)，TF1 中Ξ×F大部分?jǐn)?shù)值要大于10，而在TF2 中只有靠近火焰前鋒中心位置少部分區(qū)域的Ξ×F大于10。因此，總體上TF1對爆炸火焰的增厚效應(yīng)要強于TF2。而增厚效應(yīng)的增強，可能會導(dǎo)致火焰前鋒燃燒反應(yīng)求解的不充分，進而導(dǎo)致TF2 預(yù)測的火焰?zhèn)鞑ニ俣扰c超壓上升均稍大于TF1。

3 結(jié) 論

本文采用簡化的單步不可逆甲烷-空氣燃燒機理，基于增厚火焰模型開展了管道內(nèi)甲烷-空氣爆炸大渦模擬，研究采用了Zimont 提出的湍流火焰封閉模型評估效率因子，分析了兩種火焰探測函數(shù)對計算結(jié)果的影響，得到如下結(jié)論。

(1)基于增厚火焰模型的大渦模擬能夠準(zhǔn)確再現(xiàn)閉口管道中甲烷-空氣爆炸火焰?zhèn)鞑ミ^程，兩種火焰捕捉函數(shù)下均能觀測到郁金香形成經(jīng)歷的四個階段，且實驗和模擬結(jié)果中平板火焰形成特征時間基本吻合。

(2)火焰?zhèn)鞑ニ俣鹊拇鬁u模擬與實驗曲線基本重合，兩者火焰?zhèn)鞑ニ俣确逯底畲笙嗖顑H12%。超壓的實驗與大渦模擬曲線在火焰指數(shù)加速階段存在較大差異，但在郁金香火焰形成后超壓增長基本一致。TF2的大渦模擬預(yù)測的火焰?zhèn)鞑ニ俣扰c超壓準(zhǔn)確性要高于TF1。

(3)火焰探測函數(shù)對爆炸過程中褶皺因子變化基本上沒有影響，較大值主要分布在火焰前鋒尖端區(qū)域。TF1 捕捉到的反應(yīng)區(qū)域要明顯大于TF2，總體上前者對火焰的增厚效應(yīng)要強于后者。