帶有局部熄火現象的部分預混火焰大渦模擬研究

曾海翔， 王 平， SHROTRIYA Prshnt， 姜霖松 MURUGESAN Meentchidevi

(江蘇大學 a. 能源與動力工程學院; b. 能源研究院，江蘇 鎮江 212013)

部分預混燃燒的特點是混合氣在空間和時間上具有可變的當量比，且反應發生在若干子區域，因而同時兼具擴散火焰及預混火焰的特性[1].通過部分預混燃燒技術既可以降低燃燒過程中污染物的排放，又避免了預混火焰容易發生的回火問題[2-3].

早期工作中，對部分預混火焰主要通過實驗測量的方法加以研究.文獻[4]在高低不同雷諾數下的甲烷部分預混對沖火焰實驗中，揭示了火焰熄滅主要受大尺度的渦流控制.文獻[5]研究了正常和微重力條件對部分預混層流火焰的舉升特性及火焰結構的影響.文獻[6]在同心錐形燃燒器中產生丙烷射流部分預混湍流火焰，研究了伴流對火焰穩定性和火焰結構的影響.文獻[7]就進口速度和當量比振蕩對部分預混旋流火焰的擾動進行了實驗探究，結果表明通過改變噴油位置、噴油器阻抗和平均噴嘴速度，可以有效地控制燃燒室進口處的當量比擾動和進口速度之間的相位差，進而抑制或放大自激壓力振蕩的強度.文獻[8]采用粒子圖像測速法和平面激光誘導熒光(PLIF)技術，研究了稀薄部分預混湍流旋流火焰中流動與反應之間的非穩態相互作用，結果表明非定常流場的主導結構受剪切層中的進動渦核控制，并對火焰的穩定起著重要作用.

得益于計算機性能的提高和數值模擬技術的發展，對湍流火焰的數值模擬逐漸成為一種常用的手段.由于燃燒過程中的火焰鋒面區域發生強烈化學反應[9]，同時伴隨著流場的輸運和混合，這對火焰面的建模提出了極大的挑戰.文獻[10]利用進度變量、混合分數和標量耗散率的聯合概率密度函數(PDF)統計將非預混湍流燃燒模擬中層流火焰的概念推廣到部分預混區域，并通過湍流擴散舉升火焰和預混本生火焰的數值模擬證明了這一方法的廣泛適用性.文獻[11]將大渦模擬與條件矩封閉(CMC)方法進行耦合，通過計算湍流部分預混二甲醚射流火焰，證明了CMC方法在復雜化學條件下仍可以精確地預測湍流-化學之間的相互作用.文獻[12]在小火焰控制方程中通過引入火焰指數這一概念來識別不同區域的火焰狀態，建立了適用于湍流噴霧燃燒大渦模擬的部分預混小火焰模型.

文獻[13]首次提出“人工增厚火焰”(ATF)模型，多年以來，世界各地的專家和學者廣泛對其進行了論證和改進.文獻[14]通過引入效率函數解決了增厚后的火焰褶皺問題.文獻[15]考慮到需要增厚的火焰位置只發生在火焰的鋒面處，對于遠離火焰區的位置無需增厚處理，于是引進了動態增厚因子，這極大地拓展了增厚火焰模型的適用性.文獻[16]使用動態增厚火焰(DTF)模型結合小火焰生成流形(FGM)化學制表方法對旋流預混燃燒器進行模擬驗證，結果與實驗值吻合良好.文獻[17]為研究燃氣輪機內部的燃燒不穩定性與燃料混合狀態的關聯，分別采用完全預混和部分預混兩種假設，并利用DTF模型對PRECCINSTA燃燒器當量比為0.7和0.83兩種工況進行驗證；在部分預混方案中，針對部分預混火焰的特性對化學反應指前因子進行修正，進而調整火焰的傳播速度和厚度.結果顯示在當量比為0.7(接近貧燃極限)的工況下，采用部分預混假設方案捕捉到了自激不穩定現象和頻譜，而完全預混假設方案與實驗結果差距較大.

