基于網格化Damk?hler 數的燃燒室貧熄邊界預測方法探究

王中豪 ，張軍華 ，鄧愛明，胡 斌 ，趙慶軍,

(1.中國科學院工程熱物理研究所，北京 100190；2.中國科學院大學航空宇航學院，北京 100190；3.中國科學院輕型動力重點實驗室，北京 100190)

因貧油熄火導致的飛行安全問題一直被研究者所重視，在航發燃燒室設計階段，需要貧熄性能評估和燃燒室結構改進的多次循環，而發展一套能夠高效、高精度預測航發燃燒室貧熄邊界的方法，對于縮短燃燒室設計周期，降低設計成本具有重要意義[1].

從早期Longwell 等人提出均勻攪拌反應器模型[2]以及Ballal 等人提出回流區點火模型開始[3]，基于能量平衡和時間平衡的半經驗模型就在大量貧熄實驗的基礎上不斷被發展完善[4-9].隨著計算機技術的發展，大渦模擬方法被廣泛應用于分析熄火前的火焰行為，部分研究者將大渦模擬方法直接應用于燃燒室的貧熄邊界預測.通過大渦模擬方法對燃燒室貧熄邊界進行預測，一般是通過再現貧油熄火的瞬態過程實現的，而基于不同的燃燒類型及運行工況恰當地選擇燃燒模型，對于能否直接通過數值模擬準確再現貧熄過程有決定性作用.

Kim 等[10]應用渦破碎模型針對鈍體穩燃預混火焰進行大渦模擬，在當量比0.52 時模擬出熄火，但是在油氣比為化學恰當比時模擬的鈍體后回流區長度卻與實驗有一定偏差，Kim 等將此偏差歸因于該工況下發生的強烈的非平衡化學反應.Gokulakrishnan等[11]采用渦耗散概念模型和層流燃燒模型對鈍體穩燃預混火焰進行了大渦模擬，在當量比為0.45 時，采用渦耗散概念模型能夠模擬出全局熄火，與實驗觀察到的結果基本一致，而采用層流燃燒模型的模擬中燃燒依然維持穩定，與實驗結果不符.Gokulakrishnan等人認為這是由于渦耗散概念模型假定燃燒發生在微尺度渦團中，而層流燃燒模型在亞格子尺度下假定標量不發生脈動，前者更好地模擬了湍流和化學反應的相互作用.Black 等[12]采用概率密度輸運模型對低排放旋流霧化燃燒進行了大渦模擬.在實驗獲得的貧熄當量比下，數值模擬結果中燃燒依然穩定，未能獲得預期的貧熄結果.Black 等認為是流場中局部富油而產生的高溫區影響了預測結果的準確性，強調了精確模擬燃油軌跡及蒸發對于霧化燃燒貧熄預測的重要性.此外還有Garmory 等[13]和Zhang 等[14]采用條件矩方法對Sandia D&F 火焰進行的大渦模擬.模擬能夠再現近貧熄工況頻繁發生的局部熄火和再點火過程，Zhang 等[14]通過數值模擬獲得了不同來流流速下的貧熄邊界，與實驗結果誤差在25%～60%.Menon 等[15]采用線性渦模型對旋流預混火焰進行了大渦模擬.模擬結果展示出了火焰拉伸對火焰傳播的影響，但是由于一些尺度的渦沒有準確求解，未能捕捉到準確的貧熄邊界.

在當前計算機發展水平下，盡管能夠對燃燒室中的流場進行大渦模擬，其計算成本依然很高，并不十分適用于工程計算.在Black 等人的工作中[12]，盡管只采用了一步總包反應，每一個工況仍然需要在16 900 MHz 的計算機上計算20 d.考慮到大渦模擬仍然需要較高昂的計算成本，將貧油熄火半經驗模型與雷諾平均數方法結合，直接獲得不同工況下的貧油熄火邊界，不失為一種高效的貧熄預測方法，也有部分研究者在此方向上做出了一些嘗試.

