基于原位自適應制表與動態自適應機理的無焰燃燒氮轉化研究

成鵬飛，李鵬飛，劉 璐，胡 帆，柳朝暉，鄭楚光

(華中科技大學 煤燃燒國家重點實驗室，湖北 武漢 430074)

0 引 言

天然氣等氣體燃料是國家能源轉型期間大力發展的清潔能源，相較于其他化石燃料總體污染物生成較少，但NOx排放量亦不容忽視。我國對低氮排放的管控起步較晚但愈加嚴格，京津冀地區新建燃氣鍋爐NOx排放限值僅30 mg/m3[1]。因此，隨著我國氣體能源消費量的逐年增加，需要不斷發展低氮燃燒技術并深入研究氮轉化機理。

無焰燃燒是近年來備受關注的新型高效低氮燃燒技術，其技術關鍵在于反應混合物高速射流對爐內煙氣的強卷吸作用，被煙氣充分稀釋并加熱至超過自燃點的反應混合物將擴散至整個燃燒區且在局部低氧氛圍下緩慢燃燒，因此整體表現為彌散燃燒狀態且無明顯火焰鋒面[2-3]。相較于傳統低氮燃燒技術，如火焰冷卻技術、分級燃燒技術等[4]，無焰燃燒可從根本上抑制局部高溫產生，并在保證燃燒效率、燃燒強度和燃燒穩定性的條件下實現更低的NOx排放。

計算流體力學(CFD)是燃燒研究的重要手段之一，具有簡便快捷、成本低、數據豐富等優點，有效促進了無焰燃燒反應特征及NO生成調控的研究發展。無焰燃燒具有分布式低反應速率燃燒區和典型的中低溫燃燒特性，需建立有限反應速率模擬以考慮強烈的湍流化學反應交互[5-7]。提高NO生成的數值預測精度是無焰燃燒數值模擬的關鍵[8-10]，相較于簡化的半經驗NOx后處理模型，基于耦合含氮詳細反應機理的有限反應速率模擬可將燃燒氧化過程與氮化學反應耦合計算并考慮再循環煙氣中的初始NO對著火的加速作用，從而準確預測著火和主反應區，并得到含氮中間組分分布和NO生成路徑等氮轉化關鍵信息[11]。

由于組分數、反應數較多的含氮詳細反應導致實際數值計算過程耗時較長，通常需要耦合簡化計算方法在不顯著降低計算精度的條件下實現計算加速。原位自適應制表法(In-Situ Adaptive Tabulation, ISAT)可利用已存儲的數據代替組分反應速率微分方程的復雜求解過程以減少計算耗時，在當前燃燒模擬中得到廣泛應用[12-14]。動態自適應機理(Dynamic Adaptive Chemistry, DAC)可基于特定時間內每組當地瞬時熱化學狀態參數將詳細反應機理實時簡化為含最少組分和反應的精確子機理，實現計算加速。當前已有少數耦合DAC算法的氣體及液體燃料燃燒模擬研究[15-19]，且筆者團隊已將其拓展應用于煤粉等固體燃料[20]，但仍需進一步發展與優化該計算方法。

為提高無焰燃燒模擬的預測精度和計算速度并進一步研究分析該燃燒過程中的氮轉化機理，筆者耦合含氮詳細反應機理和簡化計算方法進行了CH4/H2熱伴流無焰燃燒模擬，發展了原位自適應制表法與動態自適應機理耦合計算方法(ISAT-DAC)并進行優化評估，最后基于經實驗驗證的模擬計算結果進行了氮化學反應動力學分析。

1 基于ISAT-DAC算法的數值模擬

1.1 原位自適應制表法

對于確定的反應動力學，最終反應狀態是初始反應狀態和時間步的函數。理論上，反應的映射關系在完成一遍計算后可制成表，當插值得出的結果滿足精度要求時，可不必進行積分運算，從而節約計算時間。但由于實際尺寸過大，包含所有溫度、壓力、組分和時間步的預制表是不允許的，須在計算過程中制表。ISAT實質上是一種利用誤差控制實時制表訪問組合空間區域的加速計算方法。

在第1個時間步Δt內，若初始反應狀態為φ0，則末狀態為φ1可按式(1)計算：

(1)

其中，S為化學源項；t為反應時間。計算式-創建表中第1項，并包含式(2)中的映射梯度矩陣A：

(2)

(3)

(4)

1.2 動態自適應機理簡化

