Analysis of knock onset based on two-dimensional direct numerical simulation and theory of explosive transition of deflagration

Engine knocking poses a challenge in the development of high-efficiency spark ignition (SI) engines, as it can lead to decreased performance, increased emissions, and even engine damage. Consequently, the knocking phenomenon has been extensively investigated for a long time in order to prevent engine knocking and improve engine performance. Starting with the understanding of various auto-ignition modes in inhomogeneous reactivity fields outlined by Zel'dovich,¹ a sequence of subsequent inquiries by Bradley's group^2–4 has substantially augmented the understanding of the knocking phenomenon.

In recent years, computational fluid dynamics (CFD) analysis employing intricate chemical reaction kinetics has been ubiquitously executed, enabling the acquisition of information on the complex interplay between fluid dynamics and chemical reactions, which is generally difficult to obtain experimentally in detail. Utilizing a 1D DNS, Ju et al. demonstrated that ignition and flame propagation are influenced by transport and low-temperature chemistry using n-heptane fuel.⁵ Subsequently, Dai et al. revealed that hot and cool spots interact with negative temperature coefficient (NTC) behavior to induce engine knocking in a genuine SI engine.⁶ Terashima et al. scrutinized the cumulative effect of pressure waves from ignition to knocking onset, grounded on 1D DNS utilizing n-heptane fuel.⁷ Apart from 1D DNS, multidimensional DNSs have been conducted, albeit a limited number of studies due to the substantial computational costs involved. In a 2D DNS for hydrogen, the multidimensional effects of pressure waves on the entire process of flame propagation, knock generation, and detonation propagation have been examined.⁸ Moreover, Luong et al. carried out 2D and 3D DNS to investigate the knock onset mechanism based on the effect of temperature inhomogeneity in the presence of turbulence.⁹ Morii et al. executed 2D DNS using n-heptane¹⁰ from ignition to knock onset, quantitatively reproducing the existing knocking experiment¹¹ for the first time. Through the implementation of such DNS, knock onset mechanisms have been progressively elucidated. In recent years, it has been reported that an auto-ignition assisted flame^12–14 engendered by low-temperature oxidation reactions in the preheat zone is observed in 1D laminar premixed flame simulations under elevated temperature and pressure conditions. The notions of the auto-ignition assisted flame are anticipated to contribute to the comprehension of the deflagration wave preceding the knock onset. Despite the progress made in understanding engine knocking, there remains a need for further investigation into the knock onset mechanisms and the role of deflagration in this process.

In this study, we aim to provide a detailed analysis of the flame structure in our previous 2D DNS including knocking phenomenon, with a particular focus on the role of deflagration in the knock onset mechanism. Initially, we analyzed the 2D DNS results¹⁰ and extracted information pertinent to the role of deflagration therein. Subsequently, a knock onset mechanism was explored based on the notion of the explosive transition of deflagration.¹⁵ 

Outline of the 2D DNS¹⁰ is given here for the readers' convenience. Target phenomenon of the 2D DNS is an experiment by Kitagawa's group¹¹ using a constant volume chamber (14.0 × 14.0 × 80.0 mm) for a stoichiometric n-C₇H₁₆/O₂/Ar mixture at the mole fraction ratio of O₂ and Ar of 21:79, an initial temperature of 480 K, and a pressure of 0.33 MPa, and the size of the constant volume chamber is approximately the same as the bore diameter of a cylinder of an SI engine. The computational domain is 40.0 × 7.0 mm, corresponding to one quarter of the experiment, which is the region surrounded by blue bold lines in Fig. 1. The left and bottom boundary conditions were set as adiabatic nonslip walls, and the right and upper boundary conditions were set as symmetric. For ignition, a hot kernel of 1400 K with a radius of 1.0 mm was set at x = 2.5 and y = 7.0 mm. A 5120 × 896 equally spaced orthogonal mesh was placed in the domain (40.0 × 7.0 mm), and the mesh size is about 8.0 μm. Fine grids were placed on the flame front by setting the AMR level as 1, and the minimum mesh size is about 4 μm. The grid convergence was examined on the knock onset timing using 1D DNS.

