Comparison Between Experimental and Simulated Knock Statistics Using an Advanced Fuel Surrogate Model

. The statistical tendency of a GDI spark-ignition engine to undergo knocking combustion as a consequence of spark timing variation is numerically investigated. In particular, attention is focused on the importance to match combustion-relevant and knock-relevant fuel properties to ensure consistency with the experimental evidence. An in-house surrogate formulation methodology is used to emulate real gasoline properties, comparing fuel models of increasing complexity. Knock is investigated using a proprietary statistical knock model (GruMo Knock Model, GK-PDF). The model can infer a log-normal distribution of knock intensity within a RANS formalism, by means of transport equations for variances and turbulence-derived probability density functions (PDFs) for physical quantities. The calculated distributions are compared to measured statistical distributions. The proposed numerical/experimental comparison constitutes an advancement in synthetic chemistry integration into 3D-CFD combustion simulations.


Introduction
Knocking combustion in spark-ignited (SI) engines is an abnormal combustion phenomenon which produces loud noise and can result in severe engine damage under certain conditions. For downsized and boosted SI engines, knock represents an extremely stringent constraint on performance and efficiency since it prevents the use of lean and stoichiometric mixtures, advanced spark timings, high compression ratios and/or high boost pressures [1]. For example, at high boost level the spark timing needs to be substantially delayed from maximum brake torque (MBT) timing to avoid self-ignition of the end-gas, preventing maximum thermal efficiency. Similarly, rich mixtures are often exploited to lower unburnt gas temperature and decrease risk of auto ignition. It is in fact generally accepted that engine knock is caused by the autoignition of a portion of the end-gas prior to the flame arrival. Extremely rapid heat release associated with the end-gas autoignition causes sharp increase in local pressure, and the propagation of pressure waves of substantial amplitude across the combustion chamber. The high-speed pressure waves exert sizeable forces on the surroundings of the combustion chamber and can cause mechanical damage to metal parts. Uncontrolled end-gas autoignition is the root cause of knock, while the temperature and pressure histories of the end-gas are, in turn, governed by the phasing and rate of development of the flame and by fuel stratification. As anticipated, knock-safe conditions are typically pursued through a combination of boosting limitation, charge enrichment and spark advance (SA) reduction. All these remedies relevantly affect engine-out performance and efficiency. Knock is a highly stochastic phenomenon that is dependent on the local flow patterns, fuel and residual gas distribution, temperature and turbulent flame history. Such nature of engine knock, related to cycle-to-cycle variability (CCV) of turbulent flow and combustion, would suggest Large-Eddy-Simulation (LES) as the most appropriate approach for CFD simulations [2,3,4]. In order to limit computational costs and times, RANS models are usually chosen to represent the average engine behaviour. Several knock models are available in literature to predict average knock onset in SI engines [5,6], though they all suffer from the intrinsic inability to account for the stochastic nature of the modelled phenomenon. Such limitations were recently overcome by the use of variance equations for knock-targeting physical variables by d'Adamo et al. [10]; a resume of the methodology is presented in a dedicated paragraph. The model results are used to reconstruct a presumed distribution function of knocking cycles thanks to turbulence-generated variance of physical fields, which in turn affect the end-gas reaction rate towards autoignition. Regardless of the approach towards knock modelling, accurate modelling of the end-gas reactivity and of the turbulent flame propagation process are key factors to successfully predict engine knock. Fuels used in internal combustion engines (ICEs) are complex mixtures of a wide spectrum of hydrocarbons. The nature and proportion of gasoline compounds are ever varying and only partially known, because they depend on the origin of the crude oil, on the oil refinery process and on regional regulations. Concerning the modelling of reacting flows in ICEs, the lack of precise knowledge is a hindrance that can be overcome through the definition of simplified recipes called fuel surrogates, i.e. mixtures of few representative chemical species (e.g. Parrafins, Iso-Paraffins, Olefins, Naphtenes, Aromatics, Oxygenates etc.) with prescribed proportions. Blending proportions aim to match a precise target property, either on the physical side (e.g. density, evaporation), or on the chemical side (e.g. laminar flame speed (LFS), ignition delay time (IDT), soot tendency etc.). Broadly speaking, iso-octane is the most popular and simple surrogate for gasoline, and it is often used for flame propagation studies; yet it is not representative for autoignition because of its high octane number (ON). Another possibility to emulate real gasoline is using a Primary Reference Fuel (PRF), i.e. a surrogate formed by iso-octane and n-heptane. Such binary mixture is seldom satisfactory because different PRF mixtures are needed to match the gasoline autoignition behaviour over different engine operating conditions [8]. Such limitation is due to the marked difference between linear and branched paraffinic fuels (such as PRFs) and other hydrocarbon classes, which account for about half the content of an actual gasoline. One of the most important parameters that distinguishes the knock attitude of a fuel is octane sensitivity (S), which is quantified as the difference between the Research Octane Number (RON) and the Motor Octane Number (MON). It is a common practice to use Toluene Reference Fuels (TRF) since they provide more accurate representation of the fuel by adding non-paraffinic compound to PRF mixtures. Measurements in shock tubes as well as in engines have proved their ability to behave like a real gasoline in term of autoignition delay. Finally, recent regulations promoting the increase of oxygenated components in transportation fuels led the authors to propose a methodology to add Ethanol to the surrogate palette [9]. The main outcomes of the methodology are recalled in a dedicated paragraph. The present paper reports a numerical study of an optically accessible research engine available at the Istituto Motori (IM-CNR) operated under soft knocking conditions [7].

