Efficient two-dimensional simulation of primary reference fuel ignition under engine-relevant thermal stratification

Despite vast research on engine knock, there remains a limited understanding of the interaction between reaction front propagation, pressure oscillations, and fuel chemistry. To explore this through computational fluid dynamics, the adoption of advanced numerical methods is necessary. In this context, the current study introduces ARCFoam , a computational framework that combines dynamic mesh balancing, chemistry balancing, and adaptive mesh refinement with an explicit, density-based solver designed for simulating high-speed flows in OpenFOAM. First, the validity and performance of the solver are assessed by simulating directly initiated detonation in a hydrogen/air mixture. Second, the study explores the one/two-dimensional (1D/2D) hotspot ignition for the primary reference fuel and illuminates the impact of transitioning to 2D simulations on the predicted combustion modes. The 2D hotspot simulations reveal a variety of 2D physical phenomena, including the appearance of converging shock/detonation fronts as a result of negative temperature coefficient (NTC) behavior and shock wave reflection-induced detonation. The main results of the paper are as follows: (1) NTC chemistry is capable of drastically changing the anticipated reaction front propagation mode by manipulating the local/global reactivity distribution inside and outside the hotspot, (2) subsonic hotspot ignition can induce detonation (superknock) through the generation of shock waves and subsequent wall reflections, and (3) while the 1D framework predicts the initial combustion mode within the hotspot, significant differences between 1D and 2D results may emerge in scenarios involving ignition-to-detonation


I. INTRODUCTION
In the pursuit of optimizing spark-ignition (SI) engines for improved efficiency and reduced emissions, understanding and mitigating the phenomenon of knock is of paramount importance.As an abnormal combustion phenomenon, knock may lead to increased engine temperatures, pressure oscillations, and potential engine damage. 1,2Typically, the chemistry involved in the combustion process within SI engines is highly complex with low-, intermediate-, and high-temperature reactions.Clearly, the interaction between fuel ignition chemistry and compressible flow physics plays a crucial role in the combustion dynamics and induced pressure levels in SI engines. 3herefore, it is essential to thoroughly investigate the complex phenomena relevant to knock in SI engines.In this study, we address this challenge through advanced computational fluid dynamics (CFD) and finite rate chemistry approaches.The primary reference fuel (PRF) 87 (13% n-heptane/87% iso-octane by volume, referred to as PRF hereafter) mixture is selected as the gasoline surrogate for our investigations.
By employing these numerical techniques and the representative mixture, we aim to gain deeper insights into the knock phenomenon.
Knock in SI engines is a result of spontaneous, abnormal combustion of the air-fuel mixture that occurs ahead of the flame front. 4ommonly, knock is classified into three categories: mild knock, strong knock, and superknock. 5,6Mild and strong knocks involve isolated pressure oscillations, DP % 5-24 bar, caused by localized hot spots, while superknock is an extreme event characterized by violent and sustained pressure oscillations, DP ' 50 bar. 5,6Furthermore, superknock is an irregular event that progresses from pre-ignition to intense knocking, often accompanied by detonation. 5It is typically triggered by the interaction of a shock wave with a propagating reaction front either on the cylinder surfaces or within a localized hotspot with favorable heat release and reactivity gradient. 7All forms of knock have adverse effects on engine performance, emissions, and durability.
Experimental investigations have traditionally been vital for understanding knock in SI engines, but they are time-consuming and resource-intensive. 1,2Computational fluid dynamics (CFD) is a powerful tool that complements experimental research by allowing analysis of knock under controlled conditions.To accurately capture the behavior of knock, it is essential to consider a finite rate chemistry that accounts for the interactions between fuel molecules, radicals, and intermediate species.The original work by Zeldovich 8 identified four modes of propagation for hotspot-induced ignition fronts in onedimensional (1D) simulations, as later expanded upon by Gu et al. 9 for H 2 -CO/air.These modes include thermal explosion (spontaneous ignition), supersonic ignition, developing detonation, and subsonic ignition and deflagration, all originating from a hotspot.Building upon this approach, Bradley et al. 10 introduced a quantitative regime diagram known as the detonation peninsula.This diagram, based on two dimensionless parameters, was initially developed for H 2 -CO/air mixtures under controlled auto-ignition (CAI) engine conditions.The detonation peninsula and the associated 1D framework have since been utilized to assess knock and superknock tendencies in SI engine scenarios. 6,11,12imulations utilizing the 1D framework have revealed complex reaction-pressure wave interactions in large hydrocarbon mixtures with low-temperature chemistry (LTC).4][15] Coolspots, multiple ignition kernels/fronts, and shock waves were identified in the presence of the fuels posing negative temperature coefficient (NTC) chemistry. 15,16Studies on n-heptane mixtures highlighted different chemical length scales, 17 low-frequency amplification of shock waves, 18 and pressure wave intensification near closed vessel walls 19 as a result of NTC.Moreover, water vapor dilution affected the NTC effects and the correlation between hotspot radius and maximum pressure in detonative cases. 20While Bates et al. 12 illustrated the detonation peninsula for iso-octane/air and n-heptane/air mixtures, a direct comparison with the original peninsula was not provided.In a recent study, Shahanaghi et al. 21investigated the prevalence of deflagrative and auto-ignitive combustion regimes for PRF and PRF-ethanol (18% n-heptane/62% iso-octane/20% ethanol by volume, PRF-E) blends.They found that NTC chemistry favored the auto-ignition regime while promoting a blended mode.In a follow-up investigation by Shahanaghi et al., 13 the hotspot auto-ignition of the aforementioned blends was studied using the 1D framework.The ignition regime diagrams were reproduced for PRF and PRF-E mixtures with air.To construct the regime diagrams, hundreds of 1D simulations were carried out for different hotspot initial conditions posing a varying level of temperature stratification.The reproduced ignition regimes were compared with the original Bradley peninsula.The results demonstrated that while NTC chemistry increases the ignitability of the mixture under specific temperature stratification levels, it may suppress detonatability (superknock) by reducing the distance traveled by a developing detonation front.
Apart from its application in simplified 1D setups, the Bradley diagram is also beneficial in multi-dimensional scenarios, offering insights into the role of chemistry in the emergence of knocking.In this context, Robert et al. 22 conducted large eddy simulations (LES) with tabulated chemistry and employed the Bradley diagram as a local detonation indicator.It was shown that the diagram predicts combustion regime transition based on Spark timing while roughly estimating the time and location of deflagration to detonation transition (DDT) within the chamber.However, it was noted that resolving the turbulent reacting flow is a necessity; otherwise, the detonation indicator fails in various cases.The turbulent time and length scale ratio can have a major role in reactivity stratification and subsequently detonation development.Two-and three-dimensional (2D and 3D) direct numerical simulations (DNS) of Luong et al. 23,24 demonstrated that reducing the length scale of temperature stratification decreases the likelihood of superknock occurrence.Moreover, LES results of Wei et al. 25 showed that increased turbulence intensity delays auto-ignition formation and reduces the ignition kernel size.Nevertheless, the simplified chemical kinetics models utilized in these studies might pose constraints for fully capturing the complexity of knock chemistry.On the other hand, utilizing detailed chemistry in such multi-dimensional simulations is computationally very expensive.Therefore, there is a need for computational methods that can balance computational efficiency and accuracy in the numerical investigation of knock using finite-rate chemistry.
The open-source toolbox OpenFOAM 26,27 has gained recognition for its versatility in handling complex reactive flow processes.With comprehensive chemical reaction mechanisms, OpenFOAM enables the integration of detailed chemistry models into CFD simulations, providing a more accurate representation of knock.To improve simulation efficiency, OpenFOAM utilizes an h-adaptivity approach, i.e., adaptive mesh refinement (AMR). 26AMR selectively increases mesh resolution in relevant regions, such as flame and shock fronts, while maintaining lower resolution in surrounding parts of the flow. 28,29dditionally, OpenFOAM version 10 introduced dynamic load balancing (DLB) as another strategy for optimizing computational efficiency in parallel simulations.DLB dynamically redistributes the computational load on local cells among multiple processors in realtime, thereby optimizing the solver's performance.
In addition to the official implementation of DLB (referred to as OF-DLB henceforth) and AMR in OpenFOAM, the community has also developed these libraries.For instance, DLB has reportedly been developed for mesh [30][31][32] and chemistry problem 33 distribution during the run-time.However, it is worth noting that to date, no direct comparison has been made between the OF-DLB and those developed by the community for DLB.Furthermore, the original 3D refinement functionality in OpenFOAM's AMR engine has been expanded to include 2D refinement (OF-AMR). 30,34,35Subsequently, the 2D AMR engine has been improved and updated to work with multi-criterion refinement and the community-made DLB technique. 32,36ased on the current literature, the hotspot ignition of large hydrocarbon mixtures with low-temperature chemistry (LTC) involves a complex interaction between reaction and pressure waves.However, conducting detailed chemistry simulations for such phenomena is computationally demanding, while simultaneously requiring a high level of spatial and temporal accuracy to capture the underlying physical processes of knock.Consequently, high-fidelity CFD simulations of gasoline surrogate hotspot ignition under engine-relevant conditions have been limited to 1D, with 2D/3D yet to be explored.In order to overcome this limitation and enable high-resolution simulations that account for finite-rate chemistry in 2D, an efficient simulation approach is needed.
In order to address the aforementioned research gap and build upon our previous studies, 13,21 the main goal of this study is to examine the propagation patterns in 2D hotspot-initiated ignition of PRF mixtures, under conditions relevant to SI engines.To achieve this goal, we have developed an effective computational solver called ARCFoam [ARC ¼ adaptive Runge-Kutta (RK) central] in the OpenFOAM framework.This solver integrates OF-DLB, AMR, and an efficient chemistry solver, namely, DLBFoam, 33,37 to facilitate 2D simulations of hotspot ignition.The solver is first tested in direct initiation of hydrogen detonation problems to demonstrate its performance and validity of the numerical schemes.Then, it is utilized in hotspot ignition for the PRF mixture in 2D cylindrical simulations with finite rate chemistry.The primary objectives of the present study are as follows: (1) Provide a solution for integrating OF-AMR with OF-DLB framework capabilities to enhance computational efficiency.(2) Develop a density-based solver based on OpenFOAM enabling 2D simulations of combustion regimes resembling knock, involving shock-chemistry interactions with finite-rate chemistry.
(3) Expand our previous research on 1D hotspot ignition simulations 13 of the PRF mixture to incorporate 2D cylindrical configuration, investigating the propagation of the ignition front beyond the borders of the hotspot region.(4) Explore and compare various ignition-front propagation patterns observed in our previous 1D study 13 within the context of 2D simulations, aiming to identify similarities and discrepancies between the 2D and 1D combustion regimes.