The governing equations are compressible Navier–Stokes equations with mass conservation of chemical species and energy equation. To perform 2D DNS with high efficiency, PeleC,^16,17 a compressible reactive flow solver with the AMReX library that can handle dynamic adaptive mesh refinement (AMR), was employed. Harten, Lax, van Leer, and Contact (HLLC)¹⁸ was used as an approximate Riemann problem solver. The piecewise parabolic method (PPM) was employed as an interpolation method to reconstruct the values at the face.¹⁹ The second-order central difference was used for the transport terms. The spectral deferred correction approach (SDC),^20,21 which is an iterative scheme, was employed as a time advance scheme, and the number of iterations was set as two. To achieve the 2D DNS for knocking experiments by further reducing computational cost, an efficient chemical equation solver, MACKS (Minimum Error Adaptation of Chemical Kinetic Solver),²² was incorporated into PeleC. Regarding chemical kinetics, the latest validated reduced chemical kinetics^23,24 developed in the SIP Engine Project in Japan was applied. Because the fuel in this study is n-heptane, reactions related to iso-octane were excluded from the original kinetics. The numbers of chemical species and elementary reactions after the above exclusion are 53 and 199, respectively. The validation cases of either full or reduced SIP kinetics²⁵ can be found in the literature.^{23,24,26–28} Additionally, since the chemical kinetics developed in the SIP needs to handle PLOG function, Zero-RK²⁹ was incorporated into PeleC.

In this study, critical information from the 2D DNS¹⁰ was extracted and juxtaposed with the outcomes of newly conducted 1D laminar premixed flame simulations to comprehend the knocking phenomena. 1D laminar premixed flame analyses were conducted utilizing Cantera 2.5.1.³⁰ The numerical conditions for the inlet pressure and temperature for the simulations were the pressure and temperature at the location illustrated in Fig. 1 for each time from the 2D DNS. The mole fractions of the chemical species were based on the stoichiometric n-C₇H₁₆/O₂/Ar mixture. The size of the computational domain for the 1D laminar premixed flame simulations varied from 1.0 to 64.0 mm, owing to the presence of an auto-ignition assisted flame, which exerts a substantial influence on the results. The employed chemical kinetics^23,24 are the same as those used in the 2D DNS.

Since comparisons of overall reaction front behavior, the first low temperature oxidation (LTO), and knock onset timings between the experiment¹¹ and 2D DNS were already made in our previous study,¹⁰ further details, in particular, three items extracted from the 2D DNS, are discussed here: the time history of pressure and temperature of the unburned mixture (Fig. 2), the state in the vicinity of the reaction front (Fig. 3), and the instantaneous flame speed (laminar burning velocity) calculated from the behavior of the reaction front and local flow velocity (Fig. 4). Note that only the 2D DNS result is indicated here, since it was found that quantitative comparison was not possible by 1D DNS due to artificially emphasized effect of pressure waves and unresolved flame front structures on the reaction front behavior.^10, Figure 2 shows the experimental pressure history¹¹ and pressure and temperature histories of 2D DNS.¹⁰ The figure shows that the 2D DNS reproduced the experimental pressure histories well. Note that the estimated timing of the knock onset in 2D DNS, which was obtained after the same post-process, i.e., the peak position of the time history of $d P / d t$⁠, with experimental knock onset was $5.76 × 10 − 3$ s, which is slightly faster than the experiment at $5.83 × 10 − 3$ s but with only 2% error. Thus, it was confirmed that the 2 D DNS could reproduce the experiment.¹¹ In the followings, 2D DNS results will be further utilized for investigating the phenomena including knocking onset.