DTF模型在保留火焰層流燃燒速度的同時，將火焰面加厚到了網格可解的尺度，通過直接求解化學反應方程即可獲得火焰的詳細信息，避免了模型封閉的困難.理論上，DTF模型廣泛適用于各種燃燒模式，然而，DTF模型對部分預混燃火焰的預測還少有報道.這是由于傳統DTF模型在使用過程中需要預設增厚因子、層流火焰傳播速度及層流火焰厚度3個參數.在處理擴散火焰時，可以認為著火只發生在化學恰當比的對沖區域，故上述3個參數可以按照化學恰當比的層流火焰進行設置；當處理著火界限內的某個預混火焰時，同樣可以根據其具體的化學當量比來設置對應的參數.然而在部分預混火焰中，可燃物的化學當量比隨時間和流場位置不斷地演化，依靠單一化學當量比的建模方式已不太適合，此時需要對著火界限內的層流火焰燃燒特性進行動態修正.

另外一個問題是對部分預混火焰的火焰位置的捕捉.在之前的研究中，通過監測CO2的生成反應[18]或是基于CO2的質量分數構建進度變量[19-20]所捕捉的火焰面都不甚理想.而由文獻[21]提出的基于熱放率捕捉火焰面的方法表現良好，故本文所使用的DTF模型同樣利用反應放熱率方法來識別火焰位置.

由于常溫常壓下，火焰的層流燃燒速度隨當量比波動較為明顯，鑒于此，所使用的DTF模型通過判斷當地可燃物的混合分數，進而自動調整層流火焰的燃燒速度和厚度.最后，利用動態k方程亞格子模型，結合DTF模型和甲烷多步化學反應機理，對悉尼部分預混燃燒裝置下FJ200-5GP-Lr75-103接近吹熄極限的工況進行了模擬驗證，下文中此工況根據其速度大小簡記為U103，證明了DTF模型在計算部分預混火焰方面具有足夠的適應性和準確性.

1 湍流燃燒模型

1.1 湍流模型

對三維納維-斯托克斯方程進行Favre濾波得到的大渦模擬控制方程參考文獻[21]，由濾波帶來的亞網格應力項及亞網格標量通量項需要模型進行處理，其中本文計算中采用文獻[22]提到的動態k方程模型對亞格子應力項加以封閉：

(1)

1.2 動態增厚火焰燃燒模型

文獻[13]考察了層流火焰的燃燒速度和火焰面厚度，發現其均勻與組分的擴散系數和化學反應有某種關聯：

(2)

(3)

文獻[14]認為增厚火焰模型降低了火焰對流場的敏感性，于是引入褶皺因子E的概念.針對部分預混火焰的特性，首先基于GRII-Mech3.0機理得到了可燃界限內甲烷常溫常壓下，層流燃燒速度sL和層流火焰面厚度δL隨當量比φ的變化情況，如圖1所示.隨后，通過MATLAB提出二者隨當量比變化規律的擬合函數，并植入DTF模型中.在所使用的DTF模型中，火焰面的位置通過多步化學反應的放熱量dq來捕捉，使得火焰的增厚系數只在燃燒火焰區內有效，修正后的火焰增厚因子Fdyn和褶皺因子Edyn分別由下式表示：

(4)

式中：ndtf為預設的增厚因子；Δe為網格尺度；a1、a2為調整流場中動態增厚曲線形態的縮放因子，此處a1=10,a2=25.褶皺因子可由下式表示：

(5)

(6)

圖1 層流火焰的燃燒特性Fig.1 Combustion characteristics of laminar flame

2 算例設置

2.1 悉尼燃燒器

