BYCFoam: An Improved Solver for Rotating Detonation Engines Based on OpenFOAM

: A rotating detonation engine (RDE) is a highly promising detonation-based propulsion system and has been widely researched in recent decades. In this study, BYCFoam, the latest gaseous version of the BYRFoam family, is developed specifically for RDE simulations. It is based on the standard compressible flow solver rhoCentralFoam in OpenFOAM and incorporates several enhancements: improved reconstruction variables and flux schemes; detailed chemistry and transport properties; the utilization of an adaptive mesh refinement (AMR) and dynamic load balancing (DLB). A series of comprehensive numerical tests are conducted, including the shock-tube problem, shock-wave diffraction, homogeneous ignition delay, premixed flame, planar detonation, detonation cellular structure and rotating detonation combustor (RDC). The results demonstrate that BYCFoam can accurately and efficiently simulate the deflagration and detonation processes. This solver enhances the capability of the BYRFoam family for the in-depth exploration of RDE in future research.


Introduction
In recent decades, detonation-based propulsion systems have received increasing interest, due to their ability to achieve higher thermal efficiency and faster energy-releasing rates compared to deflagration [1,2].The rotating detonation engine (RDE) is considered a highly promising type of detonation-based engine and has been extensively researched in experiments and numerical simulations [3,4].Due to limitations in measurement techniques, detailed flow field information cannot be obtained through experiments [5][6][7].Consequently, numerical simulations are popularly employed to gain a deeper understanding of the complicated phenomena in detonation processes.Chen et al. [8] and Liu et al. [9] investigated the effects of reversed shock waves on the operation mode in a hydrogenfueled RDE with a fifth-order WENO code.In the work by Mikhalchenko et al. [10][11][12], three-dimensional (3D) numerical modeling were conducted with different reactant mixtures and RDE configurations to study the effects of the detonation onset in the transition period.Schau et al. [13] presented a detailed analysis of the wave characteristics in the simulation of a methane-oxygen RDE by using their RAPTOR code.Li et al. [14] performed a 3D numerical simulation in ANSYS Fluent to investigate the film cooling efficiency in a hydrogen-enriched kerosene-fueled RDE.
OpenFOAM [15] is an open-source computational fluid dynamics (CFD) toolkit and has been used by many groups to simulate nonreactive and reactive flows [16,17].In the framework of OpenFOAM, a density-based solver, namely, rhoCentralFoam [18], has been developed to solve compressible high-speed flows by the central-upwind scheme [19].Based on rhoCentralFoam, Ettner et al. [20] developed ddtFoam with the HLLC scheme [21] to simulate the deflagration-to-detonation transition in inhomogeneous mixtures; Marcantoni et al. [22,23] developed rhoCentralRfFoam in predicting propagation of onedimensional and two-dimensional detonation waves; Huang et al. [24] developed a hybrid Eulerian-Lagrangian solver RYrhoCentralFoam to simulate detonation in two-phase gasliquid mixtures.
The BYRFoam family is a multiphase, multiregion, compressible and reacting flow toolbox designed for RDE modeling.It is also based on the rhoCentralFoam solver and has been continuously improved at Peking University.Xia et al. [25] first used BYRFoam to study the influence of the predetonator in RDE.After that, Shen et al. [26] considered the effects of supersonic nozzle guide vanes in a two-dimensional rotating detonation combustor (RDC) with BYRFoam.Thereafter, a two-phase Eulerian-Lagrangian version of BYRFoam was also developed.Rong et al. [27] investigated a carbon-hydrogen/oxygenrich air RDE with this two-phase gas-solid version.Recently, Hou et al. [28] conducted an unsteady conjugate heat transfer simulation of RDC with a multiregion version of BYRFoam.Generally, the BYRFoam family has achieved satisfactory results in these two-dimensional and three-dimensional RDE simulations.
In this work, an updated gaseous version of the BYRFoam family, namely BYCFoam, is developed by incorporating several new features.Firstly, we modify the original reconstruction variables in rhoCentralFoam, along with an integration of four flux schemes.Secondly, the open-source library Cantera [29] is introduced for the calculation of detailed chemistry and transport properties.Thirdly, the adaptive mesh refinement (AMR) with dynamic load balancing (DLB) [30,31] and the chemistry DLB library [32] are integrated with the solver to improve the computational efficiency.The performance of BYCFoam is validated against comprehensive numerical tests, including a zero-dimensional (0D) homogeneous autoignition problem, one-dimensional (1D) and two-dimensional (2D) shock-wave problems, 1D premixed laminar flame, 1D and 2D detonation-wave problems and a 2D rotating detonation combustor.
This paper is structured as follows.The governing equations and numerical methods are described in Sections 2 and 3, respectively.The numerical validations of BYCFoam are elaborated in Section 4. Finally, conclusions are drawn in Section 5.

