文章快速检索     高级检索
  气体物理  2020, Vol. 5 Issue (1): 56-65   DOI: 10.19527/j.cnki.2096-1642.0756
0

引用本文  

苏月涵, 沈小东, 傅嘉琛, 等. 激波诱导火焰界面失稳的并行化学加速计算[J]. 气体物理, 2020, 5(1): 56-65.
Su Y H, Shen X D, Fu J C, et al. Parallel chemistry acceleration computations of flame interface instability induced by shock wave[J]. Physics of Gases, 2020, 5(1): 56-65.

基金项目

南京理工大学本科生科研训练"百千万"计划省级项目(201810288032X);国家自然科学基金项目(11872213)

第一作者简介

苏月涵(1998-)  男, 本科, 主要研究方向为可压缩反应流的数值模拟.E-mail:852388160@qq.com

通信作者简介

董刚(1970-)  男, 研究员, 博导, 主要研究方向为可压缩反应流的数值模拟.E-mail:dgvehicle@yahoo.com

文章历史

收稿日期:2019-04-27
修回日期:2019-07-01
激波诱导火焰界面失稳的并行化学加速计算
苏月涵 , 沈小东 , 傅嘉琛 , 董刚     
南京理工大学瞬态物理重点实验室,江苏南京 210094
摘要:在带有详细化学反应机理的可压缩反应流数值模拟中,化学反应源项的计算会极大增加计算时间,基于建表技术的化学加速算法可以通过查找数据表中的数据来替代化学反应计算,从而有效提高计算效率,但数据表尺寸的过度增长会导致计算的中断.文章提出了基于两种数据表容量控制策略的并行动态存储/删除算法,并在激波诱导火焰界面失稳的数值模拟中进行了应用,以考察算法的性能.两种数据表容量控制策略分别为单表容量(Msin)控制和总表容量(Mtot)控制,当单个数据表尺寸达到Msin或总数据表尺寸达到Mtot时,对数据表进行节点删除,以保证计算的正常进行.计算结果表明,文章提出的基于表容量控制的并行加速算法,其计算准确度和计算效率之间存在关联,具有较好计算准确度算例显示了较高的计算效率.在不同的MsinMtot条件下,计算的化学加速比在2.73~3.93之间.两种表控策略的组合影响了数据表删除的频率和删除之间的同步性,当数据表删除频率小、删除同步性强时,化学加速比要更高.
关键词反应流    化学加速建表算法    化学加速比    表容量    火焰界面失稳    
Parallel Chemistry Acceleration Computations of Flame Interface Instability Induced by Shock Wave
SU Yue-han , SHEN Xiao-dong , FU Jia-chen , DONG Gang     
Key Laboratory of Transient Physics, Nanjing University of Science and Technology, Nanjing 210094, China
Abstract: In the numerical simulations of compressible reacting flow with detailed reaction mechanism, the computations of chemical source term can spend huge CPU time. A chemistry acceleration algorithm based on tabulation technique can greatly improve the computational efficiency, through replacing the direct computation of chemical reaction by searching the data in a table. However, excessive growth of the table size can give rise to a breakoff of the computations. This study presented the parallel dynamic storage/deletion algorithm based on two table size control strategies. To test the computational performance of the algorithm, the simulations of flame interface instability induced by a shock wave have been carried out. Both table size control strategies in the algorithm include the single table size (Msin) control and the total table size (Mtot) control. Once the size of single table reached Msin or the total size of the tables reached Mtot, a nodal deletion was performed for all the tables in order to maintain the table size. The results show that there exists relationship between computational accuracy and efficiency for the algorithm based on the table size control. The better computational efficiency corresponds to the higher computational accuracy. The chemical speedup ratio resulted from simulations ranges in 2.73~3.93 when Msin and Mtot vary. The combination of both table size control strategies plays an important role in affecting the deletion frequency of tables and the synchronization among the table deletion events. The chemical speedup ratio is larger when smaller deletion frequency of tables and better synchronization occur.
Key words: reacting flow    chemistry acceleration tabulation algorithm    chemical speedup ratio    table size    flame interface instability    
引言