數值模擬與半經驗公式相結合的方法主要可分3 類.第一類是將燃燒室內的流場劃分為不同的子區域，將流場簡化為由不同反應器聯接成的網絡，并對每個反應器求解詳細的化學反應過程，從而獲得熄火邊界.比較典型的有 Sturgess，Mongia 等人的工作.在Sturgess 等[16]的工作中，在基于整個燃燒室內部流動建立的反應器網絡中定義一個關鍵反應器，當此反應器溫度低于1 480 K 時認為發生貧熄.Mongia等[17-18]則是在反應器網絡中的每個反應器中應用Lefebvre 的貧熄預測公式[6]，并對每一反應器中的計算結果進行體積加權，判斷是否發生全局熄火.第二類方法將局部熄火模型與數值模擬相結合，首先預測局部熄火，再通過局部熄火發展的程度判斷是否發生全局熄火.比較典型的是Knaus 等的工作[19].Knaus等通過Damk?hler(Da)數場判斷局部熄火，當鈍體下游的局部熄火區域連成整體時，判斷發生全局熄火.此方法的局限是對局部熄火如何引發全局熄火的解釋比較粗略，沒能明確指出流場中哪些區域是控制全局貧熄的關鍵區域.第三類方法是在數值模擬的流場中先定義一個關鍵反應區，將該區域作為控制熄火的關鍵區域，并在該反應區中應用半經驗模型進行貧熄預測.比較典型的工作有胡斌等人基于燃油濃度分布定義燃燒體積并結合優化后的Lefebvre 半經驗模型發展的貧熄預測方法[20-22]；Zheng 等[23]基于速度分布定義特征面結合特征面溫度變化發展的貧熄預測方法；王中豪等[24-25]發展了基于燃燒室內OH基分布結合Da 數的貧熄預測方法.本研究發展了一種基于網格化Da 數的貧熄預測方法，將基于OH 基濃度定義的反應區均勻劃分為尺度一致的子區域，并基于子區域計算局部和整體Da 數.該方法可以降低反應區組分、溫度分布不均勻對預測結果的影響；同時基于局部Da 數場所描述的反應區燃燒組織情況，有助于進一步掌握近貧熄狀態下燃燒室內流動和化學反應的匹配情況，為燃燒室設計提供更加完善的理論支撐.

1 理論介紹

1.1 關鍵反應區

本研究中關鍵反應區是基于OH 基濃度定義的，當OH 基摩爾分數高于0.001 時，認為該區域為關鍵反應區，該邊界值是通過將數值模擬獲得的OH 基分布和可視化實驗[26]中由高速相機獲得的反應區照片定性對比獲得(實驗中的反應區定義為燃燒發光部分).圖1(a)為實驗中獲得的火焰區照片，圖1(b)為數值模擬中基于OH 基定義的關鍵反應區，將OH 摩爾分數0.001 作為邊界時，二者在形態和大小上呈現出較高的相似度.

圖1 實驗和數值模擬中的反應區[26]Fig.1 Reaction zone in the experiment and numerical simulation[26]

基于OH 基濃度定義反應區，綜合考慮了湍流、化學反應及其相互作用對貧熄的影響.OH 基在燃燒過程中產生并在燃燒結束前消耗，其存在標志著反應的劇烈進行[27].根據鏈式反應假設，自由基的消滅也是反應終結的原因[28].湍流導致的反應物濃度脈動對燃燒的影響，也可以在OH 濃度分布中得到體現[29].因此將OH 基集中分布的區域作為關鍵反應區，能夠獲得反映全場燃燒狀態的較為客觀的信息，從而對貧熄作出較為準確的預測.

1.2 反應區分區及分區數量無關性驗證

