Embedded direct numerical simulation of ignition kernel evolution and flame initiation in dual-fuel spray assisted combustion

The flame initiation process in dual-fuel spray assisted combustion is presently not fully understood. Here, diesel spray assisted combustion of premixed methane/oxidizer/EGR is explored in the post-ignition phase by scale-resolved simulations. The modified dual-fuel ECN Spray A forms the baseline configuration. An extensive local grid refinement (approaching DNS limit) around one of the first high-temperature ignition kernels is carried out in order to examine the validity of hypothesized flame initiation and deflagration. A high quality LES is used to solve the spray dynamics, while the embedded quasi-DNS (eq-DNS) region offers detailed information on the ignition kernel evolution. The finite-rate chemistry is directly integrated, utilizing 54 species and 269 reactions. Local combustion modes are investigated for the ignition kernel development toward spontaneous ignition and premixed flame propagation using various approaches, including the reaction front displacement speed, energy transport budget, and chemical explosive mode analysis. Furthermore, a new criterion based on reaction flux analysis is introduced, which is compatible with dual-fuel combustion. The spatial and temporal scales associated with the ambient methane consumption and consequent flame initiation are characterized. For the first time in dual-fuel spray assisted simulations, numerical evidence is provided on the initiation of premixed flames, and the corresponding timescale is reported. Particularly, there is a transient mixed-mode combustion phase of approximately 0.2 ms after the spray second stage ignition wherein extinction, ignition fronts, and quasi-deflagrative structures co-exist. After such a transient period, the combustion mode becomes essentially deflagrative. Finally, interactions between turbulence and premixed flame front are characterized mostly in the corrugated regime.


Introduction
Spray assisted dual-fuel (DF) combustion is considered a technological upgrade to the conventional spray diffusion flames in compression ignition engine applications.Such a technology gives an opportunity to substitute a large part of long-chain hydrocarbon fuels, such as diesel, with other fuels of lower C/H ratio, such as methane [1].In order to develop this technological approach, efforts have been increasingly dedicated to investigate the associated global and local combustion characteristics.With the help of numerical modeling and high performance computing, it is nowadays possible to investigate highly transient, highly localized combustion phenomena involving spray dynamics, spray assisted ignition, auto-ignition kernel development and flame propagation, with a high level of spatial and temporal detail.
Spray assisted DF combustion is characterized by the injection of a high reactivity fuel (HRF) spray into a hot ambient gas environment.
The inertial and aerodynamic shear forces control the liquid atomization and droplet breakup, which are typically fast processes, until spray evaporation occurs [2].Due to the multi-staged ignition process, the onset timescale of chemical reactions -denoted by τ 1 -is characterized by HRF dissociation into, e.g., alkyl radicals involving low temperature chemistry [3].After that, such radicals are further decomposed during a chemical induction phase into lighter intermediate species until thermal runaway (i.e.high-temperature ignition) is initiated [3,4], which is herein denoted by τ 2 .The resulting auto-ignition pockets are formed at regions of low scalar dissipation rate [5], also related to the most reactive mixture fraction.These sites act as an ignition source for the ambient low reactivity fuel (LRF) -air mixture.
In practical DF engines, the LRF provides 85 -98% share of the total combustion energy, whereas HRF spray is responsible for the https://doi.org/10.1016/j.combustflame.2023  remaining energy share.The pilot spray ignited pockets are chemically driven, with reaction time scales being much smaller than those for transport processes.On the other hand, turbulent structures and transport phenomena may influence the dynamics of ignition kernel development, as discussed in [6][7][8][9].The auto-ignition kernels grow, and possibly merge, while propagating downstream the spray.Excluding abnormal scenarios [10], the subsonic reaction front may propagate in auto-ignition (spontaneous), deflagration, or auto-ignitive/deflagrative mixed mode.As described by Zeldovich [11], a spontaneous front propagation is driven by the ignition delay time (IDT) gradient of the neighboring mixtures, with local exothermic reactions (i.e.zerodimensional phenomena) being responsible for the thermal runaway without much contribution from transport effects.In the deflagrating mode, however, thermochemical diffusion processes become as important as chemical reactions.Hence, the propagation can be considered to be multidimensional and balanced by the reaction-transport processes.
An illustration of the DF concept under the Engine Combustion Network (ECN) Spray A conditions is depicted in Fig. 1.Typically, LRF (methane) is port-injected early during the intake stroke, (I).Then, by the end of the compression stroke, pilot HRF (n-dodecane) spray is injected to assist the ignition and combustion of LRF (II), in which the latter is mainly consumed during the combustion/power stroke (III).By examining the spray vaporization and multi-staged ignition dynamics (red inset) corresponding to one of the HRF injector holes, the attention is herein shifted toward the evolution of a single auto-ignition kernel of n-dodecane spray (light green inset).There, a performed local 'numerical microscopy' can provide insight on the flow topology and the spatio-temporal scales affecting the kernel development.More importantly, while kernels initiate on the rich side of stoichiometric mixture fraction ( st ) [4], it is scientifically valuable to understand the modes of combustion for the kernel propagation and LRF consumption across the spray.
Various efforts have been dedicated to understand the ignition characteristics of DF spray assisted configurations using detailed threedimensional models, with little focus being devoted to the fine details in the post ignition events.Under ECN baseline conditions, Kahila et al. [4,12] conducted large-eddy simulation (LES) with finite-rate chemistry approach to investigate the ignition characteristics under methane/diesel DF mode, which was further supported via experimental findings by Srna et al. [13,14].Further LES studies have been conducted on various parameters including ambient temperature [15][16][17], chemical kinetics [18], and tri-reactivity blends [19,20], in addition to other relevant works concerning model developments [21,22].Despite recent DF studies [23,24] that also investigated reaction front structures beyond high-temperature ignition, the employed computational grid was not sufficient to capture the flame front thickness.
Experimentally, the mixed-mode combustion of flame propagation and auto-ignition has been recently discussed for various applications including staged [25], partially premixed [26], diesel/gasoline DF [27], direct DF stratification [28], and spark-initiated pilot-stabilized [29] combustion systems.In general, experimental evidence has been reported on the co-existence of deflagrative and auto-ignitive fronts.
Considering the identification of the combustion modes between spontaneous ignition and deflagration, several approaches have been adopted.Chen et al. [30] carried out two-dimensional direct numerical simulations (DNS) of thermally stratified hydrogen/air mixtures.They investigated combustion modes via displacement speed analysis of reaction wave front compared to a reference deflagrating speed.Under similar conditions, Bansal and Im [31] studied thermal and fuel stratification levels on the front propagation.They proposed a criterion based on hydroperoxyl (HO 2 ) Damköhler number to demarcate the relative importance of chemistry and diffusion processes on the heat release.In single-fuel (SF) stratified mixtures, computational singular perturbation (CSP) analysis has been applied to classify ignition regimes while considering thermal inhomogeneities [32].Further studies discussing SF combustion modes under thermal and/or compositional stratification are recalled in [33] along with the references therein.More recently, the front structure and propagation speed under auto-ignitive conditions have been fundamentally investigated using one-dimensional models [34][35][36], wherein the combustion mode is characterized according to competitive effects from residence time -ahead of the front -and the ignition delay.From the chemical perspective, Schulz et al. [37] proposed a criterion, namely the auto-ignition index, based on reaction flux analyses in [38,39] using HO 2 chemistry.Additional essential methods developed for combustion mode identification include the transport budget analysis [40], and chemical explosive mode analysis (CEMA) [41].Other tools for premixed versus non-premixed flame modes have been also noted, including the flame index [42,43], and data-driven supervised classifiers [44].
Considering LRF flame initiation in DF systems, previous studies have reported the existence of propagating flame fronts.However, the focus therein was shifted toward model validations and injection parameters effects such as the pilot HRF quantity [9,12], compression heating and injection timing [28], and the propagation modeling [45].Recently, Kannan et al. [46] investigated turbulence intensity effects on the flame formation in a shear-driven scale-resolved configuration.With relevance to constant-volume DF application, Karimkashi et al. [33,47] extended the works by Sankaran et al. [48] to develop a predictive a piori tool, namely -curve, for identifying combustion modes of stratified DF mixtures under transient conditions, while incorporating convective mixing effects.They utilized such an analytical approach for different fuel combinations as well as turbulence and fuel stratification, under DF engine relevant conditions.As a more sophisticated approach, possibly applicable for DF systems, CEMA has been widely acknowledged in the combustion community.Based on theories from dynamical systems and CSP [49], CEMA was developed as a posteriori diagnostic tool for ignition and flame structure identification [41], and it was further extended to account for diffusion [50] and recently the evaporation [51] effects.
At present, there has not been much attention on the propagation mechanisms for reacting fronts under spray assisted DF configurations.In particular, there is a lack of detailed analysis on the transition between auto-ignitive and deflagrative combustion modes.Further, most DNS studies on combustion mode identification are limited to simplified configurations, which might result in some irrelevancy with regard to real engine conditions.Such challenges are herein leveraged, to a certain extent, using the embedded quasi-DNS (eq-DNS) approach.The main idea, also noted in [52,53], is to perform local numerical microscopy within an embedded fully resolved zone -regarded as the investigation domain -while utilizing LES grid in the background.Here, LES framework is adopted to resolve the spray dynamics, under engine relevant conditions, while providing initial turbulence to the investigation domain.
In this study, we extend our previous numerical investigations on diesel spray assisted DF ignition of methane under ECN Spray A operating conditions [4] to beyond the ignition event.The main motivation is to (i) identify local combustion modes, which can be crucial for engine performance and combustion stability, and (ii) investigate the microscopic details affecting the ignition kernel development with possible flame initiation.Accordingly, starting from the secondstage ignition pocket initiated at time τ 2 , an extensive refinement of the computational grid around the most reactive mixture fraction within the investigation domain is applied.This refinement region is constructed to resolve the hypothesized deflagration front thickness and thereby shed insight toward the smallest flow scales (i.e.DNS resolution) affecting the kernel development in the post spray ignition phase.Additionally, the pilot n-dodecane spray in the background domain is resolved via LES grid (62.5 μm) to provide natural inflow/outflow boundary conditions to the refinement zone of interest.The objectives of the present study are to 1. utilize embedded quasi DNS approach in the context of ndodecane spray assisted ignition of methane.2. extend the current understanding of DF combustion mode development after spray ignition toward hypothesized deflagration.3. analyze the mode of combustion using various novel metrics.

