A split random time-stepping method for stiff and nonstiff detonation capturing

In this paper, a new operator splitting method is proposed for capturing stiff and nonstiff detonation waves. In stiff cases, an incorrect propagation of discontinuities might be observed for general shock-capturing methods due to under-resolution in space and time. Previous random projection methods have been applied successfully for stiff detonation capturing at under-resolved conditions. Not relying on random projection of the intermediate state onto two presumed equilibrium states (completely burnt or unburnt) as with the random projection method, the present approach randomly advances or interrupts the reaction process. Each one-way reaction is decoupled from the multi-reaction kinetics by operator splitting. The local temperature is compared with a random temperature within a temperature interval to control the random reaction. Random activation or deactivation in the reaction step serves to reduce the accumulated error of discontinuity propagation. Extensive numerical experiments demonstrate the ef-fectiveness and robustness of the method. For nonstiff problems, the proposed random method recovers the accuracy of general operator splitting methods by adding a drift term.


Introduction
One of the main challenges for numerical computation of chemically reacting flows are widely varying time scales of chemical kinetics, which may be orders of magnitude faster than the fluid flow time scale [1][2][3] . Such cases exhibit numerical stiffness due to the source terms representing chemical reactions [4] . When the chemical scales are not resolved numerically in time and space, a spurious solution may occur exhibiting incorrect propagation of discontinuities and nonphysical states.
This problem is well-known and has been an active area of research during the past three decades. It was first observed by Colella et al. [5] and by analysis of a scalar problem. LeVeque and Yee [6] found that the propagation error is mainly due to numerical dissipation contained in the scheme, which smears the discontinuity front and activates the source term in a nonphysical manner. To overcome this difficulty, one may reduce numerical dissipation [3,7,8] or use a sufficiently fine mesh. Front-tracking approaches [9][10][11] or local grid/timestep refinement [12,13] may obtain the correct propagation of the reactive front. However, generally full resolution of all fine scales cannot always be afforded. Since numerical dissipation is practically inevitable, another approach focuses on establishing corrected temperatures from the artificially diffused solution [14][15][16] . Tosatto and Vigevano [17] proposed a threshold method based on a variable reconstruction within bounds determined from the local cell neighbors. Difficulties with such methods are encountered in the extension to either spatially high-dimensional or multispecies/multi-reaction kinetics based reacting flows. Wang et al. [18,19] proposed a high-order finite-difference method utilizing the Harten ENO subcell resolution method for stiff source terms. In [4] , many different methods with or without operator splitting/subcell resolution/nonlinear filters are tested, showing that the degree of propagation speed mismatch of discontinuities is highly dependent on the accuracy of the numerical method, time step and grid spacing. Kotov et al. [20] further presented a realistic hypersonic non-equilibrium flow that mimics the spurious behavior and some important numerical challenges affecting the accuracy in such simulations. Zhang et al. [7] proposed the equilibrium state method where the cell average is replaced by a local two-equilibrium-state reconstruction, making its extension to high dimensions straightforward. They also extended the method to multi-reaction systems by treating the two one-way reactions independently. Methods applicable for realistic nonequilibrium chemical kinetics with multiple finite-rate reversible reactions, to our best knowledge, have not been reported in literature so far.
Bao and Jin [1][2][3]  uniformly distributed random variable. Although the random projection method cannot avoid the introduction of numerical dissipation by shock-capturing schemes, it can eliminate its effect. The method was established for scalar problems and successfully applied to model problems of 1D/2D reactive Euler equations. With the presumption of two time-independent equilibrium states of totally burnt and unburnt gases (regardless of the detailed reaction process), the method is only suitable for under-resolved stiff cases.
In this paper, we develop a split random time-stepping method for chemically reacting flows with general nonequilibrium chemistry in a unified manner, regardless of stiff or nonstiff source terms and under-or well-resolved conditions in space and time. Unlike Bao and Jin's random projection method, the activation and deactivation of chemical reactions in the reaction step is not projected onto two prescribed equilibrium states, but onto two timedependent states corresponding to advancing the reaction by one timestep forward and interrupting the reaction, respectively. The criterion to activate a reaction follows from comparison of the local computed temperature with a randomized temperature depending on the states of the forward step and its adjoint. To randomize each reaction process, the multi-reaction system is split reaction by reaction [21,22] . By adding a drift term into the random temperature sampling, the proposed method recovers the solution of a deterministic fractional step method in nonstiff cases with increasing resolution.
The paper is organized as follows. In Section 2 , we introduce the reactive Euler equations with chemical reaction source terms. A standard fractional step method is outlined by operator splitting into the convection step and reaction step. In the reaction step, a reaction-split ODE solver is developed to approximate the exact solution for general chemical kinetics, based on which random timestepping of each reaction is performed. In Section 3 , we examine the pure ODE solver and the split random time-stepping method by extensive model examples and realistic reacting flows in both 1D and 2D. Conclusions are drawn in the last section. More information about the ODE solver are provided in the appendix.