Governing Equations
The BYCFoam solver utilized in this study is based on the standard rhoCentralFoam solver in OpenFOAM [15] (version v1812) and is coupled with the open-source library Cantera [29] (version 2.5.1) for the calculation of detailed chemistry and transport properties.The governing equations for unsteady, compressible, reacting flows involving ideal gases are summarized as follows: where t is time, ρ is the density, u is the velocity, p is the pressure, T is the temperature, λ is the mixture thermal conductivity, R 0 is the universal gas constant, and W is the mean molecular weight of the mixture.h k , Y k , ωk and j k are the specific enthalpy, the mass fraction, the production rate by chemical reaction and the mass diffusive flux of the kth species, respectively.In addition, the viscous stress tensor τ and the total internal energy E are defined as: where µ is the dynamic viscosity, I is the identity tensor, and e is the internal energy.
The mixture-averaged model in Cantera is employed for the computation of transport properties, providing a beneficial balance between computational accuracy and cost.The species diffusive flux j k in Equation ( 4) becomes: where X k is the mole fraction of species k, V k is the ordinary diffusion velocity, and the correction velocity V C is introduced to ensure total mass conservation.The dynamic viscosity µ, the thermal conductivity λ and the diffusion coefficient D k are calculated by the mixture-averaged formulation [33] as follows: , and where µ k , W k and λ k are the dynamic viscosity, the molecular weight and the thermal conductivity of the pure species k, respectively.D jk are the binary diffusion coefficients for the "j-k" species pair.The default transport models in OpenFOAM are also listed as follows.The dynamic viscosity µ is calculated with Sutherland's formula [33]: where A s and T s are two fitting parameters.The thermal conductivity λ is calculated using the modified Eucken correlation [34]: where c v is the specific heat capacity at constant volume, and R = R 0 /W is the gas constant.The species diffusive flux j k is calculated with the unity Schmidt number assumption:

