# A residence time-based concept for modeling and analyzing the impact of diffusion on auto-ignition processes with thermal stratification

Zakaria Bazhar, Vincent Robin, Z. Bouali

► To cite this version:

Zakaria Bazhar, Vincent Robin, Z. Bouali. A residence time-based concept for modeling and analyzing the impact of diffusion on auto-ignition processes with thermal stratification. Physics of Fluids, 2024, 36 (10), pp.107108. 10.1063/5.0226353. hal-05016492

HAL Id: hal-05016492

https://hal.science/hal-05016492v1

Submitted on 2 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Residence Time-Based Concept for Modeling and Analyzing the Impact of Diffusion on Auto-Ignition Processes with Thermal Stratification.

Z. Bazhar,¹ V. Robin,¹ and Z. Bouali¹

Institut Pprime - UPR 3346, CNRS - ISAE-ENSMA - Université de Poitiers, Téléport 2, 1 Av. Clément Ader, Chasseneuil-du-Poitou, 86360, France

(*Electronic mail: zakaria.bazhar@ensma.fr)

(Dated: 2 September 2024)

A meticulous definition of the residence time (or fluid age) in the context of auto-ignition in a thermally stratified gas demonstrated that the fluid age can be considered as a level-set function, allowing the tracking of the a self-ignition front propagating. In order to avoid the difficulty of defining the boundary and initial conditions for this age, which can be seen as fluid birth, a normalised age is introduced. The meticulous derivation of the equation revealed additional terms related to the scalar dissipation rates and the gradient in the composition space of the auto-ignition delay. To validate this model equation, two canonical configurations were proposed: one representing a hot spot self-igniting with significant thermal diffusion and the other where a self-ignition front propagates at an almost constant speed. The analysis, based on the normalized age of the particles, reveals the impact of thermal diffusion on the acceleration or deceleration of particle aging. In certain instances, a particle rejuvenation mechanism through thermal diffusion has been identified. Furthermore, this work demonstrates the capacity of the normalized residence time field to map transitions between auto-ignition and laminar flame. Finally, the different cases studied were classified in an auto-ignition/diffusion diagram based on three non-dimensional characteristic numbers, which highlighted five combustion regimes.

## I. INTRODUCTION

The global shift towards more environmentally friendly energy sources has reignited interest in utilizing a diverse array of fuels, including sustainable options or those with low carbon footprints. Gaining a deeper insight into how these fuels behave in conventional facilities becomes fundamental for mastering combustion processes. Around two decades ago, Chiesa et al.¹ highlighted the complexities burner design, flame stability, and NO emissions involved in transitioning from natural gas to hydrogen combustion in gas turbines. More recently, Gruber et al.² performed direct numerical simulations to examine specific conditions where spontaneous ignition impacts combustion stability. Undesired ignition represents a critical phenomenon, whose repercussions on combustion efficiency, pollutant production, and safety justify extensive research efforts. This has led a numerous studies delving into diverse mixtures and conditions, see among many others²⁻⁵.

The central and most extensively investigated aspect of these ignition mechanisms is the auto-ignition delay, denoting the time required for a mixture in a given state to reach thermal runaway, driven by chemical kinetics. This delay time is used to characterize novel fuels, for instance using a rapid compression machine, see for example He et al.⁴. It also stands as the key quantity for which simplified correlations are pursued, as evidenced by the work of Zhao et al.³. Database of auto-ignition delay times also feeds machine learning methods, see for example Han et al.⁶, Huang et al.⁷, playing a central role in enhancing the comprehension and modeling strategies of complex ignition processes. Even in scenarios necessitating additional ignition characteristics such as cool flame dynamics⁸ or negative temperature coefficient behaviors⁹, the auto-ignition delay persists as a pertinent characteristic¹⁰. However, when auto-ignition processes interact with mixing processes, turbulence, compressibility effects, or heat loss, relying solely on auto-ignition time for modeling and analysis becomes insufficient. In these flows characterized by spatially and temporally varying thermochemical properties, the auto-ignition delay becomes a scalar field intricately linked to local fresh gas characteristics. It therefore requires the introduction of additional variables such as a passive scalar and a residence time to take into account of these variations in the analyses or simulations of ignition processes.

In the field of chemical engineering, Ghirelli and Leckner¹¹ established the theoretical framework underpinning the definition of the quantity of residence time and its transport equation. Residence time equations find utility in diverse applications. Shin et al.¹² employed them for statistical analysis in a turbulent round jet, while Savard et al.¹³ developed a semi-analytical model to examine the transition from spontaneous ignition to turbulent deflagration. Meanwhile, Wang et al.¹⁴ formulated a model based on decomposing the chemical rate into contributions from auto-ignition and turbulent propagation relying on the residence time concept. This approach has facilitated the independent treatment of two diffusion effects: one on passive scalar mixing and the other on the propagation of a reactive front. The presence of a diffusion term in the residence time equation remains a topic of debate. Ghirelli and Leckner¹¹ theoretically demonstrated its necessity by utilizing tracked particles as if they were experimental tracers of fluids.

Eventually, the aim of this study was to demonstrate the relevance of the concept of residence time in analyzing and/or modelling auto-ignition mechanisms during thermal mixing.

In the first section, we introduced both, a passive scalar describing the mixing between two streams of different temperatures and a residence time for which a rigorous definition of the physical meaning is provided. With these foundations, we derived
a transport equation for the normalized residence time (NRT). The aim of this theoretical development is to identify potential analogies with established concepts, such as level-set functions used for interface tracking as flame fronts¹⁵ or the Zeldovich propagation speed¹⁶. Furthermore, this mathematical approach highlights terms indicative of crucial physical mechanisms, including scalar dissipation rates (SDR). Numerous studies of combustion-related issues have highlighted the crucial importance of the SDR¹⁷⁻²¹. Bray et al.¹⁹ and Mura et al.²⁰ have even conducted detailed analyses of SDR in scenarios involving both passive and reactive scalars.

In the following three sections, we used a direct numerical simulation (DNS) tool to study three canonical scenarios: First, we examined homogeneous reactors as a validation step for the proposed NRT equation. Secondly, we studied a scenario describing hot-spot ignition in a mixture subjected to thermal mixing. Our aim was to assess the influence of the mixing mechanism on the auto-ignition delay and to verify the ability of the NRT equation to capture and analyse these effects. Finally, we explored the propagation of an auto-ignition front to demonstrate the ability of the NRT equation to track reactive interfaces. Although these scenarios may be regarded as simple problems, the corresponding mixing and chemical characteristic times yielded an ignition diagram proposal presented in the last section classifying ignition/mixing interaction regimes. Eventually, interest and limitations of the NRT equation for practical applications is discussed in conclusion.

## II. DIFFERENT FORMS OF THE RESIDENCE TIME TRANSPORT EQUATION

The main objective of this study was to evaluate the concept of residence time in a context where fresh gas simultaneously undergoes thermal mixing and auto-ignition processes. The first process was considered as binary mixing, which can be modeled using a passive scalar, while the modeling of the second process must account for auto-ignition delays. The passive scalar characterizes the thermodynamic state of the fresh gas resulting from the combination of two primary streams of identical composition but different temperatures. The auto-ignition delay times were calculated using Cantera²² as a function of the passive scalar. This section provides definitions for the passive scalar, the auto-ignition delay and the residence time, while exploring for the latter its different forms and physical meanings.

### A. Passive scalar

The use of a passive scalar in combustion problems is still widespread²⁰,²³ and rigorously defined since the work of Bilger²⁴. It significantly facilitating modeling development, albeit requiring the introduction of strong assumptions as neglecting the Dufour and Soret effects as well as the baro-diffusion effect. In the simple problems considered in this study, the composition of fresh gas was homogeneous but their temperature varied in the range corresponding to the mixing between two streams with different temperatures, $T_h$ for the highest temperature stream and $T_l$ for the lowest temperature stream. Thus, the passive scalar was defined based on the gas enthalpy $h$, including both chemical and sensible enthalpy, as follows:

$$ \xi = \frac{h - h_l}{h_h - h_l} \qquad (1) $$

Considering now low Mach number flows and unity Lewis numbers assumptions for all species, the transport equation for gas enthalpy reduces to that of a passive scalar $\xi$:

$$ \frac{\partial \xi}{\partial t} + \mathbf{u} \cdot \nabla \xi = \frac{1}{\rho} \nabla \cdot (\mathcal{F}_{\xi}) \qquad (2) $$

The density and flow velocity are respectively $\rho$ and $\mathbf{u}$, while $\mathcal{F}_{\xi}$ represents the diffusion flux of the passive scalar. For the sake of simplicity, this diffusion flux is simply modeled introducing the thermal diffusivity coefficient $\alpha$ as follow:

$$ \mathcal{F}_{\xi} = \rho \alpha \nabla \xi \qquad (3) $$

The limitations of the aforementioned hypotheses are discussed in detail in the last section of this article.

### B. Auto-ignition delay

To facilitate the exploration of various auto-ignition regimes, we selected a fuel with distinctive chemical kinetics, characterized by a sudden change in behavior at a specific fresh gas temperature. Treviño and Solorio²⁵ and Fernandez-Tarrazo et al.²⁶ highlighted that, using hydrogen, below a critical threshold, named the crossover temperature $T_c$, combustion typically manifests as
thermal runaway, characterized by gradually increasing auto-ignition times with rising fresh gas temperatures. Conversely, beyond the crossover temperature, hydrogen combustion transitions into a branched chain explosion mode, resulting in significantly reduced auto-ignition times. Results of auto-ignition calculations for stoichiometric hydrogen/air mixtures within homogeneous constant pressure reactors is presented in Figure 1-left. These computations utilized the Cantera library²² along with the Ó Conaire et al.²⁷ kinetic mechanism, which was considered the reference chemistry. This chemical kinetic scheme, comprising 9 chemical species and 40 reactions, has been subjected to comprehensive validation by Slack²⁸, across a broad range of temperatures (from 298 K to 2700 K) and pressures (from 0.05 atm to 87 atm). The section III and Figure 4 contain experimental results which validate this scheme. The temperature range investigated here spanned from 770 K to 1100 K and so includes the crossover temperature. Figure 1-right illustrates the evolution of temperature T and heat release rate $\dot{\omega}_T$ over time for a representative fresh gas temperature of 970 K. In Figure 1-left, the obtained auto-ignition delay times $\tau$ are displayed across the entire range of fresh gas temperatures $T_u$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Left: Auto-ignition delay time $\tau$ over fresh gas temperature $T_u$. Right: Temperature T and heat release rate $\dot{\omega}_T$ over time for a fresh gas temperature of 970 K</figcaption></figure>

