刚性源项处理方法在爆震模拟中的应用对比*

邓 利，马 虎，武晓松，周长省

(南京理工大学机械工程学院,江苏 南京 210094)

爆震波是激波和燃烧波的强烈耦合，涉及爆震燃烧的流动存在多个时间尺度，容易出现严重的刚性源项，使数值求解流动方程存在困难。

目前处理刚性源项问题的方法主要有分裂求解和耦合求解两类。采用显式耦合方法求解时，时间步长需要足够小以满足稳定性要求，否则可能导致计算不稳定甚至出现非物理解[1]。Bussing等[2-3]针对反应流中出现的多尺度问题，提出了点隐算法，将反应与流动的积分时间步长统一到伪时间步长上，消除了源项的刚性。Zhong[4]推广了Rosenbrock[5]从半隐出发求解刚性常微分方程组的方法，采用不同的方法分别求解方程的刚性项和非刚性项，增强了计算稳定性，提高了时间方向的精度，并在化学反应流的模拟中得到了成功应用。隐式处理刚性源项时，虽然很好地解决了刚性问题，但由于涉及矩阵求逆运算，随着反应组分增多，计算量急剧增大。

Leveque等[6]利用时间分裂法和MacCormack预测校正方法求解了含刚性源项的双曲型方程，认为分裂算法在空间分辨率足够的情况下，计算效果比耦合算法稍好。Oran等[7]对刚性源项方程的积分方法做了全面的论述。一步法虽然在求解刚性方程组时稳定性要求严格，但由于只需已知变量值，求解相对方便，并针对特征值为正且相当大的物理刚性问题时，使用一步法是比较有利的。Young等[8]采用逼近法与预测校正思想来处理化学反应源项，对刚性和普通常微分方程组都具有较好的计算性能，且方便与时间步长估算方法[9]配合使用，求解效率高，但容易出现质量分数不守恒的问题，需要在流场推进时不断修正。Mott等[10-11]在逼近法的基础上发展了αQSS方法，通过迭代校正组分值来保证质量守恒，比逼近法更加稳定和精确，并针对碳氢燃料的燃烧和爆震模拟取得了令人满意的效果。

刘瑜[12]对比了VODE(variable-coefficient ordinary differential equation solver)、梯形公式[13]和αQSS方法处理化学反应刚性源项的能力，认为在计算精度基本相同的情况下，梯形公式和αQSS方法效率更高。刘世杰等[14]利用梯形公式和αQSS方法对爆震波胞格进行了模拟，认为使用αQSS方法时计算的化学反应更完全，放热量更大。目前对刚性源项处理方法在化学反应流动模拟中的适用性，以及各方法之间的联系方面还缺少一定的认识。本文中在对各源项处理方法的稳定性分析基础上，分析各种方法之间的关系以及对化学反应源项的处理能力，并通过具体算例比较不同处理方法的计算效率。

1 计算方法介绍

1.1 控制方程

激波诱导燃烧为激波和燃烧波耦合的结果，可忽略黏性的影响，因此使用有限速率化学反应的欧拉方程描述为:

(1)

1.2 数值方法

具有二阶时间精度的Strang分裂格式表示为：

(2)

式中：算符Lf为流动求解算子，Lr为化学反应的求解算子，Δt为流场推进的时间步长，上标n为求解时间步数。

一步法为常微分方程初值问题的经典求解方法，求解变量n+1时刻的值时，只需n时刻的信息。解耦后的化学源项方程如下：

(3)

由于Strang分裂对算子求解精度的要求，选取2阶Runge-Kutta法作为化学反应求解算子。

在利用逼近法求解化学反应源项时，将组分i的净生成率写成如下形式：

(4)

式中：qi为组分i在反应中的生成项，τi为组分i的消耗特征时间。则刚性方程组中组分i在n+1时刻的密度为:

(5)

(6)

αQSS方法在逼近法的基础上进一步改进其迭代形式，在求解时将组分净生成率写为：

(7)

式中：pi为特征消耗时间的倒数，为了形式简便略去下标，如果qi与pi在时间段(t,t+Δt)之内为常数时，则可积分式(7)得到精确解：

(8)

整理式(8)得到：

(9)

(10)