根據Longwell 等[2]提出的均勻攪拌反應器理論，整個反應區可以被簡化為一溫度和組分濃度分布處處相等的反應器，該假設將反應區簡化為0 維反應器，大大降低了在反應區求解詳細的化學反應過程的計算成本，但是實際上反應區內溫度與組分分布仍然存在較大的非均勻性，在整個反應區直接應用均勻攪拌反應器假設，精確度難以得到保障.本研究通過將反應區進行分割來最大程度降低反應區內溫度和組分分布的不均勻性帶來的影響，形成長度尺度相同的子區域，再針對每個子區域進行均勻攪拌反應器假設.分區方法如圖2 所示：①統計圖2(a)所示燃燒區在x、y、z 方向的坐標范圍，依據該范圍建立一長方體如圖2(b)，其長、寬、高對應燃燒區在x、y、z 方向的跨度，此步驟相當于將不規則的燃燒區投射為規則的長方體；②對該長方體的長、寬、高進行等數目分割(例如，若長方體的每邊都分割為5 段，如圖2(c)，則長方體共被分割為53＝125 塊)；③分割后的長方體投射回不規則燃燒區，每一小長方體對應燃燒區中的一個子區域，如圖2(d).

圖2 反應區分割過程示意Fig.2 Schematic of reaction zone blocking process

子區域的長度尺度越小，其內部的溫度和組分分布就越趨近于均勻，但是考慮到計算效率等因素，不可能對反應區進行無限細分，需要通過分區無關性驗證工作，確保分區數目能夠達到計算效率和精度的最佳平衡.本研究中選取來流空氣流量1.2 kg/s 的某一鄰近貧熄工況(燃油流量略高于貧熄工況)作為驗證工況，計算不同分區數目下的整體Da 數，從1 開始，逐漸增大分區數目，當整體Da 數不再隨分區數目增大而呈趨勢性變化時，認為分區數目不再影響整體Da 數的計算.圖3 為分區無關性驗證結果，可以看出當分區數目從1 增加至216 時，Da 數單調減小，但是減小速度開始變慢，當分區數目大于216 時，Da數不再單調變化，而是呈小幅波動.因此在本研究中認為將反應區分為216 塊，能夠在保障計算精度的前提下獲得最高計算效率.

圖3 整體Da 數隨分塊數目的改變Fig.3 Variation in global Da number with the number of blocks

1.3 時間尺度

在本研究中Da 數被用于貧熄分析及預測.Da數定義為流動時間尺度(τf)與化學時間尺度(τc)的比值.該模型本身沒有修正因子，因此具有在寬工況下進行貧熄預測的潛力.

1.3.1 流動時間尺度

流動時間尺度定義為新鮮混氣在關鍵反應區的實際停留時間[24-25].該時間通過對數值模擬中的離散項顆粒進行粒子追蹤獲得，其過程如圖4 所示.考慮到實際的燃油顆粒會蒸發為氣態，設置了一個虛擬噴嘴在燃油注射點噴射低密度惰性粒子，流量在10-5kg/s 量級，以保障惰性粒子射流對燃燒不產生影響.通過追蹤粒子到達和流出某區域(反應區或其子區域)的對應時刻，確定其在該區域的實際停留時間.數值模擬中惰性粒子軌跡數量為100 條，根據每條軌跡獲得的流動時間結果，剔除明顯的問題數據之后取平均值，以消除粒子運動軌跡隨機性帶來的影響.

圖4 流動時間尺度獲得過程示意Fig.4 Schematic of obtaining the flow time scale

1.3.2 化學時間尺度

化學時間尺度定義為新鮮來流混氣能夠觸發目標區域內反應的最短停留時間[24-25]，通過商業軟件包CANTERA 中的均勻攪拌反應器模塊計算.化學時間尺度的求解方法具體為：

(1) 將均勻攪拌器初始溫度設置為子區域體平均溫度，以模擬高溫燃燒產物對新鮮來流混氣的點燃作用.

(2) 反應器的來流當量比由子區域內的混合分數反推獲得，其溫度等于燃燒室的入口溫度，模擬新鮮混氣從冷態到被加熱至反應溫度的過程.

(3) 反應器的壓力設置為子區域內的體平均壓力.

(4) 先設置較大的反應物停留時間，確保能夠成功點火，再逐步減小停留時間直至點火失敗，此時的停留時間即化學時間尺度.

采用瞬態求解器時，不同停留時間下的燃燒狀態如圖5 所示.圖5(a)為穩定燃燒狀態，反應器內溫度維持在比較高的溫度(該算例下為1 500 K 左右)，持續點燃新鮮混氣；圖5(b)為熄火狀態，反應器內溫度由初始溫度降低到來流溫度(495 K)，表明反應器內發生全局熄火.