在包含详细化学反应机理的燃烧过程的数值模拟中, 化学反应涉及的组分特征时间尺度变化极大[1], 从而使化学反应源项呈强刚性, 这极大地增加了计算花费.另外, 为刻画可压缩燃烧场中的流动间断(包括强间断和接触间断), 还需要进一步提高计算格式和/或网格的精度, 这也会增加计算化学反应源项所消耗的时间.因此, 对可压缩流场中的燃烧问题进行数值模拟时, 需要在保证计算准确度的前提下, 尽可能提高计算效率, 尤其是化学反应源项的计算效率.

传统的提高化学反应计算效率的方法主要是对化学反应机理进行简化, 由此产生了一系列的简化方法[2-7], 但是简化机理方法的不足之处在于其通用性变差, 且对某些感兴趣的痕量物质不能很好地进行预测. Pope[8]提出了一种动态自适应建表(in situ adaptive tabulation, ISAT)的方法, 以减少反应流数值模拟中化学反应源项的计算时间.该方法通过在内存中建立数据表进行数据的问询和取回等操作来替代原有的化学反应源项的直接积分(direct integration, DI)运算, 由于数据操作花费的时间远小于DI花费的时间, 因此ISAT方法可以有效提高反应源项的计算效率[9-10]. 2009年, Lu等进一步发展了ISAT方法的并行版本[11]并开发了相应的软件包x2f_mpi[12].

然而, 对热力学状态随时空变化较大的燃烧问题, 数据表的取回效率会有所下降, 这会使数据表的尺寸随计算过程迅速增长并超过计算机内存限制, 从而使计算中断.为解决这一问题, Dong等[13]在原始ISAT算法引进了节点删除机制, 得到了基于ISAT的动态存储/删除算法.该算法通过规定数据表的容量阈值来避免数据表中节点数量的过度增长, 并在气相爆轰问题中取得了3.87的化学加速比.接下来他们在动态存储/删除算法的基础上发展了并行计算的版本[14], 并将其运用到二维H2+ 2O2混合物的气相爆轰数值模拟中, 通过考察数据表尺寸对计算效率的影响, 发现并行条件下最大化学加速比可达5.52.

尽管并行版的动态存储/删除算法在气相爆轰问题中得到了检验, 但在非稳态燃烧过程中的计算性能并不清楚.就算法本身而言, 其在原理上是不依赖于燃烧问题的, 但算法的计算准确度和计算效率却依赖于问题的不同.因此构造合适的非稳态燃烧问题, 可以有效地检验并行动态存储/删除算法的有效性和通用性.为此, 本文构造了一种激波诱导的气相反应性界面(火焰界面)失稳的二维非稳态燃烧过程, 采用并行动态存储/删除算法对该过程进行了数值模拟, 考察了两种不同的控制数据表尺寸方法对算法计算效率的影响, 对影响计算效率的机理进行了分析.

1 数理模型 1.1 控制方程

对于激波诱导气相火焰界面不稳定问题, 考虑到气体来流速度比较大, 因此可以忽略其流场内分子的输运特性, 视其为理想气体.因而控制方程为二维带化学反应的多组分Euler方程, 如下

$ \frac{\partial \mathit{\boldsymbol{U}}}{\partial t}+\frac{\partial \mathit{\boldsymbol{F}}}{\partial x}+\frac{\partial \mathit{\boldsymbol{G}}}{\partial y}=\mathit{\boldsymbol{S}} $ (1)

式中

$ \begin{align} &\mathit{\boldsymbol{U=}}\left( \begin{matrix} {{\rho }_{1}} \\ \vdots \\ {{\rho }_{K}} \\ \rho u \\ \rho v \\ E \\ \end{matrix} \right),\mathit{\boldsymbol{F}}=\left( \begin{matrix} {{\rho }_{1}}u \\ \vdots \\ {{\rho }_{K}}u \\ \rho {{u}^{2}}+p \\ \rho uv \\ u\left( p+E \right) \\ \end{matrix} \right) \\ &\mathit{\boldsymbol{G=}}\left( \begin{matrix} {{\rho }_{1}}v \\ \vdots \\ {{\rho }_{K}}v \\ \rho uv \\ \rho {{v}^{2}}+p \\ v\left( p+E \right) \\ \end{matrix} \right),\mathit{\boldsymbol{S}}=\left( \begin{matrix} {{\omega }_{1}}+W \\ \vdots \\ {{\omega }_{K}}+W \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right) \\ \end{align} $

