基于多线程并行的动态自适应建表加速算法在气相爆轰模拟中的应用

吴锦涛， 董刚， 栗保明

(南京理工大学 瞬态物理国家重点实验室， 江苏 南京 210094)

0 引言

冲压加速推进技术作为一种新型发射技术，可以有效地提高弹丸的发射初速，为超远程火炮和超高速穿甲效应提供技术支持。冲压推进过程可以以亚爆轰速度、跨爆轰速度和超爆轰速度为模式工作，然而，其中涉及的可燃预混气爆轰过程十分复杂[1]，因此开展相关研究具有重要的理论意义和应用前景。

从数值模拟的角度看，要精确模拟气相爆轰过程，往往要使用包含多组分和多基元反应的详细化学反应机理。然而，在实际问题中由于组分之间的化学反应时间尺度相差悬殊[2-3]，反应源项的求解通常存在严重的刚性现象，导致计算量很大。为提高化学反应的计算效率，传统的处理方法常常涉及到对化学反应机理的简化[4-7]，但简化机理一般依赖于所研究的具体物理问题，缺乏通用性，而且对于某些重要的痕量物质往往无法预测或是严重失真，因而受到不同程度的限制。

1997年，Pope[8]提出了动态自适应建表(ISAT)方法，以减少反应流数值模拟中化学反应源项的计算时间。不同于传统的机理简化方法，该方法采用在内存中建立数据表的方式对化学反应的数据进行存储、查找和取回，通过这些操作来替代原有反应源项的直接积分(DI)运算。由于数据操作花费的时间远小于DI的时间，ISAT方法可以有效地提高化学反应源项的计算效率。在低马赫数定常燃烧问题中，流场内一个步长前后的化学反应热力学状态变化不大，大部分数据可以在表中查找成功并取回，因此ISAT方法在此类问题中具有很高的加速比。例如：在均质压燃(HCCI)发动机[9]中，其化学反应加速比相对于DI能够达到100倍以上；而在混合搅拌发生器(PSRs)[10]中，建表方法可缩短反应计算时间约1/42.

然而，国内外针对ISAT方法在可压缩流的非定常燃烧问题中的应用却报道较少。事实上在这类问题中，考虑到每个时刻流场内热力学状态相对上一时刻变化较大，大量数据可能无法从数据表取回，因此就会有大量节点不断插入表中为后续计算提供可查找的数据，这严重影响了查表效率并增加了表的尺寸，因此ISAT算法的加速效果必定会受到不同程度的限制。Dong等[11]考虑流场存在间断的非定常燃烧特点，将ISAT算法应用到二维气相爆轰波引发和传播的数值模拟中，在0阶取回判据的条件下，其化学加速比最大可达到17.88，最大总加速比可达到4.75. 在Dong等[12]后续研究中，还进一步提出了两种节点删除(全局节点删除和局部节点删除)机制来避免表的过度增长，并针对爆轰波传播的瞬态过程计算对比了两种删除机制下ISAT方法的效率。

上述研究表明，尽管ISAT方法本身通用性好，但其计算效率仍依赖于物理问题。随着近年来大规模并行计算的日益普及，ISAT方法也逐渐应用于并行计算中[13]，但针对可压缩非定常问题，目前尚未提出有效的并行算法。而以往满足可压缩非定常问题的算法则是基于串行计算的[12]，因此，发展基于ISAT的适合于可压缩非定常燃烧问题的有效并行算法十分必要。可压缩气相爆轰波精细结构的捕捉是非定常可压缩反应流计算的重要问题之一。气相爆轰波是化学反应放热支持下的激波，其化学反应特性决定了波阵面反应区的大小并因而影响了爆轰波的传播过程。另一方面，由于爆轰波阵面主要由各类强间断和接触间断所构成，其发展是不稳定的，要想精确捕捉精细结构和不稳定发展过程，必须对计算网格尺寸有十分苛刻的要求，因而计算量巨大。因此，针对这类问题来发展基于ISAT的并行算法并分析其性能(计算精度和计算效率)是非常适合的。