II. METHODOLOGY A. Theoretical background
As a theoretical background on ignition front propagation, Zeldovich 8 proposed a method to determine the propagation speed of an ignition front, S ign .Based on the theory, S ign can be correlated with the inverse of the ignition delay time (IDT) gradient of the premixed reactants.Assuming a homogeneous mixture, IDT can be expressed as a function of temperature, so that S ign can be written as where s i represents the IDT, a ¼ @s @T denotes the sensitivity of IDT to temperature variations, and rT is the temperature gradient.Note that in Eq. ( 1), the chain rule has been utilized, i.e., @s @x ¼ @s @T @T @x .Two types of waves are observed during auto-ignition: a reaction wave representing rapid chemical reactions and a compression wave generated by the thermal explosion of the mixture (reaction shocks). 38Zeldovich categorized ignition fronts into four groups based on their propagation speeds.These include deflagration, subsonic ignition, supersonic ignition, and detonation.The critical temperature gradient (rT c ) associated with this process is expressed as where U a denotes the local speed of sound.The Zeldovich theory has been further developed by Bradley and colleagues 9,10 to analyze ignition regimes within a spherical hotspot.A regime diagram, referred to as the detonation peninsula, has been established using a 1D framework that relies on two dimensionless parameters (n, e), where and In the above formula, s e is the excitation time and r h is the hotspot radius.The parameter n quantifies the coupling between the acoustic wave and the reaction front, while the parameter e represents the ratio between the chemical timescale and the acoustic timescale in the hotspot.The detonation peninsula encompasses the range of n values that lead to developing detonation inside a hotspot.Additionally, the parameter e describes the rate at which chemical energy is transferred to the pressure wave.It is important to note that the values of s i , s e , U a , and a are determined at a reference temperature (T ref ), and the choice of T ref can affect the values of n and e.In this study, T ave ¼ T (r ¼ r h /2) is chosen as T ref to maintain comparability with the original detonation peninsula presented by Bradley et al. 10 The 0D constant volume reactor model of Cantera 39 is used to calculate s i values for the PRF mixture, from which a values are estimated.Here, s e is defined as the time of reaching 5% of the maximum heat release in the constant volume reactor simulations.Finally, U a is calculated using the thermodynamic library of Cantera.