以上各式中, ρ为密度, $\rho =\sum{_{k=1}^{K}{{\rho }_{k}},{{\rho }_{k}}=\rho }{{Y}_{k}},{{Y}_{k}}$为组分k的质量分数; u, v分别代表x, y两个方向的速度分量; E为体系总能量,

$ \begin{align} &E=\rho \sum\nolimits_{k=1}^{K}{{{Y}_{k}}({{h}_{k}}-RT/{{W}_{k}})+\rho ({{u}^{2}}+{{v}^{2}})}/2 \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{h}_{k}}=\int{_{T}{{c}_{\text{p}k}}\text{d}T+h_{k}^{0}} \\ \end{align} $

cpk为组分k的定压比热, hk0为组分k的标准生成焓, wk为组分k的分子量, R为通用气体常数; p为压力; ωk是组分k的反应速率.

1.2 计算方法

针对式(1), 本文采用分裂算法进行求解, 即在每一个时间步长内, 先不考虑化学反应, 单独进行流动过程计算, 再利用更新之后的流动参数计算化学反应源项.对流动过程(式(1)中的F/∂xG/∂y), 本文采用全局Lax-Friedrich矢通量分裂格式结合9阶WENO格式[15]进行求解; 对化学反应过程(式(1)中的S), 本文则采用并行的动态存储/删除算法进行计算, 其具体描述参见第2节.此外, 为考察并行动态存储/删除算法的性能, 本文还利用基于隐式Gear格式的VODE软件包[16]进行直接积分(DI)来求解化学反应源项.最后, 式(1)中的时间导数项(U/∂t)采用了3阶Runge-Kutta法进行求解, 并取足够小的固定时间步长(Δt=1×10-8 s)确保时间推进的数值稳定性.

目前已发展了基于DI求解方法的上述算法的计算程序, 并已成功用于激波诱导的反应性Richtmyer-Meshkov(RM)不稳定问题的数值模拟[17-18], 程序代码的可靠性已得到了充分的检验.

1.3 计算构型

图 1为激波诱导火焰界面不稳定的二维计算区域示意图, 计算区域为L×H=0.1 m×0.48 m的长方形.图中灰色区域为H2+O2预混可燃气体, 白色区域为H2+0.05N2轻质气, 两种介质由于密度的不同而形成一气体界面(GI).初始条件下, 两种气体均以速度u0=2 250 m/s由下而上运动, 此时一道斜激波(图 1红线所示)与运动气体及其界面相互作用, 这一相互作用过程会使得斜激波和界面均产生偏转.本文给定界面的偏转角为δ=30°, 波前气体温度T0=298.15 K和压力p0=20 atm, 并令初始时斜激波角度为β1, 界面作用后的角度β2, 则由斜激波关系式可以分别得出β1=68.5°和β2=42.4°.当气体界面与斜激波作用后, 由于波后界面两侧气体的密度不同, 从而导致波后界面两侧气体的运动速度不同, 随着时间的进行, 界面逐渐发展为Kelvin-Helmholtz (K-H)不稳定状态, 从而进一步促进界面两侧H2+O2预混气和H2+0.05N2轻质气的混合及燃烧.波前H2+O2预混气和H2+0.05N2轻质气的充填比设置为l1:l2=1:3.

图 1 二维计算区域示意图 Fig.1 2D computational domain

本文使用了H2/O2/N2体系的详细化学反应机理[19], 该机理包含9种组分和19个基元反应, 可以很好地描述高压条件下的燃烧过程.此外, 采用均匀正方形网格覆盖整个计算域, 网格尺寸Δxy=0.4 mm, 网格数量为250×1 200(L×H).计算域的左边界为刚性无滑移壁面条件, 其余壁面均为零梯度边界条件.考虑到计算初期界面与激波相互作用时会由于RM不稳定现象产生一个不真实的大尺度涡旋结构, 因此, 为保证计算的物理真实性, 计算初期首先将左边界设置为零梯度条件并且关闭化学反应项, 待大尺度涡旋结构流出左边界, 再将其设定为壁面条件并同时打开化学反应项, 此时的流场为真实流场, 即可以展开后续计算.此外, 网格依赖性测试的结果[20]表明本文采用的计算模型和网格尺寸是可靠的.

2 并行动态存储/删除算法 2.1 加速原理

