湍流燃燒模擬中化學反應的加速算法研究進展

劉再剛, 孔文俊,3,*

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

0 引 言

燃燒是發動機中燃料能量轉換的主要形式，準確地描述燃燒過程，從而得到精確的燃燒數值模擬結果,可以高效指導發動機的設計與實驗。實際應用中，發動機中的燃燒大多是湍流的，燃燒與湍流之間的相互作用機理復雜，使得湍流燃燒數值模擬成為了重要的研究課題。

隨著計算機能力的提升，大渦模擬(Large Eddy Simulation, LES)逐漸成為湍流燃燒數值模擬的重要研究手段。在湍流燃燒LES中需要使用燃燒模型來封閉標量方程中的化學反應源項。目前精度較高的燃燒模型，如小火焰模型[1- 5]、線性渦模型[6- 7]、基于概率密度函數的模型[8- 12]、條件矩模型[13- 14]、加厚火焰面模型[5, 15]、渦耗散模型[16]和部分預混反應器模型[17]等，均需要借助詳細的化學反應機理來描述燃燒時的復雜反應過程，如燃燒場中的局部熄火和再燃、化學反應速率隨壓力的非線性變化、低溫化學反應的影響、污染物排放等。

化學反應機理是指描述化學過程的基元反應的集合。在燃燒時，燃料分解產生大量化學分子和自由基，稱為組分，組分間最基本的化學反應為基元反應。要準確模擬燃燒過程，所需的化學反應機理的規模是巨大的，可以包含數量成百上千的組分和基元反應[18]。而且化學反應組分常微分方程組(Ordinary Differential Equations，ODEs)通常是剛性的，如果使用常規的顯式格式求解，為了計算穩定必須使用極小的時間步長。常見的解決方法是使用基于隱式格式的剛性求解器，如DVODE[19]和DASSL[20]等。因此，在湍流燃燒LES中，求解化學反應方程的時間常占據主要部分[21]。

化學反應機理規模的增加將導致化學反應源項計算量的迅速增長。應用小火焰類燃燒模型，如FGM(Flamelet Generated Manifolds)或FPV(Flamelet/Progress- Variable)等，可以通過預先建表的方法回避在線求解化學反應問題，類似的還有PRISM(Piecewise Reusable Implementation of Solution Mapping)[22]、RS- HDMR(Random Sampling- High Dimensional Model Representation)[23- 24]等映射方法。相比于直接使用總包反應機理，此類方法雖然在一定程度上考慮了詳細化學反應機理的作用，但化學反應表中的變量數量嚴重受限于計算機內存空間，當燃燒的工況不處于預先存儲結果的工況范圍內時，就會缺乏靈活性和適應性，與直接在線求解化學反應相比準確性較低。為了降低直接求解化學反應的計算時間，雖然可以使用在線的建表方法，如ISAT(In Situ Adaptive Tabulation)[25- 29]，但是，ISAT方法對內存需求仍較大，并且在線建表的速度快慢也依賴于直接積分的求解速度。而提高直接求解的計算效率，可以通過使用化學反應機理簡化或采用高效的積分求解器來實現。

對于不同的燃燒問題，反應機理的組分和基元反應并不都是活躍的，可能僅僅其中的某一部分是重要的，而其他部分的反應速率則可以忽略不計。因此，研究人員發展了許多化學反應機理簡化方法，如敏感性分析法[30- 32]、準穩態假設法[33]、局部平衡法[33]、計算奇異攝動法[34]、ILDM(Intrinsic Low Dimensional Manifold)法[35]和DRG(Direct Relation Graph)法[36]及其衍生算法，如DRGEP(Directed Relation Graph Error Propagation)方法[37]、PFA(Path Flux Analysis)[38- 40]、EFA(Element Flux Analysis)[41- 43]和GPS(Global Pathway Selection)方法[44]等。

上述反應機理簡化方法需要在流場計算前完成，并在計算過程中直接使用預先簡化好的反應機理。為了保證簡化機理的準確性，這種預先簡化的機理常包含了較多的組分。然而，對于同一燃燒場中的不同時刻或當地，反應機理可以進一步簡化。因此研究人員發展了一類動態自適應化學法(Dynamic Adaptive Chemistry，DAC)[39, 45- 48]，可以在計算過程中動態簡化化學反應機理。動態簡化主要基于DRG、PFA或GPS算法，因為這些方法能以較小的額外時間花費完成機理簡化。DAC方法已應用于多種簡單模型的燃燒計算，如HCCI(Homogeneous Charge Compression Igni- tion)模型[47, 49]、部分預混攪拌反應器模型[27, 50- 52]、自著火問題和一維火焰傳播[48, 53- 54]等。