本文以二维气相爆轰波传播的过程为研究对象，在ISAT方法的基础上，建立了基于多线程信息传递接口(MPI)并行架构的化学反应数据动态存储、删除加速算法，考虑了不同的建表方法、建表容差判据和流场计算格式等因素的影响，通过数值模拟分析了并行加速算法的计算精度和计算效率。

1 计算方法

1.1 动态存储、删除的ISAT算法

在含有K种组分的反应体系中，其热力学状态函数所构成的状态空间可表达为

Φ=(φ1,…，φk, …，φK,φK+1,φK+2)=
(Y1,…，Yk,…，YK,T,ρ),

(1)

式中：Yk为组分k的质量分数；T为体系温度；ρ为体系密度。对燃烧过程，流场内每个网格(颗粒)内的热力学状态在规定步长Δt内发生的化学变化可以表达为

(2)

(2)式给出了状态函数积分的求解过程，称为DI，其中S代表化学反应速率源项，由于其具有强刚性和强非线性，因此求解将花费大量CPU时间。

在ISAT算法中，DI可被数据表的操作所取代。计算过程中首先在内存中建立一个空白数据表，并且将初次计算得到的Φ及其对应的变化量ΔΦ存入数据表的一个节点中，对后续计算中的任一变量Φ′，在数据表中先查找与之近似的Φ，并将表中满足容差范围的节点值Φ+ΔΦ取回作为Φ′的计算结果，以避免DI计算。否则，进行(2)式DI并将得到的结果存入数据表的新节点中。由于对表中节点的存储和取回操作所花费的时间远小于DI所花费的时间，这种方法可明显节约计算的CPU时间[8]。对容差范围的判定，本文采用了常数近似[11-12]以替代原始ISAT中的线性近似[8]方法，以进一步节约对数据表中节点的操作时间。同时，设定容差判据εtol以保证待查数据能够落入误差域内，否则将DI结果及其误差域作为一个数据簇存入表中新节点上。事实上，容差判据的大小决定了ISAT对节点数据的容忍范围，直接影响数据查找和取回的效率。

气相爆轰波计算中，不同节点代表的热力学状态不同，因而其被取回利用的次数也不同。为避免数据表占用内存超限，设定最大容量M，计算过程中，一旦表内节点总量超过M，则对数据表进行清理。本文采用平衡二叉树(BBT)结构来构造数据表，以确保对表中节点数据查找、取回、存储以及删除等操作能够高效完成。为了保留表中利用率高的节点，通过非递归中序遍历方法对表内节点进行查找，将表中从未取回过的节点进行删除[12]，以此来控制数据表的大小。

1.2 MPI并行计算中的建表策略

对于包含详细化学反应机理的气相爆轰数值模拟问题，由于所涉及的流场结构比较复杂且对网格分辨率要求很高，普通的串行执行方式已经无法满足计算需求，故本文采用MPI方法对整个计算过程进行多线程并行处理。

MPI并行环境下，ISAT建立热力学数据表的思路可以分为两种：第1种是各个线程分别独自建立数据表，在表所对应的分区内各线程对表中数据同时进行插入、查找和取回等操作，整个计算过程中各数据表彼此互不关联；第2种则是在整个计算区域内建立一个数据表，所有线程在计算过程中共享该表中的数据并依次对其进行操作。

第1种建表方法的优势在于各线程可同时分别对其分区内的数据进行操作，避免了数据操作的顺序等待时间，因此具有更高的计算效率；第2种建表方法通过表中数据在所有分区共享的方式避免了第一种方法中不同数据表内数据重复的不足，但共享数据表的劣势主要表现在对数据进行操作时各线程需要排队等待。

由于不同的建表策略各有优缺点，本文将分析两种建表策略在气相爆轰波传播数值计算中的性能。

2 二维气相爆轰的计算模型

2.1 控制方程和数值格式

针对二维气相爆轰波传播问题，本文采用理想气体模型，控制方程为二维带化学反应的多组分Euler方程：

(3)