并行动态存储/删除算法是基于ISAT方法[8]的并行数据表操作算法.对含有K种组分的化学热力学体系, 其热力学状态构成的状态空间可表示为

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\mathit{ = }\left({{Y_1}, \cdots, {Y_K}, T, \rho } \right) = \left({{\phi _1}, \cdots, {\phi _K}, {\phi _{K + 1}}, {\phi _{K + 2}}} \right) $ (2)

在燃烧过程中, 流场内每个网格的组分浓度在规定时间步长Δt内发生的化学变化为

$ \phi _{k}^{1}=\phi _{k}^{0}+\int{_{\Delta t}{{s}_{k}}\text{d}t\ \ \left( k=1,\cdots ,K \right)} $ (3)

式中, sk为组分k的化学反应源项; 同时, 可以根据能量守恒律采用积分方法求解温度φK+1.上述状态函数的积分求解过程被称为直接积分(DI).由于式(3)中sk具有很强的刚性和非线性, DI过程会耗费大量计算时间.相反, ISAT技术是通过建立数据表并对表中数据进行问询/取回, 以此代替式(3)的DI积分过程, 从而获得化学反应的近似解.

考虑流场中一个待问询的化学热力学状态矢量φq, 如果数据表中存在某一状态φ(j)的节点与其足够接近, 则φq所对应的化学变化R (φq)可由存储在数据表中的φ(j)和对应的R (φ(j))通过Taylor级数展开获得.通过合理控制误差范围, 可以采用如下的零阶近似方法获得热力学状态的变化

$ \mathit{\boldsymbol{R}}\left({{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^q}} \right) = \mathit{\boldsymbol{R}}\left({{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{(j)}}} \right) + o(|\Delta \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}|) $ (4)

因此, 定义一个准确的误差范围相当重要.本文定义的数据表节点φ(j)的误差范围如下

$ {e^ \pm }({\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\left(j \right)}}) = \left\{ {e_1^ \pm \left({{\phi ^{(j)}}} \right), \cdots, e_k^ \pm ({\phi ^{(j)}}), \cdots, e_{K + 2}^ \pm \left({{\phi ^{(j)}}} \right)} \right\} $ (5)

式中, ek±代表了状态参量中第k个分量的误差, 其中上标"+"和"-"分别代表φ(j)的误差上限和下限.由于式(5)表达的误差上、下限之间的误差范围共有K+2维, 因此e±的定义就构成了一个超矩形误差域, 若状态矢量φq落在φ(j)所对应的超椭误差域内, 则R (φq)通过式(4)获得R (φ(j))的线性近似值.由于数据表的问询/取回所耗费的时间远小于直接积分求解的时间, 因此可以起到加速的作用.

另外, 由于计算包含间断的瞬态可压缩反应流问题, 流场内的热力学状态随时空变化很大, 表中的大量数据无法被取回.一旦表中数据无法取回, 则流场中的化学变化将通过DI的方法进行计算(式(3)), 同时将计算得到的结果作为新的节点存入数据表中, 这样就会造成数据表内节点数量急剧上升, 甚至突破计算机内存使用上限, 导致计算中断.因此, Dong等[13]在上述零阶近似的ISAT方法中引入了数据表的节点删除机制, 有效地控制了表的快速生长过程.

2.2 并行建表策略

对于并行计算问题, 流场常常划分为多个计算子区, 如果建立数据表与子区的一一映射关系, 则需要多个数据表.若每个数据表所对应的流场状态不同, 那么每个数据表对流场状态的操作也不同, 因而直接影响了表操作之间的负载平衡.因此, 提出合适的并行建表策略是并行化学反应计算加速的关键.

一种简单直观的方法是建立计算域的每个并行子区的数据表, 化学过程的处理仅在自己的表中执行, 不同子区的数据表之间不存在数据交流, 这种方法被称为纯局部处理(purely local processing, PLP)法[11].显然, 不同分区流场变化的差异将导致数据表操作之间的差异, 从而影响并行计算的整体效率.

由于计算子区与数据表并不存在绝对的物理映射关系, 因此在计算机内存中, 可以指定某个数据表处理某个流场区域, 这就为提高表操作的负载平衡性提供了可能.能够实现的一种方法是由作者提出的转置处理(transposed processing, TP)的建表方法[21].在TP方法中, 数据表以转置的方式和不同的子区内相同位置的数据块相对应, 数据表数量与分区数量相一致.