DAC方法也在不斷地改進，并與其他簡化方法耦合應用。為了進一步降低機理簡化本身所占用的時間，Sun等[55]提出了自相關動態自適應化學(Correlated Dynamic Adaptive Chemistry，CoDAC)方法。該方法將具有相近狀態的不同空間位置和時間節點網格點進行歸類后統一處理，這與分區模型和CA(Cell Agglomeration)方法[56- 57]的思想相同，自相關技術有效地降低了化學反應機理簡化所需要的計算時間。Oluwole等[54]發展了一種DSRR(Decoupled Species and Reaction Reduction)簡化方法，相比于DRG簡化，具有同等的簡化計算速度并提供了改進的誤差控制方法。DSRR的顯著優點是不需要指定初始組分，但簡化結果可能會導致元素不守恒[54]。 Xie等[58]提出的TSRA(Time- Scale and Jacobian- aided Rate Analysis)將組分分為活躍的(Active)、直接耦合的(Directly coupled)和不重要的(Inconsequential)3類，通過迭代計算來保證元素守恒。此外，DAC方法和ISAT方法結合，構成TDAC(Tabulated Dynamic Adaptive Chemistry)[27, 51, 59- 60]方法，可以進一步提高化學反應計算速度。除了將DAC和ISAT二者結合，Xie等[60]提出的DAAM(Dynamic Adaptive Acceleration Method)可以根據DAC和ISAT的計算耗時自動選擇二者之中較快的方法進行計算。最后，分區方法或CA方法也可以和DAC及ISAT結合[52, 61- 62]，實現更高效的計算。

通過提高求解器計算效率也可以加速化學反應求解，如使用外推法[63]，點隱式求解器[64]、帶預處理的Krylov求解器[65]、剛性去除法[66]和基于多時間尺度(MTS，Multi- Timescale Method)分析的混合多時間尺度(HMTS，Hybrid Multi- Timescale Method)法[67]，HMM(Heterogeneous Multiscale Method)方法[68]和AHI(Adaptive Method for Hybrid Integra- tion)[69]方法等。

求解化學反應常微分方程組，除了使用廣泛應用的顯式或隱式格式外，還可以使用指數格式積分方法 (Exponential Integrator)[70]。這種方法構造含有Jacobian矩陣指數函數形式的積分格式，其優勢在于可以使用顯式大步長時間推進得到穩定的計算結果，適用于剛性常微分方程問題和振蕩系統問題。因此，指數格式被廣泛應用于金融與統計學[71]、正則化理論[72]、電磁場模擬[73]及淺水方程模擬[74]等領域。在指數積分方法中，需要計算Jacobian矩陣的指數函數。在多種計算矩陣指數函數的方法[75]中通用性較好、計算效率較高的方法是收縮乘方和Padé近似[76]方法。然而，這種方法所需的計算量與矩陣維數的3次方成正比，因此適用于較小規模的矩陣，對于大型矩陣，尤其是大型稀疏矩陣，其計算效率很低。一種廣泛采用的解決方法是利用Krylov子空間近似[77]，將原問題投影到其Krylov子空間，從而實現降維，減少問題規模[74, 78- 81]。Krylov子空間近似也常被應用于GMRES (Generalized Minimum Residual)算法和BiCGSTAB (Biconjugate Gradient Stabilized)算法等線性方程組求解方法中。

在燃燒數值模擬中，剛性化學反應問題也適合用指數格式求解。Bisetti[82]和Curtis等[83]的研究表明,四階指數積分格式和Krylov子空間近似方法可應用于求解自著火和HCCI的化學反應機理。然而，四階指數積分格式每步需要計算3次矩陣指數函數，計算量較大。Niesen和Wright[79]發展了一種自適應時間步長和Krylov子空間大小的算法:Phipm算法。該求解器已被應用于自著火問題數值模擬，與DVODE求解器相比，顯著提高了化學反應求解效率[84]。