Formulation
Assuming the flow is compressible, inviscid and in two dimensions for simplicity, the multi-species Euler equations coupled with reaction source terms take the form where U = ρ, ρu, ρv , ρe t , ρy 1 , ρy 2 , . . . , ρy N s −1 T , are vectors of the conserved variables, convective flux in the xor y -direction and source terms, respectively, with ˙ ω i representing rate of change of the i th species concentration in the reactive gas mixture due to the chemical kinetics consisting of N r reactions and N s species. Furthermore, e t = e + 1 2 (u 2 + v 2 ) is the specific total energy including the specific internal energy e . To close the system, an equation of state (EoS) of the form is used, with y i and W i denoting the mass fraction and molecular weight of the i th species, respectively, and R u being the universal gas constant.
The above conservation laws of mass, momentums and energy with source terms are usually solved by operator splitting. The first step is flow convection assuming no chemical reactions and passive transport of all species. The second step solves the system of ODEs of chemical kinetics S r : under adiabatic and constant-volume conditions with fixed total density and constant specific internal energy. The first-order accurate Lie-Trotter splitting scheme [23] or the second-order Strang splitting [24] can be employed to approximate the solution from the discrete time level n to n + 1 with a timestep t , i.e.
with symbol ' •' to separate each operator and to indicate that an operator is applied to the following arguments. For the convection operator S c , a shock-capturing scheme [25][26][27][28] can be adopted. For the reaction step S r , an ODE solver such as VODE [29] , CHEMEQ2 [30] and MTS/HMTS [31] can be used with or without adaptive error control. We first utilize operator splitting upon the nonequilibrium chemical kinetics so that a multi-reaction system can be decoupled into a series of single reaction steps. Then we introduce the established concept of random projection into the ODE solver in order to realize random ignition of reactions. Each reaction process is randomly advanced one timestep forward (activation) or interrupted (deactivation) instead of being projected onto two prescribed equilibrium states. In the following, we term the randomized and reaction-by-reaction ODE solver for nonequilibrium chemistry as Split Random Time-Stepping method (SPRANTS).

Split reaction-by-reaction ODE solver for chemical kinetics
For common nonequilibrium chemical kinetics, chemical production rates in Eq. (5) are derived from a reaction mechanism that consists of N s species and N r reactions where ν f ji and ν b ji are the stoichiometric coefficients of species i with description X i appearing as a reactant and as a product in reaction j . The total production rate of species i in Eqs. (2) and (5) is the sum of the production rate from each single elementary reaction as with k f j and k b j denoting the forward and backward reaction rates of each chemical reaction, and ρ l = y l ρ.
By operator splitting [21,22] , we can decouple the multireaction system, e.g., by Lie-Trotter splitting, as where the operator R j corresponds to a single reaction j and is independent of all other reactions. The reaction-by-reaction idea resembles a meso-scale model of microscopic kinetics where one molecule/atom can only experience one reaction at a time instance. This is also the case with stochastic simulation of chemical kinetics [32] . At macroscopic scale, reactions involving large numbers of species molecules/atoms are considered as simultaneously occurring processes. In [21] the second-order accurate Strang splitting is adopted, starting with the fastest reaction and ending with the slowest for half a timestep and then backwards for another half timestep. In our approach we simply take the traversal order not according to reaction rates but to the reaction-mechanism index sequence S r : where R 1 st is the reverse operator of R 1 st . Accordingly, for each R j , we have R j : We now rewrite the ODE in Eq. (11) in the following form [30] where q j i ≥ 0 is the production rate and p j i y i ≥ 0 is the loss rate for the i th species through reaction j .
Following the operator splitting of reactions, we continue to split each reaction j into a forward reaction and a backward reaction (for an irreversible reaction, it can be interpreted as a reversible reaction with zero backward reaction rate) such that the species involved will either gain mass or lose mass through the one-way forward/backward reaction from Eq. (12) , i.e. mass gain : q j i ≥ 0 , p j i y i = 0 or mass loss : with the simplified for the forward reaction of Eq. (11) . The backward reaction can be determined accordingly upon exchanging its reactants and products.
Since each elementary reaction is decoupled from the others and each reaction again is split into two opposite unidirectional reactions, finally only a single reaction equation of the type aA + bB + · · · −→ xX + yY + · · · (16) is considered in each operation. Mass conservation and positivity of mass fractions can be properly treated.
For the simple cases of Eq. (16) , one may find analytical solutions, see Appendix A . However, for the general form of Eq. (16) whose analytical solution is not explicitly known or difficult to derive, a more convenient alternative is to use quasi-steadystate (QSS) methods to obtain the approximate exact solution. QSS methods are based on the exact solution of Eq. (12) for constant p j i and q j i [33,34] , i.e.
As generally p j i and q j i depend on { y 1 , . . . , y N s } in Eq. (14) or (15) , Eq. (17) provides a linear approximation. For the QSS-based SPRANTS method, the stable timestep size is not limited to the characteristic time scales of the chemical species and thus a larger timestep implying less computational efforts is possible [21] .