燃燒火焰化學成分在不同區域差別很大，因此無須在模擬計算域全局使用統一的反應動力學機理。DAC可將詳細機理當地實時簡化為精確子機理，以實現計算加速。

簡化計算過程中首先需指定目標組分A和簡化閾值λ，目標組分通常確定為燃料主要成分及燃燒過程中的關鍵組分，簡化閾值需綜合考慮計算誤差和簡化加速效果。通過直接關系圖法[21]在每個網格單元內計算非目標組分B對目標組分A生成的貢獻rAB：

(5)

其中，ωi為基元反應i的化學反應速率；vA,i為A組分在反應i中的化學計量數。當rAB>λ時將組分B保留在動態機理中。每個非目標組分對目標組分生成的貢獻均需計算。

然后采用相同方法確定目標組分的間接貢獻組分。若組分B包含于對目標組分直接貢獻較大的組分集合中，需計算其余非目標組分K對組分B的貢獻且當rBK>λ時保留K組分。對所有組分執行上述操作以確定最終保留在簡化機理中的組分。

最后將所有與保留組分無關的反應去除，即得到每個網格單元的當地簡化機理。通過求解簡化后的低維常微分方程組，可顯著縮短計算時間。因此，當λ值較大時，可減少保留組分及反應數并加速計算，但同時導致模擬誤差的增大。

ISAT-DAC耦合算法，由CONTINO等[22]首次提出并應用于內燃機均質壓燃(Homogeneous Charge Compression Ignition, HCCI)模擬，可充分發揮2種計算方法的優勢，其計算流程如圖1所示。

圖1 ISAT-DAC計算流程Fig.1 ISAT-DAC calculation flow chart

1.3 模擬對象與計算細節

無焰燃燒模擬研究對象為同軸熱伴流射流反應器[23]。燃燒器結構如圖2所示，主要由中心燃料管和熱氧化劑伴流管組成。以內置二次燃燒器的燃燒產物作為熱氧化劑伴流，并通過控制燃燒器側面的空氣及氮氣流量調節氧濃度。中心燃料管直徑為4.25 mm, 熱伴流流管直徑為82 mm。

圖2 燃燒器結構[23]Fig.2 Burner structure[23]

數值計算域如圖3所示，將燃料噴口下游一段燃燒區域作為計算域，且由于全尺度含氮詳細反應機理模擬需較大計算量及該燃燒器結構對稱，采用軸對稱二維雷諾平均模擬來簡化計算三維圓管同軸射流[24-26]。計算域由噴口處軸向延伸500 mm，縱向延伸80 mm，其中縱向2.125 mm為燃料進口，燃料進口邊界至縱向41 mm為熱伴流進口，其余為空氣風道。實驗及模擬數據監測取中心射流與周圍空氣混合前Z=30 mm及混合后Z=120 mm截面，數據豐富且有代表性，便于進行無焰燃燒數據分析與模擬驗證[23]。采用23 420單元數的網格并進行網格獨立性計算，分析結果如圖4所示，不同網格數量的模擬結果相近且與實驗數據吻合較好，因此當前較低網格數量即可滿足高精度模擬要求并加速計算。

圖3 模擬計算域Fig.3 Computational domain of simulation

圖4 網格無關性驗證Fig.4 Grid independence validation

DALLY等[23]基于熱伴流射流反應器建立了O2質量分數3%、6%和9%實驗研究，結果表明熱伴流O2質量分數3%條件下的燃燒狀態最接近無焰燃燒。因此，筆者基于O2質量分數3%工況下的實驗數據對DAC、ISAT及ISAT-DAC耦合算法的計算效果進行評估，并分析ISAT-DAC算法在無焰燃燒模擬中的適用性。具體工況參數見表1[23]。數值計算基于FLUENT求解器。采用標準k-ε湍流模型和標準壁面函數[27]，將模型系數Cε1由1.44 修正為1.60以提升圓管射流預測精度[28]。采用高精度離散坐標(DO)輻射模型求解計算域內若干離散立體角的輻射傳遞方程，引入灰氣體加權和(WSGG)模型計算氣體輻射特性參數，總發射率的空間變化是氣體組分和溫度的函數[29]。均相燃燒與氮轉化模擬采用有限速率渦耗散概念(EDC)模型耦合GRI-Mech 3.0含氮詳細反應機理(53種組分、325步反應)[30]，以描述由基元反應和中間組分控制的瞬態燃燒過程[31]。速度壓力耦合方程采用SIMPLE算法以節約計算資源。空間離散采用高階QUICK格式以提高燃燒模擬精度。