Numerical Schemes
In the current OpenFOAM platform, a variety of spatial and temporal discretization schemes are available and flexible to use.For example, the temporal accuracy could be either first-order (e.g., the Euler schemes) or second-order (e.g., the backward and the Crank-Nicolson schemes).In the scope of this context, we chose a first-order Euler time discretization.The second-order central scheme was applied for the diffusion terms.As for the convection terms, we adopted a second-order MUSCL reconstruction [35] with a minmod limiter and implemented four schemes (Kurganov [19], HLL [36], HLLC [21] and AUSM+M [37]) to calculate the convective flux.It should be noted that the Kurganov scheme is used in the standard rhoCentralFoam solver, so we added it into our solver for comparisons.Moreover, we chose three independent primitive variables as reconstruction variables, as listed in Table 1.The improvement can be seen in Section 4.1.The time steps ∆t chem necessary for integrating the chemical source terms are typically much smaller than the time steps ∆t F needed for solving the nonreactive flows.Therefore, an operator splitting method [38] was used here.It is only first-order accurate in time, but in many practical cases, results are nearly identical compared with a second-order Strang splitting [38,39].With this operator splitting method, the source term was integrated separately in Cantera and then contributed to the governing equations in OpenFOAM.Specifically, at each fluid time step ∆t F , the Sundials CVODE library [40] in Cantera was used to integrate the chemical ordinary differential equation (ODE) systems with time steps ∆t chem in each computational cell.Then, the governing equations were solved with time step ∆t F in the whole computational domain.
Chemkin-formatted mechanisms [41] are most widely used in combustion research.They can be easily read and converted to the Cantera formatted files (cti/xml/yaml), which support different reaction types, such as Arrhenius, third-body and falloff reactions.Moreover, the pressure-dependent PLOG reaction type and analytically reduced chemistry (ARC) [42] are also included in Cantera, which are not supported in OpenFOAM.
Figure 1 provides a flow chart at each fluid time step ∆t F in BYCFoam.At the beginning of a time step, the mesh update (AMR) is carried out if current the simulation time matches with the user-defined refine interval.Multiple user-defined refinement criteria can be set to conduct mesh refinement and coarsening.The error indicator of density gradient ϵ ∇ρ defined by Zheng et al. [35] was mainly used in this work.After mesh refinement and coarsening, a mapping of field values was implemented in a straightforward and conservative manner [30].
In OpenFOAM, the MPI parallel algorithm is achieved through the process of domain decomposition, wherein the simulation domain is divided into subdomains and allocated to individual processors.When the mesh is updated, the mesh load imbalance ib max is calculated as: where C i is the number of cells on the ith processor, and C = ∑ i C i /N is the average number of cells on N processors.The mesh balancing (DLB_Mesh) is conducted when ib max exceeds the user-defined allowable imbalance, typically set at 20%.Next, convective fluxes are calculated by the user-selected scheme as listed in Table 1 and transport properties are updated from Cantera.Furthermore, the chemistry balancing (DLB_Chem) is carried out by using the DLBFoam library [32] and the chemistry ODE systems are solved by Cantera in each computational cell.At the final stage of a time step, the governing equations are solved with explicit convective terms, explicit source terms and implicit diffusion terms.
It is noteworthy that the DLB_Mesh and DLB_Chem play distinct roles in BYCFoam.The DLB_Mesh enables a balanced distribution of cells on different processors after the AMR, while the DLB_Chem helps to balance the calculation time on different processors for solving the stiff chemistry ODEs.More details on the AMR, DLB_Mesh and DLB_Chem are described in [30][31][32].

Thermodynamic Properties
The thermodynamic properties of individual species are calculated with JANAF polynomials (also known as NASA-7 polynomials) [43].The specific heat capacity at constant pressure c p,k , the specific enthalpy h k and the specific entropy s k of species k are: where R k = R 0 /W k is the gas constant of species k and a k,0 -a k,6 are two sets of polynomial coefficients for each species, corresponding to a low and a high temperature range: The thermodynamic properties of the mixture are determined by employing a weighted average of mass fractions, taking the calculation of specific heat capacity c p as an example: where A i are defined in OpenFOAM as follows: In the framework of OpenFOAM, all species should have the same intermediate temperature T Mid,k , typically set at 1000 K.Under this condition, Equation (23b) can be expressed as Equation (23c).
However, the expression of Equation ( 23c) is not applicable for some complicated negative-temperature coefficient (NTC) mechanisms, such as DME [44] and nheptane [45], in which the intermediate temperature T Mid,k of each species are not all the same.Figure 2 shows specific heat capacity c p of stoichiometric DME/air mixtures, calculated with Equation (23a,c).We can see that Equation (23c) produces discrepancies within 710-1000 K, which correspond to the intermediate temperatures of DME and air, respectively.The autoignition of homogeneous stoichiometric DME/air mixtures under constant volume conditions is calculated with different initial temperatures T 0 and initial pressures P 0 .The standard chemFoam solver in OpenFOAM is used for comparison with the results calculated by Chemkin.The NTC behavior of the ignition delay time and two-stage ignition are presented in Figure 3a,b, respectively.We can see that chemFoam underestimates the ignition delay time, with a maximum relative error around 20% when compared to Chemkin.This is partially induced by the inappropriate use of Equation (23c) in OpenFOAM, and the error decreases as the initial pressure P 0 increases.The seulex ODE solver used in chem-Foam is also responsible for such deviations in ignition delay time, as reported in Yang et al. [46].Then, in BYCFoam, we fixed it by using the original expression Equation (23a) and the CVODE solver in Cantera.The improvement can be seen in Section 4.3.