Remark 1.
The QSS approximation adopted here in SPRANTS is first-order accurate. For application to reacting flows the achievably absolute error magnitude generally is sufficient [30] .

Treatment for mass conservation
Employing QSS in Eq. (17) for all the species participating in may not necessarily be unity so that mass may be not exactly conserved. To cure this problem, one may only advance y n k to y n +1 k of a reactant k by Eq. (17) and update all other by mass conservation of a single reaction equation in Eq. (11) . This merit of knowing the exact net gain or loss of mass of other species originates from the fact that each reaction in Eq. (11) is decoupled from others. Therefore, for the reactant k , combining Eqs. (17) and (15) we have and for the other species i = k , including other reactants and all the products in reaction j , the change of mass fraction y i = y n +1 i − y n i should obey giving the update It is easy to see that N s i =1 y i = 0 , which is equivalent to N s i =1 y i = 1 for mass conservation.

Positivity-preserving treatment
Without loss of generality, we consider the forward reaction j and assume that reactant species k has ν b jk = 0 in Eq. (19) , as ν f jk > 0 is prescribed for reactants. Similarly assuming that another reactant species i also has ν f ji > 0 and ν b ji = 0 , we combine Eqs. (19) and (21) to obtain Recalling Eq. (15) for reactants i and k , we have Upon rearranging Eq. (23) and substitution into Eq. (22) we obtain  Regarding the positivity for the choosen reactant k , according to Eq. (19) , it is inherently satisfied through positivity of the exponential function. Eq. (19) implies that 0 ≤ y n +1 k < 1 due to the negative exponent such that mass fractions of all species through reaction j are bounded within [0,1] as a result of mass conservation.

Remark 2.
The present reaction-split method using analytical or approximate solutions can perform sufficiently well, as a standalone solver, for the ODE system in chemical kinetics. Its following randomization is not motivated for integrating the ODE accurately, but primarily aimed at alleviating the effect of numerical dissipation introduced by S c through shock-capturing schemes into S r .

Finite randomization of chemical reactions
Bao and Jin [1][2][3] first proposed the idea of random projection into the ODE solver in place of the deterministic projection. They also proved that the random projection method gives first-order convergence for scalar problems. For scalar problems and Euler equations with stiff source terms, the random projection method shows excellent performance in obtaining correct shocks and reacting fronts for under-resolved spatial and temporal discretizations.
Through operation splitting of the ODE system in S r , we merely need to consider the randomization of a single one-way reaction from time t n to t n +1 . In Bao and Jin's formulation, temperature is randomized and compared with a pre-set ignition temperature, T ign . Upper and lower temperature limits are needed, i.e. T u and T b (corresponding to the two equilibrium states of the initial combustible gas mixture being completely burnt and unburnt).
Here we advance the current state vector { y 1 , . . . , y N s } through a single one-way reaction with subscript j , as in Eq. (16) , where { y 1 , . . . , y N s } + represents the advance in time by operation . The change of mass fractions for the species involved in this reaction is The reverse operation from time level n is Since mass fractions of species involved are constrained in [0,1], all mass fractions have to be rescaled if necessary according to Eq. (20) . For the two states with superscripts + and −, two limit temperatures T + and T − can be implicitly obtained according to Eq. (3) with the thermodynamic relation where ρ and e are fixed during a constant-volume adiabatic reaction and h represents the specific enthalpy. If we assume that the present reaction is exothermal, we have T − < T < T + . The converse applies to endothermal reactions. T + corresponds to T b in the original random projection method while T − corresponds to T u . Given the two limit temperatures, we can assemble a local random temperature by where θ n is a uniformly distributed random real number between 0 and 1, and T * is the randomized local temperature with  the generation of random number θ n , Bao and Jin suggested the van der Corput sampling scheme [35] . Given the random temperature T * , the unidirectional reaction j is performed as The updated state solution { y 1 , . . . , y N s } j is taken as the initial state for the next reaction j + 1 .