Methodology
In the following, the governing equations and the adopted submodels are presented.The numerical techniques for combustion mode identification, namely the displacement speed, CEMA, and reaction flux analysis, are further described in Appendix A.

Governing equations and numerical framework
Considering turbulent spray combustion modeling, the gaseous phase is described by the compressible Navier-Stokes equations.The conservation of mass, momentum, species mass fractions, and enthalpy are mathematically formulated in vector notation as in the following, where ρ, ũ, p, , Ỹ , , h , h denote the implicitly filtered density, velocity, pressure, viscous stress tensor ) , th species mass fraction, mass diffusivity, sensible and total (nonchemical) enthalpy, respectively.The total enthalpy is herein defined as the summation of h and specific kinetic energy, whereas sensible enthalpy is based on the ideal gas caloric state equation, i.e. h = ∑ sp =1 Ỹ , with  0 denoting reference temperature at 298.15 K, and the isobaric specific heat capacity  , ( ) is retrieved using tabulated temperature-dependent polynomial fits.A compressible version of the Pressure Implicit with Splitting of Operators (PISO) algorithm [54] is utilized, together with ideal gas thermodynamic state equation for the pressure-velocity-density coupling.
The overbar and tilde operators in the governing equations denote unweighted and density-weighted (Favre) ensemble averaging, respectively, and the (⊗) symbol refers to the outer product.The filtered source terms   ,   ,    , and  ℎ account for the liquid spray heat and mass transfer, including phase change, in the respective LES formulations.The chemical source terms ω and ωℎ =   sp  ℎ 0  , ω denote the species overall production/consumption rate and the heat release rate (HRR), respectively, with ℎ 0  , being the th species enthalpy of formation.A Strang-type operator-splitting approach is employed to decouple combustion chemistry from fluid dynamics.The change in the local thermochemical composition vector () is described by a stiff, nonlinear initial value problem as follows, whereas the nonlinear function  ( () ) is computed using finite rate chemistry for the mass-based overall production rates of species and temperature.Temporal integration of the previous equation is performed using a stiff linearly implicit Euler extrapolation method, namely SEULEX algorithm [55].Further details on the methodology and conceptual derivations can be found in the doctoral thesis by M. Gadalla [56].
The governing equations are solved concurrently in both the LES and eq-DNS grid regions, and they are numerically discretized using the finite volume method within the OpenFOAM framework [57].Second order accurate schemes are generally utilized, with central differencing for diffusive fluxes and implicit backward differencing for time marching.Convective fluxes are treated using nonlinear limiting schemes, which are further discussed in the following section.The fluid dynamic equations are solved sequentially, and numerical stability is enforced and controlled for the global time stepping via convective and diffusive CFL criteria held below 0.3.
To enhance the computational performance, efficient strategies are adopted to accelerate the finite-rate chemistry source terms evaluation procedure [58,59].First, dynamic load balancing is applied to evenly re-distribute the computational load between parallel processes, as discussed in [58].Second, an analytical formulation of the thermochemical Jacobian matrix is employed using pyJac [60], which features minimized computational and memory operations.Third, robust linear algebra routines are utilized for matrix factorization and back substitution, as discussed in [12,59].Finally, linear scaling is achieved as noted in [56].