式中：pΔt→∞时，α→1，表示无穷快的反应；而pΔt→0时，α→1/2，表示无穷慢的反应。αQSS方法中利用预测校正的方法获得一定的精度并保证质量分数守恒，而在计算中需要不断对校正步进行迭代，直到满足迭代标准，预测和校正步的详细实施过程见文献[10-11]。

使用点隐方法时，将源项写成n+1时间步的值，通过线性化消除源项刚性的问题，实现方式为：

(11)

在与流动耦合，使用具有TVD性质的二阶Runge-Kutta法推进时，其形式可写为:

(12)

式中：I为单位质量矩阵；R为空间方程离散后的余项。源项雅各比矩阵推导见文献[2,16]。

2 源项处理方法的理论分析

2.1 源项处理方法的稳定性分析

考虑模型方程：

(13)

式中：a为常数，s(u)为反应源项。

为了简便，采用一阶迎风差分格式离散式(13)的对流项，得到如下形式:

(14)

式(14)中的源项可根据处理方法写为不同的形式。

考虑一步法，将式(14)改写为:

(15)

G=1-λΔt

(16)

其中，λ的具体形式为:

(17)

为满足求解稳定性，要求误差放大因子G的模小于或等于1，解出满足稳定要求的Δt，其形式为：

(18)

从式(17)可以看出，当τ(c)非常小时，其倒数远大于其他项的值，式(17)变为λ≈1/τ(c)。因此式(18)可简化为：

Δt≤2τ(c)

(19)

可见，在使用一步法求解刚性问题严重的化学反应流动时，化学源项积分的时间步长需要小于或等于两倍最小特征反应时间。

逼近法的分析类似，假设在Δt内q和τ为常数，则可将式(14)改写为：

(20)

该形式与梯形公式相同，其稳定性分析见文献[1]，结果表明，逼近法只要时间步长满足对流项稳定条件即可。

点隐方法的分析与一步法类似，只需将式(15)中源项的变量换为n+1时刻的值，解出误差放大因子的模为:

(21)

式中：aΔt/Δx为CFL(Courant Friedrichs Lewy)数，而Δt/τ(c)为刚性系数，可知只要满足对流项稳定条件，即CFL数小于1，误差放大因子的模是恒小于1的。

αQSS方法的稳定性分析可见文献[10]，在满足对流项稳定的条件下是绝对稳定的。

稳定性分析的结果表明:隐式方法在求解数学意义上的刚性问题时具有卓越的性能；一步法为了满足求解的稳定性，需要很小的时间步长；逼近法和αQSS方法中的特征消耗时间τi满足在求解时间步长内为常数的情况下，也具有很好的求解性能。

2.2 对化学反应流动模拟的适用性分析

化学反应流动存在多个时间尺度，但可能不完全符合刚性问题的数学定义[17]，因为化学反应源项的雅各比矩阵特征值实部可能不全为负值[7]。此时，特征值为正且相当大的情况表现为物理上的刚性问题，求解时必须保证时间步长足够小以捕捉到物理特征的变化。如果源项雅各比矩阵特征值实部全为负值，则对应的物理特征表现为组分快速变化到反应平衡状态。

基于以上分析，以点隐为首的隐式方法在求解数学意义上的刚性问题时具有优秀的性能，但针对物理意义上的刚性问题时，由于点隐对化学反应源项积分没有时间步长限制，可能会由于时间步长大而导致分辨不出物理特征变化的细节；另外，流场求解时网格数较多，使用隐式方法时，在每个网格单元都要进行矩阵的求逆运算，运算量随变量个数二次方增加[1]，导致流场推进速度很慢。一步法受稳定性要求的严格限制，在刚性很严重的时候，时间步长取值造成的计算量需求将难以接受，但在求解物理和数学上都为刚性的问题时还是具有一定的优势，至少能够保证捕捉到变化细节所需要的时间步长。逼近法和αQSS方法在求解时简化了物理问题，将组分反应源项考虑为解耦形式，从而降低了求解的难度，但造成了在计算时容易出现质量不守恒的问题，使用逼近法时需要不断修正，而αQSS方法则利用组分迭代校正到一定精度来满足质量守恒。尽管如此，αQSS在积分化学反应源项时还是具有良好的计算性能，且对于化学反应的变化适应性很好，具体分析如下。

将组分i的化学反应净生成项写为式(7)形式，并将式(5)写成式(9)形式，得到：

(22)

