湍流狀態下雙燃料混合物自燃和火焰傳播過程

周 磊，劉宗寬，袁成威，張 帆，衛海橋

(天津大學 內燃機燃燒學國家重點實驗室，天津 300350)

隨著化石能源日益枯竭，燃用石油質液體燃料帶來的環境問題日漸嚴峻，內燃機節能減排技術成為研究熱點.為控制船舶發動機對海洋環境的污染，各國政府和國際組織制定的排放法規更加嚴格.國際海事組織(IMO)出臺的TierⅢ排放法規于2016 年開始生效，規定船用發動機NOx排放較上一階段標準降低80%.

天然氣最主要的成分是甲烷，具有高氫碳比、燃燒清潔和儲量豐富等優點，有望成為發動機新型清潔替代能源的潛力.研究表明，柴油引燃天然氣船用雙燃料發動機可有效降低二氧化碳和顆粒物的排放，并可在采用稀薄燃燒技術時大幅降低NOx排放[1]，天然氣船機只需簡單地增加排放后處理設備就能滿足TierⅢ排放要求.然而由于天然氣具有比熱容大、著火極限高和火焰傳播速度較慢等理化特性，以及雙燃料混合物復雜的化學動力學原理，共同導致雙燃料低速機低負荷著火及燃燒穩定性差.同時，低速機缸內伴有湍流且空間尺度大，湍流作用下雙燃料混合氣的著火及火焰傳播過程更加復雜；此外，缸內氣流運動會造成燃料混合物濃度分層和溫度分層.因此，探究湍流對著火及火焰發展的影響有助于理解雙燃料發動機缸內自燃和燃燒模式轉變的過程.

Mastorakos 等[2]通過直接數值模擬研究湍流混合層中燃燒過程，發現湍流會縮短混合物滯燃期，而標量耗散率決定著火位置，活性最強且標量耗散率低的混合分數位置最容易發生著火.此外，Pitsch 等[3]模擬研究了正庚烷噴霧的著火過程，發現標量耗散率降低可提高燃燒速度.Jin 等[4]研究了湍流對三重火焰結構的影響，發現由于湍流的存在，沿著化學當量比混合分數對應的位置會形成多個鋒面火焰，并在傳播過程中接觸連結.Wang 等[5]通過直接模擬手段研究雙燃料層流混合層中著火到火焰穩定傳播過程的特征時間以及特征時間與壓力、混合層厚度、預混氣溫度和當量比等參數的擬合關系，發現滯燃期與溫度、壓力密切相關，峰值溫度與混合層厚度以及預混當量比相關，穩態傳播特征時間和混合層厚度以及層流火焰速度有關.

綜上所述，針對湍流狀態下雙燃料混合物著火及火焰傳播過程的研究還鮮見報道，尤其在低速機缸內高溫、高壓熱力學狀態下，對湍流流場中雙燃料著火及燃燒機制的分析不夠全面.基于此，筆者采用直接數值模擬(DNS)方法分析湍流流場對雙燃料混合層著火和火焰發展過程的影響，探究雙燃料著火特性、混合層中火焰傳播結構以及火焰發展的特征.模擬研究中分別以甲烷和正庚烷作為天然氣和柴油的替代物，相關結論可為雙燃料低速船機高效、清潔燃燒技術的開發提供一定理論依據.

1 數值模型設置

1.1 數值計算方法

雙燃料發動機內的湍流燃燒由連續方程、N-S 方程以及組分和能量的輸運方程組成，筆者使用一種高階方法數值求解控制方程[6-7].求解過程中，時間積分項采用二階對稱分裂算子方法，即在兩個半步的擴散和對流時間步之間積分一個完整的化學反應速率時間步；時間步的推進采用二階龍格-庫塔(Runge-Kutta)結合二階Adam-Bashforth(AB2)格式；空間差分采用六階中心差分格式；組分和能量方程中的對流項采用五階 WENO(weighted essentially non-oscillatory scheme)格式求解.數值方法的具體描述參見文獻[7].