表1 工況條件[23]Table 1 Operating conditions[23]

2 模擬結果與討論

2.1 簡化算法的發展評估

基于CH4/H2熱伴流無焰燃燒工況綜合分析動態自適應機理簡化法、原位自適應制表法及ISAT-DAC耦合簡化算法的預測精度和簡化加速效果，從而對無焰燃燒模擬簡化算法進行發展與評估。

首先分析比較3種簡化方法應用于無焰燃燒模擬研究中的預測精度。DAC、ISAT及ISAT-DAC算法的簡化閾值均設置為10-4。DAC和ISAT-DAC算法的初始目標組分集合確定為CH4、HO2、CO和NO，其中CH4為燃料主要組分，CO、NO為燃燒關鍵組分、HO2為著火關鍵組分。以簡化計算結果與直接積分模擬結果的相對誤差作為評估標準，并按式(6)計算：

ψφ=|φID-on-φODE|/φODE，

(6)

式中，ψφ為基于3種簡化方法的計算結果與直接積分模擬結果的相對誤差；φID-on為基于ISAT、DAC或ISAT-DAC算法模擬得到的溫度或組分濃度等燃燒關鍵參數值；φODE為直接積分計算結果。

一般認為，高精度燃燒模擬中主要參數或組分的預測誤差限值可取5%，低含量組分的預測誤差限值可取8%。選取溫度及O2、H2O、CH4、CO、OH、HO2、NO體積分數計算相對誤差值ψφ，結果如圖5所示。可知基于3種簡化方法的模擬計算結果均滿足較高精度要求，且所有參數相對誤差均未超過6%。不同機理簡化方法對溫度及O2、H2O、CH4體積分數等主要參數或組分的模擬結果影響較小，相對誤差均低于1.5%，其中CH4體積分數受不同簡化方法的影響較大。CO、OH、HO2、NO等低體積分數組分的模擬相對誤差受不同簡化方法的影響顯著。基于DAC算法的所有燃燒關鍵參數模擬結果相對誤差最低并可忽略不計，ISAT-DAC算法次之，ISAT算法的相對誤差相較其他算法顯著升高。

圖5 3種簡化方法的相對誤差統計Fig.5 Relative error statistics of three simplified methods

基于DAC和ISAT-DAC簡化算法的活躍組分分布如圖6所示。圖6(a)中單獨采用DAC算法的簡化計算域保留了更多活躍組分，因此數值預測精度略高，但過多的活躍組分會降低計算速度。圖6(b)中采用ISAT-DAC耦合算法可將不活躍反應區域內的組分基本去除，從而精確識別反應區域，且模擬計算的相對誤差亦在可接受范圍。

圖6 基于DAC和ISAT-DAC算法的活躍組分分布Fig.6 Active species distribution based on DAC and ISAT-DAC

分別統計基于直接數值積分法、動態自適應機理簡化法、原位自適應制表法和ISAT-DAC耦合簡化算法的無焰燃燒模擬計算時長，按式(7)計算加速因子并分析不同簡化算法的加速計算效果：

f=tfull/treduct，

(7)

式中，f為加速因子；tfull為基于直接數值積分法模擬的計算時長；treduct為基于簡化算法模擬的計算時長。

經歸一化計算，直接數值積分法模擬加速因子為1。不同簡化算法歸一化計算結果見表2，基于DAC、ISAT及ISAT-DAC算法的模擬加速因子分別為1.3、1.8和3.2，因此不同簡化算法均可減少計算時長。對于耦合GRI-Mech 3.0機理的氣體無焰燃燒模擬，采用ISAT算法相較于DAC算法可獲得更快的計算加速，且ISAT-DAC耦合算法可在DAC或ISAT算法基礎上大幅提高計算速度。

表2 3種簡化方法的加速因子Table 2 Acceleration factors of three simplified methods

2.2 動態自適應機理最優簡化閾值

據第2.1節不同簡化算法的評估分析可知，DAC算法可在一定程度上實現無焰燃燒計算加速且精度最高，進一步結合ISAT算法的計算優勢可在保證較高精度的條件下大幅提升計算速度。基于ISAT計算加速的方式已在燃燒模擬中得到廣泛應用，本節旨在尋找ISAT-DAC耦合算法的最佳DAC簡化閾值。

反應機理仍選用GRI-Mech 3.0[30]。ISAT閾值設定為10-4，DAC簡化閾值分別選取10-4、10-3、10-2和10-1，并按式(6)分別計算基于不同簡化閾值的相對誤差。溫度及O2、H2O、CH4、CO、OH、HO2、NO體積分數的模擬相對誤差如圖7所示。