圖5 點火成功和點火失敗時反應器溫度隨時間的變化Fig.5 Variations in reactor temperature with time under successful and failing ignition conditions

2 實驗設置

針對扇形燃燒室模型開展了貧熄性能實驗，該實驗在中科院輕型動力重點實驗室進行.該扇形燃燒室模型含三頭部，為環形燃燒室的1/4，每頭部裝有雙徑向旋流器，第一級旋流數為1.054 6，第二級旋流數為0.947 1.火焰筒壁面采用氣膜冷卻，火焰筒兩側壁各有220 個冷卻孔，防止側壁因超溫而被破壞，由于從冷卻孔流失的空氣并未實際參與燃燒，因此在計算熄火油氣比時被排除.模型燃燒室設計圖及實驗件照片如圖6 所示，圖6(a)為模型燃燒室設計圖，圖6(b)為模型燃燒室實驗件照片.燃料采用國產RP-3航空煤油，其物性參數如表1 所示.

圖6 實驗件設計圖及實驗現場照片Fig.6 Schematic and photo of test rig

表1 國產RP-3航空煤油物性Tab.1 Physical properties of Chinese RP-3 aviation kerosene

在不同空氣來流流量下獲得了模型燃燒室的貧熄邊界，不同工況的來流溫度均為495 K，來流壓力為0.23 MPa，不同來流流量為0.6 kg/s、0.7 kg/s、0.8 kg/s、0.9 kg/s、1.0 kg/s、1.1 kg/s 和1.2 kg/s.在每一運行工況下，先提供充足的燃油流量維持穩定的燃燒，然后以微小調節量梯度減少燃油供給，在每次調節燃油流量后先維持該狀態2 min，再進行下一次調節，防止因燃油流量突變引起的擾動影響獲得貧熄邊界的真實數值.通過火焰筒出口下游處熱電偶監控出口溫度，當溫度突然降低且穩定在燃燒室入口溫度(T3)附近，認為發生全局熄火.實驗反復進行2 次，以確保獲得的貧熄數據的可靠性.

3 數值模擬方法

3.1 湍流及燃燒模型

本研究通過商業軟件FLUENT 進行數值模擬，采用可實現的k-ε 方法模擬湍流流場，采用小火焰模型對燃燒進行?；?，燃燒機理采用Kundu 等[30]的煤油燃燒16 組分23 步反應機理.該機理簡化了煤油燃燒時的熱解過程，摒棄了燃燒過程中的一些冗余組分和反應，提升了計算效率，在當量比低于1.4 的燃燒反應模擬中表現尤為突出.該機理也在其他研究中被廣泛用于組分場和溫度場的獲得，其可靠性得到了多次驗證[31-32].對于湍流k-ε 方法舍棄了部分雷諾應力的數學限制，使之更符合湍流流動的真實物理現象，在模擬包含旋流和回流的流動時表現良好.小火焰模型將湍流火焰簡化為鑲嵌在湍流流場中的多個層流對沖火焰的集合，將流場內組分、溫度的分布作為混合分數、混合分數脈動以及標量耗散率的函數.由于采用查表插值的方法進行化學源項相關計算，小火焰模型耦合詳細機理進行化學計算的效率很高.數值方法驗證工作可參閱文獻[24]，限于篇幅此處不再贅述.

3.2 計算域及網格

采用循環邊界對單頭部燃燒室進行模擬，計算結果反映在全環燃燒室中每一頭部的燃燒情況.計算域如圖7 所示，機匣等與燃燒室內部流場無關的結構被簡化.

圖7 計算域及來流條件Fig.7 Computational domain and inlet conditions

網格劃分情況如圖8.考慮到該模型燃燒室結構比較復雜，除在旋流器的兩級旋流通道、入口擴壓段生成結構化網格，其他部分均生成非結構網格，網格總數約為350 萬.

圖8 模型燃燒室網格Fig.8 Mesh of model combustor

