基于簡化反應機理的甲烷-空氣爆炸增厚火焰大渦模擬

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

（1 常州大學環境與安全工程學院，江蘇常州213164； 2 重慶大學煤礦災害動力學與控制國家重點實驗室，重慶400044）

引 言

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

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

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

本文擬采用單步不可逆甲烷-空氣簡化燃燒機理開展三維大渦模擬，研究對象為出現郁金香結構[26]且基本以層流形式傳播的管道內甲烷-空氣爆炸火焰。模擬中采用Légier[13](TF1)與Durand[27](TF2)提出的兩種火焰探測函數對火焰面進行動態增厚。基于計算得到的火焰三維等值面、二維切片分析火焰動態傳播特征，將火焰傳播速度、超壓曲線與前期實驗結果進行對比以驗證模型可靠性。基于有效增厚因子分析火焰探測函數對大渦模擬結果的影響。

1 數值模型與計算步驟

1.1 增厚火焰模型

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

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

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

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

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

u'的計算公式為：

為應變率張量:

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

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

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

1.2 火焰探測函數

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

TF1為[13]:

TF2為[27]:

其中，c為反應進程變量：

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

1.3 反應機理

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

反應速率可用Arrhenius公式進行描述：

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

表1 反應機理的速率常數Table 1 Rate constants for the one step mechanism

1.4 參數設置

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

2 結果與討論

2.1 火焰傳播結構與流場

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

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

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

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

2.2 火焰傳播速度與超壓

對于管道中可燃氣體爆炸事故，火焰傳播速度和超壓是衡量爆炸危險程度的重要參數。本文中大渦模擬可靠性可通過作者前期實驗與本文模擬得到的火焰傳播速度與超壓曲線定量對比來進行驗證[32]。

圖3中給出了火焰傳播速度隨時間變化曲線的實驗和大渦模擬對比結果。爆炸發生后，實驗與大渦模擬中都可觀測到火焰傳播速度的指數形式增長。這一現象主要是由于火焰在半球形和指形階段前鋒面積的指數增長所引起[33]。火焰的加速傳播過程因前鋒接觸到管道壁面而結束。圖3 中，實驗觀測到的火焰傳播速度峰值和相應時間分別為17.5 m/s 和35 ms，TF1 和TF2 計算得到的對應值分別為15.4 m/s、35 ms 和15.7 m/s、32 ms。顯然，TF1和TF2 預測的火焰傳播速度峰值基本一致，稍小于實驗結果，但實驗與大渦計算結果最大相差僅12%。而TF1 預測的到達火焰傳播速度峰值時間與實驗結果相同，TF2 相應值則稍大。在火焰接觸到管道壁面后，鋒面面積減少，火焰傳播速度也逐漸降低，并在形成平板火焰時達到最小值。伴隨著郁金香火焰的形成，火焰再次出現加速。最終伴隨著火焰傳播速度的逐漸減小，火焰前鋒以郁金香形狀傳播直到結束。火焰傳播過程中，實驗和大渦模擬得到的火焰傳播速度曲線基本吻合，但TF1 預測結果準確度相對較高。

圖3 火焰傳播速度隨時間的變化Fig.3 The time curves of flame speed

圖4對比了實驗與大渦模擬得到的超壓隨時間變化曲線。相對于火焰傳播速度，超壓的實驗與大渦模擬曲線在“郁金香”火焰階段存在較大差異。在圖4中，與火焰傳播速度變化趨勢相一致的是，爆炸超壓在火焰傳播初期也呈指數增長。這是因為火焰傳播速度的指數增長表明燃燒反應速率呈指數增長，導致產物生成速率和放熱速率也相應提高，最終造成了超壓的指數上升。因此，實驗中觀測到的較大火焰速度導致超壓相對于大渦模擬結果要上升的更快，且兩者之間差值在火焰傳播速度達到最大值時達到最大。在火焰鋒面接觸管道壁面后，超壓曲線斜率出現下降，該過程持續到平板火焰的形成時刻。此后，在“郁金香”火焰形成后，超壓基本上呈線性增長。TF1 和TF2 預測結果均與實驗基本吻合，但后者增長要稍微快于前者。

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

2.3 有效增厚因子

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

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

由于TF1和TF2中Ξ的分布基本一致，可知Ξ×F的分布只受到F的影響。而根據有效增厚因子計算公式可知，當Fmax一定時，F是Ω的函數。因此，圖5中F與Ξ×F的分布規律均與Ω相似，且TF1和TF2預測的F最大值均小于Fmax= 10。對于TF1和TF2，受到Ξ分布的影響，Ξ×F的最大值均分布在火焰前鋒上。在反應區域內，TF1 中Ξ×F大部分數值要大于10，而在TF2 中只有靠近火焰前鋒中心位置少部分區域的Ξ×F大于10。因此，總體上TF1對爆炸火焰的增厚效應要強于TF2。而增厚效應的增強，可能會導致火焰前鋒燃燒反應求解的不充分，進而導致TF2 預測的火焰傳播速度與超壓上升均稍大于TF1。

3 結 論

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

(1)基于增厚火焰模型的大渦模擬能夠準確再現閉口管道中甲烷-空氣爆炸火焰傳播過程，兩種火焰捕捉函數下均能觀測到郁金香形成經歷的四個階段，且實驗和模擬結果中平板火焰形成特征時間基本吻合。

(2)火焰傳播速度的大渦模擬與實驗曲線基本重合，兩者火焰傳播速度峰值最大相差僅12%。超壓的實驗與大渦模擬曲線在火焰指數加速階段存在較大差異，但在郁金香火焰形成后超壓增長基本一致。TF2的大渦模擬預測的火焰傳播速度與超壓準確性要高于TF1。

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