為研究入口處燃料的混合情況對火焰的影響，文獻[23]精心設計了一套燃燒器，結構如圖2(a)所示，其中：Lr為回縮距離.燃燒器的入流管道是由兩個壁厚為0.25 mm的同心管組成，內管通入甲烷氣體，內管直徑為4 mm，外管通入空氣，外管直徑為7.5 mm.混合氣體通過外圍導管中的高溫燃燒廢氣點燃，廢氣導管的直徑為18 mm，整個裝置放置于25 mm × 25 mm的伴流風洞中.燃燒器具有凹入的內管，該內管可以從外管出口平面向上游回縮.當Lr=0時，內管與外管平齊(內管位于外管出口平面)，此時混合氣燃燒模式為完全非預混，即擴散燃燒模式；當Lr≥300 mm(40倍的外管直徑)時，此時燃氣與空氣充分混合，燃燒模式為完全預混燃燒模式.當Lr處于兩個極限距離之間時，燃燒模式為部分預混燃燒模式.為驗證DTF模型對部分預混火焰的預測能力，計算內管回縮距離選取 75 mm，對應混合氣體處于部分預混狀態，速度入流參數如表1所示，其中：U為內外管混合段的混合氣平均速度；UA為外管入流段空氣的平均速度；UF為內管入流段燃料的平均速度.為了點燃射流火焰，在射流火焰的外圍制造了一個環形的值班火焰區域，其中注入的是由C2H2、H2、CO2、9N2和空氣這5種組分按照甲烷-空氣火焰在化學恰當比時的元素配比來配置的混合氣.計算過程中值班火焰入口取為實驗測量結果，具體信息如表2所示.其中：w(N2)、w(CO2)、w(H2O)、w(O2)和w(CO)分別為對應組分的質量分數；Up為值班火焰環空內混合氣的平均速度；T為值班火焰峰值溫度.

圖2 悉尼部分預混燃燒器及計算網格示意圖[23]Fig.2 Schematic diagram of Sydney partially premixed burner and calculation grid[23]

表1 U103工況下入流速度信息Tab.1 Inflow speed information under U103 condition

表2 U103工況下值班火焰尾氣中各組分質量分數、流速和溫度信息

2.2 計算設置

本文計算基于開源的計算流體力學軟件OpenFOAM，其向來以其代碼的高效率和優良的可拓展性著稱.在原有的rhoReactingFoam 求解器基礎上，本課題組植入了DTF模型.計算采用一套分塊的包含3.82×106個單元的結構網格，如圖2(b)所示.其中，徑向、周向、軸向分別部署了76×80×628個網格單元，火焰區域內的網格最大寬度為0.2 mm.

化學反應機理對模擬的精度有著明顯的影響，文獻[24]考量了5種總包機理與文獻[25]提出的骨架機理在貧燃部分預混旋流火焰中的表現，所獲得的結論為，總包機理可以大大降低計算成本，但同時也影響結果的質量；而包含73步化學反應的Lu骨架機理表現與實驗值吻合的最好，但計算成本也遠高于總包機理.本文中化學反應機理選用文獻[26]提出的甲烷19步骨架反應機理，該機理同樣是由GRI-Mech1.2機理簡化而來.計算中燃料和空氣入口施加強度為10%的隨機擾動，時間離散為二階Backward格式，對流項采用二階有界線性格式，拉普拉斯項選用高斯線性正交.

3 結果與討論

流場在不同軸向位置處的流場信息如圖3所示.其中：r為軸向距離；x為徑向距離；d為鈍體直徑；Ux為徑向平均速度；Ux_rms為徑向均方根速度.從圖3中可以看出，可燃混合氣剛進入燃燒室時，中心軸向處的速度較高，并向外逐漸衰減，混合氣與周圍伴流氣體接觸的剪切層中的速度梯度很大.隨著混合氣沿著下游方向不斷的發展，主流的平均速度不斷降低，而周圍流體的速度卻在不斷增大.這是由于剪切層的作用，使得速度較高的可燃氣不斷地將周圍的空氣卷入主流中，以降低剪切層中的速度梯度，而周圍流體由此受益使其平均速度逐漸增大.觀察流場中的均方根脈動，可以看到在整個流場中，均方根脈動的峰值幾乎沒有發生變化，但脈動的頻譜卻越來越寬，整個流場越來越趨于一致.

圖3 不同軸向位置處的速度場信息Fig.3 Velocity field information at different axial positions