The temperature profile in Figure 1-right presents two distinctive time scales with notably disparate orders of magnitude. One corresponds to the auto-ignition time, with a value of 0.323 ms in this instance, while the other corresponds to the heat release period, with a value of 0.02 ms. The significant difference between these two characteristic times, even at such a low ignition time, validates the modeling approach utilizing the residence time concept, which treats the reaction zone as a temporal interface between fresh and burnt gas, see next sections. Consequently, if the scales associated with heat release do not necessitate resolution, the simulation's time step can be increased by approximately 50-fold. This has significant implications, particularly for 3D simulations.

The conducted numerical experiments were intentionally selected to examine intricate scenarios where thermal mixing can play a significant role during the chemical processes that are precursors to thermal runaway. Accordingly, the characteristic diffusion time $t_D$ was selected to have comparable order of magnitude to the auto-ignition delay time $\tau$. Hence, the relevant dimensionless quantity is defined here as the ratio of these characteristic times: $Ka = \tau/t_D$. This ratio is denoted as Ka because it is similar to the Karlovitz number, which is typically defined as the ratio of the chemical time to the Kolmogorov mixing time in turbulent flames. The diffusion time $t_D$ was estimated here from the diffusivity and the characteristic length of the mixing zone $l_D$, as outlined below:

$$t_D = \frac{l_D^2}{\alpha} \qquad (4)$$

For each scenario investigated, which corresponds to a range of fresh gas temperatures bounded by those of the highest and lowest temperature mixing streams, it was necessary to introduce auto-ignition delays as a function of the passive scalar $\tau(\xi)$, as illustrated in Figure 2. Furthermore, given a fixed mixing time $t_D$, the Karlovitz number, $Ka(\xi)$, exhibits the same dependence on the passive scalar as $\tau(\xi)$ does, as illustrated in Figure 2's right axis.

Both Figure 1-left and Figure 2 reveal that the behavior of auto-ignition delays varies significantly based on whether the fresh gas temperature is above or below the crossover temperature ($T_c \approx 945$ K). Thus, when considering the temperature ranges below the crossover temperature, the Karlovitz number significantly exceeds unity, whereas above it, it is notably lower, see Figure 2.

For vanishing small values of the Karlovitz number, corresponding to fresh gas temperatures higher than the crossover temperature, thermal mixing processes become insignificant. Conversely, for exceedingly large Karlovitz numbers, mixing
<figure><img src="image_3.png" /><figcaption>FIG. 2. Ignition time $\tau$ (left axis) and Ka (right axis) versus $\xi$ obtained from homogeneous reactor computations for four different ranges of initial fresh gas temperature $T_u$.</figcaption></figure>

processes are so efficient that a homogeneous mixture is formed well before the onset of ignition. Intermediate values of the Karlovitz number indicate a concurrent interplay of both mechanisms.

When considering the auto-ignition delay as a function only of the passive scalar, its transport equation behaves like that of a passive scalar, obtained by multiplying Equation 2 by the derivative $d\tau/d\xi$:

$$ \frac{\partial \tau}{\partial t} + u \cdot \nabla \tau = \frac{1}{\rho} \frac{d\tau}{d\xi} \nabla \cdot (\mathcal{F}_{\xi}) \quad (5) $$

After a few transformations, the right hand side of Equation 5 can be restated as follows:

$$ \frac{\partial \tau}{\partial t} + u \cdot \nabla \tau = \frac{1}{\rho} \nabla \cdot (\mathcal{F}_{\tau}) - \frac{d^2 \tau}{d\xi^2} \epsilon_{\xi\xi} \quad (6) $$

where $\mathcal{F}_{\tau}$ is the diffusion flux of ignition delay time, expressed as:

$$ \mathcal{F}_{\tau} = \rho \alpha \nabla \tau \quad (7) $$

As expected, the transport equation for ignition time bears a striking resemblance to that of a passive scalar, yet interestingly it introduces an extra term associated with the familiar passive Scalar Dissipation Rate (SDR)$^{17,18,20,21}$: $\epsilon_{\xi\xi} = \alpha \nabla\xi\nabla\xi$. This supplementary term also incorporates the curvature of the tabulated function $\tau(\xi)$ which can reach a very large value around the crossover temperature. Thus, Equation (6) offers the advantage of emphasizing the pivotal role of the passive SDR in influencing the auto-ignition delay field, particularly in scenarios characterized by significant changes in chemical behavior as observed with hydrogen.

### C. Residence time

The residence time $t_r$, also referred to as fluid age in the literature$^{29}$, was the pivotal quantity under study in this work for determining both the timing and location of ignition events. Defined as the product of an infinitesimal mass of fluid and the time it spends within the computational domain, the residence time serves as a measure of fluid element age. While Ghirelli and Leckner$^{11}$ elaborated on the transport equation for tracking fluid tracers, the versatility of this concept extends to various applications. Recent studies have resolved an Eulerian transport equation for fluid age in jet configurations$^{29,30}$, unveiling self-similar behaviors. In addition, Gomet et al.$^{31}$ exploited this concept to explore flame stabilization in supersonic flows, Wang et al.$^{32}$ employed it to model auto-ignition/flame-propagation interactions in RANS, and Shin et al.$^{33}$ utilized it to analyse lift flames and propose a joint passive scalar/fluid age presumed-pdf model.
The equation governing the residence time transport for fluid particles undergoing diffusion shares similarities with a passive scalar transport equation, however, it incorporates an aging source term:

$$ \frac{\partial t_r}{\partial t} + \mathbf{u} \cdot \nabla t_r = \frac{1}{\rho} \nabla \cdot (\mathcal{F}_{t_r}) + \underbrace{1}_{\text{Aging term}} \qquad (8) $$

Based on the approach of Ghirelli and Leckner¹¹ with the unity Lewis number assumption, the diffusion flux of residence time is expressed as follows:

$$ \mathcal{F}_{t_r} = \rho \alpha \nabla t_r \qquad (9) $$

As demonstrated by Shin et al.³³ through direct numerical simulation results, the distribution of reactive scalars in a turbulent jet flame is closely linked to the joint distribution of passive scalar and residence time. In addition, they asserted that the average residence time, combined with the passive scalar, can indeed serve as an effective technique for flame mapping. However, achieving this requires a precise criterion to accurately identify the location of the auto-ignition front, based on fluid age. The simplest criterion that can be employed involves directly comparing the residence time and auto-ignition time to estimate the remaining time before ignition occurrence.

## D. Remaining time to ignition: a level-set function

The remaining time until ignition, denoted as $H$, is defined in this context as the difference between the residence time and the ignition delay time:

$$ H(x, t) = t_r(x, t) - \tau(x, t) \qquad (10) $$

As a result, $H$ assumes negative values in fresh gas and positive values in burnt gas. Thus, the iso-value $H = 0$ pinpoint the location of the ignition front within the fluid domain. This level-set function $H$ holds the capability to track ignition fronts and, in this sense, bears similarities to the familiar $G$-equation employed for tracking flame fronts¹⁵. However, while $G$ denotes the distance to the flame front, $H$ represents the time remaining until ignition.

Using the definition provided in Eq. (10), a transport equation for $H$ was derived from Eqs. (8) and (6):

$$ \frac{\partial H}{\partial t} + \mathbf{u} \cdot \nabla H = \frac{1}{\rho} \nabla \cdot (\mathcal{F}_H) + 1 + \frac{d^2 \tau}{d \xi^2} \epsilon_{\xi \xi} \qquad (11) $$

It is found that the $H$-equation closely resembles the transport equation for residence time (8), albeit with the addition of an extra source term connected to the passive SDR and the curvature of ignition delay in the composition space. Consequently, the behavior of the $H$-equation mainly corresponds with that of the residence time. Given its role in determining the ignition front position, the $H$-equation is inherently tied to the behavior of reactive scalars. This straightforward demonstration complements the findings of Shin et al.³³, who employed DNS to establish the correlation between the distribution of reactive scalars and the joint distribution of passive scalar and residence time.

Another interesting aspect of Eq. (11) lies in its potential to yield a general expression for the propagation velocity of an auto-ignition front $S$, even in the presence of significant diffusion processes. Given its nature as a level-set function, the right-hand-side of Eq. (11) also writes $S \| \nabla H \|$, leading us to derive:

$$ S = \frac{1 + \frac{1}{\rho} \nabla \cdot (\mathcal{F}_H) + \frac{d^2 \tau}{d \xi^2} \epsilon_{\xi \xi}}{\| \nabla H \|} \qquad (12) $$

Equation (12) clearly highlights the significant influence of diffusion on the propagation velocity of an auto-ignition front, even without invoking the presence of a laminar flame-type propagation mechanism. Two distinct effects of diffusion are identified. Firstly, a direct effect arises from the diffusive flux $\mathcal{F}_H$ of the remaining ignition time, capable of displacing the front. Secondly, an indirect effect emerges from the thermal mixing of fresh gas characterized by the non-linear thermodynamic properties $d^2\tau/d\xi^2$. Furthermore, Equation (12) has the advantageous property of converging toward the well-known expression of Zeldovich¹⁶ when diffusion can be regarded as negligible and the domain is closed, i.e., when the residence time is homogeneous:

$$ S = \frac{1}{\| \nabla \tau \|} \qquad (13) $$

Despite successfully deriving the $H$-equation, establishing appropriate boundary and initial conditions for the $H$-field presents a significant practical challenge. A straightforward assumption involves setting $H = -\tau$ for particles upon entering the computational domain, which corresponds to the moment of particle birth, i.e., $t_r = 0$. However, this assumption implies that the gas state
at the domain inlet, determined by the $\xi$-value, must remain consistent with the initial conditions of the preliminary calculations performed to obtain $\tau(\xi)$. Such a condition is often unattainable in practical simulations, especially in cases where injected gas lack reactivity ($\tau \to \infty$), as is common when one of the inlets is cold or when fuel and oxidizer are introduced separately. To circumvent this issue, we adopt a strategy similar to the approach taken by Wang et al.³², which involves the normalization of the residence time.

### E. Normalized residence time equation

By normalizing the residence time with the ignition delay as $\hat{t}_r = t_r/\tau$, an additional level-set function is introduced for the localization of ignition fronts:

$$ \hat{H}(x, t) = \frac{H(x, t)}{\tau(x, t)} = \hat{t}_r - 1 \qquad (14) $$

Hence, akin to $H$, $\hat{H}$ exhibits negative values in fresh gas, positive values in burnt gas, and $\hat{H} = 0$ represents the iso-line demarcating the ignition front position within the fluid domain. However, its boundary and initial value remains constant at $\hat{H} = -1$, regardless of the $\xi$ value. This characteristic also holds true for non-reactive mixtures where $\tau \to \infty$.

Using the definition provided in Eq. (14), a transport equation for the Normalized Residence Time (NRT) was derived from Eqs. (5) and (8):

$$ \frac{\partial \hat{H}}{\partial t} + u \cdot \nabla \hat{H} = \frac{1}{\tau} \left( 1 + \frac{1}{\rho} \nabla \cdot (\mathcal{F}_{t_r}) - \frac{(1+\hat{H})}{\rho} \frac{d\tau}{d\xi} \nabla \cdot (\mathcal{F}_{\xi}) \right) \qquad (15) $$

The current normalized level-set transport equation retains the aging source term, now expressed as $1/\tau$, along with the diffusion terms corresponding to both residence time and passive scalar. The aging source term incorporates the local ignition delay value, thus introducing differential fluid aging within the flow. Fluid particles traversing highly reactive regions, characterized by short ignition delays, experience accelerated aging. This term poses no closure challenges, provided that the function $1/\tau(\xi)$ is well-defined. In contrast, handling diffusion terms proves to be more intricate, and their physical impacts are analytically complex to analyse.

Therefore, the diffusion terms within equation (15), namely the last two terms on the right-hand side, had been reorganized differently. After a series of manipulations, this led to the following expression:

$$ \frac{\partial \hat{H}}{\partial t} + u \cdot \nabla \hat{H} = \frac{1}{\rho} \nabla \cdot (\mathcal{F}_{\hat{H}}) + \frac{1}{\tau} \left( 1 + 2 \frac{d\tau}{d\xi} \epsilon_{\hat{H}\xi} + (\hat{H} + 1) \frac{d^2\tau}{d\xi^2} \epsilon_{\xi\xi} \right) \qquad (16) $$

where $\mathcal{F}_{\hat{H}} = \rho\alpha\nabla\hat{H}$ is the diffusion flux and $\epsilon_{\hat{H}\xi} = \alpha\nabla\hat{H}\nabla\xi$ represents the cross scalar dissipation rate. Note that the two terms related to SDR also entail the first and second-order derivatives of $\tau(\xi)$. In summation, the derived NRT equation, beyond incorporating the source term that accounts for differential particle aging ($1/\tau$), incorporates the conventional diffusion term and introduces two supplementary terms linked to SDR.

This NRT Equation (16) is intended to be used as the key transport equation for analysis (post-processing) or modeling (LES/RANS) of auto-ignition processes. However, proposing practical closure models for passive scalar and cross scalar dissipation rates may be challenging. Therefore, in this first study of NRT, we decided to focus our analyses on a modelled NRT equation, referred to as NRT-M, which is easier to use in practice:

$$ \frac{\partial \hat{H}}{\partial t} + u \cdot \nabla \hat{H} = \frac{1}{\rho} \nabla \cdot (\mathcal{F}_{\hat{H}}) + \frac{1}{\tau} \qquad (17) $$

Neglecting these SDR terms may seem unjustified, but proposing complex models to account for their effects is also difficult to justify, as the assumptions used to introduce the residence time itself may lead to biases that these models cannot compensate for. Further discussion of these assumptions can be found in the final section of this article. Thus, this study focuses on the two terms of the RHS of Eq. (17). To emphasize the effects of residence time mixing process, we also investigated Eq. (17) but excluding the diffusion term (referred to as NRT-WD):

$$ \frac{\partial \hat{H}}{\partial t} + u \cdot \nabla \hat{H} = \frac{1}{\tau} \qquad (18) $$

The results provided by the two transport equations, (17) and (18), are now analyzed in detail using canonical DNS configurations. These configurations have been meticulously devised with the explicit objective of validating the theoretical concept of the NRT model and its capacity to supplant the detailed kinetic schemes employed in CFD. Furthermore, they are intended to facilitate detailed analyses of the physical mechanisms under perfectly controlled conditions.
### III. VALIDATION WITH HOMOGENEOUS REACTOR

Numerical simulations of homogeneous isobaric reactors were used as a validation step to confirm the correct implementation of the NRT Equations in the DNS solver Asphodele³⁴⁻³⁶. This solver is a research code which solves the standard set of gaseous reacting Navier-Stokes equations and perfect gas state law within the low-Mach number approximation. The variables solved include density $ρ$, momentum components $ρu_i$, mass fractions of each chemical species $Y_k$, and sensible enthalpy $h_s$. The passive scalar $ξ$ and NRT $\hat{H}$ have also been implemented. Spatial derivatives are computed with 4th order precision finite difference schemes on regular grids, while time integration utilizes a minimal storage third-order explicit Runge-Kutta scheme. Chemical kinetics and heat release rates are determined using the SUNDIALS CVODE solver³⁷ and hydrogen kinetics from Ó Conaire et al.²⁷. A mesh sensitivity study was conducted and results gathered in Annex A. In this study, the diffusion coefficient is computed under the unity Lewis number hypothesis, which is necessary to be consistent with the assumptions required for the use of a passive scalar and residence time.

Figure 3 shows time evolution of temperature in red and $\hat{H}$ in dashed black, both computed from Asphodele, for two different initial unburnt gas temperatures, $T_u = 970$ K and $T_u = 900$ K, respectively above and below the crossover temperature.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Time evolution of temperature (red) and NRT $\hat{H}$ (black) in homogeneous reactor for $T_u = 970$ K (left) and $T_u = 900$ K (right).</figcaption></figure>

In both cases, these results demonstrate that $\hat{H}$ evolves linearly with time and approaches zero, as indicated on the right-hand scale, precisely at the moment when the temperature gradient reaches its maximum, which corresponds to the ignition instant. Therefore, this result validates the implementation of the $\hat{H}$ equation, as it accurately predicts the ignition delay.

Figure 4 compares the auto-ignition delay obtained at atmospheric pressure by solving the NRT equation with the values provided by the Cantera library²². As anticipated, this comparison yielded excellent results across the entire temperature range studied here. Figure 4-right, which corresponds to a case where the pressure is 2 atm, includes the experimental results of Slack²⁸, which serve to validate the chemical mechanism proposed by Ó Conaire et al.²⁷. This chemical scheme is a reference for hydrogen combustion and has been widely cited, see among others Varga et al. ³⁸ and Saeid et al. ³⁹. The $\hat{H}$ equation is also capable to handle a other fuels and mixtures, provided that the definition of a mixture fraction is appropriate. The results obtained with a stoichiometric methane-air mixture at atmospheric pressure and $T_u = 1500$ K using the GRI3.0⁴⁰ chemical scheme is presented in Annex B.

However, more complex configurations were also necessary to illustrate the relevance of the NRT equation in presence of diffusion.

### IV. HOT SPOT IGNITION DURING MIXING PROCESSES

#### A. Initial temperature profile

This scenario depicted a localized region of intensely heated fresh gas, commonly referred to as a 'hot spot', which undergoes mixing with the surrounding cooler gas until it spontaneously ignited.

Initially, a one-dimensional computational domain spanning 500 cells along the x-axis (with x ranging from 0 to L) was employed, featuring periodic boundary conditions. This domain was initially filled with fresh gas at two distinct temperatures:
<figure><img src="image_5.png" /><figcaption>FIG. 4. (Left): Auto-ignition delay times $\tau$ for stoichiometric H$_2$/air at 1 atm, calculated using the Cantera$^{22}$ and the NRT equation. (Right): Auto-ignition delay times $\tau$ for stoichiometric H$_2$/air at 2 atm, comparing Cantera results, NRT equation results with experimental data from Slack$^{28}$.</figcaption></figure>

$T_l$ = 770 K and $T_h$ = 900 K, as illustrated in Figure 5-left. It is worth noting that the most reactive mixture lies slightly below the crossover temperature.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Left: Temperature profile used as initial condition for the hot spot scenario with $L = 2l_D$. Right: Karlovitz number as function of passive scalar for four cases: $l_D$=1.25 mm (Ka$_h$ = 100), $l_D$=2.5 mm (Ka$_h$ = 30), $l_D$=5 mm (Ka$_h$ = 9), $l_D$=7.5 mm (Ka$_h$ = 3.64) and $l_D$=40 mm (Ka$_h$ = 0.1).</figcaption></figure>