為了研究湍流燃燒數值模擬中化學反應機理計算的加速方法，本文闡述了DAC和指數格式在湍流燃燒數值模擬中的應用情況。首先，探討了DAC在湍流燃燒數值模擬中的加速效果。然后，考察了Krylov子空間近似的指數格式求解湍流- 化學反應系統的加速效果。

1 動態化學反應機理簡化在燃燒數值模擬中的應用

Liu等[21, 85]使用53組分/325反應的GRI- Mech 3.0反應機理和基于PFA的CoDAC簡化方法對Sandia Flame D進行了LES數值模擬。Sandia Flame D是美國Sandia國家實驗室發展的經典湍流射流火焰，由主射流、值班射流和空氣伴流組成，其中主射流由25%甲烷和75%空氣組成，值班射流為甲烷- 空氣燃燒后的平衡產物。具體參數可見Barlow等[86]的研究。數值模擬結果如圖1所示，通過對比發現，使用CoDAC簡化的LES計算結果與未使用簡化得到的結果符合得很好，表明使用CoDAC可以準確模擬非預混湍流火焰特性。使用CoDAC后，LES中化學反應計算時間總體下降了29%，化學反應速率的求解仍然占據總求解時間的主要部分。Yang等[64, 87- 89]同樣使用基于PFA的CoDAC方法簡化化學反應計算，并結合CoTran方法降低組分輸運系數的計算。對湍流預混火焰的DNS模擬[64]結果表明,簡化方法的加速因子可以達到2.7，在湍流非預混火焰的DNS模擬[89]中可以達到4.0。CoDAC簡化方法也被用于模擬實驗室尺度湍流火焰Sandia Flame D[87]和Flame E[88]的LES，使用的反應機理為20組分/84反應的簡化機理，該反應機理是通過對GRI- Mech3.0機理使用GPS簡化后獲得的。該研究對比了有限速率化學模型和小火焰/進度變量方法，結果表明，2種模型都較好地反映了實驗結果的變化趨勢，但使用有限速率化學模型和CoDAC簡化相比于FPV模型，可以得到更準確的統計結果。

圖1 Sandia Flame D火焰在x/d=3、7.5、15、30、45和60截面處時均溫度的周向分布[21]

Fig.1Radialprofilesoftime-averageoftemperatureforSandiaFlameDatx/d=3,7.5,15,30,45and60[21]

DAC還被應用于超聲速燃燒數值模擬中。Wu等[90]對比了DAC、TDAC方法和骨架機理在超聲速燃燒室數值模擬中的表現。結果表明DAC和TDAC可以準確模擬燃燒室的特性和火焰結構，而骨架機理不能準確預測火焰的穩定特性。在計算效率方面，在57組分/269反應的UCSD反應機理的基礎上使用DAC可以獲得與使用30組分/143反應的骨架機理同等的加速效果，而在DAC的基礎上使用TDAC可以將計算效率再提高1倍。

在實際燃燒設備的數值模擬中，DAC方法也得到了應用，Li等[59]對比研究了TDAC方法和5種不同的化學反應機理簡化方法：DRG、DRGEP、Contino等的DAC方法[91]、PFA和EFA。研究使用非定常雷諾平均方法(URANS)數值模擬了DJHC(Delft Jet- in- Hot- Coflow)燃燒器的燃燒特性，結果表明DRGEP、DAC和EFA方法比DRG和PFA方法準確。對于小規模反應機理，TDAC方法中的制表方法貢獻了主要的計算加速，而對于大規模反應機理，DAC簡化起了主要作用。Zhou和Lu等[92- 93]使用URANS和LES方法模擬了Spray H正庚烷湍流噴霧火焰，結果表明使用DAC的加速因子可達到2.0。在化工應用方面，Huang等[94]使用基于DRG的DAC方法模擬了石英玻璃反應爐的SiCl4/H2/O2燃燒系統，加速因子為2.2。