4 結果與討論

4.1 貧熄實驗結果及分析

圖9 為不同來流流量下的熄火油氣比.由圖9可見，隨來流空氣流量增加，貧熄油氣比幾乎成線性增加，說明來流流量的增大會使燃燒室貧熄性能惡化.空氣流量的增大使反應區的熱損失增大，不利于燃油蒸發和燃燒.另外，空氣流量增大，燃燒室內氣流速度增大，導致火焰被扭曲得更為強烈，反應區更容易破碎，增大的流動速度也需要更高的火焰傳播速度與之匹配，使火焰更難穩定.綜上所述，貧熄油氣比隨來流流量的增大而增大.

圖9 貧熄油氣比隨流量的變化Fig.9 Variation in fuel-air ratio with maunder lean blow-off conditions

4.2 局部Da 數計算結果及分析

圖10 為不同工況下的Da 數場，圖10(a)～(g)均為近貧熄工況，圖10(h)為設計工況.圖中白色虛線圈起的部分為反應區，燃燒室內主要的燃燒反應發生在該區域，紅色區域為Da 數較高的區域，可以視作反應最為集中的燃燒核心.從圖中可見，在近貧熄工況，盡管隨著來流流量的增大，紅色區域所占空間及大小均有波動，但是均集中在回流區上游滯止點附近.在回流區上游滯止點，旋流引起的回流與上游來流對撞發生滯止，此區域附近的空氣流速相對較慢，利于火焰穩定.同時，上游來流帶來新鮮混氣，提供充足的反應物，高溫燃燒產物通過回流傳遞到此處，提供觸發反應需要的能量，因此該區域為反應最為集中的區域.從圖10(a)～(g)可見，來流流量逐漸增大，燃燒區逐漸擴大，流量大于1 kg/s 后，其近燃燒區尾端表面逐漸由內凹變為外凸.燃燒區形態的改變主要是由貧熄油氣比的上升導致，反應區的擴張有利于增大可燃混氣的停留時間從而維持燃燒穩定，但是由于圖10(a)～(g)均為不同工況下的近貧熄狀態，與反應區體積同時改變的還有反應區內的溫度、壓力、燃料分布，多種因素作用的效果相互耦合，使反應區擴大帶來的有利穩燃的因素被抵消.

圖10(h)顯示的是設計工況下的Da 數場.從圖中可以看出，紅色區域占據了反應區大部，表明空間上廣泛分布的劇烈的化學反應.由于燃油供給充分，反應區甚至擴展到主燃孔下游，大大增加了反應區內可燃混氣的停留時間.

圖10 不同工況下的Da 數場Fig.10 Da fields under different conditions

4.3 不同油氣參數下Da 數場分析

圖11 為來流空氣流量0.6 kg/s 時，油氣比從5.27 g/kg 逐漸增加至19.2 g/kg 的Da 數場圖，白色虛線圈起的部分為反應區.從圖11(a)和(b)可以看出隨著油氣比從5.27 g/kg 增加至9.91 g/kg，反應區顯著擴張，形狀由內凹變為外凸，但是反應區內平均Da 數沒有明顯增加，也并沒有在局部出現高Da 數的劇烈反應，表明燃油增加的初步效果是拓展燃燒反應發生的范圍，但是燃燒強度不會大幅度提升；圖11(c)表明，當油氣比提升至14.56 g/kg 時，反應區進一步擴大的同時開始出現部分表明高Da 數的紅色區域，表明進一步增加燃油會導致燃燒區域的增大和燃燒強度的提升；當油氣比增大至19.2 g/kg 時，反應區繼續小幅擴張，同時高Da 數區域范圍顯著擴大，表明此時增油帶來的主要是燃燒強度的提升.

圖11 不同油氣比下反應區Da 數分布變化Fig.11 Variation in Da fields in the reaction zone at different fuel-air ratios

4.4 整體Da 數計算結果及分析