Adopted submodels
The rationale behind the adopted modeling choices is the need for simplicity as well as reducing the parameter dependence and empirical calibration, hence avoiding further uncertainties within the sub-models.This can be achieved via utilizing (i) fine grid resolution to minimize sub-grid scale (SGS) effects, while adopting quasi laminar approaches [61,62] for the closure modeling, and (ii) implicit models that are capable of capturing SGS physics through numerical integration [63].

Spray modeling
The source terms   ,   ,    , and  ℎ , relevant to spray coupling with the gaseous phase, are evaluated using the Lagrangian particle tracking of liquid droplets.The no-breakup model approach is adopted herein, in which the droplets are introduced with a constant size corresponding to a stable Weber number (  < 12).The liquid volume fraction is herein estimated as ≈ 0.1-0.4% for the spray plume, while the maximum local value is less than 1% when sampled within a cubic control volume of side length 1 mm.Therefore, the dilute suspension assumption may be permitted for the droplet-laden flow, wherein the droplet-droplet interaction modeling becomes less important [64].Further, a cylindrical injection model is used, which extends the standard disc injection model of OpenFOAM into 3D.The implementation details of such spray sub-models are further discussed in the earlier study [2] by the authors, wherein global spray metrics under ECN conditions are thoroughly validated.Nevertheless, the present study focuses on a small spatio-temporal window of investigation wherein the vast majority of spray droplets are vaporized.

Turbulence modeling
The implicit LES (ILES) approach [63] is adopted for turbulence closure modeling of the residual SGS terms on the right side of Eqs. ( 2)-(4).Instead of additional explicit terms, the ILES approach herein imposes local numerical dissipation through a nonlinear, flux-limiting scheme termed as Gamma scheme [65], on the convective terms.This approach has been extensively used in various combustion studies [4,12,16,19,20,23,59,66].Further, in the previous ECN study [2] by the authors, a comparison between implicit and explicit LES models was performed, and it is observed that SGS effects are rather insignificant on global spray metrics when utilizing fine LES grid.Also, the reported kinetic energy spectra therein show rather similar pattern between ILES and the explicit one-equation eddy diffusivity model.Moreover, in the (DNS-like) refinement zone of interest, further illustrated in Section 3, it is expected that the governing equations reduce to the conventional form, i.e. without SGS modeling, although truncation errors due to the utilized numerical schemes are acknowledged.