Remark 3.
As the random temperature T * is uniformly distributed between the two temperature limits, the mean propagation of the reaction front recovers the physically correct position [1] . With deterministic ODE solvers, accumulation of errors may lead to nonphysical reacting front propagation, see the detailed explanation in [7,36] .
Inserting Eq. (31) into the split ODE solver in Eqs. (9) and (10) , the present SPRANTS method can be written as corresponding to the Lie-Trotter splitting or as corresponding to Strang splitting.

Remark 4.
Not requiring either the flow information at each cell and its neighbors [17] or an additional procedure to locate the reacting front in the computational domain [1] , the proposed method solves the source terms at each cell locally as a 0D problem, such that its extension to 3D reacting flows is straightforward.
For nonstiff cases when the reaction zone is well-resolved in space and time, the present SPRANTS method gradually degenerates to a deterministic ODE solver upon modification of the sampling interval in Eq. (30) as T ++ is an estimated upper bound of the temperature after N time steps (e.g., N = 5 ) and T −− corresponding to its reverse state according to Eqs. (27) and (28) , and is a small positive number. Thus f represents a dynamic measure for the time resolution of the respective reaction. One can see that, when the resolution is fine and linear approximation applies to temperature evolution, f → 1 and Eq. (34) gives for a uniformly distributed θ n in Eq. (30) . The random timestepping of reactions therefore reduces to a deterministic process according to Eq. (31) in non-stiff cases.  Remark 5. Due to the reduced randomness between activation and deactivation, the proposed SPRANTS method can also cope with nonstiff problems while the original random projection method is suitable for under-resolved stiff cases [7] .

Numerical results and discussion
In this section, we consider three types of numerical experiments. The first serves to assess the split reaction-by-reaction ODE solver based on either analytical solutions or QSS approximation  for the zero-dimensional reaction operator, ignoring fluid transport. The following two types consider the coupled fluid dynamics with chemical kinetics by using simplified model kinetics and realistic finite-rate kinetics, respectively, in 1D or 2D.

Michaelis-Menten test
The first case is the Michaelis-Menten system [37] , i.e. S 1 + where the rate constants k 1 = 10 6 , k 2 = 10 −4 and k 3 = 10 −1 . The initial concentrations are 5 × 10 −7 for S 1 and 2 × 10 −7 for S 2 with S 3 = 0 and S 4 = 0 [37,38] . For this case, analytical solutions are provided for each reaction, see Appendix A . Reactions are simulated until t = 50 . In Table 1 , the L 1 and L ∞ error norms of species S 1 and S 4 are detailed, showing the expected convergence rates, i.e. 1st order for Lie-Trotter splitting and 2nd order for Strang splitting.