為了辨別流場中的渦旋結構與位置，利用無量綱化的Q準則和渦量ω提取流場中的擬序結構，如圖4所示.Q準則云圖中可以看到在火焰的根部r/d約為0.5的位置處達到峰值，然后向下游緩緩的擴張開來，在x/d=5到x/d=20的區間內渦旋較為密集.渦量圖顯示了流體微團的旋轉規律.從圖4中可以看到，沿軸向的渦旋幾乎總是與反向渦旋成對出現，在x/d>20的下游位置渦量逐漸減弱.對于不同軸向位置，其渦旋峰值所在的位置正好與圖4流場中速度沿x方向的均方根波動峰位相吻合，這進一步表明由剪切層產生的渦旋加劇了流場的摻混和軸向速度的傳遞.

以二氧化碳的瞬時生成速率云圖RCO2及火焰指數ξ來分析火焰面的信息，兩種方案的結果如圖5所示.其中：虛線為RCO2=3 kg/(m3·s) 的等值線，用來勾勒火焰面的位置.由圖5可知，由CO2生成的速率云圖中不同時刻的火焰面的底部x/d<10的部分均沒有產生中斷，這意味著火焰底部附近的燃燒較為穩定.這是由于外部的值班火焰(根部左、右兩支角狀火焰分支)持續作用的結果.由于過量的甲烷和空氣在導管內來不及充分預混便進入到燃燒室，火焰根部的混合分數較高，只有在外層較為稀薄的地方被伴流火焰所點燃，所以根部附近的火焰面較窄，通過分析火焰指數云圖可以看到，兩種方案在此區域均為預混燃燒模式.然而在火焰的內部區域，混合氣體此時由擴散模式主導.在火焰中部1020)火焰繼續向下游發展，孔洞出現的越來越多，其范圍也逐漸增大，局部熄火現象進一步加劇，這一區域的火焰由擴散模式所主導.

圖4 燃燒室內無量綱化的Q準則及渦量ω圖Fig.4 Dimensionless Q criterion and vorticity ω in combustion chamber

圖5 CO2瞬時生成速率云圖和火焰指數云圖Fig.5 Cloud maps of CO2 instantaneous generation rate and flame index

圖6 不同軸向位置處混合分數均值及均方根波動值對比Fig.6 Comparison of mean mixture fraction and root mean square fluctuation at different axial positions

可燃氣的混合狀態決定了火焰的溫度及燃燒產物，故而計算中必須對混合分數做出精準的預測.圖6顯示了不同軸向位置處混合物的質量分數對比，其中：Zmean為混合分數均值；ZRMS為混合分數均方根波動值.由圖6可以看到，無論是火焰的根部亦或是下游區域，采用動態k方程湍流模型及結構化的382萬個網格對可燃氣的混合分數的計算結果與實驗值吻合良好.在整個流場中，混合氣的質量分數呈現出了較好的對稱性，在r/d=0的中心線上，混合氣的質量分數達到最大，然后向兩邊逐漸遞減.然而混合分數的均方根脈動峰值卻沒有出現在中心線上，而是在中心線附近并沿中心線對稱，呈現出兩個峰狀.隨著高度的增加，混合氣的質量分數峰值逐漸降低，火焰面隨之也變得越來越寬.這是由于邊緣較為稀薄的可燃氣不斷的被消耗，混合氣從濃度高的中心區域不斷向兩邊擴散，整個混合氣的濃度逐漸被稀釋的結果.

圖7為計算結果與實驗值的溫度對比，其中：Tmean為溫度均值；TRMS為溫度均方根波動值.整體上看，改進后的DTF模型計算結果較為良好的展現出了火焰中的溫度分布及波動趨勢.當x/d<10時，平均溫度峰值出現在r/d=1附近，而均方根波動的峰值出現在r/d=0.5附近；隨著下游熄火位點的增多，溫度場的均值不斷下降，而其波動也逐漸蔓延到外側區域.

圖7 不同軸向位置處的平均溫度及均方根波動值對比Fig.7 Comparison of mean value and root mean square fluctuation of temperature at different axial positions