逼近法对应αQSS方法中pΔt→0的情况，而梯形公式与逼近法本质上相同，只是缺少校正步。

对式(7)分别采用隐式和显式进行数值求解，得到如下两种形式:

(23)

(24)

分别将式(23)、(24)改写成式(9)形式，得到：

(25)

(26)

由αQSS中参数α的定义，式(25)为解耦的隐式方法，对应pΔt→∞的情况，适合快速衰减的反应，与点隐方法类似。而式(26)则为一步法，对应pΔt→-∞的情况，表明组分消耗特征时间为负值，而αQSS方法不存在这种情况。因此，由式(22)～(26)可以看出，如果组分反应可表示为解耦的形式，如式(7)，则αQSS方法几乎涵盖了其他几种方法。但是，αQSS方法中的参数α变化只能表征组分衰减到平衡状态快慢的速度，即对应特征值小于零的情况，在生成项q远大于消耗项pρ时，只依赖参数α变化并不能捕捉到快速增长的化学反应特征，此时αQSS方法通过对比生成项和消耗项的大小，自动调整时间步长以适应化学反应的变化[11]，其性质与一步法相同。从以上分析可以得到，对应物理模态快速衰减的刚性问题时，αQSS方法通过参数α的变化来适应不同的衰减速度；而对于物理模态快速增长的刚性问题时，αQSS方法通过改变积分时间步长来适应增长速度，因此对化学反应变化的适应性强。

3 数值计算

Lehr[18]进行了一系列弹道靶实验，将球形弹丸以接近爆震的速度射入一定条件下的氢气和空气预混气体中，观察到弹丸头部不稳定的燃烧现象，并拍摄了不同入射速度下弹丸头部振荡燃烧现象的阴影图，图1(a)、(b)分别对应来流马赫数为4.48和4.79的工况。McVey等[19]提出了解释激波诱导振荡燃烧的机理，其后Choi[20]对该现象做了详尽的数值模拟，认为激波诱导周期振荡燃烧现象是验证非平衡化学流计算程序的典型算例。本文中使用该算例，从计算结果、计算效率等方面比较了源项处理方法的性能。

3.1 计算结果可靠性分析

计算模型和网格选取与文献[21]相同，取球头前方2.5 mm，上方10 mm，网格点数为251×201。计算工况与Lehr[18]实验条件一致，来流温度293 K，压力42 665 Pa，预混气体来流速度分别为1 804 m/s(Ma=4.48) 和1 931 m/s (Ma=4.79)，该实验条件下对应的爆震马赫数约为4.845，CFL数取为0.85。通过记录球头轴线驻点处压力随时间变化的数据，对其进行快速傅里叶变换得到振荡频率，对比计算振荡频率和实验所得振荡频率，验证化学反应模块的正确性。为了避免由于取值区间小而造成的振荡频率计算的波动，用于计算主频的压力时间曲线选取至少100个振荡周期。

为验证源项处理方法计算结果的正确性，对来流马赫数分别为4.48和4.79的工况进行了计算，其轴线驻点压力振荡频率见表1。从表1可以看出，来流马赫数分别为4.48和4.79时，不同源项处理方法得到的振荡频率与实验值相对误差都较小。在使用251×201的网格尺度且来流马赫数4.79时的计算条件下，计算得到的组分最小反应特征时间约为6×10-11s，流动特征时间为3.5×10-9s，刚性系数Rs约为58，从稳定性分析可知，采用的化学反应时间步长应小于等于2倍特征反应时间，但文献[2]表明，在实际计算中化学反应推进步长可取为10倍最小反应特征时间。表1所列的一步法中，化学反应时间步长为6倍最小时间尺度，计算得到的振荡频率与实验值吻合，但从稳定性方面考虑，计算过程中可能存在一定的振荡。在Jachimowski化学反应模型中，组分H2O2的化学反应特征时间最小，因此出现振荡的可能性最大，图2～3所示为对应表1中Ma=4.79时，不同源项处理方法计算得到的球头轴线上组分变化曲线，其中图2为不同积分时间步长的一步法计算结果，图3(a)、(b)、(c)分别为αQSS、点隐和逼近法的计算结果。由于各组分之间的质量分数相差较大，为了对比的直观性，对图中显示的各组分质量分数进行归一化处理，归一化之后，H2的归一化质量分数为0.03，H2O的归一化质量分数采用0.25，H2O2的归一化质量分数为0.000 5，而N2的归一化质量分数为0.76。