Combustion modeling
The chemical reaction kinetics are considered through the skeletal mechanism by Yao et al. [67], developed for n-dodecane reactions under engine-relevant conditions.The mechanism involves 54 species and 269 reactions, and it has been previously utilized for DF ECN Spray A ignition of methane [4,12,15,16,18,66].In order to be consistent with the baseline configuration in [4], a unity Lewis number is assumed for all species.Furthermore, for the lean combustion of methane/air mixtures, the influence of non-equidiffusion effects, in terms of nonunity Lewis number and preferential diffusion, on the response of stretched reacting fronts are foreseen to be minor since the effective Lewis number is almost unity.
Turbulence-chemistry interaction (TCI) effects are implicitly considered via a perfectly stirred reactor approach together with a fine LES grid resolution [61,62].Such an approach considers the approximation ω ≈ ω ( p, T , Ỹ ) while neglecting SGS effects for thermal and compositional fluctuations.Various studies have considered this approach for combustion closure modeling [12,16,19,61,62,[68][69][70][71].Additionally, the adequate performance has been reported, with agreement to experimental results, which is mainly attributed to the fine LES grid resolution, as noted in [12,19,68,69,71] and the references therein.More importantly, in the post ignition phase after τ 2 which is the main investigation window of this study, the computational grid is extensively refined in the DNS limit.Such an embedded DNS-like region acts as a numerical microscope, capturing the local fine structures, and resolving the flame front thermal thickness.

Numerical setup
The present numerical study is considered a continuation to the previous works by our research group on diesel spray assisted combustion of methane, based on the ECN Spray A conditions.A summary of the thermophysical conditions and computational settings of the numerical simulation is provided in Table 1, further details of which are discussed afterwards.

Case preparation
The preparation of the numerical setup is achieved in two phases.First, the domain is discretized into a block-structured hexahedral mesh, such that the computational cells within the spray volume hold a uniform characteristic length of 62.5 μm for the innermost refinement layer.Such a grid resolution (i.e.fine LES grid) is consistent with that reported in [4] for similar DF configuration, which is rather sufficient to resolve large inertial length scales.LES is then performed to ignite the mixture, while the most reactive high temperature ignition pocket is spatially and temporally located.Second, a refinement zone is considered in the vicinity of the spray cloud for the main numerical investigation of this study.Particularly, the region surrounding the most reactive ignition kernel is refined from 62.5 μm down to 3.9 μm in the DNS limit.Such a DNS-resolution grid is embedded within the LES grid, and it is regarded as the investigation domain.
The eq-DNS grid is axially located 32 mm away from the spray injector, and it is shifted 5.5 mm and 0.75 mm away from the spray center-line axis in the other two perpendicular directions to capture the spray ignition kernel.Therefore, considering the relatively short liquid length of Spray A, the mixture within the eq-DNS region becomes entirely vaporized.Details on the embedded grid, location, and dimensions are depicted in Fig. 2. The final configuration in the whole domain results into a total of 316 million computational cells.
To initialize the new simulation, geometrical fields corresponding to the thermochemical and flow conditions at slightly before the ignition event (i.e.τ 2 ) are directly mapped to the new mesh, holding the eq-DNS grid.A short curation period (0.1 ms before τ 2 ) is herein allowed for the flow turbulence to develop and thereby to recover the fine scales of the wavenumber spectrum.Such a curation period is compromised between computational feasibility and accurate construction of the flow  physics, also considering the transient character of the spray assisted ignition process.As a remark, the eq-DNS region is initialized with an LES filtered solution [4,12,16] on 62.5 μm grid.Consequently, the eq-DNS region may not pose full history of the turbulence generation cascade starting from the injector.Indeed, there will be emerging fine scale turbulence inside the eq-DNS region, with a constructed spectrum due to the transitional grid refinement as depicted in Fig. 2.However, the boundary conditions for the eq-DNS region do not pose the finest scales, neither do we aim at artificially generating them along its boundaries.Nevertheless, due to the relatively fast turn-over time scales, such structures are expected to die out quickly, hence, having less influence on the eq-DNS results.
The choice of the eq-DNS resolution (3.9 μm) is motivated by an a priori one-dimensional analysis (not shown) of freely-propagating DF premixed flames at various mixture fractions using Cantera [72].The analysis suggests that the flame thermal thickness varies between 32 μm and 82 μm for stoichiometric ( =  st ) and ambient ( = 0) mixtures, respectively, under laminar conditions.Moreover, Pei et al. [73] estimated the Kolmogorov length scale for ECN Spray A to be within 1-5 μm.Therefore, the adopted resolution is expected to capture the reacting front within the DNS limit.
It is worth noting that the present eq-DNS setup is not necessarily identical with exact DNS, but rather a locally highly refined LES region approaching the DNS resolution limit.In essence, the concept of embedded DNS or DNS embedded LES has been previously introduced in the literature [52,74,75] in which various arguments have been made to justify the usage of well-resolved LES to initialize the solution in a DNS grid region.Such a DNS grid resolution is sufficiently high to resolve the thinnest reaction zones and Kolmogorov scale.Furthermore, while the ground truth behind DNS is the ability to resolve turbulence spectra down to the Kolmogorov scale, some works claim to achieve that using lower order accurate numerical schemes [76].In the present work, however, the ''quasi'' DNS terminology is utilized in the embedded region to acknowledge performing flame resolved simulations with (i) possible uncertainties in the full reconstruction of the turbulence spectrum across eq-DNS boundaries, and (ii) the utilization of relatively low order accurate schemes, hence the possible need for numerical dissipation rate calibration [77].

Results
In the following, simulation results are analyzed for a better understanding of the transient thermophysical phenomena associated with ignition kernel development and hypothesized flame deflagration.We highlight that the present eq-DNS is focused on post-ignition investigations.Particularly, while our previous LES studies [4,12,16,18] have been limited to ignition event, i.e.  < τ 2 , here the focus is shifted to a small spatio-temporal window in the post ignition phase, i.e.  ≥ τ 2 .For completeness, a consistency check is provided as supplementary material.

Flow topology at ignition onset
The first analysis concerns reporting the turbulence levels at the ignition onset.While the main investigation is focused inside the eq-DNS box, cf.Fig. 2, the surrounding turbulence is generated naturally by LES under the DF ECN Spray A target conditions at the background 62.5 μm grid.Fig. 3 depicts various analyses on the flow topology at τ 2 timescale.A 2D cut-plane of the instantaneous velocity magnitude in the eq-DNS region, overlaid with temperature contours at 1200 K (black) and 1500 K (white), is shown in Fig. 3(a).
The spray injector is located far upstream, away from the bottom left corner of the figure.The corresponding probability density function (PDF) plot of the 3D eq-DNS data is depicted in Fig. 3(b), while showing a mean value of 10.2 m s −1 .Accordingly, the ignition kernel is expected to be transported with a similar mean flow (convection) speed, resulting in a displacement of 1 mm within 0.1 ms.As a remark, the duration of the eq-DNS simulation is altogether 0.4 ms.Indeed, as time progresses, the reaction front propagation will further decelerate while turbulence levels decay away from the spray injector, also due to combustion, hence supporting the geometrical and temporal choices for the investigation domain.The PDF of the fluctuation velocity  ′ = √ 2∕3, with  denoting turbulent kinetic energy estimate, is depicted in Fig. 3(c).Such turbulent fluctuations exhibit a mean value of 1.47 m s −1 , with a narrow peak concentrated around 0.25 m s −1 which is rather small due to the kernel's spatial location within the spray periphery.
Fig. 3(d) depicts the 2D cut-plane of the resolved local scalar dissipation rate ( = 2‖∇‖ 2 ) field in logarithmic scale.It is well established for turbulent non-premixed combustion that most reactive mixture fractions are formed at spots of low scalar dissipation rate [5].Such spots, possibly located at vortical cores, tend to accumulate radicals during the chemical induction phase.Under pressurized and preheated engine-relevant conditions, thermal runaway emerges locally, and hence, auto-ignition kernels are formed.With relevance to the previous note, Fig. 3(e) shows temperature scatter plots in the  space, overlaid with the conditional mean (red) and standard deviation (yellow), where ignited mixtures seem to be located in regions of low .Additionally, as noted in Fig. 3(f), such most reactive mixtures occur in the richer side of  st , which is a typical scenario for diesel/methane DF ignition [4].Therefore, the observations above indicate that (i) the flow around the ignition kernel is turbulent, and the kernel will be mostly confined inside the refinement zone, (ii) regions close to the most reactive mixture hold rather small  values, which are favorable for ignition, and (iii) most reactive mixture fractions occur on the richer part of stoichiometry which is consistent with spray assisted ignition concept.) ) = 0.5 (yellow).Consistent with Section 4.1, the auto-ignition kernel first appears in the rich side of the spray, before it develops and propagates across the stoichiometric contour to burn SF methane in the premixed ambient.Here, a new timescale can also be introduced which is linked to the instance when the yellow isocontour crosses the white one.The underlying hypothesis is that the time instance of burning the spray-lean mixture is directly related to the (SF) premixed flame onset, while gradients of  below stoichiometry are rather sharp.This is further elaborated and discussed in Section 4.5.

Ignition kernel development
On the other hand, temporal evolution of OH mass fraction (lower half of the figure) demonstrates the species transient characteristics.The overlaid isocontours denote the source term ωOH at values 100 kg m −3 s −1 (green) and 1000 kg m −3 s −1 (red).As it is elaborated afterwards, it is possible to characterize the local combustion mode through OH chemistry.Particularly, while methane flame initiation is associated with a fairly steady ωOH in which OH consumption reactions are influential, during the earlier n-dodecane auto-ignition phase the OH production reactions are substantially dominant.The details of such reaction flux analysis, including the introduction of a new combustion mode identification marker for DF systems, are further discussed in Section 4.8.We note that the analysis above offers qualitative insight to the potential flame initiation, although more in-depth analysis is still required for further conclusions, as elaborated in the following.