Numerical Validations 4.1. One-Dimensional Sod Shock Tube
In the first test case, we consider the standard Sod shock-tube problem [47].The inviscid Euler equations are solved here and also in Section 4.2.A uniform grid of 200 cells was set in the one-dimensional (1D) computational domain [0, 1].The initial conditions were given by (ρ, u, p) L = (1, 0, 1) and (ρ, u, p) R = (0.125, 0, 0.1), where the discontinuity was set at x = 0.5.This discontinuity would result in the formation of a right-propagating shock wave and contact discontinuity, along with a rarefaction wave propagating to the left.The slip and adiabatic wall boundary conditions were set on the left and right sides of the computational domain.
Figure 4 shows the solutions calculated by rhoCentralFoam and BYCFoam at time t = 0.15, and the Riemann exact solution is also included for comparisons.The numerical results are almost identical, and they all show a correct agreement with the exact solution.Figure 5a is the close-up of pressure distributions behind the shock waves.Compared with BYCFoam, rhoCentralFoam yields a greater pressure fluctuation, attributed to the distinct selection of reconstruction variables as mentioned in Section 3.1.Then, we considered the numerical dissipation around the contact discontinuity.As shown in Figure 5b, AUSM+M and HLLC schemes have lower dissipation than HLL and Kurganov schemes.In general, BYCFoam demonstrates improvements over rhoCentralFoam.

Two-Dimensional Shock-Wave Diffraction
We now consider a two-dimensional (2D) problem investigated by Quirk [48] to test the shock robustness of flux schemes.A planar shock wave moving around a corner was simulated on a 400 × 400 grid.The left boundary of the computational domain was specified with an inflow condition of Mach number 5.09, while the remaining boundaries were defined by slip and adiabatic wall boundary conditions.Figure 6 shows the density contours with rich flow features, such as shock-wave diffraction, reflection and interaction.As reported in Kim et al. [49], the HLLC scheme produces shock instability at the planar moving shock.In contrast to the unphysical results of HLLC, the other three schemes work well and eliminate spurious oscillations.Considering the above simulation results in Sections 4.1 and 4.2, the AUSM+M scheme is a better choice for further validations in subsequent subsections.

Zero-Dimensional Homogeneous Autoignition
This case focuses on the homogeneous autoignition processes under constant volume conditions, to assess the performance of BYCFoam in terms of chemical reaction calculations.Two stoichiometric fuel/air mixtures with NTC behavior were considered, as shown in Table 2. Figure 7 shows the evolutions of the temperature during autoignition obtained from three different solvers.We can see that chemFoam underpredicts the ignition delay time of DME and n-heptane.Here, chemFoam used the CVODE solver in Cantera to keep consistency with BYCFoam.The absolute tolerance and relative tolerance of chemistry ODE solver were set to 10 −15 and 10 −9 , respectively.Thus, the discrepancies in Figure 7 are not attributed to the chemistry ODE solver but rather stem from the inappropriate calculation of the thermodynamic properties, as detailed in Section 3.2.On the other hand, the results of BYCFoam show excellent agreement with the calculations by Chemkin.This demonstrates that the correction in BYCFoam works well and the coupling with Cantera is successful in solving chemical kinetics.