Hydrogen-air ignition delay test
Hydrogen ignition in air considers not only temperaturedependent reversible reactions but also third-body reactions, making the approximate solution to each reaction is practically preferred. The mechanism of H 2 -air combustion follows O'Conaire et al. [39] , consisting of 9 species (including the inert N 2 ) with 23 reversible reactions (equivalent to 46 one-way reactions). This mechanism has exhibited good prediction for the ignition delay time in [40] . All temperature-dependent reaction rates are calculated using the Arrhenius law where the subscript r is f for forward reactions or b for backward reactions and T is the temperature. Parameters A, B and T ign for the forward rate of each reaction are often given in the mechanism. Backward rates often need to be calculated from the equilibrium constant K eq and k f by assuming the corresponding reaction to be in chemical equilibrium, i.e. K eq = k f /k b [41] . The third-body effect is accounted for by the summation of third-body collision efficiencies times the corresponding molar densities of species.
Initially the reactive H 2 -air mixture is at a pressure of 1 atm, and has molar ratio 2: 1: 3.76 for H 2 : O 2 : N 2 . Nitrogen is inert. All simulations end at t = 1 × 10 −3 s. First we vary the initial temperature T 0 from 950 K to 1400 K in steps of 50 K. A fixed timestep of 1 × 10 −8 s and Lie-Trotter splitting are applied. We compare the ignition delay times predicted by the present solver with the experimental data and CHEMKIN [42] results from Ref. [40] (see its Fig. 3 ) in Fig. 1 (left). The present QSS-based reaction-split method (or abbreviated as QRS) exhibits good predictions for the ignition induction of hydrogen using the present mechanism, especially in the high initial temperature range. In Fig. 1 (right), we compare the computed mass fractions with CHEMEQ2 at an initial temperature of 10 0 0 K, and good agreement is achieved especially for the ignition time. For either QRS or CHEMEQ2, there is little difference between CPU times with different initial temperatures. QRS, however, exhibits better efficiency than CHEMEQ2 for the 9-species 23-reaction mechanism at a fixed timestep, as shown in Fig. 2 . By choosing the initial temperatures at 10 0 0 K and 1200 K, respectively, we consider the mass conservation of QRS and CHEMEQ2 in Fig. 3 . It is readily to see that QRS can always preserve the mass, whereas for the CHEMEQ2 results some total mass loss or gain occurs around the ignition time.
We continue to consider the accuracy of QRS by adjusting the timestep from 5 × 10 −8 s to 8 × 10 −7 s with an amplifying factor of 2. The initial temperature is fixed at 10 0 0 K. Figure 4 (left) shows that the temperature profiles converge with decreasing timesteps. By assessing the error norms of temperature and mass fraction of H in Fig. 4 (right), it can be seen that QRS is 1st-order convergent when Lie-Trotter splitting is applied.

Reactive Euler equations with simplified model kinetics
In this part, we consider reactive Euler equations coupled with simplified model kinetics in several stiff detonation problems. In The EoS in Eq. (3) for the model problems is simplified by p = (γ − 1) ( ρe − q 1 ρy 1 − q 2 ρy 2 − · · · − q N s ρy N s ) and T = p/ρ. Numerical experiments cover single reaction to multi-reaction system in 1D and 2D detonation problems. In our computation, the AUSM+ scheme [28] is employed together with MUSCL reconstruction using a TVD Minmod limiter [43] in the Example 1 (Chapman-Jouguet (CJ) detonation) . The first case considers the simplest reacting model, which has been studied in [7] , with only one reaction and two mutually dependent species where A represents the fuel being burnt by the one-way reaction and mass fraction of the product can be given directly by y B = 1 − y A .
Note that the ignition temperature T ign is only used by the deterministic method. The initial condition to generate the detonation wave consists of two parts in only one spatial dimension, with piecewise constants given by The left part gas is at the burnt equilibrium state and moves at a speed u CJ relative to the stationary unburnt gas of the right part. The initial CJ state on the left can be obtained in theory [1,4,7] . This problem is solved on the interval [0, 30]. The left-end boundary condition is the inflow condition with fixed identical constants as the initial data on the left. Boundary condition for the right end is extrapolation from the mirror image points inside the domain. The exact solution is simply a CJ detonation wave moving to the right and we obtain the reference 'exact' solution by the deterministic method (QRS) using a well-resolved grid ( x = 0 . 0025 ) and a timestep of t = 0 . 0 0 01 . We compare the under-resolved results by SPRANTS and QRS, respectively, using two sets of grid ( x = 0 . 25 , 0 . 025 ) and timestep ( t = 0 . 01 , 0 . 001 ) with the same kinetics. Figure 5 shows the computed pressure, density, temperature and mass fraction. The proposed random method can capture the correct propagation of the detonation wave with both coarse and fine grids, while the deterministic method produces the spurious solutions in the same under-resolved conditions, i.e. a weak detonation wave propagates faster than the theoretical detonation speed of D CJ = 7 . 124 in this case. Since a coarser grid with a larger timestep renders the stiffness more severe, the deterministic method produces far more nonphysical weak detonation wave compared to the SPRANTS or the reference solution. The location of the reacting front on the coarse grid may be shifted from the exact location due to randomization, but the shift amplitude does not grow in time [1] , whereas the error accumulates with the deterministic method.
In Fig. 6 , the SPRANTS result based on a very coarse grid ( x = 0 . 6 corresponding to 50 grid points) is compared with the aforementioned under-resolved solutions by SPRANTS in terms of pressure and temperature at t = 1 . 5 . Correct location of the detonation wave is captured despite the smeared discontinuity. Convergence of pressure and temperature profiles with an increasing resolution can be seen towards the reference solution, demonstrating the accuracy of the proposed SPRANTS method in capturing the correct propagation speed of discontinuities at under-resolved conditions. We also notice that a grid of 300 nodes is employed to obtain the correct wave propagation in [7] and 50 grid points are used by a high-order finite difference scheme (WENO5/SR) in [18] for this case.
In Fig. 7 , we demonstrate that with the drift term the SPRANTS solution captures not only the correct location of reacting front but also the resolved reaction zone, in good agreement with the reference solution obtained by the deterministic method. Without the drift term, although the random method can still give the correct shock propagation, it fails to capture details of the resolved postshock reaction zone by overshooting the temperature magnitude.
Example 2 (Strong detonation) . This example considers a reacting model, which has been studied in [7] , with one reaction and three  4 , 300 , 0 , 0 , 2 , 32 , 18 ) , The initial condition of piecewise constants is given by 10 , 8 , 0 , 0 , 1 The left part gas is at the burnt equilibrium state and it is moving at a speed larger than u CJ relative to the stationary unburnt gas of the right part so that a strong detonation wave is to occur. This problem is solved on the interval [0, 50]. The exact solution consists of a detonation wave, followed by a contact discontinuity and a shock, all moving to the right. Again, we obtain the reference solution by QRS using a resolved grid ( x = 0 . 0025 ) and a very small timestep ( t = 0 . 0001 ). We compare the results by SPRANTS and the deterministic method using a coarse grid and a finer grid with stable timesteps, as explained in Fig. 8 . Note that in the deterministic method, we adopt both the Arrhenius model and the Heaviside model for the chemical kinetics. The proposed SPRANTS method can capture all discontinuities effectively, while the deterministic method produces spurious solutions at the same under-resolved conditions. In particular, using the Heaviside model, the deterministic method produces a less accurate solution due to the stronger stiffness compared to the Arrhenius model (see the right column of Fig. 8 ).