Local and statistical analysis
In this section, quantitative analysis is performed to better understand the local and statistical characteristics of the combustion mode.First, a basic energy transport budget analysis is carried out.Particularly, while auto-ignition (spontaneous) front is governed by the reaction terms, deflagration is rather controlled by diffusion and reaction processes.Second, thermal thickness and the structure of key species are statistically examined and compared against 1D reference flames.

Local energy budget
The energy budget terms, namely convection, diffusion, and reaction from Eq. ( 4) are plotted, after being normalized with the maximum absolute value of the reaction and diffusion terms, along local investigation lines across the reaction front, for two representative 2D snapshots as depicted in Fig. 5.The snapshots are selected at times τ 2 and τ 2 + 0.25 ms, as noted in the figure, to represent different combustion modes.The local sampling lines are all centered on the reaction front at  CH 4 = 0.5 (solid curves in the 2D planes, while  CH 4 = 0.3 isocontours are dashed) with a varying length of 0.75 and 0.1 mm for the top and bottom snapshots, respectively.
At τ 2 (top row) the temperature levels are elevated, yet have not reached their maximum limit above 2000 K.More importantly, diffusion effects are noted to be insignificant, and the energy transport is in balance between convection and reaction terms.This observation is aligned with the definition of a typical spray assisted ignition, which is a reaction-driven phenomenon.
As time progresses, the reaction front propagates away from the spray injector and, therefore, turbulence effects turn to be less pronounced.Moreover, diffusion processes become significant, while preheating is attained as temperature and active radicals diffuse away from the reaction zone toward the fresh mixture.Accordingly, reaction and diffusion processes become equally influential, which is characteristic to premixed flames, and the reaction peak is balanced by diffusion and convection.
The observations above are among the first numerical evidence that deflagrative fronts locally emerge after a short combustion phase which is auto-ignition driven.Moreover, it is also noted that the reaction zone denoting radical production and heat release -during the deflagrative phase -becomes narrow and thus comparable to a premixed flame structure, in contrast to the earlier time instances.Therefore, the presented analysis demonstrates locally the transition period from an ignition driven combustion to a deflagrative reaction front.Next, the flame thermal thickness for the last frame in Fig. 5 is investigated.

Flame thickness
To further verify the premixed methane flame initiation as suggested in Fig. 5 at τ 2 + 0.25 ms, here thermal thickness ( L ) in the reaction zone is investigated.By using the definition  L = ( b −  u )∕max(| ∕|), thermal thickness is estimated and compared against 1D flamelets simulated using Cantera, as presented in Fig. 6.In particular, two steady laminar unstrained premixed flames at representative mixing states ( =  st and  = 0) are considered hereafter as the canonical cases, as depicted in the left subplot, to validate the DF flamelets.The remaining thermochemical conditions in the reference cases are consistent with the presented 3D setup, cf.Table 1.While  L values of the representative canonical flames are 82.80 and 31.05μm for the respective  = 0 and  =  st , it is observed that local  L values across the reaction front of the 3D data (right subplot) are mostly within that range.Also, after noting that one investigation line of the 3D data (particularly, line-2 from Fig. 5) holds  L value close to that for canonical stoichiometry (middle subplot), we then aim at estimating thickness distribution in the investigation domain.In that regard, an automated procedure is performed to sample many investigation lines (over 500 samples) of 0.4 mm fixed length each, across the reaction front ( CH 4 = 0.5) in the outer periphery of the spray.After filtering out highly distorted samples, the scatter of local  L values along the investigation lines are plotted against their local average , as depicted in the right subplot.There,  L values vary mostly within the canonical reference values.Moreover, they correlate with , where samples of leaner mixture present higher  L .Indeed, it is expected that transient and straining effects from 3D data prevent thickness values to perfectly match with the canonical 1D flames which are unstrained, planar and at steady-state.Finally, the values of  L further confirm the employed 3.9 μm grid choice, hence the reaction fronts are well resolved.

Statistics of key species
Another analysis to verify flame initiation is reporting the species structure compared with 1D unstrained canonical flame data.Using the same reference Cantera flames at  = 0 and  =  st as previously discussed, scatter plots of key species, namely CH 2 O, HO 2 , CH 3 , and OH are herein presented in Fig. 7.The temporal evolution of the normalized mass fractions, conditioned by  ≤  st , are plotted in the methane progress variable space.The normalization constant is estimated as the maximum value between the 3D data, 1D data at  st and 1D data at  = 0.The columns denote snapshot instances τ 2 , τ 2 + , τ 2 + 3, τ 2 + 5 (left to right) with  = 0.05 ms.While the standard deviation is marked by the gray region, in the latest snapshot the structure of a single flamelet (particularly, line-2 from Fig. 5) is also overlaid in red ( ).From the figure, it is observed that during the early development of the auto-ignition kernel (i.e.τ 2 and τ 2 +0.05 ms), the structure of species data does not span the whole progress variable space, i.e. methane is not yet fully consumed.Moreover, species structure follows a pattern that is closely related to the canonical stoichiometric flame while being restricted to narrow  CH 4 ranges.On the other hand, later snapshots demonstrate fairly similar species structure as observed in the canonical flames, while spanning the entire  CH 4 space.Hence, the thermochemical states observed in 3D and 1D flamelets can be noted to be rather similar.At this point, it should be noted that while local zones could exhibit deflagrative behavior, as in the plotted flamelet ( ) previously demonstrated via transport budget, the present statistics might not be sufficient to conclude an overall deflagrative mode of the initiated fronts.Instead, details of local combustion mode analysis within the fronts are further discussed in the following sections.

Displacement speed analysis
In the previous analysis, flame initiation was investigated in terms of flame thermal thickness, statistics of key species, and energy transport budget analysis by comparing 3D data against representative simple 1D flame structures.Another criterion to distinguish the deflagrative mode is the reaction front propagation speed.In fact, previous literature has shown that the ratio of spontaneous ignition front speed to laminar deflagration front speed can vary by orders of magnitude [30,31].
The front displacement speed ( *  ), weighted by density to ensure dilatational invariance across the flame thickness [78], is herein evaluated and then compared with the laminar flame speed for a stoichiometric DF mixture,  L ( st ) = 0.79 m s −1 .Fig. 8 depicts the temporal evolution of mean (solid curve) and standard deviation (filled region) of  *  along the reaction front, divided by  L ( st ).As noted in the figure, at τ 2 the front speed is much higher than laminar deflagration speed.After a short period, the front speed rapidly conforms to a deflagration behavior, with the ratio  *  ∕ L approaching unity on average.Such an observation provides strong evidence on the emerging flame fronts shortly after spray ignition.Moreover, from the viewpoint of DF compression ignition engines, it is rather essential to understand the mode of combustion, and possible operational consequences, at different crank angles.Next, an analysis to identify the onset for initiation of ambient methane combustion is presented.