One-Dimensional Premixed Laminar Flames
This classical one-dimensional test is designed to evaluate the performance of our solver in situations where the convection, diffusion and chemical terms dominate and balance, according to Law (2010) [50].The temperature and pressure of the premixed hydrogen/air mixture were 300 K and 1 atm, respectively.The chemical mechanism was the one proposed by Li et al. [51] accounting for 9 species and 21 elementary reactions.
Figure 8 shows the laminar flame speeds obtained from three different solvers.The standard reactingFoam solver in OpenFOAM obviously has different results compared to BYCFoam and Cantera.The simplified transport models in Equations ( 14)-( 17) are responsible for such deviations.As plotted in Figures 8 and 9, the results show a great agreement between BYCFoam and Cantera.The slight variance in flame speeds can be ascribed to the difference in methodologies, as Cantera utilizes the modified damped Newton's method [52] to directly seek the converged solution of 1D steady equations, while BYCFoam achieves the quasi-steady-state flame through transient simulation.Generally, the use of the mixture-averaged model enables BYCFoam to accurately simulate one-dimensional premixed laminar flames.

One-Dimensional Planar Detonation Wave
In the following subsections, three detonation cases involving a strong coupling of shock waves and chemistry are considered.In this part, we simulated a one-dimensional planar detonation wave with kerosene/air mixtures under a grid size of 0.1 mm.A similar validation of hydrogen/air detonation was conducted by Sheng [53] with BYCFoam, which gave a satisfactory result with a maximum relative error around 0.2%.
As shown in Figure 10a, a detonation tube was filled with premixed gaseous kerosene/air mixtures at 1 atm and 373 K.The detonation wave was directly initiated by a hot spot (50 atm, 2000 K) set at the left end of the tube.The slip and adiabatic wall boundary conditions were set on the left and right sides of the tube.Due to the complex nature of the kerosene mixture, the full-scale kinetic mechanisms were too time-consuming for computations.Consequently, a reduced C 12 H 23 /air mechanism from Liu [54] was employed to model the detonation of kerosene/air.This mechanism contained 18 species and 30 reactions, as presented in Appendix A. Figure 10b illustrates the temperature and pressure history in the detonation tube with stoichiometric kerosene/air mixtures.It can be seen that a stable propagating detonation wave is captured well.Furthermore, the theoretical Chapman-Jouguet (CJ) properties calculated by SDToolbox [55] are compared with the 1D simulation results in Table 3.The maximum relative error is around 1%, which confirms that BYCFoam is capable of simulating the premixed near-stoichiometric kerosene/air detonation propagation with this simplified mechanism.Table 3. Comparisons of one-dimensional detonation properties with varying equivalence ratios ϕ.U is the detonation velocity, p and T are the pressure and temperature at the CJ points of the detonation wave, respectively.err is the maximum relative error between 1D simulations and the theoretical CJ properties of (U, p, T).It should also be noticed that the reduced mechanism used here was only verified over a range of equivalence ratios 0.8 ≤ ϕ ≤ 1.4.When we tried to further expand the range of ϕ, the error increased rapidly or the detonation initiation failed.Therefore, some more detailed mechanisms will be considered in future work.