计算时采用了分裂算法，将控制方程(3)式中的流动过程和化学反应过程解耦，针对每个步长，先不考虑化学反应，只进行流动过程的计算，再利用经过流动过程的新参数计算化学反应源项。分别采用基于ISAT的并行加速算法和DI方法来求解反应源项。流动过程采用了高精度加权本质非振荡(WENO)格式进行求解，为分析不同数值格式对并行加速算法的计算精度和计算效率的影响，分别采用了5阶WENO格式和9阶WENO格式两种数值格式，它们的具体形式可参见文献[16]. 最后，采用3阶显式Runge-Kutta方法进行时间推进，计算的定步长为1.0×10-9s.

2.2 计算模型

图1为气相爆轰波传播的初始二维计算区域示意图，计算区域Ly×Lx为9 mm×3 mm的长方形。采用均匀正交网格覆盖整个计算域，网格数量为900×300(y×x)，网格尺寸Δx=Δy=0.01 mm. 沿着y方向将流场分成15个相同大小的物理分区，各进程(进程1～进程15，见图1)在对应的分区并行计算。计算区域左边界、右边界采用周期性边界条件，上边界、下边界采用零梯度边界条件，为保证爆轰波稳定传播的距离足够长，令计算区域随着爆轰波传播一起向上方移动。计算结果表明，本文所使用的计算区域大小能够完整地捕捉爆轰波阵面的精细结构，从而达到无限长管道中二维爆轰波稳定传播的数值模拟效果。

计算初始时刻，在计算域y为0～3 mm范围内设定ZND形式的爆轰波结构，波阵面上方的区域(y为3～9 mm)内充满了初始压力p0=1 atm和初始温度T0=300 K、当量比为2∶1的H2、O2预混气，以维持二维爆轰波的自持传播。然而，由于ZND模型是理论上爆轰波传播的稳定解，它本身是一维结构解，作为本文二维计算的初始条件，需要在数值误差扰动下经过很长时间才能发展为不稳定的爆轰波结构。因此，为尽快诱导爆轰波形成不稳定的胞格结构，本文在一维ZND波阵面上方的中心位置处设定一道细长条形的扰动区(初始压力为38p0，初始温度T0不变)，以加速胞格结构的出现。计算开始后，爆轰波阵面以稳定爆速自下向上传播，受扰动作用之后发生强烈失稳并在波面上出现若干三波点，经过一段时间的衰减，爆轰波由过驱动状态逐渐转化为稳定传播状态。

2.3 计算验证

考虑到计算初期由于人为设定扰动区所导致的过驱动状态下爆轰波传播过程不具有物理上的真实性，本文从稳定传播后的某个时刻(步长10 000)开始对计算结果进行统计和处理。结果表明，达到稳定传播状态后的爆轰波传播速度vD=2 812.0 m/s，与Chapmam-Jouget(C-J)理论爆轰波传播速度vC-J=2 843.2 m/s的结果非常接近，证明了本文计算爆轰波传播状态的准确性。

图2给出了稳定传播条件下典型步长(步长20 000)时的爆轰波阵面结构和胞格宽度。从图2(a)中可以明显观察到爆轰波是由入射激波(SW)、马赫杆(MS)、横波(TW)及其它们所汇聚形成的三波点(TP)所组成。经统计，爆轰波在稳定传播过程中，阵面上约存在14个TP，相应的胞格宽度尺寸约为0.428 mm. 图2(b)给出了本文计算得到的爆轰波胞格宽度与文献[17]中实验测得的等当量比H2、O2(与本文相同)爆轰波胞格宽度随着初始压力变化的比较，从中可以发现胞格尺寸随着初始压力的增加而减小，本文在p0=1 atm下的计算结果与文献[17]的变化规律是一致的，证明了本文计算的爆轰波结构是合理的。

3 结果与讨论

本文以等当量比H2、O2气相爆轰波传播过程为研究对象，针对其中的化学反应过程，分析了建表方法、容差判据和数值格式精度等因素对基于并行计算的动态存储、删除ISAT算法性能影响，性能评价指标包括计算精度和计算效率。基于该算法的不同算例参数见表1.