Since the present 2D DNS¹⁰ includes detailed information from ignition, flame propagation and transition of deflagration into knocking, it is worth to extract information on the role and structure of deflagration during the process. Therefore, the role and structure of reaction front was examined. Here, the results at $3.958 × 10 − 3$ s where the deflagration front propagated with the distorted tulip flame shape are shown in Fig. 3. Figure 3 shows that an enlarged image of reaction front [Fig. 3(a)], the temperature and pressure profiles [Fig. 3(b)], temperature and velocity profiles [Fig. 3(c)], and velocity divergence and heat release rate profiles [Fig. 3(d)] of the 2D DNS at $3.958 × 10 − 3$ s. The gray areas in Fig. 3 indicate the region where the sign of the velocity divergence is being positive near the flame front [see Fig. 3(c)]. This gray area is a good representation of the characteristic deflagration wave as follows. For shock wave detection, a region where velocity divergence exists is widely used because the velocity divergence become infinite negative at the shock front. Since the reaction front is with a strong source term, we expect that the sign of the velocity divergence to be positive, unlike shock waves. Thus, the region where the sign of velocity divergence is positive can be regarded as the proof of the existence of the “healthy” structure of deflagration. We term this region as “hydrodynamic flame zone” and use it as a measure of the existence of healthy deflagration. Figure 3(c) shows that the hydrodynamic flame zone covers the preheat zone, reaction front with a positive heat release rate (reaction zone), and some burned gas region, which are the essential components of the healthy flame structure of 1D laminar premixed flame. Since the nature of hydrodynamic flame zone is hydrodynamic source due to exothermic reactions, behavior of this zone would be informative and it would play important role in, e.g., further discussion of turbulence–chemistry interaction³¹ using various DNS results with detailed kinetics. By the way, the present 2D DNS was for a compressible reactive fluid, and the characteristic time of the chemical reaction is much shorter than the characteristic time of the fluid near the flame front. This led to the existence of a sharp pressure peak in the reaction zone [see Fig. 3(a)]. The position of the peak is the same with the position of the peak of heat release rate. Therefore, the flame front velocity in the laboratory frame could be determined by tracking this pressure peak at the flame front, which can be easily derived from the pressure, a primitive variable.

The concept of the hydrodynamic flame zone is also of use when general definition of deflagration is required in DNS. The position where the velocity divergence becomes zero in the upstream is defined as $x unburned , DNS$ for the unburned region and that in the downstream, $x burned , DNS$ for the burned region (see Fig. 3). Here, the hydrodynamic flame zone has a thickness of $( x unburned , DNS − x burned , DNS )$⁠. The laminar burning velocity (flame speed) can be defined as the difference between the velocity at $x unburned , DNS$ and the flame front propagation velocity. Note that a low-pass filter was applied at 5000 Hz to the time history results, since the acoustic resonance mode has a very large effect on the velocity fluctuation. The procedure for extracting instantaneous and local burning velocity from the 2D DNS is summarized as follows:

1. Detecting the pressure peak in the reaction zone and evaluating its velocity.

2. Setting the edge of the unburnt region, where the velocity divergence is zero as $x unburned , DNS$⁠.

3. Subtracting the flame front propagation velocity from the velocity at $x unburned , DNS$ to obtain the instantaneous and local burning velocity.