選用由Liu 等[8]提出的正庚烷簡化機理，包含45組分、112 步反應.機理涵蓋正庚烷自燃、高溫化學以及甲烷氧化動力學，可準確預測高壓下正庚烷滯燃期[8-9]、火焰傳播速度和絕熱火焰溫度[10]等關鍵參數.此外，該機理廣泛用于預混甲烷混合氣下正庚烷噴霧、著火及燃燒的相關研究，與試驗結果一致性良好[11].

1.2 計算區域及參數定義

圖1 所示二維矩形計算區域尺寸為3.2 mm×3.2 mm，計算區域由2 0482個均勻網格劃分，大約420 萬個網格，網格解析度可達到1.56 μm[12].計算區域中正庚烷與甲烷/空氣均勻混合氣之間設置層流混合層，混合層中熱分層和組分分層根據雙曲正切函數進行設定，即

圖1 二維湍流計算區域示意Fig.1 2D schematic of turbulence-like computational domain

式中：f (x)代表組分質量分數或溫度隨空間位置x變化的通用參數；fR和 fL分別為右側值和左側值；cx 代表中間位置；δ為初始混合層厚度.

根據雙曲正切函數的特性，正庚烷質量分數為1%～90%的實際分布范圍約為初始混合層厚度的3倍.計算區域初始條件設定見表1.依據合成湍流方法[13]添加湍流場，初始湍流場使用Passot 等[14]提出的各向同性能譜方程設置.平均流速設置為0，均方根流速(RMS) u＇ ＝0.5 m/s[15].積分長度尺度lλ＝0.8 mm，為區域長度的1/4，柯爾莫哥洛夫長度尺度lk＝22μm，積分尺度雷諾數為120.

表1 二維湍流混合層參數Tab.1 Main parameters for 2D turbulence-like mixing layer

需要指出，混合層厚度δ和混合區域厚度不是同一概念，混合區域大約為3δ厚度.計算區域各組分的空間分布見圖2.計算區域左、右(垂直方向)為非反射的出流條件，上、下(水平方向)為周期性邊界條件.

圖2 初始狀態下計算區域主要組分的空間分布Fig.2 Initial distribution of major species for computational domain

研究中所用的參數混合分數ξ根據氮氣質量分數進行定義，公式為

2 模擬結果與數據分析

2.1 湍流場雙燃料著火和火焰發展

運用CHEMKIN 軟件中零維均勻反應器計算所得混合氣滯燃期在混合分數空間的分布如圖3 所示.第一階段滯燃期極小值τ1,min=0.278 ms，對應混合分數為ξMR,1=0.240.第二階段滯燃期極小值τ2,min=0.374 ms，對應的混合分數為ξMR,2=0.275.由初始當量比的分布可知，第一、二階段滯燃期極小值均發生在濃混合物中(混合物總體當量比φ＞1).第二階段著火位置對應溫度低于負溫度系數(NTC)溫度區間的下限為850 K，這是因為初始溫度繼續升高，NTC 現象將導致計算區域混合物滯燃期變長，計算區域反應活性最強的混合分數對應ξMR,2位置，此處混合物濃度合適且不受NTC 效應影響，故自燃傾向最大、滯燃期最短.

圖3 滯燃期、初始溫度和當量比在混合分數空間上分布Fig.3 Distribution of ignition delay times(IDs)，initial temperature and equivalence ratio in mixture fraction ξ space

溫度梯度的變化可表征化學反應空間位置和放熱強度的變化，圖4 展示了不同時刻下計算區域溫度梯度的二維空間分布.計算區域局部溫度和關鍵組分的變化如圖5 所示.