Fuel surrogate methodology
In [9] the authors introduced an efficient methodology to formulate gasoline fuel surrogates able to match the main chemical and physical properties, the autoignition and the flame propagation characteristics of a commercial gasoline. Three fuel surrogates of increasing complexity were formulated and validated against laminar flame speed, shock tube (ST) and rapid compression machine (RCM) experiments available in literature for oxygenated gasoline over a wide range of pressures, temperatures and equivalence ratios. Particular care is needed to select which key properties of the real fuel need to be targeted: H/C is essential to correctly describe SL, adiabatic flame temperature, heat of combustion, αs, density, molecular weight and boiling point. O/C ratio is also included as a target to extend the validity of method to oxygenated fuels [32]; finally, RON and MON are needed to properly emulate end-gas reactivity. Once these targets are known, a set of equations is set to estimate the corresponding properties of the proposed surrogate; a linear mathematical formulation (Le Chatelier's mole fraction-based mixing rule) is adopted to formulate a set of equations whose solution can be obtained following [9].
Two surrogates are here formulated and compared: a four-components Ethanol Toluene Reference Fuel (ETRF) and a six-components one (named SEN in the following). To validate the methodology focusing on end-gas reactivity, it is important to demonstrate the capability of the surrogates to match the autoignition behaviour of a RON95 gasoline with a sensitivity around S=10. Coryton and Halterman gasoline fuels ignition delays were tested in shock tubes (ST) and rapid compression machines (RCM) over a wide range of Tu= 650 -1250K, p=10-40 bar over different equivalence ratios (0.45-0.9-1.8). Ignition delay calculations were performed in DARS v2019.1 chemical solver, licensed by Siemens Digital Industries SW, using constant pressure reactors and two different chemical semi-detailed mechanisms (for the sake of brevity, only the results using PoliMi [22] mechanism will be shown in this paper). Results shown in Figure 1, suggest that a proper gasoline fuel surrogate can be used, together with a validated chemical kinetics mechanism, to adequately estimate autoignition characteristics. Since the high temperature region tends to hide the reactivity dependence on RON/MON and composition, the most interesting region to validate the proposed approach is the intermediate temperature range (~700 − 850 ). This is of crucial relevance in ICEs reacting flows as it is the temperature range typical at end-of-compression stroke and early flame development, i.e. when chemical reactivity starts to exponentially increase. As clearly visible in the figure, the two fuels are characterized by a different behaviour in the region. Coryton gasoline has a lower reactivity, consistently with its higher (97.5). Most interestingly, Haltermann gasoline exhibits a more pronounced behaviour compared to Coryton gasoline, due to its lower octane sensitivity (7.6) . The two surrogate formulations predict largely different reactivities. In particular, the surrogate overestimates the reactivity in the NTC range. Conversely, the surrogate is in very good agreement with the experimental measurements over the full range of pressures and equivalence ratios. It can be concluded that the proposed surrogate better describes the reactivity dependence on pressure and mixture composition in the region, which in turns is mainly driven by the improved estimation of the octane sensitivity. It also better matches the reactivity of both tested gasolines in the low temperature region ( < 700 ). Consistent considerations can be made for the different temperature regions: both surrogates are quantitatively able to represent the autoignition quality of the fuels in the high temperature range. Moving to the intermediate region, the surrogate better matches the fuel reactivity dependence on mixture quality, closely resembling the Coryton fuel autoignition characteristics due to the similarity of . Such outcome is of particular interest for engine knock simulations, considering that the end-gas region close to Knock-Limited Spark Advance ( ) conditions usually exhibits temperature in the = 750-900 range. In the low region, the surrogate can quantitatively match the Coryton gasoline behaviour. Conversely, the surrogate exhibits a mixture reactivity at 10 bar equal to the experimental one at 20 bar, and it overestimates gasoline reactivity over the whole mixture quality range.  to track the progress towards autoignition; such approach is usually based on a unique cell-averaged ignition delay time ̃= =(̃, ,̃,̃). The autoignition delay time is calculated either using empirical correlations or using a look-up table approach. [6]. A major limitation of this method is the use of cell-averaged values for the input variables contained in φ vector. This implies that for a generic fluid cell a perfectly uniform value of pressure, unburnt temperature, equivalence ratio and residual mass fraction is assumed, i.e. every computation fluid cell is considered as a laminar wellstirred reactor. Even in a RANS framework, ensemble averaged turbulence may affect the local thermal and mixing state, which in turn can induce variations in the reactivity of the unburnt charge. To account for mixture reactivity dispersion around the average, a statistical knock model was proposed by d'Adamo et al. [10,27,28]. The model, hereafter called GK-PDF, is based on the reconstruction of the physical states whose presence is statistically accounted for in the same fluid cell. Their dispersion around the mean value is driven by two competing factors: inhomogeneity due to field gradients and local turbulent mixing. Such competition is expressed by the transport of two additional quantities, namely temperature and mixture fraction variance Tu' and Z' (Eq.1 and 2).
The right-hand side terms are respectively a source term (i.e. variance generated by field variable spatial gradients) and a well term which describes variance dissipation due to turbulent mixing in the fluid domain. The letter term is a function of a turbulent relaxation time-scale: Where ct(Ret) -1 is a model constant calculated as a monotonic increasing function of the local turbulent Reynolds number. It is easy to note that a low turbulent level (i.e. low Ret) implies a small value of ct , this in turn resulting in a long turbulent relaxation time-scale. Consequently, the variance destruction operated by turbulent mixing is slow and the probability to find far form average in-cell states is high. The opposite situation characterizes highly turbulent conditions, for which the systems evolves towards a perfectly stirred status. Similarly, in a perfectly homogeneous isothermal field no Tu and Z gradients are present, and the variance is null. The opposite for highly heterogeneous situations, where the second moments for Tu and Z are not null. In order to model a statistical reconstructions of all the possible combinations of = ( , ) and Tu which can be present in the cell, it is possible introduce the general formulation for bi-variate Gaussian distributions (Eq. 3). This choice is motivated by its relative simplicity and its robust representation of in-cell distribution of states.
Where is the correlation coefficient which is calculated in the model at each iteration while ̃ ; ̃ are the mean values instead ; are the standard deviation of the respective mean values. The continuous PDF of the in-cell states is discretized into 9 steps for both Z and Tu (i.e. 81 possible (Z, Tu) states) covering 2 standard deviations (i.e. 95% probability) for each fluid cell. An AI delay for each (Z, Tu) combination is interpolated from the LUT, thus describing a distribution of in-cell reactivity; a 2-moments framework is then introduced to synthetize the calculated field and to account for both the average in-cell average reactivity and its faster-than-average deviation. A cell-averaged ignition delay time (τd) is calculated from all the (Z,Tu) states using the bivariate PDF as a weight. In this way higher probability states (e.g. the (Z, Tu) condition) are considered as more representative in the averaging process than lower-probability ones. As the τd value alone would not be sufficient to fully represent the reactivity distribution, a second moment is added, i.e. the calculated standard deviation τRMS. The model transports a presumed distribution of reaction rates, hence earlier-than-average knock can be predicted based on the assumed distribution function: while the use of the PDF-weighted average τd leads to the mean knock prediction, its substitution with delay times reduced by 2 or 3 standard deviations is in charge of estimating knock-proner situations. Clearly, the probability of the knocking event associated with τd'' or τd''' delay times lowers as they correspond to far from average (Z, Tu) conditions. When the AI heat release is activated, the KI associated to the knocking event is calculated and it is used to directly compare CFD and experiments. To this aim, the most reasonable choice would be to use the same knock index, e.g. MAPO. The numerical simulation of the pressure wave propagation could be used to directly estimate MAPO from CFD; however, this would result in a very challenging task: time-steps and grid density should be well resolved to adequately comply with the Courant Friedrichs Lewy (CFL) condition. Also, the pressure signal should be monitored at the very same position as the pressure transducer and the CFD signal should be filtered with the same frequency as the raw transducer trace. To avoid these issues and to retain a reasonable CPU efficiency, an alternative method to estimate the pressure rise given by knock is proposed. It is based on the energy released by AI and it can be used in RANS frameworks. The estimated MAPO function [26] (eMAPO, Eq. (4)) is defined as the pressure rise corresponding to the constant volume (instantaneous) combustion of the auto-ignited fuel. In Eq. (4) pKO is the in-cylinder pressure at KO, / is the fraction of the AI fuel over the gas mixture, Ki is the fuel lower heating value, TKO is the unburned temperature at KO, while cp and MW are the specific heat and the molecular weight of the unburned mixture and Ru the universal molar gas constant.
The eMAPO index derives from the Knock Severity Index (KSI) proposed by Klimstra in [26] relating the knock pressure wave to an isochoric auto-ignition, and it shares the same goal to be an objective and easily measurable indicator for Knock Intensity (KI). Moreover, it can be directly compared to the experimental MAPO. In fact, all the variables in Eq. (4) are available from CFD simulations while the fraction of knocking fuel is derived from the heat released by AI. This energy-based approach permits to estimate the knock pressure rise in RANS framework through the eMAPO and using the same limit as in the experiments (eMAPOmax=2 bar). The eMAPO is calculated in 3D simulation and is used to reconstruct a log-normal distribution of KI. Then, the associated CDF is used to estimate the knock probability below/above the eMAPOmax limit. The average KI is estimated using eMAPOAVER by activating the AI heat release associated with the average knock precursor (corresponding to τd), while the eMAPORMS is inferred as the difference between the eMAPO n associated to the faster-than-average knocking event and eMAPOAVER. This is formalized in Eq. (5), in which eMAPOn is the one resulting considering as source term for the knock precursor τd n divided by n.
Using the pair of eMAPO (eMAPOAVER ; eMAPOn ) it is possible to build a PDF of knocking cycles and also a Cumulative Density Function (CDF) to presume a probability of knocking cycles. In analogy to typical engine-out probability functions to describe the fraction of knocking cycles, a log-normal distribution function is chosen (Eq. 6) where μ and σ are the mean and the standard deviation of ln(x). They are calculated using the above mentioned eMAPOAVER/n. Coherently with the RANS meaning for mean and variance, the mean value M is taken as the eMAPOAVER at knock onset, while the standard deviation is the difference in burnt fraction between eMAPOn and eMAPOAVER. The linear system is defined as (Eq. 7): The two-equation system in two unknowns (μ;σ) is solved and the expression for the log-normal parameters are obtained. Finally, the fraction of knocking cycles is calculated as (Eq. 8): Methods and models Table 1 shows the main geometrical and operating characteristics of the investigated engine. A 3D model of the engine is built using a customized version of STAR-CD v4.30, licensed by SIEMENS Digital Industries SW. As shown in Figure 2, symmetry is exploited to reduce the computational effort.
As already expressed in [11] the GDI optical unit under investigation features a non-negligible crevice volume which is included in the CFD model to account for compression ratio reduction and blow-by losses. The total number of fluid cells ranges from 1.2 to 0.4 million cells at TDC and BDC, respectively. The global average mesh size is around 0.8 mm for both combustion chamber and ports.  The k-ε RNG model for compressible flows is chosen to close the set of Navier-Stokes equations. A calibrated 1D model of the engine is used to impose time varying pressure and temperature boundary conditions at both the intake and the exhaust port. An additional mass flow rate is applied at the annular area at the bottom of the crevice to model the blow-by losses. Uniform wall temperatures are applied at each engine component facing the combustion chamber and the GruMo-UniMore wall heat transfer model [12,13,14,29] is used to estimate wall heat transfer. The ECFM-3Z combustion model [15] customized with a correlation for laminar flame speed at engine-relevant conditions [16,17] is used. A relatively simple flame kernel deposition model is adopted to model spark ignition. Since the study focuses on knock prediction, the stratification of fresh charge has a key role in the simulation [23,30,31].
To reduce modelling uncertainties, a high degree of accuracy in describing the six-hole full-cone spray evolution is required. Therefore, the spray is modelled following the approach described in [18,19,20,21,24,25] for multi-hole GDI injectors. Commercial RON95 gasoline is used as a reference fuel targeting the lambda value measured at the test bench (1.11). As stated earlier, autoignition delays are pre-computed off-line using a 0D chemical solver based on a constant pressure reactor model. The CRECK mechanism with 156 species and 3465 reactions [22] is used to provide autoignition delays. To minimize errors due to extrapolation, the range for the Look-Up Tables is purposely extended to possibly cover the whole spectrum of physical and chemical states within the combustion chamber and throughout the engine cycle. Stepping and ranges are reported in Table 2, and the resulting table covers more than 50000 thermodynamic states. A demonstration of the impact of surrogate formulation on autoignition estimation is reported in Fig.  (3). Such outcome is of interest for engine knock simulations, considering that the end-gas region under incipient autoignition conditions usually exhibits temperature values in the range 750-900 K. In the low-to-mid temperature region, the differences between ETRF and SEN ignition delay times are not negligible. Furthermore, it is possible to note how the high non-linearity of the percentage difference, which is dependent on both unburnt temperature and equivalence ratio. As an example, the highest differences are met around slightly rich mixtures. Differences between the two surrogate formulations decrease for very lean mixtures and they become negligible only for very hot mixture states, above 950K. Such evidence is of high relevance, as it justifies the effort to formulate a proper surrogate fuel and it discourages from the introduction of calibration factors and constants for the ignition delays.

Results
As visible in Figure (4.a), a good agreement is visible between experimental data, 1D and 3D simulations in terms of in-cylinder pressure history and Mass Fraction Burnt (MFB). Figures (4.b-c) show temperature and equivalence ratio fields on a section plane orthogonal to the cylinder axis passing  Looking at the average ignition delay times, depicted in Fig (5.a-b), it is easy to note that the ETRF surrogate shows higher reactivity in the unburnt mixture with respect to the SENARY surrogate, while the ETRF case seems to be more prone to knock. As we explained, two simulations are needed to infer the percentage of knocking cycles, predicting the mean and the faster-than-average knock occurrence. Pressure traces for the sequence of four simulations using the two surrogates are reported in Fig 6. a-b.
To discard the effect of initial conditions, two consecutive cycles are simulated; results are reported for the second cycle. The large difference between the two surrogates in terms autoignition prediction is visible. It is worth to remind that the investigated engine shows low knock probability.  The final aim of the PDF-Knock Model is to estimate a probability of knocking cycles. In the experimental practice, a knock threshold value for KI must be defined to distinguish knocking and knock-safe cycles. KI is expressed by the estimated MAPO (eMAPO) which is proportional with the residual unburnt fraction at the knock onset. Probability and Cumulative density functions (PDF-CDF) are shown in Fig. 7. Using a value of eMAPO=2 bar as threshold (dotted black line) it can be noticed that the ETRF surrogate predicts 100% of knocking cycles. Knock occurrence is predicted not only for the faster-than-average simulation, but also for the mean cycle. This is an evident overestimation of the knock signature of the tested engine and it is attributed to the overestimation of the reactivity when using the ETRF fuel surrogate. On the other hand, a promising result is obtained by the SENARY model, which closely approaches the experimental indication of 10% of knocking cycles; late knocking is detected in the faster-than-average simulation whereas the average cycle is knock-free, in agreement with the experimental evidence. A notable improvement in the results is therefore obtained by adopting an improved methodology for fuel surrogate targeting multiple physical/chemical properties.

Conclusions
A methodology to statistically predict knock behaviour using RANS approach is applied to an optically accessible engine unit. At first, attention is dedicated to the comparison between experimental data, provided by IM-CNR, and outcomes of the CFD analyses. Experimental data show a moderate percentage (i.e. 10%) of cycles exceeding a low predefined knock threshold; such percentage is beyond the tolerable fraction of knocking events commonly accepted in modern SI engines. The numerical methodology of statistical knock prediction coupled with a recently developed methodology to formulate surrogates provides interesting results. The choice of fuel surrogate formulation plays a crucial role in knock prediction since it largely impacts the rate of increase of mixture reactivity, as shown by the dependency of ignition delay times on surrogate formulation. Two surrogates with different numbers of components are compared. The two recipes differ in the ability to match representative physical/chemical properties of the real fuel. In addition, a non-linear scaling between the delays can be observed over a range of equivalence ratios and temperatures typical of SI GDI combustion. Differences decrease for increasing unburnt temperature, but they are relevant in typical end-of-compression / early-flame-development conditions. Such observation is of interest since it demonstrates that no scaling factor can be introduced to compensate for errors and oversimplifications in the surrogate formulation. The study also shows that fuel surrogates cannot be designed only based on analogous RON/MON values, and an extended validation over the full low-medium-high temperature reactivity is necessary. Knock analysis using ETRF surrogate as reference fuel would lead to strong overestimation of the knock signature of the investigated engine, with a presumed fraction of knocking cycles approaching 100%. This is motivated by the overestimated mixture reactivity of the ETRF surrogate, showing early knock also in the average cycle. Conversely, a six-component surrogate (SEN) seems to better simulate the mixture reactivity growth as it does not predict autoignition for the average cycle simulation, while a low occurrence of knocking cycles in the faster-than-average simulation is detected.