B. Initialization and propagation patterns
The main focus of this paper is to study the formation of autoignition waves in a constant volume domain, both in 1D planar and 2D cylindrical configurations.In this setup, the hotspot is represented by a localized linear temperature gradient with specific values for the average temperature (T ave ), temperature gradient (rT), and radius (r h ).The initial pressure and mixture composition are assumed to be uniformly distributed throughout the domain.It is important to note that when simulating the combustion of gasoline surrogates, we encounter a distinctive characteristic known as the negative temperature coefficient (NTC) effect. 40A detailed analysis of these different scenarios for the PRF mixture is provided in our earlier study 13 for 1D planar configurations.Here, we target to represent these cases for 2D cylindrical configurations, which are more relevant to practical engine conditions.The NTC effect is represented by the pink region in IDT vs temperature profile of PRF at 50 bar, Fig. 1(a), and occurs under engine-relevant average temperatures and pressures.
Three distinct zones (I-III) are indicated in Fig. 1(a).Zone I includes both I h and I c , denoting a monotonous IDT variation with temperature.Within zone I, I h and I c exhibit the minimum IDT corresponding to the maximum and minimum temperature, respectively.Zone II identifies the region where the IDT vs temperature profile becomes non-linear due to NTC, with the minimum IDT value not aligning with either the minimum or maximum temperature.Zone III indicates the temperature range wherein a local IDT minimum exists while still maintaining the max/min IDT corresponding to the min/ max temperature.The presence of the NTC effect may lead to various modifications in the initiation and propagation patterns of hotspot ignition.
Figure 1(b) presents a schematic representation of the hotspot temperature profiles at the initial time (i.e., when rT < 0). Figure 1(c) portrays the expected ignition front propagation pattern when the average temperature (T ave ) and the corresponding T max min values fall within regions characterized by relatively high or low temperatures, indicated by the blue color and denoted as I h in Fig. 1(a).In this pattern, the ignition originates from the center of the hotspot, where the temperature is the highest (T max ¼ T ave À ½rT Â r h ), and then propagates radially outward toward the region where the temperature is lower (T min ¼ T ave þ ½rT Â r h ).This propagation pattern has been extensively studied in the literature.Furthermore, the term "coolspots" refers to the scenarios where the average temperature (T ave ) and the corresponding T max min values are located within the NTC region (labeled as I c in Fig. 1(a).In such cases, ignition starts from T min due to its lower local IDT.In the present work, to ensure that the ignition pattern within coolspots resembles that of normal hotspots, i.e., Pattern I in Fig. 1(c), the simulations are initialized with an inverse (positive) rT from the center.
Off-centered ignition kernels occur in scenarios where T max min are located on two sides of the NTC region, as depicted in Fig. 1(a) zone II.In such situations, the minimum IDT value does not align with either one of these temperatures.Instead, the highest reactivity corresponds to an intermediate temperature between these two points, leading to ignition initiation from the middle location inside a hot/ coolspot.As a result, two fronts, one propagating inward and the other outward, appear simultaneously, as shown in Pattern II in Fig. 1(d).This dual-front propagation pattern is a characteristic feature of offcentered ignition kernels and can lead to unique combustion behaviors compared to typical hotspot configurations.
Finally, in cases where the temperature range encompasses the NTC region [T min ( T(I c ) ( T max ], as depicted in Fig. 1(a) zone III, a secondary ignition kernel can emerge.This secondary ignition kernel appears due to the localized high reactivity within the hotspot.Initially, ignition starts from the center of the hotspot, where the temperature is the highest.However, as the combustion process evolves, a secondary ignition kernel may form ahead of the primary ignition front, as illustrated in Pattern III in Fig. 1(e).
C. ARCFoam (ARC 5 adaptive, Runge-Kutta, central) In the present study, the open-source CFD package OpenFOAM 27 is used to solve the 2D reacting simulation problems.
ARCFoam provides a fully compressible, density-based solution of the flow with explicit time stepping.Additionally, ARCFoam incorporates AMR and load balancing (mesh, chemistry) strategies in order to accelerate the finite rate chemistry solution of reacting high-speed flows.More details regarding ARCFoam's framework are provided in the following after introducing the governing equations.

Governing equations
The governing equations studied herein are the Navier-Stokes (N-S) equations for mass, momentum, species concentration, and enthalpy.The N-S equations can be written as follows: @ðquÞ @t þ r Á ðquuÞ ¼ r Á s À rp; In the above equations, q, u, p, Y k , N sp , h, s, a t , and D denote density, velocity, pressure, mass fraction of the species k, total number of species, sensible enthalpy, viscous stress tensor, and thermal and mass diffusivities, respectively.The production rates of each specie and heat release rate (HRR) are represented by _ x k and _ x h , where , in which Dh 0 f ;k is the enthalpy of formation.In this study, the unity Lewis number assumption (D ¼ a t ) has been employed for all species in the simulations.This decision is motivated by two reasons.First, the estimated effective Lewis number 41 for the PRF/air mixture in hotspot simulations, calculated using Cantera, 39 is approximately unity.This suggests that the magnitudes of diffusivity and thermal conductivity are comparable.Hence, this assumption aligns well with the characteristics of the PRF/air mixture.][44][45] Nevertheless, by adopting the unity Lewis number assumption, the effects of diffusion on species transport in H 2 /air mixture are simplified, without completely disregarding them.This approach results in a more tractable and computationally efficient representation of the detonation process in the studied mixtures, and the study on possible differences remains for future work.Additionally, viscosity and thermal conductivity are calculated using the Sutherland law and Eucken approximation, respectively. 46The ideal gas law and thermal equation of state are also employed to close the system of equations.

Numerical procedure
ARCFoam is based on the rhoCentralFoam 47 solver.As a part of OpenFOAM, rhoCentralFoam is a density-based solver utilizing second-order central, Kurganov and Tadmor (KT) 48 discretization scheme for capturing shock waves in high-speed non-reacting flows.Here, the time integration of rhoCentralFoam is modified to the second-order explicit Runge-Kutta (RK) scheme.Moreover, similar to Ref. 49, the species transport equations together with the chemical reaction source terms are included in the solver to enable the solution of reacting flow problems.Within each RK sub-step, a firstorder operator splitting scheme 50,51 is employed to advance the flow field and chemical source terms over time.The chemical source terms are calculated using direct integration of the finite rate chemistry.The Jacobian matrix for the system of kinetic ordinary differential equations (ODEs) is calculated analytically using pyJac 52 and its recent implementation to the OpenFOAM framework, DLBFoam v.1.1. 37hese numerical techniques have been successfully used for various compressible reactive flow simulations. 18,37,49Figure 2 provides an overview of the flow chart detailing the various processes involved in each time step during run-time.The diagram demonstrates essential steps, including mesh pre-processing, such as Mesh update and balance, and the solution of the flow field equations.
Each time step starts with the mesh update (AMR) followed by mesh balancing, which is discussed in Sec.II D. After that, convective fluxes (i.e., interpolants of qu; quu; quY k , and quh) are calculated using the KT central scheme.Next, the chemistry is updated.First, the system of kinetic ODEs of all the reacting cells is distributed evenly between the processors by DLBFoam.Furthermore, the chemistry problems are solved using the DLBFoam v.1.1'smodified ODE solver.Finally, the chemical source terms ( _ x k and _ x h ) are updated.For more details on the load balancing and chemistry solution methods, the reader is referred to Refs.33 and 37.In the last step of the RK iteration, the conservative variables (i.e., q; qu; qY k and qh) are updated from Eqs. ( 5)- (8).At the final stage of each time step, the primitive variables (u; h; Y k ; p) are updated.
In OpenFOAM, the effect of basic boundary conditions is applied on the cell faces by adjusting the diagonal terms within the coefficient matrix and source terms of the discretized linear system of N-S equations. 26Notably, in OpenFOAM, parallelization is achieved by geometrical domain decomposition, which creates additional boundaries in the decomposed domain, i.e., processor boundaries.Normally, the field values (i.e., velocity, density, pressure, etc.) need to be communicated between the processors across the processor boundaries.Explicit time integration of N-S equations facilitates parallel simulations using ARCFoam.The implicit formulation of N-S equations enforces additional communications since the processor boundary values are needed to converge in the solution of the linear system of equations.However, the explicit formulation does not require the additional communication.Nevertheless, utilizing AMR in parallel simulations may change the efficiency of the solver.Details on AMR and load balancing are provided in Sec.II D.

D. Mesh update and balance
In a parallel run, AMR may change the initial well-distributed (balanced) mesh to clusters of cells concentrated on a single or a few processors.Consequently, the clustered cells may increase the runtime of the respective processors compared to the other processors.Therefore, re-distributing (balancing) the mesh during the parallel run with active AMR is a crucial step to maintain the computational feasibility of the simulation.It should be noted that in ARCFoam, mesh update and balance functions are invoked in a pre-processing stage before each CFD time step, as depicted in Fig. 2.These functions might be executed separately at different intervals during the runtime.Figure 3 depicts the flow chart of AMR and balancing processes.It is important to recognize that the frequency of these function calls and the time taken by these processes can significantly affect the overall runtime, and they should be calibrated according to the total number of cells and processors.Sections II D 1 and II D 2 discuss the AMR and balancing processes in more detail.

Adaptive mesh refinement
By default, OpenFOAM supports 3D AMR limited to hexahedral meshes.However, in recent years, the default AMR engine has been extended to accommodate 2D problems by the community. 30,34,35igure 4(a) presents a schematic of a two-level refined grid in 2D problems, with two refinement regions indicated by blue/purple (w 1 =w 2 ) FIG. 2. A schematic presentation of the main processes taken within a time step in reactive compressible solution in ARCFoam.and yellow/red (w 0 1 =w 0 2 ) colors.Moreover, the flow chart of the different steps during the OF-AMR process, which is taken prior to the CFD time step, is shown in Fig. 3.The time step period in which the OF-AMR is executed is based on a user-defined parameter, i.e., refineInterval.

Physics of
In the first step of the OF-AMR process, a refinement criterion (normalized error function) is calculated for the whole domain, as shown in Fig. 3.In the present work, based on Refs.53 and 54, the ratio of the second-order to the first-order derivative terms in the Taylor expansion of density is considered as the refinement criterion.Accordingly, grid cells are marked for refinement or unrefinement if the calculated error is higher than the user-defined refineLevel or lower than the user-specified unrefineLevel, respectively.
In the second step and during the 2D refinement process, the marked (parent) cells undergo subdivision into four smaller (child) cells, each with a quarter the size.It is notable that cells are not refined in the direction normal to the empty boundaries.Moreover, the newly created child cells acquire a refinement "level" that is increased by (þ1) compared to their parent cells, with the initial level set to 0, as illustrated in Fig. 4(b).
During the third step of the OF-AMR process, the marked refined (child) cells are merged to restore a coarser cell.In the conventional method of mesh unrefinement, the cells that were added during the refinement process are removed.This approach inherently means that the OF-AMR begins with a coarse mesh and the mesh after unrefinement can only have the same (level 0) or higher level than the original mesh. 26To accomplish this, OF-AMR utilizes a "tree" data structure (refinement-tree), i.e., refinementHistory, which stores the refinement changes and allows for cell removal upon request.Depending on whether the refinement process is operating in 2D or 3D, it may also be referred to as a Quadtree or Octree refinement.An example of 2D AMR from the present paper simulations, with five levels of refinement, is depicted in Fig. 4(c).The figure illustrates the mesh and temperature distribution for the Richtmyer-Meshkov (RM) instability observed in the H 2 detonation simulation in Sec.III B.

Dynamic load balancing
Since the release of its version 10, OpenFOAM 27 has been equipped with a new framework for dynamic load-balancing, decomposition, and redistribution. 55This framework provides the capability for users to calculate the computational load on each processor (CPU load) during parallel runs, currently available for chemistry solutions.The loadBalancer then translates the CPU loads into cell-based weights, which are utilized by the distributor, zoltan 56 or FIG. 3. A schematic presentation of the preprocesses before a time step in CFD solution in ARCFoam.Adaptive mesh refinement and load balancing stages are controlled independently by user-defined refine and redistribution intervals.scotch, 57 to redistribute the mesh based on the cell CPU loads during the runtime.
As mentioned earlier, AMR typically generates local clusters of cells within the domain, as shown in Fig. 4. In parallel runs, these clusters can significantly increase the computation time on the processors responsible for that particular part of the domain.This situation can greatly reduce the efficiency of the solver, especially when dealing with a large number of cells.To address this issue, a new load-balancing framework has been developed by authors to extend OF-DLB specifically for the OF-AMR engine.This framework is provided with ARCFoam, aiming to optimize the distribution of the computational load and to improve the solver's performance.
The additional feature offered by the ARCFoam's extended DLB framework is the ability to handle the distribution of the refinementHistory in OF-AMR.In situations where child cells with the same parents are distributed across multiple processors, the OF-AMR engine faces difficulties in merging these child cells into a single parent cell.Essentially, the refinement-tree becomes incomplete on each processor because the entire parent cell is not located on any specific processor.To address this challenge, the current implementation of the extended DLB is based on a solution proposed in Ref. 31.This approach constrains the OF-DLB to create new balanced cell distributions that encompass the entire parent cells (level zero).By ensuring that the parent cells are fully present in one processor, the implementation enables the correct merging of child cells and maintains the integrity of the refinement-tree during the parallel execution.
Similar to OF-AMR, the frequency at which the load balancing process occurs is determined by the user through the redistributeInterval parameter, specifying the number of time steps between each load balancing operation.It is important to note that the mesh redistribution process involves parallel communication between the processors, and if used excessively, it can increase the overall runtime of the simulation.To control this, a secondary parameter, defining the maximum imbalance between the processors denoted as maxImbalance (i max ), needs to be set by the user.The maximum imbalance across all processors (n proc ) is calculated using the following equation: where (n local ) is the number of cells per processor and n ideal ¼ n total = n proc represents the ideal number of cells per processor when the total cells (n total ) are evenly distributed across all processors.The flow chart illustrating the various steps of the load balancing process can be seen in Fig. 3.
In the first step, the visible cells (which include both child cells and unrefined parent cells) for each processor are counted, and the number of child cells is stored as weights on the parent cell maps.If a parent cell remains unrefined, its corresponding weight is equal to one.Figure 4(a) visually represents the cell clusters, visible cells, and the corresponding weights for each cluster.Specifically, the cell cluster in the lower-left corner represents the visible cells that contain level 1 and 2 child cells.The weights w 1 and w 2 indicate the numbers of level 1 and level 2 cells, respectively, as denoted by the blue and purple colors in Fig. 4(a).The same comprehension is applicable to upper-right cluster where the weights are indicated by w 0 1 and w 0 2 in Fig. 4(a).
Subsequently, these weights are projected onto their respective coarse cell maps, as depicted in Fig. 4(b).
In the second step, the coarse cell maps, along with their corresponding weights for each processor, are sent to the balancer (distributor), such as zoltan 56 or scotch, 57 to generate new cell maps.This process optimizes the distribution of cells across processors based on the assigned weights.In the final step, the fvMesh's distribute function is called to exchange new cells, including their associated fields, among the processors.This function also handles the creation of new boundaries based on the updated cell maps.Additionally, the refinement-tree is distributed among the new host processors, and the refinementHistory is reconstructed for each processor, ensuring the integrity of the refinement information on all processors.
It is worth emphasizing that in ARCFoam, the DLBFoam and extended OF-DLB have distinct roles.DLBFoam addresses the imbalance caused by stiff chemistry in thin reaction fronts, while the extended OF-DLB ensures a balanced distribution of cell clusters created by OF-AMR across the processors.These libraries are applied at different stages of the solution procedure, see Fig. 2, with DLBFoam being used during the chemistry update phase and the extended OF-DLB being employed during the pre-processing.As a result of their distinct purposes and execution in different solution phases, these libraries are considered independent components.
As a final remark, the number of processors in the 2D cylindrical parallel simulations was increased as the simulation progressed.This adjustment was necessitated by growing cell counts from expanding reaction fronts and the addition of new cells by AMR engine, see Sec.II E. To maintain efficiency, the process required pausing the simulation, performing "re-decomposition," and then resuming from the last time step.Similar to dynamic cell distribution in standard domain decomposition algorithms within OpenFOAM, the current implementation does not preserve the integrity of the refinementHistory on each processor.To address this, a customized utility named redecomposePar was created.This utility controls domain decomposition by ensuring the integrity of refinementHistory, following the approach proposed by Voskuilen. 31

E. Case setup and initialization
In Secs.II E 1 and II E 2, we present the computational setups and initiation methods employed for both the 1D planar and 2D cylindrical simulations of H 2 /air detonation and PRF/air hotspot ignition cases.

Direct initiation of hydrogen cylindrical detonation
The simulation involves a stoichiometric H 2 /air mixture using the chemical kinetics mechanism developed by Marinov et al., 58 which consists of ten species and 27 reactions to address the finite rate chemistry.For the 2D simulations, a rectangular domain representing a quarter area of the cylinder is considered.Symmetric boundary conditions (BCs) are applied at the bottom (y ¼ 0) and left (x ¼ 0) sides, while the right and top boundaries are treated as walls.Homogeneous temperature and pressure are initially considered along the domain with an energy source located on the bottom left corner.On the other hand, for the 1D simulation, the left boundary is set to symmetry, while a wall condition is applied to the right boundary.
The specific values for the initial pressure and temperature (P 0 , T 0 ), as well as the pressure, temperature, and radius of the energy  38 which is calculated using the SDToolbox. 59Furthermore, to induce transverse waves in 2D simulations, a sine wave with an amplitude of A ¼ 0.3 r s and a frequency of x ¼ 12 is superimposed on the initial pressure and temperature profiles of the energy source.It is notable that the initial perturbation has been demonstrated to have no discernible impact on the decay of the average detonation velocity. 60he selected source pressures (P s ) for the present study, as listed in Table I, ensure that both configurations fall within the supercritical regime. 61Notably, there is an order of magnitude difference between the P s values for 1D and 2D configurations, consistent with the observed variation in critical initiation energy between planar and cylindrical setups. 61Additionally, according to the conducted mesh sensitivity study in Appendix A, a minimum grid spacing corresponding to 12 cells per induction length is sufficient to achieve mesh convergence.Furthermore, in Appendix A, the results have been validated against prior numerical studies, 43,62,63 which include studies conducted using OpenFOAM, both with static mesh 43 and AMR 63 configurations.
In both 1D and 2D setups, initially, a refinement region is considered for the interior section containing the energy source (i.e., r < r s ) with five levels of refinement (corresponds to a total number of N cells ¼ 0:3 M in 2D settings).However, as the blast wave is generated and the subsequent detonation propagates, the front location undergoes refinement using the OF-AMR.This progressive refinement results in a gradual increase in the mesh cell count, in the case of 2D simulations, reaching N cells % 9 M as the front expands radially.Consequently, the number of cores was adapted in response to the escalating mesh count during the simulations using the redecomposePar utility, rendering the complete simulation runtime inaccessible.The simulations were executed utilizing up to 100 CPUs, with an approximate allocation of 70 000 cells per processor (N cpp ).

2D hotspot ignition of the PRF mixture
The simulations are conducted using a stoichiometric mixture of PRF/air as the reactants.To address the finite rate chemistry, we utilize the 115 species and 856 reactions kinetic mechanism developed by Stagni et al. 64 It is notable that in our previous study, 21 we validated this mechanism against experimental data and a detailed kinetic mechanism using high-pressure 0D constant volume reactor and 1D laminar flame simulations in Cantera. 39The simulation case is a 2D rectangular domain that corresponds to a quarter of the cylinder's area with symmetric boundary conditions on the bottom (y ¼ 0) and left (x ¼ 0) boundaries, while the upper and right boundaries are modeled as solid walls.The 2D hotspot is initiated using a linear temperature profile similar to the 1D framework, details in Sec.II B. Furthermore, uniform temperature and pressure conditions are initially assumed throughout the rest of the domain.Table II provides the configuration details for the 2D cases, including the domain length (L D ), minimum mesh size (Dx min ), number of cells (N cells ), number of processors (N proc ), and characteristics of the detonation front, such as the von Neumann (ZND) pressure (P z ), equilibrium (CJ) speed (Ucj), and induction length (D i ).
As indicated in Table II, the estimated induction length (D i ) using the ZND model is 1 Â 10 -5 m.In our previous work, 13 it was demonstrated that utilizing a solution method with a grid resolution of 12 points per D i ensures grid-independent results and captures the ZND detonation structure and its peak pressure.It is notable that 1D simulations are conducted for the selected points on the map, utilizing the same mesh and domain sizes.The additional 1D simulations aim to extend our previous investigation and explore the propagation and potential transition of the ignition regimes beyond the boundaries of the hotspot in the 1D context.
The 2D cylindrical simulations are initiated by refining the interior of the hotspot with approximately 2 Â 10 6 cells, while the surrounding areas have a lower cell count.As rapid ignition occurs, pressure waves emerge, prompting the OF-AMR engine to refine the fronts down to the minimum mesh size (Dx min ).Subsequently, the total number of cells gradually increases to (%) 90 Â 10 6 cells as pressure/density waves form and the front propagates radially.Similar to the previous case, the number of cores was scaled up using redecomposePar utility as the mesh count expanded during the simulations.Therefore, the availability of a unique CPU time for the entire simulation is not feasible.The 2D simulations demanded substantial computational resources, utilizing up to 5000 processors, which translated to approximately N cpp % 18 000, ensuring a reasonable runtime for each simulation case.

III. RESULTS AND DISCUSSION
In Secs.III A-III C, we first report and analyze the performance of the ARCFoam framework through two distinct test cases.Next, we

A. ARCFoam performance report
In this section, we report the performance of ARCFoam through 2D simulations of both H 2 /air detonation (Sec.III B) and PRF/air hotspot detonation (Sec.III C 3) cases.The primary objective is to present the distinct execution times of individual stages within the CFD solution under specified case conditions.Furthermore, we compare the ARCFoam's performance between moderately (H 2 case) and highly (PRF case) computationally demanding simulations.Finally, we showcase the solver's scalability for the massively parallelized simulations of PRF hotspot ignition in Sec.III C.
The performance is reported in the late stages of simulation for both H 2 and PRF mixtures, obtaining % 9 and 90 Â 10 6 cell counts, respectively.It is important to highlight that the simulation parameters for the H 2 case, including refineInterval, maxImbalance, and N cpp , were adjusted to align with the PRF case for the comparison.Consequently, the performance depicted here does not correspond to that of the H 2 /air simulations presented in the results, with details in Sec.II E 1.All the simulations in this study were conducted on Mahti supercomputer, provided by the CSC-the Finnish IT Center for Science.Mahti features 1404 nodes, each equipped with two 64-core AMD EPYC (Rome) processors running at 2.6 GHz.The nodes are interconnected with a high-speed 200 Gbps infiniband HDR network, and the system operates on the RHEL 7.8 Linux operating system. 65able III presents the configurations and results of the test simulations for both the H 2 /air direct detonation (H 2 test case) and the PRF/air hotspot detonation (PRF test case).The provided configurations include refineInterval, redistributeInterval, maxImbalance, the dynamic distribution algorithm (distributor), number of iterations (N itter ), N cpp , number of processors (N cores ), approximate number of cells (N cells ), and number of species (N sp ).Based on the N cores , the H 2 and PRF test cases can be classified as medium-scale and massively parallel simulations, respectively. 66he H 2 case exhibits significantly shorter execution times compared to the PRF case, as indicated by the test results in Table III.8][69] Therefore, the varying sizes of the utilized H 2 (Ref.58) (ten species) and PRF 64 (115 species) mechanisms result in significant differences in the chemistry solution time per computational cell for these cases.The efficiency of the utilized chemistry solver, DLBFoam v.1.1, 37has been demonstrated in terms of scaling performance by Tekgul et al. 33 Its effectiveness and speed-up have been documented across a range of 2D/3D reacting simulation cases, employing diverse levels of parallelization with processor counts ranging from 8 to 1920. 37he test results reveal a noteworthy difference in the number of refinement (N refine ) and balancing (N balance ) executions.However, in both examined test cases, the respective durations for AMR and mesh balancing processes are similar.It is important to note that the AMRdriven mesh addition is influenced by the pace at which the underlying physical phenomenon and its corresponding refinement region evolve within the computational domain.The similarity indicates an even division of mesh pre-processing time between these two stages.Figure 5 presents pie charts illustrating the cumulative percentage of various operations throughout the total execution time, as detailed in Table III.
The pie chart for the PRF test case illustrates that nearly all (95%) computational time is allocated to the chemistry update operation throughout the simulation.Notably, in this case, the extra AMR and mesh balancing stages do not result in substantial processor intercommunication within the massively parallelized setup.As a result, the speedup achieved in the PRF test case is comparable with the previously reported effective performance of DLBfoam. 33,37Nevertheless, the remaining 5% of the computation is distributed across other processes, with 60% attributed to flow solution tasks (N-S solve and flux calc), 14% dedicated to boundary condition corrections, and the remaining 26% divided between the AMR and mesh balancing processes.
Given the same N cpp for both cases, the H 2 case exhibits the highest percentage of computations (35.1%) allocated to BC corrections, surpassing the flow solution time (N-S solve þ flux calc ¼ 32.5%).This significant contribution underscores the overhead generated by excessive processor boundaries. 70Notably, the chemistry update accounts for only 7.8% of the total execution time, implying that the optimal number of processors for the H 2 test case is less than 500.Furthermore, it is worth noting that excluding the chemistry time, the cumulative percentage of AMR and mesh balancing processes constitutes 26.6% of the remaining processes.This closely aligns with the value reported for the PRF test case, illustrating the scalability of these processes within the framework of ARCFoam.

B. Direct initiation of hydrogen detonation in a cylindrical configuration
Next, the ARCFoam performance is demonstrated by simulating the direct initiation of detonation in a hydrogen/air mixture using both 1D planar and 2D cylindrical configurations.These simulations serve to showcase the solver capabilities in scenarios resembling 2D hotspot simulations, where the ignition front exhibits curvature, influencing the front propagation.Prior to the PRF simulations, we establish the effectiveness of the numerical approach by first demonstrating credible outcomes for a simpler fuel system comprising a smaller number of species (H 2 /air).
The simulation starts with the initiation of a blast wave caused by a high-energy deposit.Next, transverse waves (TWs) emerge as a result of the sinusoidal perturbation superimposed on the initial energy source.Such perturbation is applied to accelerate the development of cellular instability in the present case.Shortly after the blast wave is generated, a decaying over-driven detonation starts to propagate in the radial direction at speeds higher than U cj .Figure 6(a) illustrates the trajectories of the normalized pressure gradient, which represents the evolution of the detonation cellular structure in the current simulation.According to Fig. 6(a), the detonation front propagation consists of two main stages: 60,71 (I) the no-cell stage, where the overdriven detonation initially decays, and (II) the cell-growing stage (x/D i > 260), where the front speed decays to CJ values.In this stage, TWs, normal to the direction of propagation, appear as a result of the cellular instability, leading to the growth of cells as the detonation front expands radially.
Figures 6(b) and 6(c) superimpose the numerical Schlieren (density gradient) and temperature profiles at different time instances: t 1 ¼ 0.000 155 s, t 2 ¼ 0.000 49 s, t 3 ¼ 0.000 825 s, t 4 ¼ 0.001 165 s.Furthermore, the development of RM instabilities can be observed in the lower-left corner of Figs.6(b) and 6(c) at time instances t 1 and t 2 .These instabilities are commonly observed at the interface separating two fluids of different densities subjected to impulsive acceleration as occurs when a shock wave travels through the interface. 72In the present work, RM instabilities emerge at the interface of the initial high-pressure and temperature deposit with the ambient.The three distinct mushroom-shaped structures observed in the figures are a result of the sinusoidal perturbation applied to the initial energy source.
The zoomed-in structure of the detonation front at t ¼ t 4 is illustrated in Figs.6(d The scatterplot of P Àq À1 for the cylindrical detonation and the 1D planar profile at the front location (t ¼ t 4 ) are presented in Fig. 7(a).4][75] Notably, the scatter data for the 2D cylindrical detonation are collected from points near the detonation front, while in the 1D case, a linear profile shows the P Àq À1 variations across the detonation front.As depicted in Fig. 7(a), while the 1D profile intersects the CJ and ZND (von Neumann) states, signifying stable propagation, certain scatter points from the 2D simulations are positioned above the von Neumann state and below the CJ state.These disparities suggest the presence of over-drive and active particles downstream of the detonation front, respectively.The former observation can be attributed to particles located near strong MSs and TPs.Meanwhile, the latter observation may be attributed to delayed reactions occurring behind weaker incident shocks on the detonation front.In such scenarios, the establishment of equilibrium is delayed compared to more intense reactions behind the MSs.
The front propagation speeds against the propagation distance normalized by D i are displayed in Fig. 7(b) for both 1D and 2D setups.In the 2D cylindrical configuration, the front speed is derived by monitoring the maximum pressure at various time intervals.These measurements are acquired from the bottom (symmetric) boundary, revealing fluctuations attributed to cellular instability.Figure 7(b) depicts that the initial phase of both 1D planar and 2D cylindrical detonations exhibits a substantial over-drive caused by the high energy deposition at the start.Subsequently, it is noteworthy that the cylindrical case decays rapidly toward the CJ value, maintaining an average propagation speed below the CJ threshold.In contrast, the one- dimensional planar detonation decays at a comparatively slower rate and stabilizes at the average CJ propagation speed, consistent with previous observations. 60As a remark, this difference may be due to the fact that the 1D simulations represent a planar configuration restricting the radial decay of pressure intensity.
The observed deficit of the average propagation speed of the cylindrical detonation from the CJ value during later stages of propagation (x/D i > 400) can be attributed to the curvature effects, 60,61,76 which enhance instability and prolong the induction times at the front location.It should be noted that achieving the theoretical CJ value in the cylindrical case would require an infinitely large radius.Additionally, the scatter data in Fig. 7(b) demonstrate irregular fluctuations in the front speed, which can be correlated with the inhomogeneity of the cellular structure as shown in Fig. 6(a).Nevertheless, it is evident from Fig. 7(b) that the majority of scatter points beyond the normalized distance of x/D i ¼ 400 lie within the bounds of 0.7-1.6 times the CJ speed, indicating sustained propagation of the detonation front. 74

C. 2D hotspot ignition of PRF mixture
Section III C provides the results of 2D simulations investigating hotspot ignition with different propagation regimes.As mentioned, "hotspot ignition" refers to simulations that are initiated by a linear temperature gradient (rT) and radius of the hotspot (r h ).Building upon our previous work, 13 which focused on 1D ignition regime diagrams (Bradley peninsula) for gasoline surrogates (PRF and PRF-E), here we examine distinct propagation patterns arising from the NTC behavior of the PRF mixture in 2D context.These patterns include off-centered ignition and the emergence of secondary ignition kernels ahead of the main front, as detailed in Sec.II B.

Revisiting hotspot ignition regime diagram (Bradley peninsula)
A reproduction of various ignition regimes obtained from our 1D simulations 13 is illustrated in Fig. 8 at a relatively high initial pressure of P ¼ 50 bar and an average temperature within the NTC range, i.e., T ave ¼ 900 K.The blue and pale blue zones correspond to subsupersonic ignition, while the red zone denotes the detonation regime.These zones are outlined based on hotspot pressure levels, dividing the map into regions of relatively high and low pressures using a predefined threshold.In our earlier study, 13 two different approaches were tested: one comparing the average hotspot pressure to 1.5 times the constant-volume reactor pressure (P cv ), while the other comparing the maximum pressure to the Chapman-Jouguet pressure (P cj ).Remarkably, both methodologies yielded similar regime diagrams.Subsequently, the areas surpassing the threshold are identified as detonation zones, while the remaining areas are divided between supersonic and subsonic ignition regions.In Fig. 8, the bottom-left blue zone represents supersonic ignition, while the top-right blue zone corresponds to subsonic propagation.Notably, the reproduced detonation zone (red region) does not align precisely with the original Bradley peninsula indicated by the black dashed line. 10This discrepancy can be primarily attributed to the impact of NTC on the propagation of the ignition front, as discussed in our previous study. 13Additionally, the green and purple dashdotted lines mark regions characterized by off-centered ignition and coolspot ignition.However, it is important to acknowledge that the boundaries of the different regimes within the reproduced detonation peninsula are sensitive to the choice of fuel 15,[77][78][79][80] and kinetic mechanism, 81,82 as indicated by previous research.
To investigate different propagation patterns (see Sec. II B) in a 2D cylindrical context, three cases were selected based on the ignition regime diagram, as depicted in Fig. 8. Point 1 is positioned in the subsonic ignition propagation region, close to the off-centered ignition and detonation zone predicted from the 1D simulations.Point 2 corresponds to the detonation peninsula, where off-centered ignition is expected due to NTC chemistry.Finally, point 3 is located within the supersonic ignition propagation region, where ignition is initiated by a coolspot.Sections III C 2-III C 4 will provide a detailed discussion of the 1D planar and 2D cylindrical simulation results with hotspot initialization corresponding to points 1-3.

Analysis of subsonic ignition front propagation (point 1)
In this section, we present the results of both 1D planar and 2D cylindrical simulations conducted at point 1, as indicated in Fig. 8.The simulation considers a hotspot with a radius of r h ¼ 0:012 62 m (e % 25) and an initial temperature gradient inside the hotspot of rT ¼ 15:8 K/mm (n % 40).Based on Fig. 8, this point is located in the subsonic ignition region, close to the off-centered ignition and the detonation zone boundaries.
The transient profiles of temperature and normalized pressure (P/P z ) from 1D planar simulation are depicted in Figs.9(a) and 9(b), respectively.The profiles are displayed at different time intervals, with time instances 1-6 plotted in 0.05-ms intervals and time instances 7-15 in 0.1-ls intervals.It can be observed that after the primary ignition at times 0-1 in Fig. 9(a), a subsonic ignition front (S ign % 21:2 m/s) emerges from the center of the hotspot (times 2-6).However, as it progresses beyond r =r h % 0:55, a rapid secondary ignition occurs at time 7, leading to the appearance of secondary ignition fronts propagating to the left and right.
The rapid secondary ignition leads to an increase in the pressure level up to von Neumann pressure (P z ) inside the hotspot, Fig. 9(b).The left ignition front collapses with the primary subsonic front, while the right propagating front evolves into a detonation front outside the borders of the hotspot.These results indicate that although the primary ignition follows the expected behavior based on the point 1 location on the map, Fig. 8, the presence of NTC can alter the dynamics of front propagation both inside and outside the hotspot, resulting in elevated pressure.Regarding 2D simulations, the temporal evolution of temperature (top row) and pressure (bottom row) in the cylindrical simulation of point 1 is depicted in Fig. 10.Initially, the primary ignition takes place at the center of the hotspot [Fig.10(a)], giving rise to an ignition front that propagates at a subsonic speed (primary ignition front, hereafter called PIF).As the burnt gas expands behind the ignition front, it induces outward propagating pressure waves [Fig.10(b)], which accelerate the low-density burnt products toward the unburnt gas and trigger RM instabilities on the edge of the PIF [Figs.10(a) and 10(c)].Furthermore, at t ¼ 0.001 44 s, a rapid secondary ignition occurs, resulting in the formation of converging and diverging ignition fronts, (SIF 2 and SIF 1 , where SIF stands for secondary ignition front) and a subsequent increase in pressure at the same location (indicated by PW 2 and PW 1 ), as shown in Figs.10(c) and 10(d).
At t ¼ 0.001 44 s, as depicted in Figs.10(e) and 10(f), the pressure waves (PW1 and PW2) can be observed to decouple from the SIFs.PW1 propagates outward in the radial direction, while PW2 interacts with PIF and fastens the growth of RM instabilities on the front, as shown in Fig. 11(e).Moreover, Figs.10(e) and 10(f) demonstrate the presence of relatively high pressures (140 bar) and temperatures (3300 K) behind the imploding pressure wave (PW2), which arises from the transverse reflected shock waves propagating in tangential direction caused by the collapse of PW2. 83,84ubsequently, the converging SIF 2 is merged with the PIF, while SIF 1 continues to propagate with a subsonic speed.Furthermore, PW 2 reflects from the center, called RPW 2 , and follows PW 1 , as shown in Fig. 10(h).In the later stages, when PW 1 and RPW 2 reach the top and right wall boundaries, they reflect and cause an increase in the ambient temperature to 1200 K, as illustrated in Fig. 10(i).However, a collision occurs between the reflection of PW 1 from the right and top boundaries (RPW 1 ) and PW 1 .This is followed by another collision between the reflection of RPW 2 from the same boundaries (R 2 PW 2 ) and RPW 2 at the upper-right corner, leading to a temperature increase to  approximately 1400 K.This elevated temperature level then triggers a third ignition kernel at the upper-right corner.
The temperature (top row) and pressure (bottom row) profiles at the following stage of combustion for point 1 are depicted in Fig. 11.It is notable that the maximum limit of pressure contours is increased to an order of magnitude higher values than in Fig. 10, and hence, the details of the complex system of shock waves are not observable in these sub-figures.
Following the occurrence of the third ignition, a detonation front emerges and further develops from the top-right corner, shown in Figs.11(a demonstrate the subsequent fourth and fifth ignitions at the bottomright and top-left corners, resulting in increased pressure in these regions.Eventually, the newly emerged ignition points, along with the subsonic ignition and the detonation front, lead to the complete consumption of the remaining unreacted gas inside the chamber through a volumetric ignition process. In summary, comparing the outcomes of the 1D and 2D simulations, the following observations can be made.First, both scenarios exhibit a primary subsonic ignition front propagation, aligning with the anticipated combustion mode indicated in the reconstructed regime diagram (Fig. 8) at point 1.Second, both simulations exhibit a secondary ignition event occurring ahead of the primary ignition kernel due to the NTC effect (Pattern III in Sec.II B).However, notable distinctions emerge after the secondary ignition between the 1D and 2D results.In the 1D planar case, the strong secondary ignition develops into a detonation front.In contrast, within the 2D cylindrical case, the secondary ignition transforms into a shock front accompanied by a subsonic ignition front.This deviation is attributed to the curvature of the shock front in the cylindrical setup, which induces additional gas expansion behind it.Consequently, this expansion prompts cooling and decay of chemical reactions, ultimately preventing the sustained development of detonation. 61

Analysis of detonation front propagation (point 2)
The current case presents the results of 1D and 2D simulations conducted for point 2 on the detonation peninsula, as shown in Fig. 8.At this specific location, the hotspot radius is r h ¼ 0:015 144 m (with e % 30), and the initial temperature gradient is rT ¼ 6:33 K/mm (with n % 16).It is noteworthy that this point is located in the coolspot and off-centered ignition region of the regime diagram.
Transient profiles of temperature and pressure from the 1D planar simulation, plotted against the normalized radius (r/r h ) of the hotspot, are shown in Fig. 13.The profiles are presented at various time intervals, where each interval from times 3 to 13 corresponds to a duration of 0.1 ls, illustrating the evolution of detonation front propagation.As observed in Fig. 13(a), the initial ignition kernel is located off-center at r =r h % 0:7, indicating the off-centered ignition effect of NTC chemistry (times 0-2).The off-centered ignition leads to left and right propagating ignition fronts, as depicted in Fig. 13(a).In Fig. 13(b), it can be observed that the rapid energy release from the ignition leads to an increase in pressure.Subsequently, the ignition front develops into a left/right propagating detonation fronts, and the pressure reaches the von Neumann pressure.It is noteworthy that another ignition kernel is observed at the temperature profiles at the center of the hotspot (r/r h ¼ 0); however, it does not significantly affect the dynamics of the ignition fronts and is not discussed further in this case.
Similarly, in the 2D cylindrical detonation case, the primary ignition occurs at a certain distance (r/r h % 0.7) from the center of the hotspot, leading to inward/outward propagating, i.e., converging/ diverging detonation fronts (CDF/DDF).Figure 14 illustrates the temperature profiles obtained from the cylindrical simulation, at t ¼ 0.001 503 s, showcasing the DDF located at r=r h % 1.18 and the CDF located at r=r h % 0.13.As the detonation fronts advance, the curvature of the DDF decreases, while the curvature of the CDF increases.The change in curvature is associated with the variation in the leading shock area for both CDF and DDF.The intensified compression resulting from the reduced shock area accelerates the CDF, leading to additional amplification of the post-detonation pressure and temperature (Fig. 14).
Providing a closer look at the temperature profiles for the CDF and DDF, Figs.14(b) and 14(c) depict a zoomed view of these profiles, respectively.These profiles reveal the presence of high local temperatures at the locations of the TPs.Characterized by the temperature gradients, Figs.Magnified profiles of the logarithm (log 10 ) of heat release rate [log 10 (HRR)] and the pressure gradient normalized with its maximum value (rP n ) are plotted in Fig. 15.The color bands are limited between the maximum (10 18 ) and minimum (10 13 ) HRR for better clarity and comparison.Moreover, the zones with negative heat release are colored in white in Figs.15(b) and 15(d).
Compared to DDF, the CDF exhibits a higher number of transverse waves per unit area, Fig. 15 bottom row, leading to a decrease in the scale of front cellular instability.On the other hand, DDF enters the cell-growing stage, where the detonation cell irregularity gradually increases due to reduced curvature.A closer examination of the front structure of the zoomed-in heat release rate and pressure gradient profiles in Fig. 15 reveals the prompt occurrence of negative HRR behind the MSs, ISs, and TWs.Based on the reactants' species profiles of our simulations (not shown here for brevity), the negative HRR values behind the MSs and ISs are attributed to the decomposition of reactants (iso-octane and n-heptane).However, in the absence of the reactants, the negative HRR values at the location of TWs primarily may arise from the decomposition of intermediate species.This finding is in alignment with the outcomes of shock tube experiments and kinetics simulations conducted by He et al., 85 which investigated various gasoline surrogates including n-heptane and iso-octane.These studies demonstrated the occurrence of endothermic reactions during the decomposition of reactants and potential intermediate species behind the reflected shock wave.
Furthermore, a gradual decay in chemical reactivity [13 < log 10 (HRR) < 18] behind both CDF and DDF is observed in Fig. 15.However, a comparison between Figs. 15(a) and 15(c) indicates that the chemical relaxation extends further downstream for the CDF compared to the DDF.This observation suggests that subsequently the mechanical and thermal relaxation (dissipation of hydrodynamic perturbations 73 ) behind the CDF may be influenced by the curvature.To further investigate this phenomenon, the Mach number in the detonation frame of reference [Mach re ¼ ðu À DÞ=U a , where D is the detonation speed] has been calculated and plotted together with O 2 mass fraction profiles (Y O2 ) in Fig. 16.
The sonic plane (Mach re ¼ 1) is outlined by the white contour lines in Figs.16(f) and 16(h).It can be observed that the sonic plane is shifted further downstream for the CDF compared to the DDF due to the influence of the increasing front curvature.This shift indicates a prolonged mechanical and thermal relaxation caused by the presence of TWs.A similar trend was observed by Han et al. 60 when comparing planar and cylindrical detonation fronts.Furthermore, the white contour lines in Figs.16(b) and 16(d) represent the boundary extending from the edge of the detonation fronts to one percent of the O2 massfraction in the reactants (Y O2 ¼ 0.022), highlighting the reaction thickness of the detonation fronts.The reaction zone of the CDF is visibly thinner than that of the DDF, which contrasts the trends observed in log 10 (HRR) profiles.The difference in reaction zone thickness, in Figs.16(b) and 16(d), can be linked to the elongated induction region behind the ISs within the DDF.These findings indicate that the CDF is characterized by a short induction region featuring a high rate of reactions followed by a long tail of gradually decaying reactions.In contrast, the DDF exhibits an elongated induction zone, albeit with relatively rapid chemical relaxation.To gain further insights into the structure of detonation fronts, scattered data of log 10 (HRR) at locations where 14 < log 10 (HRR) < 18 are plotted against temperature in Fig. 17.
The scatter plots in Figs. 17  CDF and DDF's conditional averages indicates a relatively higher HRR at the low-temperature peak for the CDF corresponding to a narrower distribution band and a smaller number of points.Additionally, it can be observed that the scatter points cover the entire temperature range in Fig. 17(a), but a more clear distinction between low-temperature and high-temperature regions is evident in Fig. 17(b).This can be attributed to the regular cellular structure of the converging front, while the outward wave exhibits a broader range of cellular instabilities.Scatter plots of OH mass fractions (Y OH ) vs temperature for DDF and CDF are displayed in Fig. 18(a).It is observed that the scatter points for the CDF are concentrated at higher temperatures, indicating rapid combustion of reactants with shorter IDTs.On the other hand, the scatter points for the DDF exhibit slower decay rates and cover a wider temperature range, suggesting slower reaction progress.Moreover, Fig. 18(b) presents a scatterplot of p À q À1 for the cylindrical detonation, alongside the 1D planar profiles extracted from the results at the same normalized radial positions (r =r h ) as the 2D simulations, corresponding to the DDF and CDF.It is observed that both 2D cylindrical and 1D planar DDF and CDF profiles pass through the CJ point indicating a sustained front propagation.However, upon examining the 2D profiles, it is evident that the scatter points for the DDF are extended to relatively higher pressures and lower densities, indicating a larger amplitude of fluctuations (higher instability) for the DDF front.
Finally, the spatial evolution of the front speed for the 1D planar and 2D cylindrical CDF and DDF is depicted in Figs.19(a) and 19(b), respectively.The front speed is determined by tracking the maximum pressure at different time instances.In the 2D simulations, these measurements are taken on the bottom (symmetric) boundary, which demonstrates fluctuations due to the cellular instability.The profiles show an initial burst of relatively high speed, followed by a developing stage with speeds below the CJ value.After traveling a run-up distance of approximately 0.25 times the hotspot radius (Dr/r h % 0.25), both the DDF and CDF speeds approach the CJ value.Furthermore, while the 1D planar profiles do not exhibit notable fluctuations, the scatter plots for both the 2D CDF and DDF show low-and high-amplitude fluctuations, respectively.This difference in fluctuations aligns with our previous observations on the instability levels of the CDF and DDF.Nonetheless, the spatial average speed for both fronts converges to the CJ value, indicating sustained propagation.
The averaged CDF speed reveals that the CDF does not become over-driven until the later stages of propagation near the center (r/r h % 0.05), Fig. 19(b).Initially, the detonation is primarily sustained by the chemical heat release, as shown in Fig. 19(b).However, as the convergence progresses, the dominant mechanism for pressure rise shifts from chemical reaction heat to compression due to the shock area convergence.In this stage, the pressure rise resembles that of a shock wave. 86,87t is notable that the chemical mechanism used in this study 64 and its associated thermodynamic properties (such as c p , enthalpy, and entropy) are validated and applicable for temperatures below 5000 K. Additionally, when solving the ordinary differential equations (ODEs) for reactions near the hotspot center, extremely high temperatures were encountered, as previously reported in a study by Jiang et al. 87 To prevent numerical issues, the reactions were suppressed when the temperature exceeded 4500 K.It is important to note that such high-temperature values are only reached in the final stages of the CDF convergence and they are not expected to significantly impact the DDF propagation.
Furthermore, the CDF undergoes reflection from the center (bottom-left corner), resulting in the propagation of a reflected nonreacting shock wave (RSW) behind the DDF.The DDF, on the other hand, continues to propagate with an average equilibrium (CJ) speed until reaching r/r h % 2, Fig. 19(b).At this point, the unburnt mixture ahead of the DDF, which initially exhibits relatively high temperatures [% 1200 K, as shown in Fig. 14(a)], undergoes spontaneous ignition and begins to burn the remaining mixture ahead of the DDF.More detailed information regarding the temperature contours, mass fractions, relative Mach number, and pressure gradient, which provide insights into the characteristics of the DDF at this stage, can be found in Appendix B.
In summary, the comparison of 1D planar and 2D cylindrical simulations consistently demonstrated off-centered ignition (Pattern II in Sec.II B) leading to detonation for both configurations.This phenomenon, attributed to the NTC effect, aligns with the predicted combustion mode in the reconstructed regime diagram (Fig. 8) at point 2.Moreover, apart from the inherent differences between the 1D and 2D setups, including curvature effects, it was observed that the dynamics of CDF were notably influenced by the reduced shock area.This reduction led to higher propagation speeds, pressures, and temperatures, underscoring the pivotal role of curvature considerations in predicting temperature/pressure levels, especially at point 2.

Analysis of supersonic ignition front propagation (point 3)
Section III C 4 presents the results of 1D and 2D simulations conducted at point 3 on the detonation peninsula, as indicated in Fig. 8.At this specific location, the hotspot radius is r h ¼ 0:012 62 m (with e % 25), and the initial temperature gradient is rT ¼ 0:4 K/mm (with n % 1).It is noteworthy that this point is initiated as a coolspot within supersonic ignition zone based on the 1D simulation results as depicted in the regime diagram.
Figure 20(a) illustrates temperature profiles obtained from the 1D planar simulation at point 3.The profiles are presented in 0.2-ls time intervals, showcasing three distinct stages of the ignition process.The initial volumetric ignition (0-1) occurs with a rapid temperature rise and an approximate propagation speed of 20 000 m/s.During the front formation stage (1-4), the temperature profile exhibits a steeper gradient, and the ignition front begins to form with an average speed of approximately 10 000 m/s.Finally, the process enters the final volumetric ignition stage (4-9), where the ignition occurs uniformly throughout the unburnt mixture with an average speed of approximately 27 000 m/s.The temperature profiles indicate that the minimum temperature reached during the ignition process is approximately 1500 K, inducing rapid ignition propagation.Moreover, the pressure contours depicted in Fig. 20(b) demonstrate that the pressure at the burnt product is approximately equal to the constant volume pressure (P cv % 187 bar).This observation suggests a volumetric ignition process occurring without significant pressure amplification, as there is no pronounced coupling between the ignition front and the pressure perturbations induced by the expansion of burnt gases.For the present case, both the 1D planar and 2D cylindrical simulations are initiated as coolspots due to the influence of the NTC chemistry, a behavior consistent with the predicted combustion mode indicated in the reconstructed regime diagram (Fig. 8) at point 3.The coolspot ignition, in both the 1D and 2D simulations, exhibited a conventional propagation pattern (Pattern I in Fig. 1).Furthermore, both the 1D and 2D configurations demonstrated three primary stages: initial volumetric ignition, formation of a supersonic front, and final volumetric ignition.The comparison with the predictions from the regime diagram suggests that, when complex shock reaction front systems and their interactions due to NTC behavior (i.e., off-centered ignition or multiple ignition kernels) are absent, the Bradley peninsula approach performs well in predicting ignition front propagation and associated front characteristics.

IV. SUMMARY AND DISCUSSION
The results from the 1D planar and 2D cylindrical hotspot simulations for points 1-3 are summarized in Table IV.This table encompasses the initial ignition regime, the impact of the NTC chemistry on the outcome and initiation of the simulations, along with the main observations drawn from these simulations.
In accordance with the results provided in Secs.III C 2-III C 4, as well as the summarized data in Table IV, the following observations can be made:  (1) Both 1D planar and 2D cylindrical hotspot simulations demonstrate that the theoretical regime diagram can successfully predict the dynamics of ignition front propagation inside the hotspot, i.e., the initial phases of the chemically reactive flow.However, certain events, such as the induced pressure waves in the subsonic mode and the NTC behavior, have the potential to significantly alter the thermochemical conditions of the mixture, and consequently, impact the characteristics of the ignition fronts.
(2) The study reveals that in both 1D planar and 2D cylindrical hotspot scenarios involving detonation, the process initiates with an initial burst (spontaneous ignition).This is followed by a developing stage with reaction front speeds lower than the CJ value and eventually transitions to sustained detonation propagation.Nevertheless, notable differences are observed between these configurations.These disparities include the distinct critical ignition energies required for initiating the detonation, the presence of cellular instabilities, and the effect of curvature on cellular instability and the reactions occurring within the front in the 2D configuration.These factors contribute to the dynamic behavior and spatial distribution of the detonation front in the 2D cylindrical configuration, leading to variations in the detonation propagation compared to the 1D planar scenario.The reflection of these diverging pressure waves can lead to localized temperature and pressure increases close to the walls, resulting in subsequent strong ignition and detonation.Such events in realistic engine scenarios are known as shock wave reflection-induced detonation (SWRID), which is the most commonly occurring detonation initiation in spark ignition engines. 1 This finding highlights the critical importance of the subsonic regime on the detonation peninsula, as it may significantly influence the occurrence of SWRID and, consequently promote the deflagration to detonation transition phenomenon.

V. CONCLUSIONS
The present work continues the authors' previous research on the effect of thermal stratification on secondary combustion modes initiated in end gas under SI engine-relevant conditions.Here, we study the different ignition regimes and their impact on pressure oscillation amplitudes in the end-gas region.The work introduces ARCFoam, a 2D numerical framework based on OpenFOAM to capture the interaction between fuel chemistry and compressible flow physics.The framework is first validated and assessed in the 2D simulation of directly initiated detonation of H 2 /air mixture problem.Second, the 1D Bradley regime diagram from our previous work 13 is revisited at initial pressure (50 bar) and average temperature (900 K) for the PRF mixture.Furthermore, three points were selected for 2D cylindrical hotspot simulations in order to better understand the differences between 2D and 1D simulations.The following conclusions are made: (1) ARCFoam's methodologies, including adaptive mesh refinement, dynamic mesh balancing, and chemistry balancing, are crucial for effective 2D simulations.These techniques reduce cell counts through localized refinement and balance the computational load across the processors in parallel reacting CFD simulations.For less computationally demanding simulations (e.g., H 2 case, 9 M cells, ten-species chemistry), ARCFoam improves the computational efficiency.However, for resourceintensive simulations with intricate chemistry (e.g., PRF cases, 90 M cells, 115-species mechanism), ARCFoam is considered to be a key methodological enabler.(2) Notable differences exist between 1D and 2D hotspot ignition scenarios, such as the emergence of the Richtmyer-Meshkov instability due to shock wave-induced impulsive acceleration at the ignition front interface.Moreover, differences in critical initiation energies and the lack of curvature effects in the 1D setting result in distinct detonation initiation and development characteristics.Although average detonation speeds may align, the occurrence of 2D cellular instabilities introduces considerable fluctuations in the transient detonation front structures.(3) The 1D ignition regime diagram is proficient in predicting initial ignition propagation within the hotspot, yet the influence of factors like curvature, confinement, and NTC effects can modify the behavior of the primary front and induce secondary ignition modes both within and outside the hotspot.(4) NTC chemistry plays a key role in initiating and transitioning ignition modes.Both 1D and 2D simulations highlight NTC chemistry's impact on secondary ignitions within and outside the hotspot, leading to changes in combustion regimes and the emergence of phenomena like shock wave reflections and converging shock/detonation fronts, which may significantly alter engine cylinder combustion dynamics.(5) In confined hotspot ignition scenarios, subsonic ignition is pivotal due to potential shock wave reflections.In cases where ignition energy is insufficient for direct detonation initiation a decoupled shock wave reflects off walls, elevating pressures and temperatures, leading to a new ignition front near the wall that may result in shock wave reflection-induced detonation.
experimental and numerical studies. 43,62,63,88The present simulation serves as a showcase to validate the ARCFoam numerical methods, which are coupled with the load-balanced AMR and DLBFoam techniques.In addition, a mesh sensitivity study is conducted to find the required mesh resolution for the simulation in Sec.III B.
The schematic representation of the 2D channel, along with the specified boundary conditions (BCs), is depicted in Fig. 22(a).Here, the detonation is initiated through the use of three equidistant highenergy sources located at the left end of the channel.The essential dimensions of the channel, namely, its length (L) and height (h), are reported in Table V.Moreover, the source and ambient temperatures and pressures denoted as T s , P s , T 0 , and P 0 , respectively, are also listed in Table V.
For the chemistry solution, the Marinov et al. 58 H 2 mechanism is utilized in the simulations.To ensure accurate results, the grid resolution is adjusted in relation to the ZND induction length (D i ¼ 3 Â 10 -3 m), detailed in Sec.II E 1.The base grid size is set at Dx; Dy % 0.0015 m, and during the run-time, the AMR is applied with additional 2, 3, and 4 levels of refinement.This choice of refinement levels corresponds to having 6, 12, and 24 cells per D i , respectively.Furthermore, to maintain numerical stability, an acoustic Courant number [Co ¼ ðjUjþjajÞDt Dx ] of approximately Co ¼ 0.02 is observed to be necessary.slightly over-predicts the detonation cell length, while the results obtained with 12 and 24 cells/D i meshes show good agreement with each other.This finding indicates that the resolution of D i =12 is sufficient to capture the H 2 cellular instability using the current chemical mechanism and under the considered conditions of P 0 and T 0 .
Furthermore, to validate the results, the detonation cell aspect ratio (k=a) is estimated for detonation cell I, which is indicated by the red-dashed box in Fig. 22(b).The detonation cell aspect ratio is the ratio between the height (k) and the cell length (a), as shown in Fig. 22(b).This aspect ratio is then compared with the results from numerical studies that used a similar case configuration. 43,62,63The comparison shows that the present cell aspect ratio (approximately 0.43) is in good agreement with the reported values (ranging from 0.42 to 0.43) from these previous studies.The consistency between the present results and the literature values further supports the validity and capability of the utilized numerical methods in accurately simulating 2D detonation in the hydrogen/air mixture.

APPENDIX B: THE PRF POINT 2 LATE-STAGE IMAGES
The late-stage profiles of point 2 in the cylindrical hotspot ignition simulation are shown in Figs.23-25.The temperature and pressure profiles in Fig. 23 reveal the state of the diverging detonation front (DDF), which reaches a radial distance of r ¼ 0.03 m (r =r h % 2), while the reflection of the converging detonation front from the center (RSW) propagates behind the DDF.In Figs.23(c) and 23(f), a zoomed-in view of the temperature and pressure profiles is provided, emphasizing the presence of strong and weak triple points (TPs) with non-uniform distance, resembling irregular cellular structures.Additionally, the temperature distribution illustrates a relatively wide induction zone behind the longer IS in the top left, while the induction zones behind the other ISs in the lower right are shortened due to the presence of multiple TPs.
The OH and O 2 mass fraction profiles (Y O2 and Y OH ) are presented in Fig. 24.The reaction thickness, indicated by the location where Y O2 ¼ 0.022, corresponds to approximately 1 Â 10 -4 m.A comparison between the Y O2 and Y OH profiles, shown in Figs.24(c) and 24(f), reveals a sharp decrease in O 2 mass fraction and a prominent OH peak at the location of the MS, indicating the hightemperature region of active combustion.On the other hand, behind the IS, a gradual decrease Y O2 is observed, with delayed OH production in the reaction zone.Furthermore, both Y O2 and Y OH profiles demonstrate that the mixture shocked at the front is almost entirely consumed by the diverging detonation front (DDF), leaving only a small amount of unreacted gas in the slipping layers behind the front.These regions are further burnt with slower rates behind the TWs.
The Mach re and log 10 (HRR) profiles are presented, in Fig. 25.The sonic plane, indicated by white contour lines in Fig. 25(c), represents the hydrodynamic distance of approximately 1 Â 10 -4 m, which corresponds to the estimated thickness of the reaction layer.In Fig. 25(f), the zones with negative heat release are colored in white.Examining the log 10 (HRR) profile, it is observed that negative values occur behind the TWs, indicating the occurrence of endothermic reactions (cooling effect) involving the intermediate species breakup.

FIG. 1 .
FIG. 1. IDT vs temperature diagram of PRF at 900 K and 50 bar, and schematic presentation of the propagation patterns indicating: (a) different temperature zones (I-III) associated with various propagation patterns (I-III), (b) hotspot/coolspot initial temperature profile, (c) centered ignition and regular propagation inside the hotspot/coolspot (zones I h / I c ), (d) off-centered primary ignition due to NTC and double ignition front propagation (zone II), and (e) occurrence of primary ignition and subsequent formation of the secondary ignition kernel due to NTC (zone III).

FIG. 4 .
FIG. 4. Schematic demonstration of (a) visible cells in the two refinement regions, (w 1 , w 2 ) and (w 0 1 ; w 0 2 ) indicate local weights (number of cells) of each refinement region, (b) different cell levels (1 and 2) on each refinement region and their projection on parent cells (level 0), (c) mesh and temperature distribution of the Richtmyer-Meshkov instability in H 2 detonation with five levels of refinement.
), 6(e), and 6(f), which show the distribution of OH mass fraction, pressure, and temperature within the detonation front, respectively.In particular, incident shock (IS), Mach shock (MS), triple point (TP), and TWs are highlighted in Fig. 6(e).Furthermore, Fig. 6(f) displays the longer induction zone behind the weaker incident shocks, as well as the shear layers (SLs) formed by the slipstream between TWs and MSs.Notably, the negligible OH mass fraction, Fig. 6(d), in the induction zone indicates the extended ignition delay time behind the ISs, while the relatively high values of OH mass fraction behind TWs signify the completion of reactions facilitated by the presence of TWs.

FIG. 5 .
FIG. 5. Pie chart demonstrating the cumulative percentage of operations during the total execution time for H 2 and PRF test cases.This includes processes such as chemistry update, AMR, mesh balancing, flux calculation (flux calc), N-S solution (N-S solve), and boundary condition correction (BC correct).

FIG. 6 .
FIG. 6. 2D contour plots illustrating the characteristics of H 2 cylindrical detonation: (a) cellular structure showing the presence of a no-cell region (I) and a cell growing region (II).(b) and (c) Density gradient and temperature profiles at different time instances: t 1 ¼ 0.000 155 s, t 2 ¼ 0.000 49 s, t 3 ¼ 0.000 825 s, t 4 ¼ 0.001 165 s. (d)-(f) Zoomed-in profiles depicting the OH mass fraction, pressure, and temperature at the specified location within the box.Moreover, Mach shock (MS), incident shock (IS), triple point (TP), transverse waves (TWs), shear layer (SL), and induction zone are annotated on panels (e) and (f).

FIG. 7 .
FIG. 7.Comparison between one/two-dimensional simulations: (a) p À q À1 profile (represented by the red line in 1D) and scatterplot (in 2D) at time t ¼ t 4 , illustrating variations in particle states along the detonation front, along with highlighted CJ and ZND post-shock states, (b) front propagation speeds vs normalized propagation distance, showcasing the dynamics of the ignition front propagation.The Green dashdotted and the gray dashed lines depict U cj and 1.6-0.7Ucj , respectively.

FIG. 9 .
FIG. 9. 1D transient profiles of temperature and normalized pressure at point 1 for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, demonstrating subsonic ignition front propagation and the occurrence of secondary ignition.The profiles are presented at different time intervals: times 1-6 correspond to 0.05-ms duration, and times 7-15 correspond to a duration of 0.1 ls in each interval.Dotted lines continue the profiles intersecting with other time instances for improved clarity.

FIG. 10 .
FIG. 10.Temporal evolution of 2D temperature (top row) and pressure (bottom row) profiles of point 1, [(a) and (b)] before the occurrence of secondary ignition (0.001 435 s), [(c) and (d)] after secondary ignition and creation of outward (PW 1 ) and inward (PW 2 ) propagating pressure waves, [(e) and (f)] detachment of PW 1 and PW 2 from the reaction front, [(g) and (h)] reflection of PW 2 from the bottom left corner (RPW 2 ), and [(i) and (j)] reflection of PW 1 and RPW 2 from the top and right walls, i.e., RPW 1 and R 2 PW 2 waves, respectively.PIF and SIF 1,2 indicate the primary and secondary ignition fronts.
)-11(d).The front propagation speed vs normalized radius is depicted in Fig. 12.It can be observed that the process initiates with spontaneous ignition at r =r h ¼ 3, leading to the development of the detonation and the convergence of the front speed to the CJ speed [approximately 1870 (m/s)] at r =r h % 2. In addition, Figs.11(e)-11(j)