圖4a 中，t＝0.313 ms 時刻白色虛線代表混合物分數ξMR,1=0.240 等值線.圖4b 所示t＝0.383 ms 時刻白色實線代表ξMR,2=0.275 等值線.第一階段著火發生后的t＝0.313 ms 時刻、沿ξMR,1等值線左側更濃混合物位置逐漸發生反應放熱，形成多個連接的放熱區域.結合圖5a 中t＝0.313 ms 時刻線條A 上溫度和組分分布特征，發現溫度梯度峰值位置和低溫化學表征組分過氧化氫酮KET(OC7OOH)峰值重合，且H2O2和HO2質量分數逐漸升高.文獻[16]指出第一階段著火后，H2O2逐漸積累增多，HO2在第一階段和第二階段著火后快速上升形成峰值.綜合可知，此處火焰為正庚烷低溫化學反應放熱主導形成的冷焰.文獻[17]研究混合層中低溫化學反應時指出，冷焰逐漸形成于比ξMR,1更濃的混合分數位置，與筆者分析結論一致.

第二階段著火(τig=0.372 ms)發生后，根據圖4b中t＝0.383 ms 時刻的空間放熱情況分布可知，沿著ξMR,2等值線右側稀混合物中生成多個膨脹核心，正庚烷側的濃混合物中形成薄條狀的火焰結構.由圖5b 中線條B 上的溫度、溫度梯度和組分分布可知，膨脹核心為著火后形成的高溫膨脹核心.膨脹核心邊緣放熱形成兩個溫度梯度峰值，核心內部低溫反應組分KET 快速消耗，高溫反應區H2O2和HO2也在第二階段著火后開始反應消耗.薄條狀的放熱區域KET 快速生成，定義KET 組分質量分數峰值超過20%區域為低溫化學反應區域(綠色區域)，與窄條狀區域基本重合，可知其本質是低溫化學反應形成的冷焰區域.

隨著反應進行至圖4c 所示t＝0.403 ms 時刻，多個高溫膨脹核心逐漸連結，形成高溫的火焰傳播前沿.正庚烷側和甲烷側高溫火焰前沿區域為高溫已燃區.在正庚烷側(濃混合物側)火焰下游仍存在低溫化學形成的冷焰區，甲烷側火焰前沿反應更強，溫度也快速升高.該時刻的突出特點是出現同時存在濃預混冷焰和高溫火焰的雙火焰結構，如圖5c 所示.

圖5 標志線條上溫度、溫度梯度和關鍵組分質量分數分布Fig.5 Spatial distribution of temperature，temperature gradient and key species along lines

反應繼續進行如圖4d 和圖4e 所示，正庚烷側傳播方向高溫火焰前沿和低溫冷焰之間的距離逐漸縮短，更濃混合物中的火焰前沿處當量比高、氧氣含量少；而濃引燃油側反應逐漸變緩，發生中溫化學反應.相較于層流混合層，湍流使反應前沿變得更加褶皺，反應持續進行，甲烷側火焰鋒面褶皺程度低于正庚烷側.

圖4 不同時刻下溫度梯度在計算區域的空間分布Fig.4 Distribution of temperature gradient for a sequence of time instants in the computational domain

第二階段著火后在ξMR,2右側生成多個高溫膨脹核心，為探究湍流對自燃點和膨脹核心出現位置的影響，檢測了著火時刻、膨脹核心在混合分數梯度場及ξMR,2等值線上混合分數梯度的分布特征.其中，混合分數梯度表示混合分數在物理空間上的梯度變化，即式中：ξ表示混合物分數.混合分數梯度是標量耗散率的組成部分，可用于解釋湍流等因素導致的混合不均對熱量和組分的耗散作用[2,18].圖6 為湍流混合層中高溫膨脹核心分布.其中，標記區域Z1～Z4用于現象解釋，Z2和Z4為高溫膨脹核心區域.

由圖6a 可知，在ξMR,2等值線左側，混合分數梯度變動很大，而右側相對平緩，只在Z1和Z3兩個區域存在大的梯度波動.而該區域并不存在高溫膨脹核心，即該區域不適合著火和化學反應發展.Z2和Z4兩個高溫膨脹核心生成在當量比梯度平緩區域.分析可知，混合分數梯度大的區域，基團和熱量快速耗散，不利于著火核心的形成，即使著火核心形成，熱量也很快被耗散掉，無法繼續燃燒膨脹，故Z1和Z3區域無法形成高溫膨脹核心.初始火核更容易生成在ξMR等值線上標量耗散率小的區域.