表1 算例及对应参数

表1中，建表方法I代表每个线程各自独立建表的方式，方法II代表所有线程共享一个数据表的建表方式，参见1.2节。其中，方法I中设定独立建表的数据最大容量M=3×106，而方法II中设定表中数据最大容量M=4.5×107，从而保证方法I中所有数据表(15个)的总容量和方法II中的单个数据表容量相等。

3.1 建表方法

图3给出了步长20 000时刻爆轰波阵面附近的流场压力及温度分布图，以比较两种不同建表方法(算例1和算例2)的计算精度，为了便于比较，图3中还给出了DI对应的计算结果。从图3可以看出，两种建表方法的计算结果和DI相比，阵面上的压力和温度分布以及TP个数均保持一致，证明建表方法在本文气相爆轰波计算中的准确度，而波系结构显示出的轻微不规则性是由9阶WENO格式较低的耗散特性引起的，与加速算法本身没有关系。

为了定量比较两种建表方法的计算精度，图4给出了步长20 000时刻加速算法及DI计算所得爆轰波阵面附近的组分OH和H2O2(痕量物质)分布曲线。图4结果显示，两种建表方法(算例1、算例2)对应的计算结果均展现出与DI相似的分布规律，峰值大小及所处的位置十分接近。从痕量物质H2O2的分布曲线可以看出，建表加速算法对小浓度自由基(10-5量级)的捕捉仍具有较高的准确性，这是以简化机理为基础的加速方法所不具备的。

MPI环境下ISAT加速算法的核心在于通过建立和查找数据表的途径取代热力学状态的DI，从而能有效地缩短计算化学反应过程所需的CPU时间，因此，计算效率是检验算法性能的一个重要指标。图5给出了不同建表方法(算例1和算例2)和DI计算化学反应所累计消耗的CPU时间tc曲线及其对应的加速比Rt(定义为DI方法与ISAT方法计算化学反应的CPU时间比)。从图5(a)可以看出：各线程独自建表的ISAT方法(算例1)耗时较DI明显缩短，表明该方法在气相爆轰数值模拟中的优越性；而各线程共享数据表的ISAT方法(算例2)在计算前期(步长16 230之前)整体效率低于第1种建表思路，但总体上节约源项计算的CPU时间依然比较可观。之后，该曲线的斜率突然上升并不断接近DI，甚至在步长17 560之后超过DI的时间，这是因为16 230时刻表中节点总数达到最大容量M后ISAT方法对数据表进行了节点删除操作，接着需要通过DI对表中节点进行补充。该过程在共享数据表的方法中属于顺序操作，相比于算例1同步操作所耗费的时间大幅上升，之后每经历一次表的节点删除，时间曲线都会出现一个对应的拐点(如步长17 490和18 810时刻)，因此该方法在提高计算效率上存在一定的局限性。

由图5(b)可知：加速比曲线定量表明了算例1的优势，经过一小段时间的波动之后最终稳定下来的源项计算时间大约只占DI的40%；算例2由于数据表结构相对于算例1中小容量数据表更复杂，因此其加速比一开始就低于算例1，同时，查找、取回和存储的顺序操作导致其加速比一直呈现不断下降的趋势并在计算后期低于算例1. 考虑到两种并行ISAT建表方法的计算精度相当，显然，各线程独自建表的ISAT方法(算例1)无疑是更优的选择。

3.2 容差判据

容差判据在整个并行ISAT计算方法中，直接决定了数据表中节点的取回精度，分析不同容差对ISAT方法的影响可以兼顾其在气相爆轰波传播问题中的计算精度和计算效率。从图3不同容差加速方法(算例1和算例3)的压力和温度分布可以看出，相比于DI，在容差较大(算例3)的情况下，加速算法对数据的近似取回更多，这种近似作为数值扰动不断累加，导致其不同数据表之间的差异性变得更明显，因而导致了波阵面下方出现大小与分区尺寸相同的低压区(见图3(a)中算例3)，这说明其相对于小容差加速算法(算例1)会损失一部分计算精度。

