劉再剛, 孔文俊,3,*

(1. 中國科學(xué)院工程熱物理研究所 中國科學(xué)院輕型動(dòng)力重點(diǎn)實(shí)驗(yàn)室, 北京 100190; 2. 中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院, 北京 100049; 3. 北京航空航天大學(xué) 宇航學(xué)院, 北京 100191)

0 引 言

燃燒是發(fā)動(dòng)機(jī)中燃料能量轉(zhuǎn)換的主要形式，準(zhǔn)確地描述燃燒過程，從而得到精確的燃燒數(shù)值模擬結(jié)果,可以高效指導(dǎo)發(fā)動(dòng)機(jī)的設(shè)計(jì)與實(shí)驗(yàn)。實(shí)際應(yīng)用中，發(fā)動(dòng)機(jī)中的燃燒大多是湍流的，燃燒與湍流之間的相互作用機(jī)理復(fù)雜，使得湍流燃燒數(shù)值模擬成為了重要的研究課題。

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

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

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

對(duì)于不同的燃燒問題，反應(yīng)機(jī)理的組分和基元反應(yīng)并不都是活躍的，可能僅僅其中的某一部分是重要的，而其他部分的反應(yīng)速率則可以忽略不計(jì)。因此，研究人員發(fā)展了許多化學(xué)反應(yīng)機(jī)理簡(jiǎn)化方法，如敏感性分析法[30- 32]、準(zhǔn)穩(wěn)態(tài)假設(shè)法[33]、局部平衡法[33]、計(jì)算奇異攝動(dòng)法[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]等。

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

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

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

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

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

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

1 動(dòng)態(tài)化學(xué)反應(yīng)機(jī)理簡(jiǎn)化在燃燒數(shù)值模擬中的應(yīng)用

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

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

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

DAC還被應(yīng)用于超聲速燃燒數(shù)值模擬中。Wu等[90]對(duì)比了DAC、TDAC方法和骨架機(jī)理在超聲速燃燒室數(shù)值模擬中的表現(xiàn)。結(jié)果表明DAC和TDAC可以準(zhǔn)確模擬燃燒室的特性和火焰結(jié)構(gòu)，而骨架機(jī)理不能準(zhǔn)確預(yù)測(cè)火焰的穩(wěn)定特性。在計(jì)算效率方面，在57組分/269反應(yīng)的UCSD反應(yīng)機(jī)理的基礎(chǔ)上使用DAC可以獲得與使用30組分/143反應(yīng)的骨架機(jī)理同等的加速效果，而在DAC的基礎(chǔ)上使用TDAC可以將計(jì)算效率再提高1倍。

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

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

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

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

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

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

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

圖2 在各處理器核心上的(a) ODE時(shí)間(tODE,p) 以及(b) ODE時(shí)間、(c)溫度T、(d) PFA占比(φPFA,p) 和(e)簡(jiǎn)化機(jī)理中組分?jǐn)?shù)(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和簡(jiǎn)化閾值εr=0.1時(shí)相對(duì)計(jì)算時(shí)間和溫度隨時(shí)間的變化[21]

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

2 Krylov子空間近似的指數(shù)積分格式

對(duì)于如下的燃燒化學(xué)反應(yīng)常微分方程問題：

(1)

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

(2)

式中φ表示包括溫度和組分質(zhì)量分?jǐn)?shù)在內(nèi)的微分方程變量，ωk表示各個(gè)變量的化學(xué)反應(yīng)速率，h為時(shí)間步長(zhǎng)，當(dāng)時(shí)間步長(zhǎng)較短時(shí)可認(rèn)為Jacobian矩陣Jk=?ωk/?φk是近似恒定的，矩陣的φ函數(shù)定義為：

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

(3)

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

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

(4)

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

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

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

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

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

(a)

(b)

圖4 使用反應(yīng)機(jī)理(a)GRI- Mech 3.0[102]和(b)USC Mech II[103]得到的表1中算例的加速因子γ與最大平均相對(duì)誤差ε的關(guān)系。Xsc(εr)，Xsm(α)和Xscm(εr,α)表示表1所列工況算例，括號(hào)中的數(shù)值εr和α分別表示CoDAC簡(jiǎn)化閾值和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]時(shí)，表1中的算例中用于計(jì)算快方程、慢方程和CoDAC簡(jiǎn)化的壁面時(shí)間[99]

Fig.5Walltimespentonfastequationintegration,slowequationintegrationandCoDACreductionforthecasesinTable1withUSCMechII[99, 103]

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

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

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

3 結(jié) 論

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

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

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