Two-Dimensional Detonation Cellular Structure
In contrast to the one-dimensional planar results from Section 4.5, detonation waves actually exhibit complex multidimensional structures.Here, the two-dimensional detonation cellular structure was simulated to further verify the reliability of our solver.A 6 cm height channel with a H 2 /O 2 /Ar mixture (2/1/7 by volume) was set with the same initial conditions and chemical mechanism as Oran et al. [56].As shown in Figure 11, the one-dimensional CJ detonation was placed on a two-dimensional domain, and the transverse perturbation was created by placing an unreacted gas pocket (10 mm × 14 mm) behind the leading shock.The temperature and pressure of the gas pocket were set to seven times the initial conditions.The slip and adiabatic wall boundary conditions were applied to all sides of the computational domain.The base grid size was 200 µm, and the finest grid size was 25 µm after a three-level mesh refinement.A combination of error indicator ϵ ∇ρ , mass fraction Y H 2 and temperature T were used as adaptation criteria, listed in Table 4.The default threshold value of the error indicator was set to 0.04.The computation could fastly converge to a well-developed detonation structure but lost some details in transverse waves and slip lines.Then, we continued the calculations by adjusting the threshold to 0.01, which enabled a better resolution of transverse waves and slip lines.Figure 12 shows the final flow fields of a mode-four detonation, which is consistent with that of Oran et al. [56].Here, the mode means the number of triple points or transverse waves.Figure 13a gives a detailed look around the detonation front, and the triple wave structure is captured well in this simulation.The triple-point trajectories were also recorded by the history of the maximum pressure, as shown in Figure 13b.The detonation cell size was approximately 53 mm × 30 mm, which was close to the values reported in [38,56,57], as listed in Table 5.In general, the current simulation confirmed the capability of BYCFoam to resolve complex detonation wave systems.BYCFoam and its early version BYRFoam have been widely used in the RDC simulations of our group [25][26][27][28]53].In this part, we performed a two-dimensional test to assess the computational efficiency improved by the AMR and DLB algorithms.Figure 14 illustrates the schematic of an RDC in a 2D computational domain (12 cm × 4 cm).The premixed stoichiometric hydrogen/air mixtures were injected into the combustor through the bottom boundary with Laval nozzle inlet conditions [8,9].The area ratio of nozzle was 1.7.The stagnation pressure and temperature of the fresh gas were fixed at 8 atm and 540 K, respectively.Under this condition, the inlet mass flow rate in the simulation was about 486 kg/(m 2 • s).The wave transmissive boundary condition [58] was adopted at the outlet to prevent the shock reflection from the top boundary.Periodic boundary conditions were set on the left and right sides of the computational domain.The time step in the simulation was about 10 −9 s, and the simulation time was 800 µs.Based on the error estimation method proposed by Smirnov et al. [59,60], the maximum allowable integration steps in the present study was about 3.8 × 10 11 when the total error was controlled within 1%.The integration step in this simulation was around 8 × 10 5 .Therefore, the accumulation of numerical errors was well controlled, and the simulation was reliable.To initialize the RDC flow field, a hot spot (30 atm, 3000 K) was set at the bottom of the computational domain, as shown in Figure 15.Stoichiometric hydrogen/air mixtures with 1 atm and 300 K were located at the right of the hot spot to form the detonation.The rest of the flow field was filled with air at ambient pressure and temperature.The chemical mechanism proposed by Ó Conaire et al. [61] was used in the present study.As listed in Table 6, six cases were considered to assess the computational time cost with different acceleration algorithms.Case 0 was a benchmark with static grids of ∆x = ∆y = 0.1 mm.Cases 2-5 were calculated with the AMR.The base grid size was 0.4 mm, and the finest grid size was 0.1 mm after a two-level mesh refinement.An error indicator ϵ ∇ρ with a threshold value of 0.04 was used as adaptation criteria.Dynamic mesh balancing was used in Cases 4-5 and chemistry balancing was used in Case 1, Case 3 and Case 5.The relative CPU time in the computations is also recorded in Table 6. Cmpared with Case 0, a maximum reduction in CPU time of 11% was achieved by using the AMR and DLB in Case 5. Figure 16 plots the history of mass flow rates in Case 0 and Case 5.The mass flow rates of the inlet and outlet converged to a value around 486 kg/(m 2 • s) after 400 µs.The results of the AMR in Case 5 were close to those of the static mesh in Case 0.Moreover, the well-developed flow fields of the RDC are shown in Figure 17.Almost identical rotating detonation wave structures were obtained, which further verified the reliability of the AMR and DLB.Overall, the introduction of the AMR and DLB algorithms in the RDC simulation yielded a substantial enhancement in computational efficiency with a small compromise in the accuracy of the results.