圖6 湍流混合層中高溫膨脹核心分布Fig.6 Distribution of high-temperature expanding kernels in the turbulence-like mixing layer

為進一步分析湍流混合層中燃料著火和火焰發展過程，圖7 展示著火后溫度梯度、OH 質量分數、溫度峰值和CH4質量分數在混合分數空間的分布特征.由溫度梯度變化和高溫化學表征組分OH 質量分布可知，高溫化學發生在ξMR,2=0.275 位置(濃混合物且當量比大于1)，然后向兩側傳播，分別進入正庚烷側和甲烷側，甲烷側的反應強度遠大于正庚烷側火焰前沿的反應強度.稀引燃油側傳播的火焰前沿逐漸向化學當量比混合分數(ξst=0.024 9)位置傳播.2D 湍流混合層算例由于混合層厚度較大以及著火后膨脹作用，ξst位置超出計算區域外，未發現邊緣火焰火焰結構和預混甲烷的引燃過程[4].圖7c 顯示著火后計算區域溫度峰值在混合物空間的變化，發現溫度峰值逐漸向化學當量比混合分數位置傳播.圖7d 展示了甲烷在濃混合物逐漸生成、在稀混合物中快速消耗的過程.第二階段著火后，濃混合物側火焰前沿由于氧含量和其他活性基團不足，甲烷和正庚烷無法完全燃燒，正庚烷裂解反應生成甲烷.稀混合物傳播分支，氧氣含量增多，反應溫度逐漸升高，甲烷成分逐漸燃燒消耗.

圖7 不同時刻下溫度梯度、溫度峰值、羥基和甲烷質量分數隨混合分數的變化Fig.7 Conditional averaged temperature gradient，peak temperature，mass fraction of CH4 and OH against mixture fraction at different time instants

2.2 湍流場化學爆炸模式分析

化學爆炸模式分析(chemical explosion model analysis，CEMA)是由Lu 等[19]提出的一種定量檢測方法，可以用來識別臨界火焰特征(熄火極限和火焰前沿)、火焰和著火結構等，已在相關研究中廣泛應用[19-22].CEMA 基于化學源項雅可比矩陣 Jw的特征值確定混合物自燃傾向.

對反應系統離散守恒由源項和混合項構成，即

式中：y 為因變量向量，包括組分濃度和其他穩態變量，在CEMA 中，組分摩爾濃度和溫度包括在其中；為物質導數算子；w 為源項；S 為混合項.

對控制方程的兩邊取雅可比矩陣，主要包括化學源項的雅可比矩陣以及混合項的雅可比矩陣，即

式中：Ja為因變量向量物質導數的雅克比矩陣；Js代表混合項的雅克比矩陣.CEMA 分析中，通過化學雅可比矩陣確定臨界的火焰特征.通過求解雅可比矩陣特征值λe，并利用特征值分析化學反應特征[19,22-23]，其表達式為

進而可得CEMA 特征值Λe，且當Λe＞0 對應化學爆炸模式(CEM)，Λe＜0 對應非化學爆炸模式(Non-CEM)，Λe計算公式為

式中：be和αe分別是雅可比矩陣的左、右特征向量；eλ為雅可比矩陣的特征值；sgn()為符號函數；Re()表示取特征值的實數部分.

化學爆炸模式(CEM)下的反應物活性較高，自燃傾向大，且Λe越大表示反應物活性越高，自燃傾向越大，著火后燃燒劇烈，通常對應未燃反應區；非化學爆炸模式(Non-CEM)對應的區域反應活性較低，無自燃傾向.研究表明，CEM 向Non-CEM 的轉變過程與臨界火焰結構、著火以及預混火焰前沿位置密切相關，可用來確定著火和火焰前沿的位置.同時，Λe的零臨界點可用于區分已燃區域和未燃區域.筆者采用CEMA 分析混合層中自燃和火焰傳播特征，關于CEMA 方法的詳細介紹可參見文獻[21].CEMA 分析主要針對混合物中單點的組分濃度以及溫度，未考慮混合物中的擴散混合作用，但并不意味混合層中組分和熱量擴散作用不重要.