The characteristic length for the diffusion process was approximated as half of the domain length ($l_D = L/2$). Consequently, the Karlovitz number is influenced by both the passive scalar, as illustrated in Figure 2, and the physical domain's length, which in turn depends on the diffusion time.

In this context, we conducted an analysis of four distinct cases with identical initial temperature profiles but varying domain lengths ($l_D$ = 1.25 mm, $l_D$ = 2.5 mm, $l_D$ = 5 mm, and $l_D$ = 40 mm). The range of variation in the Karlovitz number for each case is presented in Figure 5-right.

For this hot spot configuration, the Karlovitz number of utmost significance is associated with the most reactive mixture, denoted as Ka$_h$ and corresponding to $\xi = 1$. It is crucial to keep in mind that heterogeneity relates to temperature, while composition consistently remains homogeneous. Based on this value, three distinct regimes emerge:

* **Homogeneous auto-ignition** characterized by rapid diffusion ($Ka_h > 1$) leading to a homogeneous mixture well before ignition. It corresponds to case Ka$_h$ = 100 (green line), see Figure 5-right, and case Ka$_h$ = 30 (black line).
* **Heterogeneous auto-ignition** where ignition occurs significantly before the onset of mixing ($Ka_h < 1$). It corresponds to case $Ka_h = 0.1$ (orange line).

* **Diffusion-altered homogeneous auto-ignition** is observed for $Ka_h$ values closer to unity but still above it. It corresponds to case $Ka_h = 9$ (blue line), and case $Ka_h = 3.64$ (purple line).

## B. Homogeneous auto-ignition: $Ka_h = 100, Ka_h = 30$

<figure><img src="image_7.png" /><figcaption>FIG. 6. Temperature profiles (left) at three instant ($t$={0 s; 0.015 s; 0.962 s}) and temporal evolution (right) of temperature and NRT $\hat{H}$ at three specific locations ($x_1$=1 mm, $x_2$=3 mm, $x_3$=5 mm) for the case $Ka_h = 30$. (Multimedia available online)</figcaption></figure>

Figure 6-left (Multimedia available online) presents temperature profiles at three different time points and (right) the temporal evolution of temperature and the NRT-M $\hat{H}$ at three specific locations ($x_1$ = 1 mm, $x_2$ = 3 mm, $x_3$ = 5 mm) obtained for the case $Ka_h = 30$. The temperature profile became completely flat at approximately $t \approx 0.05$ s, whereas thermal runaway was observed at around $t \approx 0.95$ s. The NRT exhibited linear growth with time, and thermal runaway coincided with $\hat{H} = 0$. The swift achievement of a homogeneous mixture renders this case highly similar to the homogeneous reactor scenario discussed in the previous section. Similar results were obtained for the case $Ka_h = 100$.

Nevertheless, even in cases where the mixing process is rapid, a mixing phase exists, see figure 6-right. Thus, it is intriguing to explore the impact of the diffusion term within the NRT-M Eq. (17) by conducting a comparative analysis with results derived from the NRT-WD Eq. (18).

Figure 7 illustrates the temporal evolution of the NRT at the center of the domain ($X/L = 0.5$) obtained using Eq. (17) in comparison to that obtained using Eq. (18). The graph on the left corresponds to the case with $Ka_h = 30$, while the graph on the right corresponds to the case with $Ka_h = 100$.

As anticipated, for extremely high values of the Karlovitz number, the diffusion term in Eq. (17) has no discernible effect. However, in the case with $Ka_h = 30$, a significant difference between the two solutions becomes apparent. This outcome underscores the non-negligible influence of the diffusion term, which improves the solution, even in situations where the mixing mechanism seems to be negligible. As a result, the auto-ignition time obtained using Eq. (17) closely aligns with the time obtained using the reference chemical mechanism.

For a more comprehensive and precise analysis of the mechanisms involved during the initial mixing process, we also explored the cases with Karlovitz number closer to unity.

## C. Diffusion-altered homogeneous auto-ignition: $Ka_h = 9, Ka_h = 3.64$

For a better understanding of the effects of thermal diffusion, Figure 8 shows two different representations of the temporal evolution of the NRT $\hat{H}$: the t-x diagram (top) and the $(\hat{H}, \xi)$ point cloud (bottom). These representations are for the case $Ka_h = 9$ and are obtained using either the NRT-M Eq. (17) (left) or the NRT-WD Eq. (18) (right). These representations highlight the very different behavior depending on whether diffusion is taken into account in the NRT $\hat{H}$ equation. Thermal diffusion leads to rapid homogenization of the mixture fraction field ($\xi = 0.58$ throughout the entire volume at around 0.3 seconds), whereas
<figure><img src="image_8.png" /><figcaption>FIG. 7. Temporal evolution of temperature and NRT $\hat{H}$ at the center of the domain ($X/L = 0.5$) obtained using both modeled Equation NRT-M, Eq. (17) and Equation NRT-WD, Eq. (18) with $K_{a_h} = 30$ (left) and $K_{a_h} = 100$ (right)</figcaption></figure>

the $\hat{H}$ field remains heterogeneous if the diffusion term is not included. Thus, when the NRT-M Eq. (17) is used, all points ($\hat{H}, \xi$) in the cloud converge to a single point before auto-ignition ($\hat{H} = 0$); auto-ignition is homogeneous within the volume in this case. Conversely, when using the NRT-WD Eq. (18), the $\hat{H}$ field remains heterogeneous when the first points self-ignite, resulting in a self-ignition front propagating. Therefore, the physical behavior of the gas obtained with or without a diffusion term is completely different.

Figure 9 (Multimedia available online) compares the temporal evolution of NRT $\hat{H}$ for the case with $K_{a_h} = 9$ (top) and $K_{a_h} = 3.64$ (bottom) obtained using both the modeled Equation NRT-M (left graph, Eq. (17)) and modeled Equation without diffusion NRT-WD (right graph, Eq. (18)). Once again, auto-ignition occurs homogeneously throughout the entire domain in both scenarios ($K_{a_h} = 9$ and $K_{a_h} = 3.64$) if the diffusion term is included. The initial temperature profiles are identical in both cases, but the mixing phase is prolonged in the second scenario. Since the mixture retains its heterogeneity for a more extended period in the second case, the resulting auto-ignition delay is slightly shorter. This straightforward test case underscores the influence of heterogeneity on auto-ignition delay. The duration during which fluid particles exist in the most reactive state significantly impacts the delay. This characteristic is one that the NRT field can replicate. Consequently, the discrepancy between the solutions from the two models — i.e., with or without diffusion terms — becomes quite significant. By omitting the diffusion term in the NRT equation, the heterogeneity of the $\hat{H}$ field is maintained over time, illustrating a propagating auto-ignition front rather than homogeneous auto-ignition. On the other hand, when the diffusion term is included, NRT-M model consistently produces excellent results.

Furthermore, the behavior of NRT-M model is notable: at the center of the domain, where the mixture is most reactive, $\hat{H}$ initially increases and reaches a maximum before decreasing to a minimum just before the mixture becomes well-mixed and increases again. This maximum indicates that particle age can decrease as it mixes with cooler gas, suggesting the presence of a diffusion-induced rejuvenation mechanism. In the absence of the diffusion term, this mechanism cannot be taken into account, as $\hat{H}$ can only increase, as dictated solely by the aging term.

To highlight this rejuvenation mechanism, Figure 10 compares the temporal evolution of the chemical aging source term $1/\tau$ with the diffusion term $1/\rho\nabla\cdot(F_{\hat{H}})$ at the center of the domain ($X/L = 0.5, \xi_{ini} = 1$, left graph) and near the left-side of the domain ($X/L = 0.1, \xi_{ini} \approx 0$, right graph). Both terms have absolute values higher at the center of the domain where the hottest region is located. As the chemical aging source term is significantly higher in this region at the beginning of the mixing process ($\xi \approx 1$), $\hat{H}$ increases more rapidly at the center than on the cold sides ($\xi \approx 0$). Thus, as the $\hat{H}$ gradient increases, the diffusion process is set in motion, leading to negative diffusion term at the center (so contributing to decrease $\hat{H}$) and positive diffusion term on the sides (so contributing to increase $\hat{H}$). These results show that the diffusion term (in absolute values) can be greater than the chemical term (ageing term) itself, thus contributing to accelerate or decelerate ageing of fluid particles. The diffusion-induced decelerated ageing (negative diffusion term) can even lead to rejuvenation if the RHS of Eq. (17) becomes negative.

To highlight this rejuvenation mechanism, Figure 10 compares the temporal evolution of the chemical aging source term $1/\tau$ with the diffusion term $1/\rho\nabla\cdot(F_{\hat{H}})$ at the center of the domain ($X/L = 0.5, \xi_{ini} = 1$, left graph) and near the left side of the domain ($X/L = 0.1, \xi_{ini} \approx 0$, right graph). Both terms have higher absolute values at the center of the domain, where the hottest region is located. As the chemical aging source term is significantly higher in this region at the beginning of the mixing process ($\xi \approx 1$), $\hat{H}$ increases more rapidly at the center than on the cold sides ($\xi \approx 0$). Thus, as the $\hat{H}$ gradient increases, the diffusion process is set in motion, leading to a negative diffusion term at the center (thus contributing to a decrease in $\hat{H}$) and a positive
<figure><img src="image_9.png" /><figcaption>FIG. 8. Temporal evolution of the NRT $\hat{H}$, presented as a t-x diagram (top) and $(\hat{H}, \xi)$ point cloud (bottom), obtained with the NRT-M Eq. (17) (left) and NRT-WD Eq. (18) (right) for case $Ka_h = 9$</figcaption></figure>