图4给出了两种不同容差对应的加速算法中组分OH和H2O2(痕量物质)的分布曲线。由图4可以看出，相对于小容差(算例1)算法，大容差(算例3)算法与DI计算结果相比，虽然仍具有相似的规律，但流场内组分含量值特别是峰值大小存在一定的差异，尤其在痕量物质H2O2的分布曲线上，大容差的结果呈现出比较明显的震荡，表明其对波系结构的刻画相对于小容差计算而言不够准确。

为了分析容差放大判据对计算效率的影响，图6给出了算例1和算例3的节点取回率Rn，这里Rn定义为数据表中数据取回次数和数据查找总次数的比值，该值反映了并行ISAT方法的数据查找效率。从图6可以看出：两种误差放大判据对应的加速算法中，线程1～线程14的数据取回率都高于85%，部分线程的取回率甚至在95%以上，证明并行ISAT方法在源项计算中的数据查找效率令人满意；线程15因其所在分区包含爆轰波阵面在内，TP周期性碰撞使得该分区热力学状态变化非常频繁，进而导致二者取回率均有所下降，其中，算例1取回率下降到约60%，算例3取回率下降到约74%，总体上大容差(算例3)算法相对于小容差(算例1)算法的查找效率更高。

为了进一步分析不同容差判据对线程15查找效率的影响，图7给出了该线程一次完整建表过程中节点取回次数的分布曲线。考虑到不同容差下并行ISAT建立热力学数据表的规律不同，为了保证一致性，二者均选取时间步步长在14 100～14 940之间的典型建表过程。结果显示，在包含M=3×106个节点的数据表中，算例1中约有700个节点的取回次数不小于1，而算例3中约有7 000个节点的取回次数不小于1. 剩余节点没有在后续计算过程中被取回，因而视为“无效”节点被删除，这部分调用次数高的节点在经过平衡操作之后仍保留在表中，并在下一次建表过程继续使用，从而保证节点高利用率的同时有效地节约内存空间。对比二者可以发现，节点取回次数的分布规律相同，但算例3由于放大了容差范围使其能够被取回重复利用的节点数量和节点调用总次数均有所增加，再次证明了对容差区间的放大可以有效地提高并行ISAT方法的查找效率。

根据图6和图7对数据查找效率的对比，直观上看，算例3的计算效率应当高于算例1，然而，通过图5 CPU时间和加速比曲线的对比可以看出，除计算初期(步长12 500之前)算例3具有较高的加速比外，后期不同容差的加速算法具有基本相同的加速比。显然，在建立、删除数据表的过程中，部分操作“拖累”了算例3的计算效率。

为探讨是何种操作“拖累”了算例3的计算效率，图8给出了不同容差条件下几个典型线程对应的数据表建立、删除过程。从图8(a)结果可以看出，小容差(算例1)计算中，线程15表的建立、删除过程十分频繁，这是由于爆轰波阵面位于此线程，其化学反应和热力学状态变化剧烈，表的取回率较低(见图6(a))，因此插入表中的节点迅速增加，从而导致了频繁的表建立、删除过程。相比之下，其他线程所处理的流场热力学状态变化较小，表中节点取回率高，表的删除频次较小(线程11)甚至不删除(线程12～线程14)。图8(b)结果则表明，容差判据放大后，线程15的表中节点取回率提高，表的建立、删除频率下降，但由于大容差导致的差异性变大，从上游线程(线程15)计算得到的结果传播到下游线程(线程11～线程14)后，会导致下游线程取回率下降(见图6(b))，因此这些线程的建立、删除频率有所提高，进而使得这些线程中表的建立、删除同步性显著下降。因此，从计算效率来看，容差的区别并不能补偿同步性变差所带来的影响，故算例1和算例3的加速比之间并没有显著差别。

图8计算结果也进一步说明，各线程数据操作的同步性也是影响并行算法计算效率的重要原因之一，发展具有同步性的高效算法也是今后需要进一步开展的工作。