圖8～圖10 分別展示不同時刻下計算區域內Λe、O H 質量分數和溫度的分布情況.圖8a 所示t＝0.313 ms 時刻計算區域中央呈現窄條狀的青色區域(Λe＜0)，存在非化學爆炸模式(Non-CEM)，其形成主要與冷焰反應有關，主要集中在第一階段著火后濃混合物中緩慢發生低溫化學反應的區域.分析該時刻下OH 組分質量分數和溫度分布，如圖9a 和圖10a 所示，發現Non-CEM 區域右側區域溫度逐漸上升，形成多個高活性核心，但處于化學反應初期，化學反應活性仍很強，故對應的位置為化學爆炸模式(CEM).

圖8 不同時刻下計算區域Λe 空間分布Fig.8 Spatial distribution ofΛe in the computational domain at different times

圖9 不同時刻下計算區域OH空間分布Fig.9 Spatial distribution of OH concentration in the computational domain at different times

圖10 不同時刻下計算區域溫度空間分布Fig.10 Spatial distribution of temperature in the computational domain at different times

圖8b 所示t＝0.383 ms 時刻第二階段著火發生后，著火核心燃燒后發熱膨脹，形成高溫膨脹核心，因而離散的高溫膨脹核心著火后形成藍色(Λe＜0)的Non-CEM 區域.實際上，難以找到合適的OH 組分質量分數限值來確定反應前沿，在圖9b 中該時刻下OH 質量分數較高處和已燃區并非完全匹配，而使用CEMA 分析方法可以量化方式準確地確定火焰前沿和反應區域.

反應繼續進行至t＝0.403 ms 時刻，多個位置相繼進入第二階段著火時期，已存在的高溫膨脹核心膨脹連結，開始在空間上形成兩個獨立的區域.同時，正庚烷側下游存在冷焰形成的條狀Non-CEM 區域.高溫已燃區和冷焰區中間的中溫反應區仍然保持較高的自燃傾向.t＝0.430 ms 時刻低溫冷焰和高溫已燃Non-CEM 區域已經連結，繼續向兩側傳播.

3 結論

(1) 第一階段著火后，湍流作用下在ξMR,1等值線邊緣偏濃區域逐漸生成冷焰；第二階段著火后，在ξMR,2等值線偏甲烷側生成多個高溫的膨脹核心；高溫核心在混合分數平緩的區域更易生成；混合分數梯度過大會增大標量耗散率，加強熱量和活性基團的耗散，不利于燃燒反應的穩定進行.

(2) 不同于層流混合層中火焰前沿，湍流作用下火焰前沿形成褶皺向兩側傳播；位于熱膨脹核心位置的火焰前沿傳播較快，混合分數梯度大的區域反應進行緩慢，反應前沿傳播滯后；甲烷側氧含量高，前沿傳播加快，火焰前沿褶皺程度逐漸降低；正庚烷側火焰前沿在傳播下游存在冷焰反應區域，呈現“雙火焰”結構，冷焰位置低溫化學組分KET 大量生成，火焰前沿處溫度高，高溫化學反應組分OH 大量生成；隨著反應進行，火焰前沿傳播至更濃混合物中，雙火焰面之間的距離逐漸縮短.

(3) 由混合分數空間上溫度、組分和溫度梯度時間變化特征可知，甲烷側的化學反應強度更大，反應前沿逐漸向ξst位置傳播；另外，甲烷組分在混合分數空間上同時存在生成和消耗反應區，甲烷側的溫度高且氧氣充分，甲烷分子氧化消耗，而正庚烷側氧含量低，正庚烷無法完全氧化，部分裂解反應生成甲烷分子；CEMA 分析發現特征值可以準確地區分冷焰和高溫膨脹核心區，相較于表征組分等值線的方式，可以量化確定反應區、火焰前沿和雙燃料NTC 效應影響的中溫反應區.