SF combustion onset
In this section, an attempt is made to identify the temporal onsetherein introduced as τ 3 -for SF methane combustion in the premixed ambient.Such a timescale can be directly connected to the premixed flame initiation analysis, which is the main topic of the present work.The performed analyses in the previous sections already suggest that local deflagrative fronts may emerge after about 0.2 -0.3 ms after ignition.
To estimate the value of τ 3 , two approaches are herein considered, noting that  gradients below stoichiometry are rather sharp.In the first approach, τ 3 is defined as the timescale at which the front's mean mixture fraction falls below stoichiometry.The PDFs of ( CH 4 = 0.5) at various time instances are depicted in Fig. 9 (left).The corresponding mean and standard deviation values during the entire temporal window of investigation are realized in Fig. 9 (right) with blue and cyan colors, respectively.Using this approach, τ 3 is estimated as τ 2 + 0.22 ms when the front crosses  st , on average.Also, regarding the observed trend of the PDF plots, the physical interpretation is rather straightforward: the most reactive  (at τ 2 ) develops to a front that gets leaner while expanding and spanning a wider distribution of  as time progresses.Again, the wide distribution of  further highlights the locality of the mixing, hence the thermochemical state.Therefore, initiated fronts may pose a mixed-mode combustion in which local zones only become deflagrative according to their local thermochemistry, which is discussed in the following section.
In the second approach, discretized domain statistics are used to define τ 3 .First, the data are partitioned according to their burning state.Particularly, each computational cell is considered either burned or unburned according to its local  CH 4 , whether above (and equal) or below 0.5, respectively.Then, considering the unburned data inside the spray (called set-I), they are defined by those cells with ( CH 4 < 0.5) ∩ ( ≥  st ).Similarly, the burned data outside the spray stoichiometry (i.e.spray-lean or ambient, called set-II) are conditioned by ( CH 4 ≥ 0.5) ∩ ( <  st ).As time progresses after τ 2 , it is expected that the cell count of set-I will decay, whereas set-II (being empty at τ 2 ) will grow.Such a note is precisely depicted by the red curves in Fig. 9 (right).Accordingly, an estimate for τ 3 here could be the intersection of both curves, where a global minimum is realized with regard to both the unburned zone inside the spray and the burned zone outside the spray.Using this approach, τ 3 is estimated as τ 2 + 0.19 ms, which is rather close to the resulting value from the previous approach.Here, τ 3 is also considered as the onset at which the reaction front crosses stoichiometric mixture fraction, cf.Fig. 4. Further, it is observed that a steady minimum in set-I does not entirely diminish at later instances, which implies that some zones inside the spray have not reached a full ignition state, possibly due to the entrainment of freshly vaporized spray mixtures.

Chemical explosive mode analysis
In this section, chemical explosive mode analysis (CEMA) is utilized to characterize local combustion modes in terms of auto-ignition ( ), deflagration ( ), and diffusion-based extinction ( ).As discussed in Appendix A.2, auto-ignition and deflagration are defined through the local marker  =  s ∕  , in which the combustion mode can be considered the former when  magnitude is less than unity, and the latter when  is greater than unity.Additionally, local extinction is defined by  that is less than negative unity, for dominating inhibition effects from diffusion against ignition.Moreover, here a temperature threshold of value 1000 K is chosen to distinguish fresh mixture from the combusted one.Fig. 10 depicts the local combustion mode development starting from  2 .The plots correspond to the reaction zone within  CH 4 = 0.3 and 0.7, while the black dashed curve represents  CH 4 = 0.5 isocontour.According to CEMA, it is noted that, first, throughout ≈ 0.1 ms after τ 2 auto-ignition is the main mode in the reaction zone.During this phase, the front is highly wrinkled, and various local extinction zones are observed.Then, at the outer periphery of the spray, turbulent levels decay and deflagrative zones emerge, initially from zones with  CH 4 ≥ 0.5 (i.e.inside the black dashed contour).Overall, there is a transitional period of approximately 0.2 ms after τ 2 wherein complex combustion structures co-exist, including deflagration, auto-ignition, and local extinction.After that, slow propagation of the premixed flame front becomes the main mode.This analysis further confirms the timescale for flame initiation, i.e. τ 3 , previously discussed in Sections 4.4 and 4.5.Indeed, due to dynamic boundaries of the investigation domain, newly wrinkled auto-ignition zones may appear in the far bottom of the later snapshots, also reflected in Fig. 14.Further, in the bottom subplot of the figure,  is used to demarcate the various modes along a representative flamelet which is sampled at τ 2 + 0.24 ms.Such a sample flamelet exhibits a typical pattern of combustion modes as reported for classical premixed flames in the literature [50].
After demonstrating the combustion mode development through CEMA on 2D cutplanes, the overall 3D statistics are presented in Fig. 11.The pre-ignition reactive zone, herein defined by ( e > 0) ∩ ( CH 4 > 0.3), is inspected for auto-ignition, deflagration, and diffusionbased extinction.The auto-ignition regime (i.e.reaction dominating) is further decomposed into diffusion-promoting (0 <  < 1) and diffusion-inhibiting (−1 <  < 0) modes, although diffusion effects are insignificant during τ 2 time scale.Consistent with previous results, first, auto-ignition mode dominates the entire statistics at τ 2 .Shortly after that, deflagration zones emerge and start to grow while diffusion-inhibiting auto-ignition regions decay, i.e. diffusion mainly assists ignition.It is not before 0.2 ms when ≈ 50% of the normalized statistics become in deflagrating mode.According to the present analysis, even during the deflagrating fronts, reaction-dominating zones will be still evident, as it is further elaborated by the flamelet plot in Fig. 10.