3.3 数值格式精度

计算的数值格式可以影响流场分布，进而会改变反应体系的热力学状态分布，这也会影响加速算法的计算性能。本节采用了5阶WENO计算格式替代前面的9阶WENO计算格式(算例4)，以此分析数值格式的精度对加速算法影响。

图9给出了步长20 000时刻爆轰波阵面附近流场的压力和温度分布，其中，为了便于比较，DI的计算结果也是采用5阶WENO格式计算的。结果显示，采用5阶WENO格式，所有计算结果中三波结构分布得更整齐，爆轰波稳定传播之后阵面上的TP结构为12个，相比于9阶WENO格式计算的TP个数(图3中为14个)有所减少，这说明数值格式对计算结果有明显影响。从图9的5阶WENO格式下算例4与DI计算结果的比较来看，并行ISAT加速算法的结果仍能与DI的结果保持一致，尤其是波阵面及波系结构的形态与DI的结果吻合得均很好。

图10给出了5阶WENO格式下ISAT和DI计算的组分OH和H2O2分布曲线。从图10可以看出，由于爆轰波阵面分布和波系结构更加规则，并行ISAT算法得到的结果与DI吻合很好，即便对浓度很小的痕量物质H2O2，二者在分布规律和组分含量的峰值大小、位置上均吻合得很好(见图10(b))，再次表明了并行ISAT加速算法在计算精度上的优势。

图11为5阶WENO格式下算例4的加速比曲线。图11结果显示，随着爆轰波的传播，加速比曲线经过一段时间的波动后，最终稳定在约2.2倍左右，这与9阶WENO格式下算例1稳定后的加速比非常接近，说明并行ISAT加速算法的计算效率对数值格式的依赖性不强。

综合以上讨论可知，采用9阶WENO格式的算例1和5阶WENO格式的算例4在计算精度和计算效率上均表现出较好的性能。表2总结了计算截止时刻(步长20 000)两个算例各项计算时间以及加速比的情况，并与相应计算格式下的DI计算结果进行了比较，其中，总时间ttot定义为整个计算花费的总CPU时间，总加速比Rtot定义为DI计算总CPU时间与加速算法计算总CPU时间的比值，φc表示化学反应计算花费的CPU时间在总CPU时间中的占比。

表2 不同数值格式下的计算效率统计

从表2中的数据可以看出，两种数值格式下，加速算法均能有效地提高整体计算效率，其总加速比都维持在1.5倍以上。相对于5阶WENO格式而言，9阶WENO下的并行加速算法可以将整个计算过程中化学反应计算时间的占比从DI方法中的58.5%降低至36.7%，表明本文提出的并行算法在高精度格式下的优势更明显。

4 结论

本文在并行计算的基础上，针对非定常可压缩反应流中化学反应源项计算量巨大的特点，提出了各线程独自建表和所有线程共有单表的两种基于MPI环境的多线程并行策略，并成功地将其与带节点删除机制的ISAT方法相结合，取代原始的DI运算，以在不损失计算精度的同时提高计算效率。为分析不同因素对计算结果的影响，本文从建表方式、容差范围和数值计算格式3个方面给出了并行加速算法在反应性瞬态流场中的计算精度和计算效率。所得结论如下：

1) 相比于DI方法，各线程独自建表的并行加速方法在保证计算精度的同时，化学反应的计算加速比可达到2.17～2.43.

2) 本文提出的两种加速方法计算精度与DI方法相比几乎一致，但各线程独自建表的算法具有更高的计算效率。

3) 不同容差范围下并行ISAT算法的数据查找效率不同，但小容差下的加速算法在保证更高精度的同时，仍具有和大容差条件下相当的计算效率。

4) 不同数值格式下，并行加速算法在保证计算精度的基础上，具有大致相当的加速比，表明本文算法的计算效率对数值格式并不敏感。