FIG. 11 .
FIG. 11.Continued 2D temperature profiles (top row) and pressure profiles (bottom row) of Point 1 for the PRF mixture at (T ¼ 900 K and P ¼ 50 bar) from the previous time step shown in Fig. 10.The profiles illustrate the formation and propagation of the third, fourth, and fifth ignition fronts.
14(b) and 14(c) exhibit a relatively short induction zone behind the MSs and a longer induction zone behind the ISs.This observation indicates the delay in the occurrence of reactions and heat release behind the IS, resulting in lower temperatures compared to the TP location.

FIG. 13 .
FIG. 13. 1D transient profiles of temperature and normalized pressure at point 2 for the PRF mixture at (T ¼ 900 K and P ¼ 50 bar), demonstrating the appearance of off-centered ignition and subsequent converging/diverging detonation.Time intervals from 3 to 13 represent 0.1 ls each, showcasing detonation front propagation.Dotted lines continue profiles intersecting with other time instances for enhanced clarity.
FIG. 15.Zoomed-in 2D profiles of log 10 (HRR) (top row) and the pressure gradient normalized by its maximum value (rPn) (bottom row) at point 2 for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, illustrating converging and diverging detonation fronts (CDF and DDF) at normalized radial positions of r =r h % 0.13 and 1.18, respectively.The white color in panels (b) and (d) indicates negative HRR values.