diffusion term on the sides (thus contributing to an increase in $\hat{H}$). These results show that the diffusion term (in absolute values) can be greater than the chemical term (aging term) itself, thereby contributing to either accelerate or decelerate the aging of fluid particles. The diffusion-induced decelerated aging (negative diffusion term) can even lead to rejuvenation if the RHS of Eq. (17) becomes negative. This RHS is denoted here as the aging factor $\Theta = (\tau^{-1} + \rho^{-1}\nabla \cdot (F_{\hat{H}}))$. Figure 11 shows the temporal evolution of the Favre average $\tilde{\Theta}$ (black), minimum $\Theta_{min}$ (red) and maximum $\Theta_{max}$ (blue) values of the aging factor with (left) and without (right) diffusion terms. Thus, these curves give a representation of the evolution of the aging factor distribution.

Without the age-diffusion mechanism, the aging factor distribution is initially large, indicating that the mixture is thermally heterogeneous: $\Theta_{min}$ and $\Theta_{max}$ are significantly different. This heterogeneity progressively disappears as the mixture becomes homogeneous ($\Theta_{min}$ and $\Theta_{max}$ converges). Eventually, the aging factor reaches a constant value throughout the entire domain.

On the other hand, accounting for the age-diffusion mechanism provides a more rapid convergence of $\Theta_{min}$ and $\Theta_{max}$. However, just before full convergence is obtained, at $t \approx 0.015s$, these extreme values of the aging factor reverse their behavior, causing $\Theta_{min}$ and $\Theta_{max}$ to begin diverging. The minimum $\Theta_{min}$ even reaches negative values, indicating a rejuvenation process for some gas particles. All curves eventually converge to the mean value.

Eventually, age-diffusion mechanism leads to accelerated particle aging in some region of the flow, decelerated aging in some other and even particle rejuvenation processes.

### D. Heterogeneous auto-ignition:

If the Karlovitz number of the most reactive mixture ($Ka_h$) is reduced to below unity, it will ignite first before the mixing process becomes significant. However, in this case, different scenarios may arise based on the Karlovitz number of the less reactive mixture ($Ka_l$). Consequently, analyzing $Ka_h$ alone is no longer sufficient to identify the auto-ignition/diffusion interaction regime. When the Karlovitz function $Ka(\xi)$ intersects the line $Ka = 1$, as illustrated in Figure 5 for $Ka_h = 0.1$ ($l_D = 40$ mm), the mixing process predominates after ignition of the most reactive mixture. Hence, two scenarios may unfold: either the Karlovitz number of the less reactive mixture is significantly larger than unity ($Ka_l >> 1$), implying that the gas is too cold for ignition,
<figure><img src="image_10.png" /><figcaption>FIG. 9. Temporal evolution of temperature and NRT $\hat{H}$ at three specific locations ($x_1$ = 1 mm, $x_2$ = 3 mm, $x_3$ = 5 mm) for the cases $Ka_h$ = 9 (Top) and $Ka_h$ = 3.64 (Bottom) obtained using both modeled Equation NRT-M (left) and NRT-WD (right). (Multimedia available online)</figcaption></figure>

allowing sufficient time for the stabilization of a laminar flame structure (refer to Section V B), or $Ka_l$ remains close to unity, resulting in the propagation of an auto-ignition front altered by diffusion (refer to Section V A).

Nevertheless, the initial hot-spot canonical configuration is inadequate for examining scenarios involving the propagation of auto-ignition fronts. Consequently, we have opted to introduce a second canonical configuration, characterized by an initial temperature profile that induces the propagation of an auto-ignition front at a constant speed.

# V. CONSTANT SPEED AUTO-IGNITION FRONTS

## A. Auto-ignition front propagation altered by diffusion

Defining an initial condition to uphold a constant front propagation velocity in the Zeldovich sense presents a challenge, as hydrodynamic interactions during the initial stages of propagation can significantly alter the initially imposed temperature gradient in the remaining fresh gas. Nevertheless, this approach remains the most pertinent for identifying the specific impacts of diffusion on front velocity.

Recently, Taileb et al.⁴¹ investigated a canonical initial condition, akin to the one examined here, with a specific focus on compressibility effects. They explored scenarios characterized by constant and high propagation speeds of the fronts, opting to neglect diffusion mechanisms. The primary objective in that study was to examine the mechanisms potentially leading to transitions towards detonation. In the present context, the focus shifts to the influence of diffusion mechanisms.