Local turbulent combustion regime
After providing numerical evidence regarding the existence of emerging deflagrating fronts, the local premixed turbulent combustion regime is investigated.Using the classical Peters-Borghi diagram, scatter density plots for the local integral length scales (  ) and the turbulent fluctuations, normalized by stoichiometric canonical flame thickness and flame speed respectively, are depicted in Fig. 12.The plots are considered for the snapshot τ 2 + 0.32 ms in which the combustion mode is essentially deflagrative.As noted earlier, the formulation  ′ = √ 2∕3 is utilized for turbulent fluctuations while assuming isotropic turbulence.The integral length scales ( defined as   = ( ′ ) 3 ∕ ) span a wide range due to the large variation of the total kinetic energy dissipation .The black symbol denotes a mean value of the joint density function for both  ′ and   .
The figure indicates that turbulent combustion of the eq-DNS region can be mostly characterized by the corrugated flamelet regime, with local zones occurring in thin reaction and wrinkled flamelet regimes.When Karlovitz number Ka = ( ′ ∕ L ) 3∕2 (  ∕ L ) −1∕2 is below unity, deflagrative fronts are thinner than the smallest turbulent scales.Therefore, the inner flame structure is not affected by turbulence and the present combustion occurs in the flamelet limit even though the fronts themselves are wrinkled by turbulent motions on scales much larger than the flame thickness.In addition, laminar combustion regime is also recognized, which could be associated to local zones away from the spray.More importantly, the normalized scales are mainly concentrated about the Klimov-Williams criterion (i.e.Ka = 1 line), which is depicted by the island peak in blue and also by the circled-cross symbol.

Reaction flux analysis
One of the benefits behind detailed chemistry modeling herein considered is the ability to track particular reactions contributing to the methane flame initiation.This motivates the development of a new combustion mode transition index, which is purely based on reaction flux analysis.In this regard, the criterion previously introduced by Schulz et al. [37], namely the auto-ignition index, was established for SF methane flames.Nonetheless, in DF configuration involving ndodecane the validity of this metric becomes ambiguous.Particularly, while the aim is to investigate transition of n-dodecane ignition toward the ambient methane flame propagation, it is unclear how HO 2 elementary reactions would perform, as being involved in both the HRF and the ambient methane oxidation mechanisms.Accordingly, the analysis herein was decided to be based on the production/consumption reactions of OH, as the radical of choice, to distinguish deflagration from auto-ignition.
During auto-ignition (i.e.0D phenomenon), high OH production rates are expected.On the other hand, for a deflagrative flame where transport effects play an important role on radical diffusion, the net production rate will be in quasi-steady state along a particular progress variable, as consumption and production of OH become competitive.Such a hypothesis on the expected OH production rates at different instances is further supported by observations in Fig. 4, where the overall production rate ωOH peaks shortly after τ 2 , before which it gets reduced, possibly implying significant effects of particular consumption reactions.
A reaction flux analysis is performed to identify the important OH reactions along representative investigation lines for various snapshots.After ranking all OH reactions (totaling 67) based on their rates, the most important production (Eqs.( 6) and ( 7)) and consumption (Eqs.( 8) and ( 9)) reactions are noted to be (R0, R17) and (R25, R88), respectively.
The production/consumption rates of the aforementioned reactions are plotted along investigation lines, in Fig. 13, for two representative snapshots denoting auto-ignition (at τ 2 ) and deflagration (at τ 2 + 0.25 ms).The location of the sampling lines are chosen consistently with those in Fig. 5. First, it is noted that R17 is a high temperature, high pressure chain-branching reaction.Accordingly, it is associated with the hot reaction front, breaking up the meta stable H 2 O 2 and producing abundance of OH radical pool throughout the entire combustion phases, which is reflected in the figure.Second, regarding consumption reactions, R25 is considered a combustion completion reaction.Moreover, it is rather dependent on R88 products, such that methyl radicals (CH 3 ) initiate a pathway for CO formation, which is required for R25 to be triggered as elucidated in the following.
During τ 2 , reactions are dominated by R17.Temperature levels (cf.black curves in Fig. 13) are not sufficiently high, yet, to properly initiate R88 which has high activation energy required to break the stable CH 4 molecule.Moreover, at this timescale the abundance of CH 4 and OH may be affected due to spray evaporation, hence the dilution of CH 4 and OH, which affects R88 activation.At later time scales, i.e. after τ 3 , temperature levels are elevated and R88 peaks.The shift between R25 and R88 peaks is possibly due to their dependence, as previously mentioned.Once the pools of H radicals and CO molecules are built up sufficiently, R0 and R25 will peak.
The previous analysis inspires the introduction of a new flame initiation index ( ) which is entirely based on detailed chemistry, also can be considered an extension to Schulz's auto-ignition index for DF combustion.The combustion mode transition from n-dodecane autoignition to a premixed flame deflagration could be inspected through relative OH consumption and production rate magnitudes.A simple metric is defined as wherein the value of  = 2 is considered for normalization in this case, since R17 rate magnitude peak is noted to be almost twice as that for R88.Moreover, the choice of these particular rates in Eq. ( 10) can be justified where R17 and R88 are linked to the auto-ignition and the LRF consumption initiation, respectively.By making use of the newly proposed metric in Eq. ( 10), local combustion modes at different snapshots are distinguished for autoignition ( ) and deflagration ( ), as shown in Fig. 14.The metric  is herein used, with a threshold value of 0.35 to distinguish deflagration ( ≈ 1) from auto-ignition ( ≈ 0).Interestingly, the observed trends in the figure are consistent with results from CEMA, hence suggesting this detailed chemistry-based criterion to be a successful marker.
It is worth mentioning that the flame initiation index performance might be sensitive to the threshold choice.Moreover, it has been proposed and tested for a single DF configuration involving n-dodecane and methane.Following the same workflow, the reaction flux analysis based on OH chemistry can be argued as an efficient tool that is fuelagnostic.Therefore, it would be interesting to verify the feasibility of such a new criterion for different fuels at different DF operating conditions, while adapting the particular OH production and consumption Onset when over 90% of 2D data within reaction zone becomes deflagrative reactions according to the employed fuels.As a final note, a summary of the flame initiation timescale according to the various employed metrics is depicted in Table 2.

