Step 4

Description of Step 4

In this last benchmark, Step 3 is repeated while now activating chemical reactions from the start. To this aim, the kinetic scheme of Boivin et al. involving 9 species and 12 reversible reactions should be used in order to enable direct comparisons. The initial profiles of temperature and of mass fractions show the initial profiles for temperature, heat release, mass fractions of $H_{2}$ and of $O_{2}$ (????). As it can be observed from this figure, the species profiles are now different along the fuel/oxidizer interface compared to those in the non-reacting case (Step 3), shown in the initial profiles of temperature and of mass fractions. Now, the local species concentration corresponds to the equilibrium values obtained after initializing composition with Eqs (3) and (4) of Step 3 at initial temperature $T_{0}=300K$ and atmospheric pressure; in this manner, the initial conditions already mimic a real flame front, with local composition and temperature in agreement with each other.

Again, periodic boundary conditions are used in each direction. There is thus no net fluxes of mass, momentum or energy in any direction: the simulation domain can thus be considered to be closed. As a consequence, and contrary to Step 3, the thermodynamic pressure evolves in time during the present simulation because of the heat released by the chemical reactions and of the resulting dilatation. To enable meaningful comparisons, it is necessary to take into account this effect in the energy or temperature equation:

${\frac {\partial T}{\partial t}}+u_{j}{\frac {\partial T}{\partial x_{j}}}={\frac {1}{\rho C_{p}}}\left[{\frac {dP_{0}}{dt}}+{\frac {\partial }{\partial x_{j}}}\left(\lambda {\frac {\partial T}{\partial x_{j}}}\right)-{\frac {\partial T}{\partial x_{j}}}\sum _{k=1}^{N_{s}}\rho C_{p,k}Y_{k}\theta _{k,j}-\sum _{k=1}^{N_{s}}h_{k}{\dot {\omega }}_{k}\right]\,$ (1)

All symbols used in this equation are defined below. Usually, most simulations of low-Mach flow conditions assume that the thermodynamic pressure $P_{0}$ does not vary in time in case of open boundaries, considering that $P_{0}$ corresponds to the surrounding (usually atmospheric) pressure; therefore, $dP_{0}/dt$ vanishes in Eq.(1). On the other hand, for a closed computational domain -- as assumed here --, $P_{0}$ may change in time due to dilatation and must be evaluated at each timestep as follows:

$P_{0}={\frac {M_{0}R}{\int _{\Omega }\left(T\sum _{k=1}^{N_{s}}Y_{k}/W_{k}\right)^{-1}dV}}$ (2)

This equation is obtained by integrating the equation of state for perfect gases on the whole computational domain while assuming that the total mass is fixed and that the thermodynamic pressure is homogeneous inside the domain. In Eqs. (1) and (2), $\rho$, $u_{j}$, $T$, $Y_{k}$, $W_{k}$, $N_{s}$, $C_{p,k}$, $h_{k}$, $\theta _{k,j}$ and $\mathrm {R}$ are the density of the gas mixture, $j$-th component of flow velocity, gas temperature, $k$-th species mass fraction, $k$-th species molecular weight, number of species, $k$-th species specific heat capacity at constant pressure, $k$-th species enthalpy, $j$-th component of the species molecular diffusion velocity for each species $k$, and ideal gas constant, respectively. Finally, if the system is closed, the total mass $M_{0}$ of the mixture inside the computational domain $\Omega$ can be computed as:

$M_{0}=\int _{\Omega }\rho \,dV\,.$ (3)

or an additional ODE on mass can be integrated if some mass fluxes are present at boundaries. All other physical and numerical parameters are identical for Step 3 and Step 4, the only difference between the two benchmarks being the activation of the reaction terms. As in Step 3, the simulation should be performed for a physical time of at least $t=0.5\;\mathrm {ms} =2\tau _{\mathrm {ref} }$ and can be continued up to $t=2.5\;\mathrm {ms} =10\tau _{\mathrm {ref} }$ for a full comparison with the data available on the benchmark repository.

Because of the various time-integration techniques used in each code it is difficult to give recommendations regarding the timestep here. For example, DINO being fully explicit, it is always limited by the Fourier Condition and thus uses small timesteps. On the other hand, YALES2 can treat diffusion implicitly and thus uses convective timestep as a limiter. Nek5000 has no operator splitting and uses implicit schemes: its timestep constraints are again different. In all cases, the chemical source terms are treated by a stiff integrator which is the only solution to deal with the extremly small timescales present in flames. The conservative values $CFL=0.3$ and $Fo=0.2$ that were proposed in Step 3 can be used as a first attempt. In any case, it is recommended to check that the timestep is small enough to not influence the final result.

Aside suggestions

Should we prove that it is the case for our codes? It could be done easily, for each mesh, on a small number of iterations, say simulate 0.1ms with 2 or 3 different timesteps. This could also be done for Step 2 and 3 actually

References

TBD