Next, the status of flame propagation in 2D DNS was investigated based on the discussion in Sec. III B. Figure 4 shows the time history of the instantaneous flame speed of 2D DNS and the 1D laminar burning velocity using pressure and temperature of unburned mixture of 2D DNS with different computational domain sizes. Here, the four computational domain sizes of 1.0, 4.0, 16.0, and 64.0 mm were used. Note that the flame front of 1D laminar premixed flame in Cantera is set to one-third of the computational domain by default, i.e., about 0.33 mm when the computational domain is 1.0 mm. When the computational domain is large, i.e., over 4.0 mm, it is no longer possible to reproduce the results of the 2D DNS. Note that from the time zero to $0.5 × 10 − 3$ s, the burning velocity of the 2D DNS does not agree with the 1D laminar burning velocities. This is just due to the fact that the flame front is affected by the ignition energy and the initial pressure wave generated from the hot kernel. The results showed that the flame structure within the hydrodynamic flame zone obtained in the 2D DNS can be reproduced using the 1D laminar flame simulations using pressure and temperature of the unburned mixture of 2D DNS as an inlet condition. This is corroborated by the fact that the computational domain size of a 1D laminar premixed flame, capable of reproducing a 2D DNS, aligns well with the thickness of the hydrodynamic flame zone as defined in Fig. 3 (approximately 1.0 mm). Furthermore, a more detailed examination of Fig. 3 reveals that the front of the hydrodynamic flame zone is situated just ahead of the flame front. As a result, we propose that a 1D laminar flame simulation, with the inlet boundary placed as close as possible to the flame front, would yield a more accurate replication of the 2D DNS. Therefore, we conducted a 1D laminar premixed flame simulation adhering to this setup.

The threshold conditions are determined as the condition when we could not get the converged solution. As a result, the threshold condition was $5.73 × 10 − 3$ s, where the threshold temperature for the calculated burning speed was 1287.2 K, and the pressure was 2.18 MPa. The time at which the 1D laminar premixed flame cannot reach the converged solution is very close to the knock onset timing of $5.76 × 10 − 3$ s in 2D DNS. Thus, there is a strong correlation between the knock onset and the limit of the existence of the 1D laminar premixed flame structure. The transition to knocking is considered to be caused by the inability of the flame structure to hold. To support this discussion, a newly proposed concept of explosive transition of deflagration¹⁵ can be applied to 2D DNS. The outline of the theory is explained in Sec. IV.

In this section, we introduce the theory of explosive transition of deflagration. This theory, which is based on the governing equations of 0D homogeneous ignition and 1D laminar premixed flame, provides a framework for understanding the transition from deflagration to overall ignition, which is equivalent to detonation. We discuss the implications of this theory for our understanding of the knocking phenomenon. For more comprehensive information on the theory discussed in this section, along with additional figures, please refer to our earlier work.¹⁵ 

Figure 5 shows the relation between $T ̃$ and $Y ̃$ for a 1D laminar premixed flame under an initial pressure and temperature of 0.1 MPa and 300 K. The mixtures are stoichiometric H₂/O₂/N₂, H₂/O₂/Ar, n-C₇H₁₆/O₂/N₂, and n-C₇H₁₆/O₂/Ar mixtures, and their Lewis numbers of fuels are 0.41, 0.40, 2.64, and 2.60 at 0.1 MPa and 300 K, respectively. The multi-step chemical reaction model for hydrogen fuel was the UT-JAXA model,³² and the model for n-heptane fuel was the reduced SIP model.^23,24 The mole fraction ratios of O₂:N₂ and O₂:Ar are 21:79.

As discussed in the theory subsection for $L e f = 1$⁠, 0D homogeneous ignition using normalized fuel mass fraction and temperature is equivalent to 1D laminar premixed flame using normalized fuel mass fraction and temperature after spatiotemporal transformation. As described in the paper by Aspden et al.,^33, Fig. 5 shows that for $L e f < 1 , T ̃$- $Y f ̃$ curve for a 1D laminar premixed flame is convex on lower left, which is below the $T ̃$- $Y f ̃$ line of the $L e f = 1$ case (i.e., $T ̃$- $Y f ̃$ line of 0D homogeneous ignition), while for $L e f > 1 , T ̃$- $Y f ̃$ curve is convex on upper right, which is above that of the $L e f = 1$ case (i.e., $T ̃$- $Y f ̃$ line of 0D homogeneous ignition) and may intersect $T ̃$- $Y f ̃$ line of the $L e f = 1$ case (i.e., $T ̃$- $Y f ̃$ line of 0D homogeneous ignition) in the high temperature region of the flame (i.e., reaction zone) where the heat release is high. This is due to the fact that for $L e f < 1$⁠, the diffusion coefficient of the fuel is greater than the thermal conductivity, and the fuel is preferentially diffused rather than heated, while for $L e f > 1$⁠, the thermal conductivity of the fuel is greater than the diffusion coefficient and the fuel is preferentially heated rather than diffused.³³ The influence of the Lewis number on the relationship between temperature and fuel mass fraction (⁠ $T ̃$- $Y f ̃$ curve) is discussed theoretically in Buckmaster's book,³⁴ and the results in Fig. 1 are consistent with his theory. From now on, the discussion will be in the relatively low temperature region (i.e., the preheat zone), where the influence of the Lewis number is clear in Fig. 5. Thus, in the construction of the present theory, it can be assumed that the Soret and Dufour effects are negligible.