Example 3 (Strong detonation)
. This case considers a multi-step reaction mechanism with two one-way reactions and five species 1)   4 , 0 , 0 , −20 , −100 , 0 ) , 32 , 17 , 18 , 28 ) , The initial condition of piecewise constants is given by 20 , 10 , 0 , 0 , 0 . 17 , 0 . 63 , 0 . 2 ) , x < 2 . 5 , ( 1 , 1 , 0 , 0 . 08 , 0 . 72 , 0 , 0 , 0 . 2 ) , This problem is solved on the interval [0, 50].  Figure 9 presents different computational conditions and results obtained accordingly. All waves are captured with the correct speeds by the SPRANTS method, in good agreement with the reference solution. However, the deterministic method obviously cannot handle the Heaviside model with the same under-resolved grids and timesteps. Although the slower propagation of the reacting front is captured by the Arrhenius model, the deterministic method still results in spurious weak detonation, refer to the transit points around x ≈ 40 especially in the profiles at the right column of Fig. 9 .

Example 4 (CJ detonation in 2D
) . This 2D case extends EXAM-PLE 1 to a radially symmetric point-source explosion, where A in Eq. (38) is amplified by a factor of 10 , 0 0 0 to approximate the infinitely fast reaction with extreme stiffness. Similar tests have been studied in [3,16] .
A quarter domain is considered exploiting sectorial symmetry on [0, 50] × [0, 50]. The hot-spot area of the initial hightemperature high-pressure burnt gas is a circle with radius 10 and the reactive unburnt gas takes the outside. Initial condition is the same as in Example 1 except the initial velocity of the circle area is adjusted to along the radial direction, i.e.
where r = x 2 + y 2 . In our computations, a coarse grid (200 × 200) and a finer grid (20 0 0 × 20 0 0) are employed referring to Example 1. Corresponding timesteps are t = 1 × 10 −2 and 1 × 10 −3 , respectively. With the finer grid, the deterministic method still produces a spurious solution at t = 1 . 5 , see the left column of Fig. 10 , in that a nonphysical weak detonation wave is generated and the reacting front is no more circular. In contrast, the SPRANTS method can capture shape and location of the CJ detonation front accurately, see the right column of the figure, by observing the radial velocity vector in the pressure contour even in the low resolution and the self-similarly circular outwards-developing detonation fronts in black/white lines of two resolutions at different times. The line-marked locations calculated by the random method in two resolutions agree excellently with each other and thus a grid convergence to the exact solution is reasonable to expect for the proposed SPRANTS method. With negligible curvature effects [44,45] and the under-resolved reaction zone being infinitesimal, the calculated speed of the detonation front approaches the 1D theoretical speed of D CJ = 7 . 1247 as in Example 1. 2D) . The present case considers the same multi-step reaction mechanism as in Example 3 except that q OH in Eq. (38) changes into −50 . This is a 2D case used to prove the dimension-independent nature of the proposed method, unlike the original random projection method which requires a dimension-by-dimension scanning for local projection. Geometry and initial condition of piecewise constants in the 2D domain can be referred to [7] .