Reactants
It should be noticed that the 2D simulation in this subsection is an approximation of an RDE.The actual RDE problems are essentially 3D, and the modeling of an RDE poses a series of challenges, involving complex interactions among shock waves, chemical reactions, turbulence and so on.For example, the existing turbulence modeling methods for detonation are relatively constrained, as the characteristics of turbulence in detonation waves are different from classical Kolmogorov turbulence due to the compressibility effects and chemical kinetics [62].The turbulence model parameters and grid resolution require careful considerations; otherwise, it may lead to inaccurate results in numerical simulations [13].The purpose of this 2D simulation was to test the computational acceleration techniques, as mentioned above.The turbulence model was beyond our scope in this study, and more realistic 3D RDE applications will be presented with BYCFoam in future work.

Conclusions
In this paper, we presented an improved rhoCentralFoam-based solver for a rotating detonation engine (RDE), BYCFoam, which is the latest version of the BYRFoam family.The following achievements were obtained: 1.
The calculation of convective fluxes was improved by choosing three independent primitive variables as reconstruction variables, along with a wide selection of flux schemes.The AUSM+M scheme demonstrated the best performance among the chosen flux schemes.

2.
The disadvantage of OpenFOAM when calculating the thermodynamic properties was fixed.The coupling with Cantera enabled us to deal with more complex reaction mechanisms, efficiently solve chemical kinetics and accurately evaluate the detailed transport properties.

3.
The adaptive mesh refinement (AMR) and the dynamic load balancing (DLB) algorithms were integrated with the solver, resulting in a substantial enhancement in computational efficiency.

Appendix A. Reaction Mechanism
The reduced C 12 H 23 /air mechanism from Liu [54] contains 18 species and 30 reactions, as listed in Table A1.

Figure 2 .
Figure 2. Specific heat capacity c p of stoichiometric DME/air mixtures.

Figure 3 .
Figure 3. Homogeneous autoignition of stoichiometric DME/air mixtures: (a) ignition delay time with different initial temperatures T 0 and initial pressures P 0 , (b) temperature evolution during autoignition with T 0 = 625 K and P 0 = 0.5 atm.

Figure 8 .
Figure 8. Laminar flame speeds of premixed hydrogen/air mixture with varying equivalence ratios.

Figure 9 .
Figure 9. One-dimensional planar flame structure of stoichiometric hydrogen/air mixture: (a) temperature and density, (b) mass fraction of OH and H 2 O 2 .

Figure 10 .
Figure 10.Numerical simulation of one-dimensional detonation tube with kerosene/air mixture at 1 atm and 373 K: (a) schematic of the detonation tube, (b) temperature and pressure distributions at different simulation times.

Figure 11 .
Figure 11.The computational domain and initial setup in the simulation of a detonation cellular structure.

Figure 14 .
Figure 14.Schematic of an RDC in a two-dimensional computational domain.

Figure 15 .
Figure15.Schematic of the initial conditions.A hot spot is located at the bottom of the computational domain to ignite the reactants.The rest of the mixture is air at 1 atm and 300 K.

Figure 16 .
Figure 16.The mass flow rates of the inlet and outlet.Solid lines: Case 0 with a static mesh.Dashed lines: Case 5 with AMR.

Figure 17 .
Figure 17.The well-developed flow fields of an RDC.(a) Temperature and (b) mass fraction of OH in Case 0. (c) Temperature and (d) mass fraction of OH in Case 5. (e) Temperature and pressure distributions in the circumferential direction at y = 0.2 cm.Solid lines: Case 0 with static mesh.lines: Case 5 with AMR.
* c is sound speed.
Flow chart within a time step ∆t F in BYCFoam.Calculations carried out in Cantera are denoted by the blue font.

Table 2 .
Initial conditions of auto-ignition.

Table 4 .
Refinement criteria used in the simulation of the detonation cellular structure.

Table 5 .
Comparisons of the detonation cell size with other numerical results.L and λ are the length and the width of the detonation cells, respectively.

Table 6 .
Test cases for the AMR and DLB algorithms.Case 0 is a benchmark without using any acceleration methods.Cases 0-1 are calculated with uniform static grids.Cases 2-5 are calculated with AMR.The CPU time costs in the computations are listed with a relative form.