针对本文斜激波与反应性界面相互作用下的PLP和TP两种建表方式可如图 2所示.本文所有算例划分的计算子区均为25个, 图 2(a)中, 相同的矩形子区(S1-S25)沿y方向划分, 每个子区分别建立自己的数据表(T1-T25), 计算过程中, 数据表仅处理自己对应的子区内的化学反应; 图 2(b)中, 流场子区与数据表的排布方式互为转置, 每个数据表仅处理所有子区相同位置的数据簇, 处理结果返回相应的子区位置.

图 2 并行建表策略示意图 Fig.2 Parallel tabulation strategy diagram

在建表策略中, 每个数据表均采用了平衡二叉树的方式来组织数据结构, 对数据表的操作主要包括了问询、取回、插入等过程.

2.3 并行删表策略

如前所述, 在包含流场间断的瞬态可压缩反应流问题中, 如果流场内的热力学状态随时空变化很大, 则数据表能够被取回的数据的次数将减少, 从而导致需要插入表中的节点数量急剧上升, 甚至超出计算机的内存上限导致计算中断.因此, 对数据表中的数据节点做适当的删除十分必要.对单表单一计算域(非并行)问题而言, 可设定数据表的节点最大容量M, 当数据表节点数量超过M时, 则对表中某些节点予以删除, 从而达到对数据表尺寸进行控制的目的.

对并行建表方法, 删表策略则比较复杂.除了设定每个数据表的节点最大数量(容量)Msin外, 还可以设定所有数据表总的最大数量Mtot.这两个容量的不同可以显著影响表操作过程, 进而对化学反应计算效率产生影响.为此, 本文设定了不同的MtotMsin来考察删表策略对计算效率的影响, 由此形成的算例见表 1.

下载CSV 表 1 算例及表容量 Tab.1 Cases and table sizes

在删表策略中, 每个数据表均采用了节点删除的操作方式, 即对从未被取回过的节点予以删除, 而表中其他节点则予以保留.

3 结果与讨论 3.1 建表策略计算性能

算法的计算性能主要包括计算效率和计算准确度.本节以表 1中的case10为基本算例, 首先讨论了不同建表策略(PLP和TP)下的算法计算效率, 然后分析了算法的计算准确度.为便于比较, 本文还采用DI方法对算例进行了数值模拟和结果分析.

3.1.1 计算效率

为评估算法的计算效率, 本文采用了如下化学反应加速比Rt来定量表征

$ {{R}_{\text{t}}}=\frac{{{t}_{\text{DI}}}}{{{t}_{\text{algorithm}}}} $ (6)

式中, tDI为采用DI算法计算化学反应所耗费的CPU时间, talgorithm为采用本文算法计算化学反应所耗费的CPU时间.由此可知, 采用算法耗费的时间越少, 则加速比越大, 计算效率越高. 图 3给出了两种并行建表策略下的PLP算法和TP算法计算得到的Rt随计算过程的变化.可以发现, 两种算法下Rt都经历了快速上升、达到最大值后缓慢下降并趋于平缓的过程, 两种算法的Rt变化差异不大.在计算终点时刻(时间步长为353 000), TP和PLP算法的加速比分别是3.73和3.39, 这说明两种策略均能有效提高计算效率, TP策略较PLP策略显示了更好的加速效果.

图 3 TP和PLP加速比曲线 Fig.3 Acceleration ratio curves of TP and PLP
3.1.2 计算准确度

图 4给出了时间步长为353 000时, 两种算法(PLP和TP)以及DI方法计算得到的流场分布.结果表明, 燃烧场存在两个明显特征:首先, 在预混气体区(靠近左边界处), 燃烧是以爆轰波(DW)的形态发展的, 此时, 反应阵面(以图 4(a)温度分布表征)和激波阵面(以图 4(b)密度梯度表征)强烈耦合; 其次, 气体界面处的燃烧是以非预混燃烧的方式进行的, 火焰阵面(non-premixed flame, NPF)具有多尺度的涡旋结构.从图中的结果来看, 两种并行建表策略下的算法与DI方法得到的流场温度、波系结构和涡旋结构相差不大.需要说明的是, 在分层气体的界面附近存在沿下游逐渐增大的涡系结构, 不同算法得到的这些涡结构在细节上仍然存在着一定差异, 这种差异是由界面处K-H不稳定的内在性质所决定的, 不可能完全消除.