Example 5 (Strong detonation in
A uniformly distributed coarse grid (300 × 100) and a refined grid (30 0 0 × 10 0 0) are employed. Corresponding timesteps are t = 5 × 10 −4 and 5 × 10 −5 , respectively. The reference solution is obtained by the deterministic method using the fine grid and tiny timestep. The comparison of the SPRANTS method and deterministic method on capturing stiff detonation waves is based on the under-resolved grid and timestep. In Fig. 11 , at t = 0 . 1 the spurious solution given by the deterministic method on the coarse grid contains a too fast weak detonation wave, which has passed half of the domain. However, the correct detonation waves from the SPRANTS method on the same resolution and the deterministic method on a fine grid agree with each other excellently. Good agreement of the self-similar propagation of the detonation wave from t = 0 . 1 to 0.3 also can be seen in the mass fraction contour given by the reference solution and the under-resolved SPRANTS solution, respectively. The slight difference between the two solutions lies in some small near-shock statistical fluctuations due to the random nature of the method [1] .

Reactive Euler equations with realistic nonequilibrium kinetics
In this subsection, we validate the SPRANTS method for capturing stiff detonation waves governed by the reactive Euler equations coupled with realistic chemical nonequilibrium kinetics which introduces multiple temperature-dependent reactions with distinct timescales. To our knowledge, both the two test cases below are reported for the first time, taking into account the detailed hydrogen-air combustion mechanism as in Section 3.1.2 . Two different scenarios with a CJ detonation and strong detonation wave, respectively, are simulated in 1D or 2D domain. The convection operator adopts an ordinary shock capturing scheme as in the former subsection, and the reaction step is solved by the proposed SPRANTS method and CHEMEQ2 as the deterministic method to make a comparison. Reaction splitting in the SPRANTS method employs the 2nd-order Strang scheme to reduce splitting errors.
Example 6 (Realistic CJ detonation) . The setup consists of two parts divided by a shock moving to the right in a 1D domain of length L = 4 m, as in Table 2 . The theoretical CJ detonation states for the unburnt gas can be generated using the NASA Chemical Equilibrium Analysis (CEA) program [46] , and according to the CJ condition [1,4,7] , i.e.
we adopt u b = 800 m/s ≈ u CJ for the initial velocity of the burnt gas, to generate a CJ detonation wave sweeping the stationary unburnt gas. The shock is initially located at x = 0 . 4 m. Boundary condition for the left/right end is simply extrapolation from the mirror image points inside the domain. All simulations are terminated at t = 1 . 2 × 10 −3 s and use the same mechanism [39] .
The exact solution is a steady self-similar CJ detonation wave traveling from left to right. We obtain the reference exact solution by the deterministic method using a very fine grid with 10,0 0 0 points and a fixed tiny timestep of t = 5 × 10 −8 s. Two sets of under-resolved grid and timestep are considered, i.e.
In Fig. 12 at the given time: although the resolution of the grid and timestep is far lower than the resolved solution, the SPRANTS method predicts the properties of the flowfield in quite good agreement with the reference solution, including the location of the detonation wave and profiles of the mixture pressure and density. The obtained profiles tend to converge to the reference solution with increasing resolution (and decreasing stiffness). In contrast, using the same under-resolved grid and timestep, the deterministic method yields the spurious nonphysical weak detonation ahead of the shock and the flowfield profiles are totally changed in an incorrect way. In Fig. 13 , wave propagation at different times is presented by looking into the pressure distribution. Despite the deviation by few grid points, SPRANTS can always capture the correct wave location while the error in the location of reaction front by the deterministic method deteriorates by showing a too fast weak detonation wave. Note that the von Neumann spike inside the reaction zone of the reference solution can be calculated only by very fine resolution both in space and time.
In Fig. 14 (left), we additionally obtain several solutions by the deterministic method with CHEMEQ2 using N c = 40 0 , 80 0 up to 6400 grids with linearly decreasing global timesteps (lower to t = 5 × 10 −8 s as the reference solution with N c = 10 , 0 0 0 ). It can be seen that the pressure profiles converge to the reference solution (with 10,0 0 0 grid points) including the spurious weak detonation waves with N c = 400 to 1600. When the number of grid points increases to 3200 or more, the weak detonation wave disappears and the correct location of the reacting front is captured. We compare the CPU times for the reaction step of two methods based on different grids, listed in Table 3 and plotted in Fig. 14 (right), as the computational cost of integrating the ODE system dominates in reacting flow simulations. With the same mechanism, SPRANTS consumes more CPU time in the reaction step than the deterministic method using the same resolution, since each random reaction needs to assume a forward state or backward state to determine the random temperature, invoking a costly iterative root-finding operation. The deterministic method requires a much higher resolution in both space and time to reach the same prediction accuracy so that its overall computational efficiency dramatically decreases.
and the unburnt gas occupies the remaining domain in front of the initial shock. Initial states are identical with those in Example 6 except for the x -velocity of the post-shock part being increased to u b = 20 0 0 m/s > u CJ , to create a strong detonation wave. The boundary condition for the left/right end is simply extrapolation from the mirror image points inside the domain and the top/bottom boundary is considered as a slip wall. All simulations are finished at t = 1 × 10 −3 s and still use the 9-species 23reaction mechanism [39] .
With the previous 1D example, it was shown that the deterministic solution based on a grid of 3200 points in the 4 m long domain recovers the correct shock position in Fig. 14 (left). Therefore, we generate a reference solution in 2D by the deterministic method using a fine grid with 30 0 0 × 10 0 0 points and a fixed tiny timestep of t = 2 . 5 × 10 −8 s. A set of under-resolved uniform grid and timestep is also considered, i.e. 150 × 50 , t = 2 . 5 × 10 −7 s. Figure 15 displays the density distributions along with locations of the detonation wave at different times in three solutions. In comparison with the reference solution, the SPRANTS method gives reasonable locations of the reacting front at all times. Due to the low resolution used in SPRANTS, detailed characteristics presented in the reference solution such as the triple points, slip lines, small vortices and peak values of density are diffused while the overall flowfield including the profile of reacting front has been correctly captured. For the deterministic method with the same resolution, a spurious weak detonation wave can be observed with a maximum error of nearly 10% of the domain length within only 1 ms.

Conclusions
A new operator splitting method for simulating chemically reacting flows, especially for capturing stiff detonation waves in under-resolved conditions has been developed. Two procedures based on operator splitting are included: for the convection step, any shock-capturing scheme can be used; for the reaction step, the multi-species multi-reaction ODE system in the source terms is further split in a reaction-by-reaction manner. Each reaction either proceeds a timestep forward or is interrupted according to a local random temperature rather than a deterministic process with growing error accumulation. A wide range of numerical experiments including not only simple model kinetics but also realistic nonequilibrium chemistry such as the temperature-dependent finite-rate hydrogen-air combustion are considered in 1D and 2D flows, demonstrating the following properties: 1. Mass conservation and positivity of species concentration can be guaranteed by the reaction-split ODE solver, which is almost unconditionally stable due to its using either analytical or approximate exact solutions. 2. The proposed SPRANTS method can effectively predict the correct propagation of discontinuities as well as the overall flowfield information in under-resolved conditions, for both model kinetics and realistic finite-rate nonequilibrium kinetics. 3. Compared with the deterministic method using CHEMEQ2, the present SPRANTS method exhibits better computational efficiency as it can correctly capture the detonation wave with a larger timestep on coarse grids for nonequilibrium reactive flows. 4. By adding a drift term into the random temperature sampling, SPRANTS can recover the deterministic solution as the resolution improves with decreasing stiffness. 5. The dimension-independent algorithm for the source terms makes further 3D extension of the proposed method straightforward.
Employing high-order low-dissipation schemes for the present method is a subject of future research.   .
After the determination of the new state of the reactant species [ A ], states of the remaining species including all the products and other reactants can be updated by mass conservation in Eq. (20) .