表1 实验数据与理论模型结果对比Table 1 Comparison of experimental and theoretical data

图2(a)表示使用一步法时积分步长为6倍最小化学反应特征时间的计算结果，与图2(b)使用两倍化学反应特征时间的计算结果对比可发现，在燃烧阵面之后，H2O2变化曲线存在非物理振荡现象，与稳定性分析的预测结果一致，即积分时间步长超过稳定性要求的最大时间步长后，计算结果将出现不稳定或者发散现象，说明在使用一步法积分化学反应刚性源项时，要想获得准确的组分质量分数变化曲线，时间步长必须满足稳定性要求。积分时间步长取值满足稳定要求时，几种源项处理方法的结果相同，如图3(a)～(c)所示，因此，在计算化学反应流动时，逼近法、αQSS和点隐在稳定性方面要优于一步法。

观察图2(a)，发现振荡并没有随时间和空间变化而加剧，这是因为化学反应具有强非线性与随时间变化的特点，源项的雅各比矩阵特征值将在反应过程中不断变化，使得引入的小扰动在增长模态和衰

减模态之间反复转变，因此扰动不会进一步增强。另外，从表1可知，在使用一步法计算时，积分时间步长不满足稳定性要求时，也能获得正确的燃烧振荡频率。这是由于振荡频率主要由化学反应放热量决定的，H2O2组分质量分数很小，对于整个化学反应放热量的影响很微弱，因此其质量分数的小波动对于振荡频率几乎没有影响。此外，图3中组分变化曲线表明这几种源项处理方法对本算例中化学反应特征的捕捉效果基本相同，其原因在于化学反应是一个动态的非线性系统，某组分在这一时刻呈快速的增长模态，在下一时刻变为快速衰减的模态或者趋于平衡；并且本算例的刚性系数不大，化学反应特征在单个流动时间步内的变化也很小，因此各个源项处理方法对于该算例的捕捉效果基本相同。

3.2 计算效率分析

在计算硬件条件相同且CFL数为0.85的情况下，以来流马赫数4.79且完全发展的振荡燃烧流场作为初始计算流场，比较几种源项处理方法的计算效率。表2所示为流场推进物理时间473 μs时，不同源项处理方法所用的CPU时间，其中一步法的化学推进时间步长为6倍最小反应特征时间。

表2 不同源项处理方法消耗CPU时间Table 2 The CPU time in differentstiff source term methods

从表2中可以看到，点隐与一步法消耗的CPU时间最多，采用源项完全求逆的点隐所需时间约为αQSS的2倍，点隐方法需要计算源项的雅各比矩阵，并涉及复杂的矩阵求逆运算，求解时间随组分的增加呈指数增长，在使用基元反应模型求解爆震燃烧问题时，计算量很大；而使用一步法时，时间步长需取足够小以满足稳定性要求，极大降低了计算效率；αQSS计算时间与组分求解迭代时指定的误差ε有关，ε越小，所需计算时间越长；逼近法相对αQSS来说，缺少迭代过程，因此所用时间最少。由此可见，在计算结果精度基本相同的情况下，逼近法和αQSS的计算效率较其他两种处理方法高。

4 结 论

本文中从稳定性分析和数值计算方面对比一步法、逼近法、αQSS和点隐方法计算化学反应刚性源项问题的性能，得到结论如下：(1)一步法在积分刚性源项时，积分步长需小于或等于2倍最小时间尺度，而其他3种源项计算方法对时间步长取值没有影响；(2)一步法、逼近法实质上是αQSS的特例，αQSS方法可根据化学反应特征自动调整系数α和时间步长，适用范围广；(3)一步法在积分化学反应刚性源项时，可适当放宽稳定性要求的时间步长限制；(4)在计算结果相近以及和精度基本相同的前提下，逼近法和αQSS在求解刚性源项时具有较高的效率，并且αQSS能够适应化学反应的变化，是采用分裂法求解刚性源项时综合性能最好的方法。

[1] 刘君,周松柏,徐春光.超声速流动中燃烧现象的数值模拟方法及应用[M].长沙:国防科技大学出版社,2008:76-79.

[2] BUSSING T R A. A finite volume method for the Navier-Stokes equations with finite rate chemistry[D]. Cambridge: Massachusetts Institute of Technology, 1985.