圖7 ISAT-DAC算法中不同DAC閾值相對誤差統計Fig.7 Relative error statistics of different DAC thresholds in ISAT-DAC

由圖7(a)可知，溫度及O2、H2O、CH4體積分數模擬結果受不同DAC簡化閾值影響較小，相對誤差均不超過4%。DAC簡化閾值10-2條件下的模擬相對誤差仍較小，進一步增大至10-1時相對誤差顯著升高。圖7(b)中CO、NO、HO2、OH體積分數對不同DAC簡化閾值更加敏感，且DAC簡化閾值10-2條件下的低含量組分模擬相對誤差未超過8%。DAC簡化閾值10-1條件下的NO和HO2體積分數的模擬相對誤差超過10%，且CO和OH體積分數的模擬相對誤差亦趨近10%。綜上，DAC簡化閾值10-2條件下仍可保證較高的模擬精度。

基于ISAT-DAC耦合算法和不同DAC簡化閾值的活躍組分分布如圖8所示。隨著DAC簡化閾值的增大，ISAT-DAC耦合算法的計算簡化能力增強。DAC簡化閾值10-4時整體反應區域基本包含GRI-Mech 3.0詳細機理的所有組分，DAC簡化閾值10-3時整體反應區域組分數略有減少，但大部分區域保留了較多組分。當DAC簡化閾值增至10-2時，反應區域中保留的活躍組分數進一步減至約30種。簡化閾值10-1時的反應區域活躍組分數目僅余20種左右，因此CO、NO、HO2、OH等體積分數較低的組分在此閾值條件下的模擬相對誤差較大。

圖8 ISAT-DAC算法中不同DAC閾值活躍組分分布Fig.8 Active species distribution of different DAC thresholds in ISAT-DAC

統計基于ISAT-DAC算法和不同DAC簡化閾值的模擬計算時長，歸一化計算后加速因子見表3。DAC簡化閾值10-2條件下的模擬計算中，計算速度相較于10-4時顯著提升。繼續增大簡化閾值對加速因子影響較小，而燃燒關鍵參數，尤其是CO、NO、HO2、OH體積分數的模擬精度顯著降低。因此，綜合相對誤差與簡化加速效果，ISAT-DAC耦合算法中DAC最佳簡化閾值為10-2。

表3 ISAT-DAC算法中不同DAC閾值加速因子Table 3 Acceleration factors of different DAC thresholds in ISAT-DAC

2.3 不同復雜度機理的簡化效果

當前數值計算中一般無須選用高復雜度詳細機理即可滿足所需精度要求，但基于不同組分數和反應數的詳細機理可對ISAT-DAC耦合算法在無焰燃燒模擬中的適用性和優越性進行更加全面的評估。采用3種詳細機理，機理中所包含的組分數和反應數由低到高分別為DRM22機理(22種組分、83步反應)、GRI-Mech 3.0機理(53種組分、325步反應)[27]及PG2018機理(151種組分，1 397步反應)[32]。3種詳細反應機理經ISAT-DAC算法簡化后的活躍組分分布如圖9所示，由于無焰燃燒過程中存在復雜的湍流化學相互作用，3種不同詳細機理模擬的燃燒反應區域均保留了大多數組分和反應。由于DRM22機理的組分數和反應數較少，簡化空間較小，燃燒反應區域幾乎保留了該機理中的全部組分，非反應區域的組分數基本為0以節省計算耗時。GRI-Mech 3.0機理和PG2018機理的組分數及反應數較多，ISAT-DAC算法簡化效果顯著。基于GRI-Mech 3.0機理的燃燒反應區域中約保留原機理中75%的組分，且反應區域上游的活躍組分數進一步減至20種左右。PG2018機理簡化空間最大，燃燒反應區域中消除了原機理中約50%的組分。

圖9 基于ISAT-DAC算法的不同復雜度機理活躍組分分布Fig.9 Active species distribution of different complexity mechanisms based on ISAT-DAC