FIG. 16 .
FIG. 16.Zoomed-in 2D profiles of O 2 mass fraction (Y O2 ) (top row) and relative Mach number (Mach re ) in the detonation frame of reference (bottom row) at point 2 for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, depicting converging and diverging detonation fronts (CDF and DDF) at normalized radial positions of r =r h % 0.13 and 1.18, respectively.The white contour lines in panels (b) and (d) indicate the reaction zone (Y O2 ¼ 0.022), and in panels (f) and (h) indicate the sonic plane (Mach re ¼ 1).The values shown in panels (b), (d), (f), and (h) serve as representative length scales in the respective zoomed-in profiles, offering a reference for facilitating the interpretation of the plots.

FIG. 17 .
FIG. 17. Scatter plots of log 10 (HRR) vs temperature for (a) diverging and (b) converging detonation fronts.The yellow and red lines depict the conditional average of scatter plots and 1D planar results, respectively.

Figure 21
displays temperature (top row), log 10 (HRR) (middle row), and pressure (bottom row) contours obtained from the 2D cylindrical simulation of point 3. Notably, the temperature contours [Figs.21(a), 21(d), and 21(g)] show the front propagation within the hotspot (r/r h < 1) obtaining a relatively large temperature gradient and a distinct log 10 (HRR) peak [Figs.21(b), 21(e), and 21(h)] at the front location.The pressure profiles in Figs.21(c), 21(f), and 21(i) demonstrate that the pressure behind the front reaches approximately P cv .As the reaction front approaches the hotspot border, the mixture ahead of the front starts reacting, resulting in a decrease in the temperature gradient [Fig.21(j)] and a more distributed heat release pattern [Fig.21(k)].In later stages, the unburnt mixture ahead of the front continues to react, resembling a volumetric reaction process [Figs.21(m)-21(o)].It should be noted that in both the 1D and 2D cases, the reaction front speed generally exceeds the local speed of sound, preventing compression waves from catching up and impacting the mixture ahead of the ignition front and hence, influencing the dynamics of ignition front propagation.