基于并行ISAT方法的加速算法不依赖于具体的物理问题，有很好的可移植性，为数值研究冲压加速推进技术中的气相爆轰过程提供了一种高效的计算手段。需要指出的是，针对其他不同类型的非定常可压缩反应流问题，本文所提出的并行加速算法性能仍有待于进一步考察。同时，后续工作将拓展现有的建表策略，以进一步提高不同线程表操作之间的平衡性和同步性。此外，还要进一步研究并行加速算法在非均匀网格或非均匀物理分区下的计算精度和计算效率。

参考文献(References)

[1] 金志明. 高速推进内弹道学[M]. 北京: 国防工业出版社, 2001.

JIN Zhi-ming. High speed internal ballistics [M]. Beijing: National Defense Industry Press, 2001. (in Chinese)

[2] Blasi J M, Kee R J. In situ adaptive tabulation (ISAT) to accelerate transient computational fluid dynamics with complex heterogeneous chemical kinetics[J]. Computers & Chemical Engineering, 2016, 84(1): 36-42.

[3] 胡宗民, 牟乾辉, 张德良,等. 爆轰波在弯管内传播过程数值分析[J]. 计算物理, 2004, 21(5):408-414.

HU Zong-min, MU Qian-hui, ZHANG De-liang, et al. Numerical simulation of gaseous detonation wave propagation through bends with a detailed chemical reaction model [J]. Chinese Journal of Computational Physics, 2004, 21(5): 408-414. (in Chinese)

[4] Oran E S, Gamezo V N. Origin of the deflagration-to-detonation transition in gas-phase combustion[J]. Combustion and Flame, 2007, 148(1): 4-47.

[5] Law C K. Comprehensive description of chemistry in combustion modeling [J]. Combustion Science and Technology, 2005, 177(5): 845-870.

[6] Massias A, Diamantis D, Mastorakos E, et al. An algorithm for the construction of global reduced mechanisms with CSP data[J]. Combustion and Flame, 1999, 117(4): 685-708.

[7] Nafe J, Maas U. Hierarchical generation of ILDMs of higher hydrocarbons[J]. Combustion and Flame, 2006, 135(12): 17-26.

[8] Pope S B. Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation[J]. Combustion Theory and Modelling, 1997, 1(1): 41-63.

[9] Zhou L, Wei H. Chemistry acceleration with tabulated dynamic adaptive chemistry in a realistic engine with a primary reference fuel [J]. Fuel, 2016, 171(1): 186-194.

[10] Adhikari S, Sayre A, Chandy A J. In situ adaptive tabulation (ISAT) for combustion chemistry in a network of perfectly stirred reactors (PSRs) [J]. Computers & Chemical Engineering, 2017, 97(1): 124-134.

[11] Dong G, Fan B C, Chen Y L. Acceleration of chemistry computations in two-dimensional detonation induced by shock focusing using reduced ISAT [J]. Combustion Theory and Modeling, 2007, 11(5): 823-837.

[12] Dong G, Fan B C. Chemistry acceleration modeling of detonation based on the dynamical storage/deletion algorithm [J]. Combustion Science and Technology, 2009, 181(9): 1207-1216.

[13] Lu L, Lantz S R, Ren Z, et al, Computationally efficient implementation of combustion chemistry in parallel PDF calculations[J]. Journal of Computational Physics, 2009, 228(15): 5490-5525.

[14] Mcbride B J, Gordon S. Computer program for calculating and fitting thermodynamic functions [M]. Washington, US: NASA， 1992.

[15] Kee R J, Grcar J F, Smooke M D, et al. PREMIX:a Fortran program for modeling steady laminar one-dimensional premixed flames, SAND85-8240 [M].Albuquerque, NM, US: Sandia National Laboratories, 1985.

[16] Shi J, Zhang Y T, Shu C W. Resolution of high order WENO schemes for complicated flow structures[J]. Journal of Computational Physics, 2003, 186(2): 690-696.

[17] Voitsekhovskii B V, Mitrofanov V V, Topchiyan M E. Structure of the detonation front in gases (survey) [J]. Combustion, Explosion and Shock Waves, 1969, 5(3): 267-273.