图 4 不同算法得到的流场温度和密度梯度分布 Fig.4 Distributions of temperature and density gradient in flow field resulted from different algorithms

为定量考察算法的计算准确度, 图 5给出了两种算法分别在左边界附近(x=0.004 m)和H2+O2预混气充填区中间位置(x=0.037 5 m, 即l2/2处), 流场温度、压力和典型活性组分O的质量分数沿y方向的变化曲线. 图 5(a)中, 由于爆轰波会在左边壁面发生反射, 导致左边界附近热力学状态变化剧烈, 因此, 两种算法计算结果与DI结果表现出一定差异, 其中TP算法较PLP算法高估了计算的温度和低估了计算的O质量分数; 图 5(b)中, 由于反应区远离壁面, 其热力学状态的变化变缓, 因此两种算法与DI方法的计算结果更加接近.

图 5 不同位置处不同算法计算的温度、压力和O质量分数曲线 Fig.5 Profiles of temperature, pressure and mass fraction of O at different x locations

图 4图 5的结果表明, 尽管本文给出的两种算法与DI方法在热力学状态变化较为剧烈的区域存在一定差异, 但计算结果的幅值变化规律仍然是合理的.

3.2 删表策略计算性能

由上节可知, 由于终了时刻TP算法的计算效率略高于PLP算法的计算效率, 因此本节首先以TP算法为基础考察表 1中不同并行删表策略下的化学反应计算效率, 然后再分析不同删表策略对计算准确度的影响.

3.2.1 计算效率

为综合考虑不同的表容量阈值(MsinMtot)对化学加速比的影响, 本文定义了一个表容量因子, 表达如下

$ {{p}_{\text{M}}}=\frac{n\times {{M}_{\text{sin}}}}{{{M}_{\text{tot}}}} $ (7)

式中, n代表分区数量. pM代表了实际数据表容量总和(n×Msin)与总数据表容量(Mtot)的比值, 该值越大, 说明单个数据表容量表越大, 此时表的删除主要受总表容量(Mtot)控制; 反之该值越小, 则删除受单表容量(Msin)的控制.

图 6给出了表 1所有算例的pMRt之间的关系, 图中数字代表了表 1中典型算例序号.通过观察可以发现, 化学加速比Rt随着pM的增加呈现出迅速增加然后再缓慢减小的趋势. Rt较高的算例主要分布在pM处于2~6之间的范围, 其中最高化学加速比出现在算例11(pM=4.2, Rt=3.93);而化学加速比最小的为算例8(pM=1.0, Rt=2.73).上述结果表明, 设定不同的数据表容量, 对加速比有显著的影响.

图 6 RtpM的关系曲线(图中数字代表表 1中的算例序号) Fig.6 Relationship between Rt and pM (the number in figure represents the sequence number of the cases in Table 1)

为进一步探讨数据表容量对加速性能的影响机制, 本文选择了4个典型的算例(算例8, 7, 2和11, 见图 6中带编号的算例), 分别给出了每个算例所有数据表尺寸随计算的变化过程, 见图 7(a)~(d).图中的结果显示, 对每一个数据表, 在数值模拟过程中其尺寸都经历了锯齿形的变化, 这意味着表的尺寸在计算过程中会不断增长, 当达到MsinMtot容量阈值后, 数据表发生节点删除, 表尺寸迅速降低然后开始新的增长, 这一周期性的删除/增长行为导致了数据表尺寸的锯齿形变化.

图 7 不同算例的数据表尺寸增长/删除变化规律 Fig.7 Variations of growth/deletion of table size for different cases