FIG. 19 .
FIG. 19.Front speed vs normalized distance profiles of 2D and 1D simulations indicating: (a) diverging fronts and (b) converging fronts, from ignition instance to the fully developed detonation.Scatterplots illustrate the 2D cylindrical simulation results and black markers depict the 1D results.The green dash-dotted lines mark the location of ignition for 2D and 1D simulations, and the black dashed line indicates the CJ speed.

FIG. 20 .
FIG. 20. 1D transient profiles of (a) temperature and (b) normalized pressure at point 3 for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, illustrating the supersonic ignition front propagation.The dashed line indicates the constant volume pressure value P cv .

( 3 )
FIG. 21. 2D profiles of temperature (top row), log 10 (HRR) (middle row), and pressure (bottom row) for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, illustrating the appearance of supersonic ignition front propagation.

Figure 22 (
FIG. 22. 2D planar detonation simulation results, (a) schematic presentation of the simulation setup, (b) detonation cellular structures inside the channel, (c) sensitivity of the triple point trajectories to the grid resolution.

TFIG. 23 .
FIG.23.2D profiles of the pressure gradient normalized by its maximum value (rPn) (top row) and temperature (bottom row) and at point 2 for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, illustrating diverging detonation and reflected shock (RSW) fronts.

FIG. 24 .
FIG. 24.2D profiles of OH mass fraction (Y OH ) (top row) and O 2 mass fraction (Y O2 ) (bottom row) at point 2 for the PRF mixture at T ¼ 900 K and P ¼ 50 bar, The white contour lines in panels (d)-(f) indicate the reaction zone (Y O2 ¼ 0.022).
, T s , r s ), are provided in TableI.Additionally, the ZND (Zeldovich-von Neumann-D€ oring) induction length (D i ), Chapman-Jouguet (CJ) speed (U cj ), and the domain length (L D ) are indicated in the table.In this study, D i is defined as the distance between the leading shock and the location of maximum heat release, Phys.Fluids35, 126102 (2023); doi: 10.1063/5.017477835, 126102-7 V C Author(s) 2023 source (P s

TABLE I .
Conditions for the three cases: initial pressure and temperature, von Neumann pressure, CJ speed, source pressure, and temperature, ZND induction length, domain length, source radius, and minimum mesh size.

TABLE II .
64nditions for the three cases: pressure, average temperature, von Neumann pressure, CJ velocity, induction length, domain length, minimum mesh size, number of processors, and number of cells.'s capabilities by presenting results from simulations involving the direct initiation of detonation in an H 2 /air mixture at relatively low temperatures and pressures.This simulation employs a ten-species chemical mechanism58to handle finite rate chemistry.Finally, we delve into more computationally intensive scenarios, focusing on finite rate chemistry in low to intermediate temperature reactions of the PRF mixture, utilizing a 115-species kinetic mechanism.64Thesecases aim to demonstrate the auto-ignition behavior of a common gasoline surrogate (PRF mixture) under high pressure and intermediate temperature conditions relevant to SI engines.It is notable that an additional 2D planar direct initiation of H 2 /air detonation simulation is conducted and presented in Appendix A to explore mesh sensitivity and validation.Throughout Secs.III B and III C, 1D planar simulations are also performed alongside 2D cylindrical simulations to facilitate comparison and verification.

TABLE III .
Configuration and results of the PRF and H 2 test cases.

TABLE IV .
Summary of initial ignition regime, the corresponding NTC effect, and the key observations from the 1D planar and 2D cylindrical simulations of points 1-3.

TABLE V .
The 2D planar detonation case configuration.