文獻[2]在層流部分預混射流火焰的直接數值模擬中指出，部分預混燃燒火焰中包含兩個反應層，內層中的甲烷和氧氣不斷消耗產生CO、H2、OH等中間產物，然后在外層的氧化層中被反應為CO2和H2O.這里選擇CO作為中間產物的代表，并與實驗值進行對比，結果如圖8所示，其中:w(CO)mean為一氧化碳質量分數均值;w(CO)RMS為一氧化碳質量分數均方根波動值.對于中上游中間產物的計算結果，DTF模型結果較好的展現了中間產物CO的含量及波動特征.同時可以看到在火焰的下游區域，計算結果較實驗值出現了一定的偏差，此區域的CO質量分數均方根波動也變得較為劇烈，這可能是由于采用簡化機理無法較為全面地描述CO的反應過程.從圖8左圖可以看到，不同平面高度上CO的峰值均出現在r/d=0.75附近，然后向兩邊緩慢擴散；到了下游區域，曲線逐漸變得平滑.與圖7中溫度曲線的峰值位置進行對比，可以發現放熱峰值出現在r/d=0.9附近，稍遠于CO峰值面，這表明Suresh的假設是比較合理的.由于避免了環境因素的影響，計算值比實驗值呈現出了更好的對稱性.

將混合氣的溫度隨質量分數分布制成散點圖如圖9所示，其中：T為溫度；Z為混合分數；藍色虛線為化學恰當比所在線，將火焰分割為左、右兩支.左半支的混合分數較低，而對應的溫度較高，故其為高溫燃燒廢氣與伴流空氣的混合支，右側的混合分數較高但對應的溫度比較低，對應可燃氣與高溫燃燒廢氣的混合支.由于入流速度過高(接近吹熄極限)，導致火焰根部x/d=1附近大量的可燃氣體來不及引燃便被帶到下游.進一步在x/d=5的位置觀察火焰燃燒狀態，可以發現此時混合氣的燃燒狀態較為理想.從x/d>10，左、右兩個混合分支間開始出現熄火點，越到下游熄火點出現的越為密集.由圖3中的火焰指數云圖可以看出，在流場的中下游，擴散燃燒模式逐漸占據主導地位.對比實驗和計算的散點圖可以發現，采用改進后的動態層流火焰燃燒速度及層流火焰厚度的DTF模型整體非常良好的捕捉到了這些變化趨勢.散點圖中火焰的根部x/d=1處，無論是實驗亦或是計算模擬，散點數據都比較少；且大部分散點都集中在化學當量比分割線附近.當混合分數Z超過0.1，300 K附近出現一條水平直線.這是由于混合氣剛進入燃燒室時，底部主流區域的可燃氣混合分數過高且流速過快，來不及與周圍的燃燒廢氣進行熱量交換，這部分散點全部集中在300 K的帶狀分布附近.由于這部分的散點比例過大，導致整個散點圖中的點顯得比較稀少.

圖8 不同軸向位置處混合物的CO質量分數平均值及均方根波動值對比Fig.8 Comparison of mean mass fraction and root mean square fluctuation of CO at different axial positions

圖9 不同軸向位置處混合物的溫度隨混合分數分布散點圖Fig.9 Scatter plot of temperature with mixture fraction distribution at different axial positions

4 結論

為了應用動態增厚火焰模型計算當量比不均勻分布的湍流部分預混火焰，需要考慮到當量比變化對火焰褶皺函數中出現的層流火焰速度、厚度的影響.本文通過計算一系列的層流火焰構造了兩個擬合函數，其可以在DTF模型中自適應調整火焰參數，很好地解決了上述問題.結合動態k方程亞格子湍流模型，對接近吹熄極限的悉尼部分預混火焰進行了大渦模擬研究.本文研究所獲得的結論如下.

(1) 通過流場中混合分數均值及均方根脈動對比，顯示了動態k方程亞格子湍流模型可以很好地預測非均勻入流的射流流動問題.

(2) 觀察CO2生成速率云圖和OH基渲染的三維瞬時火焰面，不難發現采用改進后的DTF模型對此類帶有熄火現象的部分預混火焰具有較好的預測能力，溫度-混合分數散點圖分布也很好的說明了這一點.

(3) 本文使用的甲烷多步簡化機理對某些組分如CO的質量分數預測較實驗值偏低，若采用更為詳細的化學反應機理或結合化學制表方法，這種現象可能會有一些改善.