图 7(a)对应的是加速比最低的算例8(Rt=2.73)的表删除/增长行为, 可以发现由于该算例单表容量过小(Msin=3×106, 见表 1), 即使总表容量很大(Mtot=72×106), 每个数据表还是很容易增长到单表容量上限, 因此使得数据表的删除变得频繁; 另外, 由于每个表增长快慢也不相同, 因此表的删除也不同步.这种频繁和非同步删除数据表的行为导致了最低的加速比. 图 7(b)(c)分别对应了两种中等加速比的算例7和算例2(Rt=3.08和3.39).在图 7(b)中, 算例7单表容量和总表容量都很大(Msin=9×106, Mtot=126×106), 尤其是后者, 因此表的删除/增长行为主要受单表容量控制, 因此较大的单表容量使得数据表的删除频率明显下降, 这是导致算例7比算例8加速比提高的主要原因, 但由于表增长快慢的不同, 删表并不同步, 因此其加速比提高得并不显著. 图 7(c)代表了与图 7(b)相反的另一种情况, 由于此时总表容量下降至Mtot=36×106, 因此, 数据表的删除主要受总表容量控制, 当数据表尺寸的总和达到总表容量时, 所有数据表都被删除, 这意味着表删除的同步性得到增强, 但由于总表容量小, 表删除仍过于频繁, 因此该算例的加速比提高得也不显著. 图 7(d)的表删除/增长行为代表了计算加速比最好的情况, 即较大的单表容量(Msin=12×106)和适当的总表容量(Mtot=72×106)保证了较低的删表频率和较强的同步性, 因此化学加速比达到最大.

3.2.2 计算准确度

为进一步考察不同数据表容量对计算效率的影响, 仍取算例8, 7, 2和11的结果进行分析. 图 8给出了时间步长为353 000时, 4个算例在x=0.004 m处压力的相对误差(er)的变化曲线.可以发现, 4个算例的压力相对误差在y=0.14~0.22 m范围内均存在明显波动.由于该范围对应了爆轰波在壁面反射的波后区(见图 4), 其热力学状态沿空间变化较大, 从而导致算法的预测出现较大误差.进一步可以发现, 算例7和11的相对误差要小于算例2和8.前两个算例对应的是大尺寸数据表的情况, 而后两个算例对应的则是小尺寸数据表的情况, 如图 7所示, 这意味着算法的计算准确度和计算效率之间存在一定关系.

图 8 基于TP算法的典型算例的压力相对误差, x=0.004 m Fig.8 Relative errors of pressure in typical cases based on TP algorithm, x=0.004 m

图 9进一步给出了本文所有算例在时间步长为353 000时, 位于x=0.004 m处波后沿y变化的压力相对误差绝对值的平均量(${{{\bar{e}}}_{\text{r}}}$)与加速比(Rt)之间的关系.可以发现, 计算准确度和计算效率之间存在显著的关联, 较小的平均相对误差对应了较高的计算加速比.加速比最小的算例8(Rt=2.73), 其平均相对误差为9.02%;而加速比最高的算例11(Rt=3.93), 其平均误差为4.04%.结合图 7进行分析可知, 相对于通过DI计算插入新节点(相当于被问询节点没有取回), 对数据节点的取回(相当于被问询节点成功取回)具有更高的精度.因此当数据表尺寸越大, 表中数据节点就越多, 其成功取回的概率也就越大, 这导致了大数据表的算例具有更高的计算准确度.

图 9 平均相对误差与加速比的关系 Fig.9 Average relative error vs. speedup ratio

图 9的结果还表明, 本文使用的大多数算例, 在x=0.004 m处的相对误差平均值均小于5%, 由于该处在整个流场中热力学状态变化最为剧烈, 故对算法性能的要求也最为苛刻, 图 9的结果表明本文算法的计算准确度是满足要求的.

4 结论

本文在动态自适应建表(ISAT)的基础上, 采用提出的PLP和TP两种并行建立数据表的策略, 和两种基于表容量控制(MsinMtot)的并行删除数据表的策略, 形成了并行动态存储/删除算法.数值模拟了激波诱导的反应性界面失稳过程, 以考察算法的计算性能.得到的主要结论如下:

(1) 针对激波诱导的燃烧不稳定问题, 本文采用的两种并行化学加速建表算法(PLP, TP)均能够在满足计算准确度的同时, 提高计算效率.计算准确度和计算效率存在关联, 较好的计算准确度对应了较高的计算效率.

(2) 在TP算法下, 不同的数据表容量尺寸导致计算的化学加速比明显不同, 其中, 较大单表容量和适当的总表容量条件下化学反应加速比最高, 为3.93;而过小的单表容量条件下的计算加速比最小, 为2.73.

(3) 数据表的删除频率和删除的同步性决定了算法的计算效率.当数据表的删除由大容量单表尺寸控制时, 所有数据表的删除频率会减小、同步性增强, 因此计算的化学加速比显著提高; 而当数据表的删除由小容量单表尺寸控制时, 数据表的删除频率会明显增加、同步性也会变差, 因而使得化学加速比降低.