DAC在內燃機數值模擬中受到了廣泛關注，DAC方法及其他多種化學反應機理簡化方法被應用于加速計算。Shi等[95]使用基于DRGEP的DAC進行了DI(Direct- Injection)發動機的二維數值模擬，并改進了確定初始組分的方法。計算速度的對比表明，在多維發動機數值模擬中的DAC簡化加速效果不如在單區絕熱發動機模型模擬中顯著，但仍有效加速了化學反應計算。對于34組分和61組分的正庚烷反應機理，計算時間分別節省了30%和50%。Contino等[91]使用TDAC模擬了HCCI發動機的燃燒過程，加速因子可高達500.0，而對于柴油發動機，加速因子在9.0左右。Xie等[60]將DAAM方法用于PCCI(Premixed Charge Compression Ignition)發動機的數值模擬，模擬結果顯示DAAM相比于ISAT方法使計算加速了50%，其計算效率也高于ISAT- DAC方法，且數值模擬的準確性處于同一水平。此外，DAC和ISAT方法也可以和CA方法聯合[61- 62]，在分區的基礎上再進行DAC簡化或ISAT建表，可以進步一提高計算速度。 Zhou等[61]的研究表明，對于約40個組分的反應機理，CA聯合DAC的方法可以使計算加速一倍，而對于約140個組分的反應機理，加速因子可達到9.0左右。

DAC在湍流燃燒數值模擬中的加速效果不盡相同，其影響因素是多方面的：

(1) DAC簡化的加速效果與數值模擬結果誤差是相關的，用戶所能容忍的誤差越大，DAC簡化能實現的加速就越大。誤差的調節方式與DAC中具體應用的簡化方法有關，如DRG、DRGEP、PFA等通過簡化閾值調節，合理閾值的大小還根據所使用簡化機理的不同而不同，最終的誤差需要通過結果來驗證；而DSRR則通過給定容許誤差來調節。

(2) 加速效果與燃燒特性有關。Liu等[21]的研究表明，對于燃燒反應不活躍的情況，如自著火過程發生前的慢速反應區域和自著火完成后的平衡區域，以及一維預混火焰面上游的未燃區和火焰面下游的已燃區，燃燒過程中活躍的組分數較小，使用DAC簡化可以顯著降低局部簡化機理中的組分，從而大大提高計算速度。而對于反應較為劇烈的區域(如自著火時刻或火焰鋒面處)，化學反應過程復雜，涉及組分多，因此DAC簡化產生的局部簡化機理規模較大，甚至與輸入的反應機理規模相當，此時DAC的加速效果不明顯。

(3) 加速效果與輸入反應機理有關，應用規模越大的反應機理通常能獲得越明顯的加速效果。例如Wu等[90]的研究表明，同等網格數量的情況下，對于57組分、75組分和111組分的反應機理，DAC的加速因子分別為1.9、2.2和2.5。

在實際應用中，DAC如果與其他加速手段耦合，加速效果也會受到影響。例如DAC方法與ISAT方法耦合時[51]，DAC的加速效果隨ISAT誤差閾值的增大而略有降低。在湍流燃燒的并行計算中，化學反應計算的并行方式也對DAC的加速效果有較大影響。Liu等[21]的研究表明，由于并行LES中化學反應的計算效率受火焰面附近處理器核心的計算效率控制，在并行計算中，化學反應求解完成后，核與核之間需要進行一次同步并通信，然后再進行下一步計算，因此總體的化學反應計算時間由各個核心中最慢的決定。例如，圖2所示是不同核心上的化學反應計算時間tODE,p，可以看到不同計算區域中的tODE,p差別很大，火焰反應區附近的tODE,p顯著高于其他反應活性較低的區域。放大圖2(a),顯示出最高的tODE,p出現在火焰中r=d, 10d

圖2 在各處理器核心上的(a) ODE時間(tODE,p) 以及(b) ODE時間、(c)溫度T、(d) PFA占比(φPFA,p) 和(e)簡化機理中組分數(NS)的局部放大視圖[21]

Fig.2Contourplotof(a)ODEtimeoneachprocessor(tODE,p)togetherwithzoom-inviewof(b)ODEtimeoneachprocessor, (c)temperatureT, (d)percentageofcellsperformingPFA(φPFA,p)and(e)numberofselectedspeciesinlocallyreducedmechanism(NS)[21]

圖3 使用DAC和CoDAC在初始條件為φ=1.0,T0=1300K和簡化閾值εr=0.1時相對計算時間和溫度隨時間的變化[21]

Fig.3Evolutionofrelativecomputationaltimetrandtemperatureforauto-ignitionsimulationwithDACandCoDACatφ=1.0,T0=1300Kandreductionthresholdεr=0.1[21]

2 Krylov子空間近似的指數積分格式

對于如下的燃燒化學反應常微分方程問題：

(1)

二階的指數積分格式可以表示為[70]：

(2)