[3] BUSSING T R A, Murman E M. Finite-volume method for the calculation of compressible chemically reacting flows[J]. AIAA Journal, 1988,26(9):1070-1078.

[4] ZHONG X. Additive semi-implicit Runge-Kutta methods for computing high-speed non-equilibrium reactive flows[J]. Journal of Computational Physics, 1996,128(1):19-31.

[5] ROSENBROCK H H. Some general implicit processes for the numerical solution of differential equations[J]. The Computer Journal, 1963,5(4):329-330.

[6] LEVEQUE R J, YEE H C. A study of numerical methods for hyperbolic conservation laws with stiff source terms[J]. Journal of Computational Physics, 1990,86(1):187-210.

[7] ORAN E S, BORIS J P. Numerical simulation of reactive flow[M]. Cambridge: Cambridge University Press, 2005:114-158.

[8] YOUNG T R, BORIS J P. A numerical technique for solving stiff ordinary differential equations associated with the chemical kinetics of reactive-flow problems[J]. The Journal of Physical Chemistry, 1977,81(25):2424-2427.

[9] CHIANG T, HOFFMANN K. Determination of computational time step for chemically reacting flows[C]∥Proceedings of AIAA 20th Fluid Dynamics, Plasma Dynamics and Laser Conference. Buffalo, New York, USA, 1989.

[10] MOTT D R, ORAN E S, VAN L B. A quasi-steady-state solver for the stiff ordinary differential equations of reaction kinetics[J]. bJournal of Computational Physics, 2000,164(2):407-428.

[11] MOTT D R, ORAN E S. CHEMEQ2: A solver for the stiff ordinary differential equations of chemical kinetics[R]. Naval Research Lab, Washington D C, 2001.

[12] 刘瑜.化学非平衡流的计算方法研究及其在激波诱导燃烧现象模拟中的应用[D].长沙:国防科技大学,2008.

LIU Yu. Investigations into numerical methods of chemical non-equilibrium flow and its application to simulation of shock-induced combustion phenomenon[D]. Changsha: National University of Defense Technology, 2008.

[13] 刘君,张涵信,高树椿.一种新型的计算化学非平衡流动的解耦方法[J].国防科技大学学报,2000,22(5):19-23.

LIU Jun, ZHANG Hanxin, GAO Shuchun. A new uncoupled method for numerical simulation of non-equilibrium flow[J]. Journal of National University of Defense Technology, 2000,22(5):19-23.

[14] 刘世杰,林志勇,孙明波,等.采用不同化学反应源项处理方法的胞格爆轰数值研究[J].国防科技大学学报,2010,32(5):1-6.

LIU Shijie, LIN Zhiyong, SUN Mingbo, et al. Numerical simulation of cellular detonation using different chemical reaction source term methods[J]. Journal of National University of Defense Technology, 2010,32(5):1-6.

[15] JACHIMOWSKI C J. An analytical study of the hydrogen-air reaction mechanism with application to scramjet combustion[R]. NASA TechnicaI Paper 2791, 1988.

[16] 潘沙.高超声速气动热数值模拟方法及大规模并行计算研究[D].长沙:国防科技大学,2010.

PAN Sha. Hypersonic aerothermal numerical simulation method and massive parallel computation research[D]. Changsha: National University of Defense Technology, 2010.

[17] TORO E F. Riemann solvers and numerical methods for fluid dynamics: A practical introduction[M]. Springer Science & Business Media, 2009:531-542.

[18] LEHR H F. Experiments on shock-induced combustion[J]. Astronautica Acta, 1972,17:589-597.

[19] MCVEY J B, TOONG T Y. Mechanism of instabilities of exothermic hypersonic blunt-body flows[J]. Combustion Science and Technology, 1971,3(2):63-76.

[20] CHOI J Y, JEUNG I S, YOON Y. Computational fluid dynamics algorithms for unsteady shock-induced combustion: Part 1: Validation[J]. AIAA Journal, 2000,38(7):1179-1187.

[21] 刘世杰,孙明波,林志勇,等.钝头体激波诱导振荡燃烧现象的数值模拟[J].力学学报,2010,42(4):597-606.

LIU Shijie J, SUN Mingbo, LIN Zhiyong, et al. Numerical research on blunt body shock-induced oscillation combustion phenomenon[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010,42(4):597-606.