参考文献
[1]
Warnatz J, Mass U, Dibble R W. Combustion: physical and chemical fundamentals, modeling and simulation, experiments, pollutant formation[M]. 4th ed. Berlin: Springer, 2006.
[2]
Law C K. Comprehensive description of chemistry in combustion modeling[J]. Combustion Science and Technolo-gy, 2005, 177(5/6): 845-870.
[3]
Lu T F, Law C K. A directed relation graph method for mechanism reduction[J]. Proceedings of the Combustion Institute, 2005, 30(1): 1333-1341.
[4]
Jones W P, Rigopoulos S. Reduced chemistry for hydro-gen and methanol premixed flames via RCCE[J]. Combustion Theory and Modelling, 2007, 11(5): 755-780.
[5]
König K, Maas U. On-demand generation of reduced mechanisms based on hierarchically extended intrinsic low-dimensional manifolds in generalized coordinates[J]. Proceedings of the Combustion Institute, 2009, 32(1): 553-560.
[6]
Lam S H, Goussis D A. The CSP method for simplifying kinetics[J]. International Journal of Chemical Kinetics, 2010, 26(4): 461-486.
[7]
Wang P, Hou T Z, Wang C J, et al. Large eddy simulations of the darmstadt turbulent stratified flames with REDIM reduced kinetics[J]. Flow, Turbulence and Combustion, 2018, 101(1): 219-245.
[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]
Contino F, Jeanmart H, Lucchini T, et al. Coupling of in situ adaptive tabulation and dynamic adaptive chemistry:An effective method for solving combustion in engine simulations[J]. Proceedings of the Combustion Institute, 2011, 33(2): 3057-3064.
[10]
Kumar A, Mazumder S. Adaptation and application of the in situ Adaptive Tabulation (ISAT) procedure to reacting flow calculations with complex surface chemistry[J]. Computers & Chemical Engineering, 2011, 35(7): 1317-1327.
[11]
Lu L Y, Lantz S R, Ren Z Y, et al. Computationally efficient implementation of combustion chemistry in parallel PDF calculations[J]. Journal of Computational Physics, 2009, 228(15): 5490-5525.
[12]
Lu L Y, Ren Z Y, Lantz S R, et al. Investigation of strategies for the parallel implementation of ISAT in LES/FDF/ISAT computations[C]. Fourth Joint Meeting of the U.S. Sections of the Combustion Institute, Philadelphia, PA: Drexel University, 2005.
[13]
Dong G, Fan B C. Chemistry acceleration modeling of detonation based on the dynamical storage/deletion algori-thm[J]. Combustion Science and Technology, 2009, 181(9): 1207-1216.
[14]
Wu J, Dong G, Li Y. Parallel chemistry acceleration algorithm with ISAT table-size control in the application of gaseous detonation[J]. Shock Waves, 2019, 29(4): 523-535.
[15]
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.
[16]
Brown P N, Byrne G D, Hindmarsh A C. VODE:A variable-coefficient ODE solver[J]. SIAM Journal on Sci-entific and Statistical Computing, 1989, 10(5): 1038-1051.
[17]
Jiang H, Dong G, Chen X, et al. A parameterization of the Richtmyer-Meshkov instability on a premixed flame interface induced by the successive passages of shock waves[J]. Combustion and Flame, 2016, 169: 229-241.
[18]
Chen X, Dong G, Li B M. Numerical study of three-dimensional developments of premixed flame induced by multiple shock waves[J]. Acta Mechanica Sinica, 2018, 34(6): 1035-1047.
[19]
Burke M P, Chaos M, Ju Y G, et al. Comprehensive H2/O2 kinetic model for high-pressure combustion[J]. International Journal of Chemical Kinetics, 2012, 44(7): 444-474.
[20]
胡明媛, 董刚, 陈耀慧, 等. 斜激波与反应性气体界面作用的数值模拟[J]. 弹道学报, 2019, 31(1): 60-67.
[21]
Wu J T, Dong G, Li B M. Parallel chemistry acceleration algorithms based on ISAT method in gaseous detonation computations[J]. Computers & Fluids, 2018, 167: 265-284.