Setting such a specific initial condition required a preliminary procedure. First, following the approach outlined by Zeldovich ¹⁶,
<figure><img src="image_11.png" /><figcaption>FIG. 10. Temporal evolution of the chemical aging source term $1/\tau$ with the diffusion term $1/\rho \nabla \cdot (\mathcal{F}_H)$ at $X/L = 0.5$ (left) and at $X/L = 0.1$ (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11. Time evolution of the Favre average $\tilde{\Theta}$ (black), minimum $\Theta_{min}$ (red) and maximum $\Theta_{max}$ (blue) values of the aging factor $\Theta$ corresponding to the RHS of the NRT-M Eq.(17) (left) and NRT-WD Eq.(18) (right)</figcaption></figure>

we deduced from Eq. (13) the spontaneous ignition front velocity in the normal direction as:

$$S_n = \left( \frac{d\tau}{dn} \right)^{-1} \qquad (19)$$

Integrating this equation provided the profile of the ignition delay time along the normal direction. Assuming a constant velocity, we derived:

$$\tau(n) = \frac{n - n_h}{S_n} + \tau_h. \qquad (20)$$

Here, $n_h$ represents the location of the first ignited point in the normal coordinate. By inverting the known tabulated function $\tau(\xi)$ and combining it with Eq. (20), we obtained the passive scalar profile $\xi(n)$ and, consequently, the fresh gas temperature profile $T_u(n)$ along the normal direction of the auto-ignition front, for which the speed $S_n$ remains constant.

Figure 12 illustrates the initial temperature and ignition delay profiles corresponding to three different auto-ignition front speeds: $S_n = \{25.82; 15; 11.36\} \text{ m s}^{-1}$. The domain length is fixed at $L = 10$ mm, and the temperature of the less reactive mixture, denoted as $T_l = 970$ K, is set higher than the cross-over temperature. Thus, the last point of ignition occurs precisely at the domain boundaries and the corresponding Karlovitz number approaches $Ka_l \approx 0.12$.
<figure><img src="image_13.png" /><figcaption>FIG. 12. Initial profiles of temperature (left) and ignition delay time (right) corresponding to three different auto-ignition front speeds $S_n$ = {25.82 ; 15 ; 11.36} m s⁻¹</figcaption></figure>

As anticipated, the initial ignition delay time profile consists of straight lines, a characteristic indicative of a constant Zeldovich's speed, Eq. (13). The scale on the right-hand side of Figure 12's right graph highlights the corresponding Karlovitz range.

<figure><img src="image_14.png" /><figcaption>FIG. 13. Temperature profiles (left) at three different times $t = \{0 ; 0.395 ; 0.474\}$ ms and time evolution (right) of temperature and NRT $\hat{H}$, using NRT-M model Eq. (17), for three different locations $x_1, x_2, x_3$, for $S_n = 25.82$ m s⁻¹. (Multimedia available online)</figcaption></figure>

The left graph of Figure 13 (Multimedia available online) presents the profile of temperature at three different time points and the right graph the time evolution of temperature and NRT $\hat{H}$ at three different locations for the case where the initial profile of temperature is set to obtain $S_n = 25.82$ m s⁻¹. For these very low values of the Karlovitz number, the mixing term of NRT is negligible until the front reach the domain's boundaries and the same results are obtained using either the modeled equations NRT-M Eq. (17) or NRT-WD Eq. (18). The time evolution of $\hat{H}$ is observed to be quasi-linear regardless of the x-position, and the agreement between the value $\hat{H} = 0$ and the front position, based on the temperature profile, is excellent, see Figure 13-right. A small bump appears on the profile due to expansion-induced convection.

Figure 14 shows the $(\hat{H}, \xi)$ point cloud obtained with the NRT-M Eq. (17) for the three initial temperature profiles considered ($S_n$ = {25.82 ; 15 ; 11.36} m s⁻¹). This behavior significantly contrasts with the diffusion-altered homogeneous auto-ignition observed in Figure 8. In these cases, diffusion in the fresh gas (characterized by the passive scalar) is negligible: $\xi$ remains between 0 and 1, or nearly so, throughout the calculation. On the other hand, the effect of the aging differential is very significant on the slope of the isochrones when the propagation velocity is low.

Figure 15 depicts x-t diagrams illustrating the displacement of the front based on the maximum heat release position in red, along with the movement of the $\hat{H} = 0$ isovalue, in dashed black, for three different initial profiles. The blue line represents the
<figure><img src="image_15.png" /><figcaption>FIG. 14. ($\hat{H}, \xi$) points cloud obtained with the NRT-M Eq. (17) for initial temperature profiles associated to $S_n = 25.82 \text{ m s}^{-1}$, $S_n = 15 \text{ m s}^{-1}$, and $S_n = 11.36 \text{ m s}^{-1}$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 15. Temporal evolution of position at the critical points ($\hat{H} = 0$ and maximum of heat release rate $\dot{\omega}_T$) depicted on the x-t diagram for initial temperature profiles associated to $S_n = 25.82 \text{ m s}^{-1}$, $S_n = 15 \text{ m s}^{-1}$, and $S_n = 11.36 \text{ m s}^{-1}$, along with the initial profile from the slope (Sn)</figcaption></figure>

displacement that would be observed if the front strictly followed the profile imposed in the initial condition.

While this canonical numerical experiment emphasizes that we are still far from capturing the complexity of a practical flow, its analysis highlights the intricate interplay between thermal diffusion and chemical reactions.

Given that auto-ignition initiates at the center of the domain, the propagation consistently commences at $X/L = 0.5$ for all cases. However, the front initiates sooner for lower propagation speeds, as illustrated in Figure 15. This phenomenon can be attributed to the initial temperature profile, where the temperature at the center is higher for cases with lower propagation speeds, as depicted in Figure 12.

Secondly, Figure 12 demonstrates that the initial temperature gradients are particularly pronounced at the center of the domain, resulting in localized mixing prior to the onset of self-ignition. This pre-ignition mixing induces a slight modification to the initially set temperature profile, resulting in a front speed during the initial stages of propagation that is faster than expected. This acceleration is evident in the divergence between the blue and red lines in Figure 15, with a more significant effect observed at

<figure><img src="image_17.png" /><figcaption>FIG. 16. Temporal evolution of the maximum value of the heat release rate $\max(\dot{\omega}_T)$ and the corresponding diffusion term $\text{Diff}_T$ of temperature, taken at the same location, for initial temperature profiles $S_n = 25.82 \text{ m s}^{-1}$, $S_n = 15 \text{ m s}^{-1}$, and $S_n = 11.36 \text{ m s}^{-1}$</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 17. Temporal evolution of position at the critical points ($\hat{H} = 0$ and maximum of heat release rate $\dot{\omega}_T$) depicted on the x-t diagram for initial temperature profile associated to $S_n = 161.5 \text{ m s}^{-1}$, along with the initial profile from the slope (S$_n$)</figcaption></figure>

lower propagation speeds when the initial temperature gradient at the center is larger. While diffusion plays a non-negligible role in this central zone, the $\hat{H} = 0$ isovalue consistently aligns with the front position, regardless of its propagation speed, showcasing its remarkable ability to track auto-ignition processes even in the presence of diffusion mechanisms.

Nonetheless, outside the narrow central zone, the red and blue lines run parallel, indicating that the front's velocity aligns with that of the initial condition. Consequently, the front speed requires some time to stabilize at the value dictated by the initial profile. However, at low propagation speeds, we observe an interesting behavior: the position of the $\hat{H} = 0$ isovalue (dashed back line) tends to return completely to the initial condition (blue line), while the front (red line) continues its trajectory and deviates from the position of the $\hat{H} = 0$ isovalue. This mechanism is characteristic of the early stages of the transition to a laminar flame, where propagation is no longer solely driven by the temperature profile in fresh gas, but by diffusion between the burnt and fresh gas.

This characteristic is clearly demonstrated when comparing the maximum value of the heat release rate of the temperature equation with the corresponding diffusion term at the same position, see Figure 16. The greater the velocity of the front, the smaller the initial jump in the value of the diffusion term, which occurs when the fresh gas/burnt gas front is established, just after ignition. As expected, this diffusion term becomes increasingly negligible as the velocity of the front increases. Nevertheless, it is not zero, and the mere comparison of these terms does not permit a clear differentiation between a self-ignition front and a laminar flame. The divergence of the $\hat{H} = 0$ position ($X_{\hat{H}=0}$) from the flame position, see Figure 15, on the other hand, provides a precise quantification of the alteration of the self-ignition front by diffusion between fresh and burnt gases.

Figure 16 also shows that the values of chemical source term and the diffusion term do not stabilize because the temperature of the fresh gases changes as the propagation progresses. The last case presented in the Figure 16, i.e. $S_n = 11.36 \text{ m s}^{-1}$, is a case close to the limit of the transition between an auto-ignition front and a laminar flame whose propagation speed in these fresh gas conditions is close to $8.5 \text{ m s}^{-1}$ (see next section).

In circumstances where both Karlovitz numbers are markedly less than one ($K_{a_h} = 2.9 \times 10^{-3}$ and $K_{a_l} = 1.01 \times 10^{-2}$), the influence of diffusion on the auto-ignition front is negligible, see Fig. 17-left. In this scenario, the position of the maximum heat release rate coincides with that of the $\hat{H} = 0$ position, as illustrated in Figure 17-right.

The lower the propagation speed of the auto-ignition front, the more significant the role of diffusion mechanisms characteristic of laminar flames becomes. In this manner, the NRT field precisely detects the onset of the transition to a laminar flame. However, due to the short distances and times examined in these cases, achieving a stabilized laminar flame is impossible. Accordingly, the regime emphasized here is commonly referred to as a 'diffusion-altered auto-ignition front propagation' characterized by Karlovitz values below but close to unity.

Finally, the last regime yet to be studied is termed the 'heterogeneous ignition and flame propagation', wherein the transition to a laminar flame occurs rapidly in the computational domain.

## B. Heterogeneous ignition and flame propagation

Figure 18-left (Multimedia available online) illustrates the temporal evolution of temperature and $\hat{H}$ at several positions within the domain for the scenario where the initial temperature profile is set to obtain an auto-ignition front of $1.3 \text{ m s}^{-1}$, see the x-t
<figure><img src="image_19.png" /><figcaption>FIG. 18. Temporal evolution of temperature and NRT $\hat{H}$ at three specific locations using NRT-M model (left), temporal evolution of x-t for $S_n$ = 1.307 m s$^{-1}$. (Multimedia available online)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIG. 19. Temporal evolution of the maximum value of the heat release rate max($\dot{\omega}_T$) and the corresponding diffusion term Diff_T of temperature, taken at the same location, for initial temperature profile $S_n$ = 1.307 m s$^{-1}$.</figcaption></figure>

diagrams of Figure 18-right. Such a velocity is significantly lower than the laminar flame speed that can be obtained within the fresh gas, i.e. a premixed mixture near 900 K, i.e. $S_l \approx 8.5$ m s$^{-1}$. As already proved in previous section, the $\hat{H}$-field accurately predicts the ignition time at the domain's center but quickly diverges from the flame front position whose structure is a laminar flame structure with displacement speed near 10 m s$^{-1}$.

Figure 18-right illustrates the propagation of the fronts defined by the maximum heat release rate (red line) and the propagation of the isovalue $\hat{H} = 0$ (black dashed line). The divergence between the two curves, precisely indicating the time and position of the transition from self-ignition to laminar flame propagation, occurs in this case right at the initiation of propagation when auto-ignition takes place. The red curve reflects the propagation of a laminar flame, with its slight curvature attributed to mixing in the fresh gases that is not yet completely complete, as shown in Figure 19-right. Figure 19-left shows essentially the same result as Figure 16, but the values of the maximum heat release and diffusion terms becomes stable. This is a characteristic of quasi-stabilized laminar flame.

## VI. DIAGRAM SUMMARISING THE AUTO-IGNITION/DIFFUSION REGIMES AND DISCUSSIONS

In Paragraph IV, we illustrated the utility of the NRT-M transport equation, Eq. (17), in accounting for particle rejuvenation effects during auto-ignition in a hot gas mixing with a cooler gas. Moving on to Paragraph V, we demonstrated its applicability in tracking propagating auto-ignition fronts. And eventually, in situations where laminar flames can propagate after auto-ignition, the
$\dot{H}$-field is capable of tracking the front's evolution until a laminar flame structure emerges. Consequently, an auto-ignition/flame map can be derived from this field. Therefore, the suggested canonical numerical experiments have enabled us to assess the relevance of the NRT concept and its diffusive transport equation. Furthermore, the analysis has led to the identification of various regimes governing the interaction between thermal diffusion mechanisms and auto-ignition fronts.

To summarize these regimes, we have introduced a diagram based on three characteristic interaction times: diffusion time, hot gas ignition time, and cold gas ignition time. From these characteristic times, three dimensionless numbers are derived, including the two Karlovitz numbers $Ka_h$ and $Ka_l$, and a parameter $R_{Ka}$ indicative of the heterogeneity of the gaseous mixture:

$$Ka_h = \frac{\tau_h}{t_D}, \quad Ka_l = \frac{\tau_l}{t_D}, \quad R_{Ka} = \frac{Ka_l}{Ka_h} = \frac{\tau_l}{\tau_h} \qquad (21)$$

The diffusion time is evaluated from Eq. (4). The heterogeneity number is the ratio of the two Karlovitz numbers.

<figure><img src="image_21.png" /><figcaption>FIG. 20. Combustion diagram for auto-ignition with thermal stratification. Each red dot represents a simulation case</figcaption></figure>

<table><caption>TABLE I. Summary of simulations: $Ka_h$, $Ka_l$ and $R_{Ka}$ values with corresponding regimes</caption><thead><tr><th>Ka<sub>h</sub></th><th>Ka<sub>l</sub></th><th>R<sub>Ka</sub></th><th>Regime</th></tr></thead><tbody><tr><td>1.31 x 10<sup>2</sup></td><td>9.84 x 10<sup>4</sup></td><td>7.53 x 10<sup>2</sup></td><td>Homogeneous auto-ignition</td></tr><tr><td>3.27 x 10<sup>1</sup></td><td>2.46 x 10<sup>4</sup></td><td>7.53 x 10<sup>2</sup></td><td>Homogeneous auto-ignition</td></tr><tr><td>8.17</td><td>6.15 x 10<sup>3</sup></td><td>7.53 x 10<sup>2</sup></td><td>Diffusion-altered homogeneous auto-ignition</td></tr><tr><td>3.64</td><td>2.73 x 10<sup>3</sup></td><td>7.53 x 10<sup>2</sup></td><td>Diffusion-altered homogeneous auto-ignition</td></tr><tr><td>7.70 x 10<sup>-2</sup></td><td>1.22 x 10<sup>-1</sup></td><td>1.59</td><td>Auto-ignition front propagation altered by diffusion</td></tr><tr><td>4.83 x 10<sup>-2</sup></td><td>1.22 x 10<sup>-1</sup></td><td>2.52</td><td>Auto-ignition front propagation altered by diffusion</td></tr><tr><td>1.89 x 10<sup>-2</sup></td><td>1.21 x 10<sup>-1</sup></td><td>6.42</td><td>Auto-ignition front propagation altered by diffusion</td></tr><tr><td>2.9 x 10<sup>-3</sup></td><td>1.0 x 10<sup>-3</sup></td><td>3.46</td><td>Auto-ignition front propagation</td></tr><tr><td>1.89 x 10<sup>-2</sup></td><td>6.15 x 10<sup>3</sup></td><td>3.25 x 10<sup>5</sup></td><td>Heterogeneous ignition and flame propagation</td></tr></tbody></table>

In Fig. 20, the diagram depicts the heterogeneity number $R_{Ka}$ against the Karlovitz number of the most reactive mixture, $Ka_h$. This diagram reveals five distinct auto-ignition/diffusion regimes. The table I summarizes the key simulations along with their $Ka_h$, $R_{Ka}$ values, and the corresponding regime.
The first regime, located on the right side of the diagram in Fig. 20 and along the bottom axis, corresponds to 'homogeneous auto-ignition'. In this regime, either the diffusion time is small compared to the two auto-ignition times, or the heterogeneity number vanishes. In the former case, both Karlovitz numbers exceed unity. It is characterized by a well-mixed mixture that undergoes auto-ignition in the whole volume simultaneously without significant influence from the rapid mixing processes.

The second regime, situated just to the left of the first, corresponds to 'diffusion-altered homogeneous ignition'. Both Karlovitz numbers remain greater than one but approach closely. In this scenario, the most reactive mixture is partially diluted by the least reactive mixture before ignition takes place. The mixing processes are rapid enough to achieve a well-mixed mixture before ignition, but the ignition time is delayed compared to the equivalent fully homogeneous mixture.

On the opposite side of the diagram, when the diffusion time exceeds the two auto-ignition times, the gas ignites successively from the most to the least reactive mixture without being influenced by the mixing mechanism. Both Karlovitz numbers are significantly below unity. This corresponds to the 'ignition front propagation' regime.

If the Karlovitz number of the least reactive mixture alone reaches the order of unity, the mixing processes may influence the ignition front during its propagation toward the least reactive mixture. This regime is referred to as 'diffusion-altered ignition front propagation'.

Eventually, if the Karlovitz number of the least reactive mixture becomes significantly larger than unity, while the Karlovitz number of the most reactive mixture remains considerably lower than unity, ignition takes place in the most reactive mixture. Subsequently, a laminar flame propagates toward the least reactive mixture. This corresponds to the 'heterogeneous ignition and flame propagation' regime.

All these regimes have been identified within the canonical numerical experiments proposed in this study. Each red dot on the diagram in Fig. 20 represents a simulation case. The most remarkable physical mechanism identified is particle rejuvenation, which can be more easily explained through a molecular description of the gas as follows. Consider two distinct fluid particles at disparate temperatures, comprising molecules each equipped with a chronometer measuring their respective age. The chronometers are activated simultaneously. However, the normalisation employed results in a age that is relative to the local ignition delay. Consequently, the molecules in the higher temperature state age at a faster rate, see Figure 21. Subsequently, the two fluid particles are mixed after a designated time interval, resulting in the formation of a new fluid particle composed of molecules with varying ages. However, by analogy with the temperature concept at thermodynamic equilibrium, it is considered that a given fluid element must have a unique age. It can thus be concluded that either the fluid age is the average of the molecular ages, or that the molecular ages may be transmitted from one molecule to another in a manner analogous to the diffusion of temperature. The consequence is that the age of the mixed fluid particle is less than that of the hottest fluid particle prior to mixing and greater than that which would have been expected for the cold fluid particle. As a result, the older fluid particle has become younger and the younger has aged prematurely by mixing with each other. The Figure 21 provides a visual representation of the aforementioned process. In conclusion, the rejuvenation mechanism can be considered as a direct consequence of the application

<figure><img src="image_22.png" /><figcaption>FIG. 21. Rejuvenation and accelerated aging in auto-ignition processes: a consequence of fluid age diffusion.</figcaption></figure>
of the conservation principle to the age of a fluid particle. From a molecular perspective, this implies that two molecules of disparate ages, situated in close proximity to one another, will undergo a process of age equilibration. Although this concept may appear to be contrary to intuitive reasoning, the results presented in this article demonstrate its efficacy in modelling auto-ignition processes in the presence of thermal diffusion.

However, the NRT concept is based on a number of simplifying assumptions, including the assumption of a unity Lewis number for all species and a low Mach number. These assumptions are a direct consequence of the utilisation of a unique mixture fraction to describe the mixing between two streams. It is therefore to be expected that the results will be biased in configurations where these assumptions are no longer valid. Moreover, the NRT exhibits similarities with a chemical reaction progress variable. Therefore, the application of a single residence time to describe the evolution of all chemical reactions is insufficient to account for specific physical mechanisms where several characteristic times are involved, such as for cold flames or for the evaluation of pollutant production. Nevertheless, as with any modelling work, improvements or corrections to the model will enable it to be used more widely, but will necessitate an increase in its complexity. In particular, the use of additional variables, either to describe the mixing in a more refined way, as proposed by Gomet et al.³¹, who use several mixture fractions, or by using additional progress variables to describe different chemical characteristic times, will facilitate the model’s application in a broader range of scenarios. Eventually, it has been demonstrated that two scalar dissipation terms in the NRT equation have been overlooked despite their potential to play a significant role in auto-ignition in regions where scalar mixing is particularly strong.

## VII. CONCLUSIONS

This paper presents a comprehensive investigation that integrates theoretical advancements with original canonical numerical experiments to explore the auto-ignition in temperature heterogeneous mixtures. The primary breakthroughs include a detailed analysis of the normalized residence time concept and the validation of its modeled transport equation. Additionally, the paper identifies auto-ignition regimes within heterogeneous mixtures, categorizing them within an auto-ignition/diffusion interaction diagram.

Two specific numerical experiments, each with unique initial conditions, were meticulously designed to explore these various regimes. The first involves a step function representative of a hot spot, while the second incorporates a temperature profile providing a constant propagation speed of ignition fronts. The initial thermodynamic properties of the gas were selected with consideration for Karlovitz numbers and the crossover temperature.

Analytical developments, which involve representing the residence time as a temporal level-set to track an auto-ignition front, have provided justification for the inclusion of a diffusive term. Additionally, attention has been drawn to supplementary source terms, related to scalar dissipation rates, along with the slope and curvature of the auto-ignition delay in passive scalar space.

The results have shown the relevance of the modeled equation for residence time in the 'homogeneous ignition,' 'diffusion-altered homogeneous ignition,' and 'ignition front propagation' regimes. Additionally, the findings highlight the capability of the residence time field to map transitions between auto-ignition and laminar flame in the 'heterogeneous auto-ignition' regime.

Finally, a unique mechanism, described as particle rejuvenation during the inflammation phase under the influence of mixing with cold gas, has been identified and accounted for through the diffusion term in the normalized residence time equation.

In summary, this article demonstrates that the concept of residence time, when associated with a precise definition and suitable equation, holds promise for studying auto-ignition in a heterogeneous mixture.

## VIII. ACKNOWLEDGEMENTS, SUPPLEMENTARY MATERIAL, AND REFERENCES

This work was performed using HPC resources from GENCI-[IDRIS, CCRT, CINES] (Grant 2022-A013B07456). This work was supported by the CPER FEDER Project of Région Nouvelle Aquitaine. This work relates to the French Government program "Investissements d'Avenir" (EUR INTREE, reference ANR-18-EURE-0010)

<footnote n="1">¹P. Chiesa, G. Lozza, L. Mazzocchi, Using hydrogen as gas turbine fuel, J. Eng. Gas Turbines Power 127 (2005) 73-80.</footnote>

<footnote n="2">²A. Gruber, M. R. Bothien, A. Ciani, K. Aditya, J. H. Chen, F. A. Williams, Direct numerical simulation of hydrogen combustion at auto-ignitive conditions: Ignition, stability and turbulent reaction-front velocity, Combust. Flame 229 (2021).</footnote>

<footnote n="3">³Z. Zhao, Z. Chen, S. Chen, Correlations for the ignition delay times of hydrogen/air mixtures, Chin. Sci. Bull. 56 (2011) 215-221.</footnote>

<footnote n="4">⁴X. He, B. Shu, D. Nascimento, K. Moshammer, M. Costa, R. X. Fernandes, Auto-ignition kinetics of ammonia and ammonia/hydrogen mixtures at intermediate temperatures and high pressures, Combust. Flame 206 (2019) 189-200.</footnote>

<footnote n="5">⁵W. Yang, K. K. J. Ranga Dinesh, K. H. Luo, D. Thevenin, Direct numerical simulations of auto-igniting mixing layers in ammonia and ammonia hydrogen combustion under engine-relevant conditions, Int. J. Hydrogen Energy 47 (2022) 38055-36074.</footnote>

<footnote n="6">⁶W. Han, Z. Sun, A. Scholtissek, C. Hasse, Machine learning of ignition delay times under dual-fuel engine conditions, Fuel 288 (2021).</footnote>

<footnote n="7">⁷Y. Huang, C. Jiang, K. Wan, Z. Gao, L. Vervisch, P. Domingo, Y. He, Z. Wang, C.-H. Lee, Q. Cai, et al., Prediction of ignition delay times of jet a-1/hydrogen fuel mixture using machine learning, Aerosp. Sci. Technol. 127 (2022).</footnote>

<footnote n="8">⁸Y. Ju, C. B. Reuter, O. R. Yehia, T. I. Farouk, S. H. Won, Dynamics of cool flames, Prog. Energy Combust. Sci. 75 (2019).</footnote>
⁹P. Zhao, C. K. Law, The role of global and detailed kinetics in the first-stage ignition delay in ntc-affected phenomena, Combust. Flame 160 (2013) 2352-2358.

¹⁰S. S. Goldsborough, A chemical kinetically based ignition delay correlation for iso-octane covering a wide range of conditions including the ntc region, Combust. Flame 156 (2009) 1248-1262.

¹¹F. Ghirelli, B. Leckner, Transport equation for the local residence time of a fluid, Chem. Eng. Sci. 59 (2004) 513-523.

¹²D.-h. Shin, R. D. Sandberg, E. S. Richardson, Self-similarity of fluid residence time statistics in a turbulent round jet, J. Fluid Mech. 823 (2017) 1-25.

¹³B. Savard, E. R. Hawkes, K. Aditya, H. Wang, J. H. Chen, Regimes of premixed turbulent spontaneous ignition and deflagration under gas-turbine reheat combustion conditions, Combust. Flame 208 (2019) 402-419.

¹⁴X. Wang, V. Robin, A. Mura, A normalised residence time transport equation for the numerical simulation of combustion with high-temperature air, Combust. Theor. Model. 23 (2019) 821-853.

¹⁵F. A. Williams, Turbulent combustion, in: J. D. Buckmaster (Ed.), The Mathematics of Combustion, SIAM Philadelphia, 1985, pp. 115-147.

¹⁶Y. Zeldovich, Regime classification of an exothermic reaction with nonuniform initial conditions, Combust. Flame 39 (1980) 211-214.

¹⁷T. Mantel, R. Borghi, A new model of premixed wrinkled flame propagation based on a scalar dissipation equation, Combust. Flame 96 (1994) 443-457.

¹⁸H. Pitsch, H. Steiner, Scalar mixing and dissipation rate in large-eddy simulations of non-premixed turbulent combustion, Proc. Combust. Inst. 28 (2000) 41-49.

¹⁹K. Bray, P. Domingo, L. Vervisch, Role of the progress variable in models for partially premixed turbulent combustion, Combust. Flame 141 (2005) 431-437.

²⁰A. Mura, V. Robin, M. Champion, Modeling of scalar dissipation in partially premixed turbulent flames, Combust. Flame 149 (2007) 217-224.

²¹S. Yao, B. Wang, A. Kronenburg, O. Stein, Conditional scalar dissipation rate modeling for turbulent spray flames using artificial neural networks, Proc. Combust. Inst. 38 (2021) 3371-3378.

²²D. G. Goodwin, H. K. Moffat, I. Schoegl, R. L. Speth, B. W. Weber, Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, https://www.cantera.org, 2022. Version 2.6.0.

²³N. Swaminathan, N. Chakraborty, Scalar fluctuation and its dissipation in turbulent reacting flows, Physics of Fluids 33 (2021).

²⁴R. Bilger, Turbulent flows with nonpremixed reactants, Turbulent reacting flows (2005) 65-113.

²⁵C. Treviño, F. Solorio, Asymptotic analysis of the high-temperature ignition of CO/H₂/O₂ mixtures, Combust. Flame 86 (1991) 285-295.

²⁶E. Fernandez-Tarrazo, A. L. Sanchez, F. A. Williams, Hydrogen-air mixing-layer ignition at temperatures below crossover, Combust. Flame 160 (2013) 1981-1989.

²⁷M. Ó Conaire, H. J. Curran, J. M. Simmie, W. J. Pitz, William, C. K. Westbrook, A comprehensive modeling study of hydrogen oxidation, Int. J. Chem. Kinet. 36 (2004) 603-622.

²⁸M. Slack, Rate coefficient for h + O₂ + m = hO₂ + m evaluated from shock tub measurements of induction times, Combust. Flame 28 (1977) 241-249.

²⁹D.-h. Shin, R. D. Sandberg, E. S. Richardson, Self-similarity of fluid residence time statistics in a turbulent round jet, J. Fluid Mech. 823 (2017) 1-25.

³⁰D.-h. Shin, A. J. Aspden, V. Aparece-Scutariu, E. S. Richardson, Unsteady self-similarity of jet fluid age and mass fraction, Phys. Fluids 35 (2023).

³¹L. Gomet, V. Robin, A. Mura, Influence of residence and scalar mixing time scales in non-premixed combustion in supersonic turbulent flows, Combust. Sci. Technol. 184 (2012) 1471-1501.

³²X. Wang, V. Robin, A. Mura, A normalised residence time transport equation for the numerical simulation of combustion with high-temperature air, Combust. Theor. Model. 23 (2019) 821-853.

³³D.-H. Shin, E. S. Richardson, V. Aparece-Scutariu, Y. Minamoto, J. H. Chen, Fluid age-based analysis of a lifted turbulent dme jet flame dns, Proc. Combust. Inst. 37 (2019) 2215-2222.

³⁴J. Réveillon, C. Pera, Z. Bouali, Examples of the potential of dns for the understanding of reactive multiphase flows, Int. J. Spray Combust. Dyn. 3 (2011) 63-92.

³⁵S. Zhao, Z. Bouali, A. Mura, Computational investigation of weakly turbulent flame kernel growths in iso-octane droplet clouds in cvc conditions, Flow Turbul. Combust. (2019) –.

³⁶Z. Bouali, J. Reveillon, C. Pera, Development of an extended reactor configuration to analyze preferential segregation impact on spray autoignition, Fuel 302 (2021) 120869.

³⁷A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, C. S. Woodward, Sundials: Suite of nonlinear and differential/algebraic equation solvers, ACM Trans. Math. Software 31 (2005) 363-396.

³⁸T. Varga, T. Nagy, C. Olm, I. Zsély, R. Palvölgyi, E. Valko, G. Vincze, M. Cserhati, H. Curran, T. Turanyi, Optimization of a hydrogen combustion mechanism using both direct and indirect measurements, Proc. Combust. Inst. 35 (2015) 589-596.

³⁹M. H. Saeid, J. Khadem, S. Emami, Numerical investigation of the mechanism behind the deflagration to detonation transition in homogeneous and inhomogeneous mixtures of h₂-air in an obstructed channel, Int. J. Hydrogen Energy 46 (2021) 21657-21671.

⁴⁰G. Smith, D. Golden, M. Frenklach, N. Moriarty, B. Eiteneer, M. Goldenberg, C. Bowman, R. Hanson, S. Song, W. Gardiner, J. Lissianski, Z. Qin, Gri-mech website, http://www.me.berkeley.edu/gri_mech/.

⁴¹S. Taileb, G. Farag, V. Robin, A. Chinnayya, A canonical numerical experiment to study detonation initiation from colliding subsonic auto-ignition waves, Phys. Fluids 35 (2023) 76-101.

## Appendix A: Mesh convergence study

A mesh sensitivity study was conducted with 100, 500, and 1000 cells ($N_x$) for a one-dimensional case within a temperature range of 950 K to 970 K. This corresponds to a regime of auto-ignition front propagation altered by diffusion. The corresponding extreme values of auto-ignition delays were found to be $\tau_{950} = 5.2 \times 10^{-4}$ s and $\tau_{970} = 3.2 \times 10^{-4}$ s and the characteristic diffusion time was $t_D = 4.3 \times 10^{-3}$s.

The Figure 22 illustrates the temporal evolution of temperature and $\dot{H}$ at the midpoint ($L_x/2$) for these three meshes ($N_x$ = 100, $N_x$ = 500, and $N_x$ = 1000). The auto-ignition delays, temperature and age evolution obtained with the three meshes are largely comparable, with a minor discrepancy observed for the coarsest mesh. In contrast, the outcomes for the two most refined meshes are nearly indistinguishable, which suggests that the spatial discretisation employed in this paper is adequate.

Regarding the time step, Asphodel employs the CVODE solver from the SUNDIALS suite, which utilizes an adaptive time-stepping method. Among the various CFL conditions, the one defined by the chemical reaction rates is the most limiting and is the one that controls the time step. The minimum and maximum values of the time step obtained for the three meshes are summarized in Table II. This result demonstrates that, irrespective of the meshes under consideration, the numerical time step is several orders of magnitude smaller than both the chemical and mixing time scales. This ensures the precise resolution of the rapid dynamics associated with auto-ignition while maintaining a reasonable overall computational time.
<figure><img src="image_23.png" /><figcaption>FIG. 22. Time evolution of temperature and NRT $\hat{H}$ at the midpoint of the domain ($L_x/2$) for three mesh resolutions, $Nx = \{100, 500, 1000\}$, in the auto-ignition front propagation regime within the temperature range of 950 K to 970 K.</figcaption></figure>

<table><caption>TABLE II. Minimum and maximum time step sizes for different spatial discretizations.</caption><thead><tr><th>Nx</th><th>Δt<sub>min</sub> [s]</th><th>Δt<sub>max</sub> [s]</th></tr></thead><tbody><tr><td>100</td><td>5.14 × 10<sup>-8</sup></td><td>1.58 × 10<sup>-6</sup></td></tr><tr><td>500</td><td>1.40 × 10<sup>-8</sup></td><td>4.06 × 10<sup>-7</sup></td></tr><tr><td>1000</td><td>3.30 × 10<sup>-8</sup></td><td>1.01 × 10<sup>-7</sup></td></tr></tbody></table>

## Appendix B: Validation of the NRT model for methane-air auto-ignition

Figure 23 illustrates the results obtained with a stoichiometric methane-air mixture at atmospheric pressure and at $T_u = 1500$ K using the GR13.0⁴⁰ chemical scheme. The results obtained using the $\hat{H}$ equation are excellent for both fuels, hydrogen and methane, thereby demonstrating the equation's capacity to handle a broad spectrum of fuels and mixtures for the prediction of auto-ignition delay times.

<figure><img src="image_24.png" /><figcaption>FIG. 23. Time evolution of temperature (red) and NRT $\hat{H}$ (black) for a homogeneous reactor initial state of unburned stoichiometric methane-air mixture at $T_u = 1500$ K.</figcaption></figure>