In Fig. 5, at a given normalized temperature, the reduction of the normalized fuel mass fraction, i.e., fuel consumption, can be considered as the progress of the reaction. Thus, the progress of reaction is of the order of $L e f < 1 , L e f = 1$ (i.e., 0D homogeneous ignition), and $L e f > 1$⁠. Therefore, if $L e f < 1$⁠, a 1D laminar premixed flame structure can exist even if the inlet temperature and pressure are very high. However, if $L e f > 1$⁠, depending on the inlet temperature and pressure conditions, there may be conditions where a 1D laminar premixed flame structure cannot exist. For example, since the Lewis number of gasoline fuel is greater than one, vanishment of premixed flame structure that corresponds to knocking can be expected to occur in SI engines whenever the unburned gas region rises to temperature and pressure conditions that make flame propagation impossible. The Lewis number should also affect the transition condition from deflagration to detonation, since the characteristic time of the flame behind the shock wave is of the same order of magnitude as that of ignition. It is interesting to note that the $T ̃$– $Y f ̃$ curve discussed here is also used with respect to the discussion of Type Ia supernovae.³³ 

Figure 6 shows the calculated burning velocity for different inlet temperatures at the pressure of 0.1 MPa. A 1D laminar premixed flame simulation was performed by varying the inlet temperature from 300 to 3000 K in 100 K increments, at a pressure of 0.1 MPa, and an equivalence ratio of 1.0. Note that the computational domain was kept as short as possible because $T ̃$ and $Y ̃ f$ must be convex functions to preserve the Legendre transform even with a multi-step chemical reaction model. After confirming that the simulation fails no matter how much the computational domain is reduced, the next step is to run the simulation varying the inlet temperature in 10 K increments to determine the temperature limit at which the simulation can obtain the converged solution. From Fig. 6, for hydrogen fuel (⁠ $L e f < 1$⁠), the burning velocity increases with increasing inlet and initial temperatures, indicating that flame propagation is possible even at 3000 K. In the case of n-heptane fuel (⁠ $L e f > 1$⁠), however, the burning velocities have certain upper limit values for the existence of flame structure. For instance, it is 1180.0 K for n-C₇H₁₆/O₂/Ar and it is 1270.0 K for n-C₇H₁₆/O₂/N₂. This fact is consistent with our discussions in the last paragraph, which shows that for $L e f > 1$⁠, there are certain conditions where flame structures cease to exist. Note also that the threshold temperature is apparently lower for n-C₇H₁₆/O₂/Ar than for n-C₇H₁₆/O₂/N₂, despite the small difference of the Lewis number. This is because, as discussed in the theory section, it is not strictly the Lewis number that determines the threshold temperature but $d d x ( ρ D f d Y ̃ f d x )$ and $1 c p d d x ( λ d T ̃ d x )$ that determines the threshold temperature. The adiabatic flame temperature varied with the specific heat ratio of the diluted gas and was 2345.4 K for n-C₇H₁₆/O₂/Ar and 2133.8 K for n-C₇H₁₆/O₂/N₂. In other words, the argon dilution has a higher adiabatic temperature, and it leads to a steeper temperature gradient, resulting in a lower temperature for the explosive transition. This is also corroborated by Fig. 5, which shows a faster reduction of normalized fuel mass fraction for n-C₇H₁₆/O₂/Ar as compared to n-C₇H₁₆/O₂/N₂.