圖12 為不同來流工況下的時間尺度，圖12(a)為整體流動時間，圖12(b)為整體化學反應時間.整體流動時間尺度通過追蹤惰性粒子在整個反應區的停留時間獲得，整體化學反應時間為子區域化學反應時間的體積加權平均值.從圖12(a)和(b)可以看出，整體流動時間和整體化學反應時間都隨來流空氣流量的增大而整體呈現上升趨勢.整體流動時間主要決定于反應區的體積和反應區內氣流的平均流速.來流空氣流量的增大及燃油流量的增大，拓展了反應區的范圍，增大了反應區體積，而反應區體積的增大會增大混氣在燃燒區的停留時間；與此同時，來流流量的增大直接增大了燃燒室內的氣流流速，又會減少混氣在燃燒區的停留時間.在本研究中，隨空氣流量的增大，燃燒區體積的增大對流動時間尺度的影響占據主導，但是燃燒區體積的增大會因為受到諸多因素限制而漸漸減慢，氣流速度卻與來流流量呈一定的線性關系，因此流動時間尺度的增加漸漸減慢.化學時間尺度受到各子區域體平均溫度、壓力、當量比的影響.體平均溫度、壓力、當量比的提升，都有利于反應的進行，減小化學反應時間.隨著來流流量的增大，燃燒區內整體體平均溫度、當量比都呈上升趨勢，但是由于壓力損失的增大，體平均壓力呈下降趨勢，最終化學時間受到體平均壓力的主導，隨來流流量的增大而增大.

圖12 近貧熄工況整體流動時間尺度和整體化學反應時間尺度Fig.12 Global flow time scales and chemical reaction time scales under different lean blow-off conditions

圖13 為近貧熄工況整體Da 數的計算結果，整體Da 數將整個反應區的燃燒狀況用一個Da 數表示，Da 數越大表示燃燒進行得越穩定.圖13 所統計的不同來流流量下的Da 數，平均值在7.9 左右，并圍繞該值上下輕微波動(平均波動范圍在5%左右，最大波動不超過14%).這表明由整體Da 數刻畫的貧熄工況對于來流流量的變化不敏感，通過整體Da數判斷貧熄具有對不同來流工況的普適性.來流工況的變化會引起燃燒室內溫度、組分場變化，而本文中的流動時間尺度和化學反應時間尺度均基于流場獲得，因此二者可以在一定程度上反映出由于來流工況變化引起的流場變化，從而使得全局Da 數能夠在不同來流流量的近貧熄工況下保持穩定.

按照最理想情況，貧熄工況臨界Da 數應等于1，實際計算結果與該值有偏離.初步分析認為，這是因為在本文所涉及的預測方法中對不同狀態下的霧化質量未充分區分，未能體現貧熄狀態下單油路旋流噴嘴霧化質量惡化對燃燒區的影響，未來需要補充霧化相關的實驗，將不同狀態下的霧化質量信息整合到預測模型中.另外，需要對分塊后的每個反應器入口參數進行更精確的計算.當前方法中每個反應器的入口當量比均近似等于體平均當量比，此種計算來流當量比的方法和實際情況存在一定差距，會造成一定的誤差.

圖13 近貧熄工況整體Da 數Fig.13 Global Da numbers under different lean blow-off conditions

5 結論

在本研究中，發展了基于網格化Da 數的貧熄預測方法，基于貧熄實驗結果通過局部Da 數和整體Da 數進行了熄火機理分析.

(1) 整體流動時間尺度和化學反應時間尺度均隨來流流量增大而增大，流動時間的增大主要是因為燃燒區體積隨來流流量的增大而增大，而化學反應時間的增大是因為隨來流流量增大，反應區內平均壓力減小.

(2) 從近貧熄狀態逐漸增加燃油供給，反應區向下游擴張，反應區內的反應強度也隨之提升.

(3) 基于本文整體Da 數計算結果可知，對于空氣流量從0.6 kg/s 增加到1.2 kg/s 的貧熄工況，整體Da 數在7.9 上下波動，平均波動范圍在5%左右.該結果表明本文的預測方法可以在一定程度上捕捉到貧油熄火的關鍵信息，對不同的來流工況具有一定的普適性，可為未來燃燒室貧油熄火邊界預測及工程設計提供理論支持.