式中φ表示包括溫度和組分質量分數在內的微分方程變量，ωk表示各個變量的化學反應速率，h為時間步長，當時間步長較短時可認為Jacobian矩陣Jk=?ωk/?φk是近似恒定的，矩陣的φ函數定義為：

φ0(z)=exp(z),φl(z)=zφl+1(z)+(1/l!)

(3)

此時，常微分方程的求解變成了對一般形式為φ1(τA)v的向量的求解，其中A為在n×n維實空間Rn×n中的矩陣，v為在n維實空間Rn中的向量。

當n很大時，φ函數的計算需要巨大的計算量，利用Krylov子空間的性質，組合φ1(τA)v可以方便地在更小的m維Krylov子空間中得到近似。Krylov子空間是由向量序列v,Av,A2v,…,Am-1v張成的空間，構造Krylov子空間的算法稱為Arnoldi算法[77]。原向量φ1(τA)v在Krylov子空間中的近似可以表示為[77]：

(4)

式中β=‖v‖，Vm=(v1,v2, …,vm)為Krylov子空間的基向量，Hm為上Hessenberg矩陣，是原矩陣A在Krylov子空間中的投影，e1和em分別表示第1維和第m維單位向量。當矩陣A在Krylov子空間中投影為Hm時，Hm的特征值稱為Ritz值，隨著m的增加，收斂于原矩陣A模最大的特征值。這一重要性質使Krylov子空間廣泛應用于近似求解矩陣特征值問題[97]。通過近似式(4)，原n維問題降為m維問題，使指數格式求解所需的計算量顯著降低。

Bisetti[82]使用四階Rosenbrock指數格式exp4[98]計算了火焰的自著火過程，結果表明,盡管使用大時間步長時非線性系統剛性很大，但指數格式方法仍是穩定的。Curtis等[83]對比了CPU和GPU計算時exp4和expbr43指數積分格式的計算效率，結果表明2種指數積分格式在CPU和GPU上都慢于隱式求解器CVODE和Radau- IIA。然而，四階指數積分格式每推進1個時間步長需要計算3次矩陣指數函數。而矩陣指數函數的計算十分耗時，為了實現加速計算，應盡量減少計算次數。因此，劉再剛等[84]使用二階指數格式，如式(2)，并引入Schur分解方法有效控制了算法的舍入誤差[99]，發展了用于加速燃燒化學反應計算的指數積分求解器。

指數格式算法的計算效率可以使用簡單的燃燒模型來估計。帶配對混合(Pairwise mixing)的局部攪拌反應器(Partially- Stirred Reactor, PaSR)是一個假設的反應器[100]。這種反應器可以用于模擬湍流燃燒數值模擬中的小尺度混合和反應過程，常用于檢驗燃燒化學反應計算方法和ISAT方法的準確性和計算效率。模擬使用的燃料和參數與Liu等[21]的研究相同，PaSR入流由3部分組成，分別為甲烷、甲烷- 空氣燃燒后的平衡產物以及空氣，并分別對應Sandia Flame系列火焰的主射流、值班射流和空氣伴流。在Liu等[21]的研究中，CoDAC在PaSR模擬中得到了與LES模擬相近的加速結果，表明了PaSR能合理估計LES的計算效率。對于燃燒化學反應問題，還可以同時使用一些加速算法。CoDAC簡化應用于指數格式可以減少算法中矩陣A和v的維數，SpeedCHEM(SC)庫[101]中的稀疏矩陣算法和解析Jacobian矩陣算法可以降低矩陣存儲所需的內存空間，并加速Jacobian矩陣和Arnoldi算法中矩陣- 向量乘法Av的計算。MTS方法[69]可以分離剛性問題中的剛性部分和非剛性部分，從而降低剛性方程組的數量。在指數格式中，作用與CoDAC相同，可以減少算法中矩陣A和v的維數。Liu等[99]將指數積分格式與SC庫、CoDAC簡化和MTS方法結合，并與DVODE求解器的計算結果作了詳細對比，不同的組合算例列于表1。

表1 使用不同加速方法的算例設置Table 1 Simulated cases for different acceleration methods