Next, the flame structure is discussed in terms of the temporal evolution of temperature in 0D homogeneous ignition and 1D laminar premixed flames, as we have described in the theory subsection as well as our previous work.¹⁰  Figure 7 shows the temperature time histories of 0D homogeneous ignition and 1D laminar premixed flame for H₂/O₂/N₂ and n-C₇H₁₆/O₂/N₂ cases. The 0D homogeneous ignition with the multi-step chemical reaction model was computed using Cantera with constant pressure and enthalpy.³⁰ Note that the temperature profile of the 1D laminar premixed flame was obtained in the time direction by a spatiotemporal transformation using the residence time. From Fig. 7(a), it can be seen that the converged solutions of 1D laminar premixed flame can be obtained, although the temperature profile of the 0D homogeneous ignition intersects that of the 1D laminar premixed flame under conditions where the initial temperature is above 1500 K. However, Fig. 7(b) shows that the converged solutions of 1D laminar premixed flame cannot be obtained if the characteristic times of the 0D homogeneous ignition and the 1D laminar premixed flame are of the same order of magnitude.

Finally, we delve into the relationship between the 2D DNS and the explosive transition of deflagration. In 1D laminar premixed flame computations, if the inlet temperature and pressure are set to match those of the unburned gas region and no converged solution is found, auto-ignition will invariably occur in the unburned gas region. This happens because flame propagation in the unburned gas region is unachievable. Figure 8 illustrates the relationship between pressure and temperature in the unburned gas region, as extracted from the 2D DNS, and the boundary of the explosive transition of deflagration when pressure varies from 0.3 to 3.0 MPa in increments of 0.3 MPa. This pressure–temperature relationship is frequently employed to assess knock tendency using the Livengood–Wu integral.³⁵ When the unburned gas region conditions pass the explosive transition boundary depicted in Fig. 8, signifying the explosive transition of deflagration, auto-ignition occurs within the unburned gas region as the presence of a flame structure becomes untenable. Essentially, this translates into the initiation of a knock onset. Interestingly, these observations align with the timing of knock onset as illustrated in the 2D DNS (see Fig. 8). The aforementioned findings suggest that a reliable prediction of knocking could be feasible if we portray the explosive transition of deflagration in terms of temperature and pressure, provided the trajectory followed by the unburned gas under variable temperature and pressure conditions is known.

Data analysis of a 2D DNS that reproduced the knocking experiment was performed, and it was confirmed that the reaction front behavior in 2D DNS can be reproduced as a 1D laminar premixed flame simulation at extreme conditions. By investigating knocking with a 1D laminar premixed flame, we found that there is a strong correlation between the timing of knock onset and the conditions under which a 1D laminar premixed flame structure cannot exist. From the theory of explosive transition of deflagration, it was shown that when the Lewis number is unity, the time evolution of the normalized fuel mass fraction and temperature in a 0D homogeneous ignition is equal to the temporal evolution observed in a 1D laminar premixed flame, if the spatiotemporal transformation is properly applied. Furthermore, the rate at which the normalized fuel mass fraction decreases in the preheat zone was found to depend on the Lewis number, and when the Lewis number is greater than unity, no flame structure exists above a certain threshold temperature. Finally, by considering the theory of explosive transition of deflagration, the mechanism of knock onset was explained, and it was found that the knock onset occurs when the temperature and pressure conditions in the unburned gas region cross the explosive transition boundary.

This work was partially supported by JSPS KAKENHI (Grant No. 19KK0097).

The authors have no conflicts to disclose.

Youhi Morii: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Akira Tsunoda: Investigation (equal); Methodology (equal); Writing – review & editing (supporting). Ajit Kumar Dubey: Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Kaoru Maruta: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.