分別統計基于3種不同復雜度的詳細反應機理和不同簡化算法的模擬計算時長并按式(6)計算加速因子，結果見表4。隨著機理復雜度的提高，模擬簡化空間增大，基于不同簡化算法計算的加速因子均增大。耦合DRM22機理或GRI-Mech 3.0機理時的DAC簡化加速因子略低于ISAT算法，而耦合PG2018機理時的結果相反，因此DAC算法所實現的簡化加速效果相較ISAT算法對反應機理的復雜程度更加敏感，即DAC算法對高復雜度反應機理表現出更加優異的簡化計算潛力。不同簡化算法對于耦合低復雜度DRM22骨架反應機理模擬的加速因子較小，基于ISAT-DAC算法模擬得到的加速因子僅為1.7，基于DAC算法模擬的簡化加速效果可忽略不計。相較于DRM22機理，不同簡化算法對于耦合GRI-Mech 3.0詳細機理模擬的加速因子略有增大，且ISAT-DAC耦合算法的簡化加速效果較為可觀。相較于DRM22機理或GRI-Mech 3.0機理，對于耦合多組分、多反應的PG2018大型詳細反應機理模擬，基于不同簡化算法的加速效果顯著增加，ISAT-DAC耦合算法所對應的加速因子高達10.8，從而大幅縮減計算時長。

表4 不同復雜度機理的簡化加速因子Table 4 Acceleration factors of different complexity mechanisms

2.4 實驗驗證與氮轉化路徑分析

將基于GRI-Mech 3.0機理及ISAT-DAC耦合算法且DAC簡化閾值10-2條件下的數值計算結果與O2質量分數(y(O2))3%工況下的CH4/H2熱伴流無焰燃燒實驗數據進行對比分析以驗證所選簡化算法及誤差閾值的準確性。

軸向30、120 mm截面處分析計算結果如圖10所示(x表示物質的體積分數)，測量不同截面的溫度及O2、H2O、OH、CO的質量分數及NO體積分數，相應模擬計算結果與實驗數據基本吻合，可較準確預測上述燃燒關鍵參數，包括廣為關注的低濃度CO和NO組分。

圖10 模擬與實驗的數據驗證Fig.10 Data validation between simulations and experiments

NNH中間體由H與N2在火焰條件下快速生成，并進一步結合O原子生成NO：

(8)

(9)

N2O中間體由O與N2及火焰中的其他組分M通過復合反應生成：

(10)

部分N2O中間體通過下列反應被重新還原為N2：

(11)

(12)

(13)

其余N2O中間體進一步結合O原子生成NO：

(14)

當前燃燒模擬研究大多基于ISAT算法實現詳細機理模擬的計算加速，因此本文對比分析了采用ISAT算法模擬時的氮轉化路徑，如圖11(b)所示。基于單獨ISAT算法模擬時的氮轉化路徑和反應速率與ISAT-DAC簡化算法有所差別，但主要反應路徑基本一致。大多數反應路徑的反應速率在同一數量級，少數反應速率相差一個數量級。由于DAC算法對含氮反應路徑的簡化作用，基于ISAT-DAC耦合算法的模擬反應路徑相較ISAT算法更加精簡。

圖11 CH4/H2熱伴流無焰燃燒氮氧化物轉化路徑Fig.11 NOx conversion path of CH4/H2 flameless combustion with a hot co-flow

3 結 論

1)不同簡化計算方法對CO、OH、HO2、NO等低濃度組分的模擬相對誤差影響顯著。采用單一DAC算法精度最高，但計算加速效果并不明顯；ISAT算法精度最低，而計算速度相較于DAC算法提升較小；ISAT-DAC算法可保證模擬所需的較高精度并大幅實現計算加速，是綜合預測精度與簡化加速效果的最優算法。

2)基于ISAT-DAC耦合算法且DAC簡化閾值設定為10-2的模擬計算中，相較于DAC簡化閾值10-4條件下可保留主要燃燒組分及反應，并在保證較高模擬精度的條件下實現顯著計算加速。當簡化閾值增至10-1，燃燒關鍵參數的模擬相對誤差顯著升高，而計算速度提升較小。因此DAC簡化閾值10-2時綜合效果最優。

3)多組分、多反應的復雜機理具有更高的簡化計算空間。基于ISAT-DAC算法和PG2018大型含氮詳細機理模擬的燃燒反應區域僅保留約50%的原機理組分，可在保證所需模擬精度的條件下實現約10.8倍的計算加速。DAC算法相較于ISAT算法對高復雜度機理表現出更優異的簡化計算潛力。

4)基于實驗驗證與氮轉化路徑分析，CH4/H2熱伴流無焰燃燒NO生成主要取決于NNH路徑和N2O中間體路徑，熱力型NO及快速型NO可忽略不計。相較于ISAT算法，基于ISAT-DAC耦合算法的氮轉化路徑可保留主要NO反應并適當精簡。