圖4所示是分別使用GRI- Mech 3.0[102]和USC Mech II[103](111組分/784反應)計算得到的結果。圖中的加速因子和誤差均以表1中的參考工況R為基準，加速因子和誤差的詳細定義可參考文獻[99]。從圖中可以看出，在同等精度并同時采用了SC的情況下，使用指數格式的計算速度均顯著快于使用隱式格式求解器DVODE，并且加速效果顯著優于使用DVODE結合3種加速算法。但指數格式結合CoDAC和MTS加速算法時，不僅計算精度下降，計算效率也有所下降。從圖5可以看出，效率下降主要是由CoDAC和MTS帶來的額外開銷導致的。因此只使用含SC的指數格式可以得到最佳的加速效果。

(a)

(b)

圖4 使用反應機理(a)GRI- Mech 3.0[102]和(b)USC Mech II[103]得到的表1中算例的加速因子γ與最大平均相對誤差ε的關系。Xsc(εr)，Xsm(α)和Xscm(εr,α)表示表1所列工況算例，括號中的數值εr和α分別表示CoDAC簡化閾值和MTS安全因子[99]

Fig.4ThespeedupfactorγasfunctionsofmaximumaveragedrelativeerrorεforthecasesinTable1with(a)GRI-Mech3.0[102]and(b)USCMechII[103].Xsc(εr),Xsm(α)andXscm(εr,α)representthecasesinTable1andnumbersεrandαintheparenthesesrepresentthereductionthresholdofCoDACandthesafetyfactorofMTS,respectively[99]

圖5 使用USC Mech II[103]時，表1中的算例中用于計算快方程、慢方程和CoDAC簡化的壁面時間[99]

Fig.5Walltimespentonfastequationintegration,slowequationintegrationandCoDACreductionforthecasesinTable1withUSCMechII[99, 103]

指數格式的加速效果一方面來源于大步長顯式時間推進的優勢，另一方面是由于與Krylov子空間近似的降維效果。從圖6中可以看到，Krylov子空間的大小增長速度遠低于剛性方程組規模的增速。由此可知，在化學反應計算中應用指數格式在保證計算精度的同時可以顯著提高計算速度。因此，其在燃燒LES中有很好的應用前景。指數格式求解器直接求解化學反應方程組，如果與查表類加速方法，如ISAT相結合，可以進一步降低化學反應計算時間，從而實現快速高精度的湍流燃燒數值模擬，從而使得詳細化學反應機理用于實際發動機燃燒室湍流燃燒數值模擬成為可能，并以此指導湍流燃燒研究和燃燒器設計與實驗。

圖6 使用GRI- Mech 3.0[102]和USC Mech II[103]機理時表1中算例Ascm和Bscm的Krylov子空間維數(mKrylov)隨快方程數量(nfast)的變化[99]

Fig.6ThedimensionoftheKrylovsubspace(mKrylovasafunctionofthenumberoffastequations(nfast)forcasesAscmandBscmshowninTable1withthemechanismsofGRI-Mech3.0[102]andUSCMechII[99, 103]

3 結 論

在湍流燃燒數值模擬中使用詳細化學反應機理可以顯著提高燃燒反應計算精度，但是由此導致的計算量激增是阻礙大規模高精度燃燒數值模擬的主要因素之一。在計算過程中動態、自適應地進行化學反應機理簡化可以獲得一定的加速效果。然而，在并行燃燒數值模擬中，由于燃燒場中火焰反應區域的分布不均勻，處理器核心的負載極度不平衡，且核心之間需要進行必要的同步通信，加速效果在很大程度上不能充分發揮。相比之下，使用剛性求解器得到的加速將直接作用于每個處理器核心，更有利于整體計算效率提高。相比于常見的隱式格式以及CoDAC和MTS加速方法，在同等精度下，Krylov子空間近似的指數積分格式對化學反應計算的加速效果更為顯著，加速主要來源于指數格式的顯式大時間步推進能力和Krylov子空間近似的降維效果。

為了進一步提高化學反應的求解速度，首先，可行的方法是以高效積分求解器為基礎，緊密結合查表類方法，互相彌補2種方法的不足：用查表代替重復性計算，用高效積分補充表中的元素。其次，使用自適應網格以及自適應平衡處理器核心負載的算法，可以有效避免計算資源的空閑，提高整體計算效率。最后，燃燒化學反應的求解涉及的浮點數運算密集，較適合使用圖形處理器架構進行運算，這也是燃燒數值模擬加速的一個發展方向。

致謝:感謝國家自然科學基金項目(91441131和U1738113)的支持。