Conclusions
In this work, the combustion mode development of n-dodecane spray assisted dual-fuel combustion of methane under the modified ECN Spray A baseline is numerically investigated.A highly parallelized scale-resolved simulation is performed while using direct integration of finite-rate chemistry.Local numerical microscopy is conducted through extensive refinement of the large-eddy simulation (LES, 62.5 μm) in the resolution limit of direct numerical simulation (DNS, 3.9 μm).Such an embedded DNS/LES resolution approach enabled the study of the ignition kernel evolution toward the hypothesized deflagration.
Regarding the main results, first, turbulence levels at the ignition onset are reported for dual-fuel configuration.Then, combustion mode development after spray auto-ignition is investigated through various metrics.The performed analyses confirm the existence of emerging deflagrating fronts shortly after the spray second-stage ignition (τ 2 ).In addition, a new timescale (τ 3 ) is defined and estimated for the premixed ambient methane combustion, which can be also linked to the flame initiation onset.Moreover, the turbulent premixed combustion regime is investigated using the classical Peters-Borghi diagram.Finally, a new metric based on reaction flux analysis is proposed and verified against CEMA results.The overall findings are highlighted in the following.
1. Combustion mode transition from spray auto-ignition to flame deflagration is numerically confirmed from different perspectives, including the (i) diffusion-reaction balance (CEMA and energy transport budget), (ii) front structure statistics and thermal thickness, (iii) front displacement speed, and (iv) detailed chemistry and reaction flux analysis.2. The flame initiation timescale τ 3 is estimated to be approximately 0.2 ms after τ 2 , which is about 1.25-1.3τ 2 .
3. OH chemistry can be used as a marker to distinguish flame deflagration from spray auto-ignition.The newly proposed metric, based on reaction flux analysis, provides results that are consistent with CEMA. 4. Using the Peters-Borghi diagram, the premixed turbulent combustion in the isolated investigation domain is mainly characterized by the corrugated flamelet regime, in which the deflagrative fronts are thinner than all turbulent scales.Moreover, additional zones occur locally in the thin-reaction and wrinkled regimes, while laminar combustion is recognized far away from the spray.5.The combustion mode is a highly localized phenomenon, and deflagrative fronts can emerge, locally, as rapid as 0.05 ms after τ 2 .There is a transient period of ≈ 0.2 ms after τ 2 wherein mixedmode combustion occurs, including extinction, ignition fronts, and quasi-deflagrative structures.Later, after τ 3 = τ 2 + 0.2 ms, the combustion mode becomes essentially deflagrative.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Volume-rendered artistic illustration of the DF concept and underlying processes: (I) LRF port injection and premixing with in-cylinder oxidizer, (II) HRF pilot spray assisted multi-staged ignition, and (III) LRF combustion providing main (85-98%) energy.Focus is herein shifted toward small spatio-temporal window for auto-ignition kernel development and LRF consumption.

Fig. 2 .
Fig. 2. Domain discretization, resulting into total of 316 million cells.The LES grid (62.5 μm) is considered, to resolve spray dynamics and provide turbulence to the embedded DNS grid (3.9 μm) which is the main investigation domain.The eq-DNS grid is axially located 32 mm away from the spray injector, and it is shifted 5.5 mm and 0.75 mm away from the spray centerline axis in the other two perpendicular directions.

Fig. 3 .
Fig. 3. Flow topology at ignition onset (i.e.τ 2 ) within eq-DNS region of size 3 × 3.5 × 1 mm 3 : (a) 2D cut-plane of velocity magnitude, (b) PDF of velocity magnitude, (c) PDF of turbulent fluctuations, (d) 2D cut-plane of logarithmic  with ε denoting machine epsilon, (e) scatter plot of temperature in  space, and (f) scatter plot of temperature in  space.Both scatter plots are overlaid with conditional mean (red) and standard deviation (yellow) values.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 4 )
Fig. 4 depicts temporal evolution of the temperature (top half) and hydroxyl radical (OH) mass fraction (bottom half) fields in the eq-DNS region.The snapshots are presented for the fields starting from τ 2 with a temporal increment of  = 0.04 ms.Additionally, the upper snapshots are overlaid with isocontours of  st (white) and methane progress variable  CH 4 ≡ ( ( CH 4 −  u CH 4 )∕( b CH 4 −  u CH 4

Fig. 4 .
Fig. 4. Top: Cut-plane field plots of temperature overlaid by isocontours of stoichiometric mixture fraction ( st , white) and methane progress variable ( CH 4 = 0.5, yellow).Bottom: cut-plane field plots of OH radical mass fraction overlaid by isocontours of source terms ωOH at values 100 kg m −3 s −1 (green) and 1000 kg m −3 s −1 (red).Snapshots start from τ 2 with a temporal increment of  = 0.04 ms.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. Local energy budget along sampled investigation lines, normalized by maximum absolute value of reaction and diffusion terms.Rows denote snapshots at τ 2 and τ 2 +0.25 ms to represent reaction-driven and deflagrative modes, respectively.Solid and dashed curves in the 2D planes denote methane progress variable at 0.5 and 0.3, respectively.

Fig. 6 .
Fig. 6.Thermal thickness values for two representative 1D Cantera simulations at  = 0 and  st (left), single sampling line of the 3D data (middle), and scatter of many samples plotted against their mean mixture fraction values (right).

Fig. 8 .
Fig. 8. Temporal evolution of density weighted displacement speed, averaged along the reaction front defined by  CH 4 = 0.5, and divided by canonical laminar premixed flame speed value under stoichiometric conditions.Filled region denotes the corresponding standard deviation.Front speed becomes deflagrative, on average, after τ 2 + 0.2 ms.

Fig. 9 .
Fig.9.Left: PDFs of  along reaction front at various time instances.Right: Two approaches to locate ambient burn onset, first one defined when mean() falls below stoichiometry, and second one is defined when statistics of the unburned data inside spray and burned data outside spray both intersect with global minimum.

Fig. 10 .Fig. 11 .
Fig.10.CEMA markers applied to 2D planes for combustion mode identification.Overall, there is a transitional period of ≈ 0.2 ms after τ 2 wherein mixed-mode combustion structures co-exist.After that, slow deflagrating fronts become the main mode.Last subplot denotes combustion modes along a representative flamelet sampled at τ 2 + 0.24 ms.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Turbulent combustion regime using classical Peters-Borghi diagram for snapshot τ 2 +0.32 ms.Scatter density plot denotes local integral length (  ) and turbulent fluctuations ( ′ ) normalized by stoichiometric 1D flame thickness and flame speed, respectively.Symbol denotes mean of the joint density function for  ′ and   .(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 13 .
Fig.13.Reaction flux analysis for OH production/consumption rates at two representative snapshots denoting auto-ignition vs. deflagration.After ranking all reactions, top production () and consumption ( ) rates are plotted along sampling lines.The shift between peaks is due to reactions dependence.The plots motivate a new metric for combustion mode transition based on relative consumption to production rate magnitudes.

Fig. 14 . 2
Fig. 14.Flame initiation index based on reaction flux analysis.Red zones indicate auto-ignition, while green zones suggest a premixed flame.Results are consistent with CEMA.

Table 1
Summary of thermophysical operating conditions, together with associated computational settings, for the DF ECN Spray A numerical simulation.