A RANS Knock Model to Predict the Statistical Occurrence of Engine 1 Knock 2

16 In the recent past engine knock emerged as one of the main limiting aspects for the achievement of higher 17 efficiency targets in modern spark-ignition (SI) engines. To attain these requirements, engine operating 18 points must be moved as close as possible to the onset of abnormal combustions, although the turbulent 19 nature of flow field and SI combustion leads to possibly ample fluctuations between consecutive engine 20 cycles. This forces engine designers to distance the target condition from its theoretical optimum in order to 21 prevent abnormal combustion, which can potentially damage engine components because of few individual 22 heavy-knocking cycles. 23 A statistically based RANS knock model is presented in this study, whose aim is the prediction not only of 24 the ensemble average knock occurrence, poorly meaningful in such a stochastic event, but also of a knock 25 probability. The model is based on look-up tables of autoignition times from detailed chemistry, coupled 26 with transport equations for the variance of mixture fraction and enthalpy. The transported perturbations 27 around the ensemble average value are based on variable gradients and on a local turbulent time scale. A 28 multi-variate cell-based Gaussian-PDF model is proposed for the unburnt mixture, resulting in a statistical 29 distribution for the in-cell reaction rate. An average knock precursor and its variance are independently 30 calculated and transported; this results in the prediction of an earliest knock probability preceding the 31 ensemble average knock onset, as confirmed by the experimental evidence. The proposed model estimates 32 not only the regions where the average knock is promoted, but also where and when the first knock is more


Introduction
In the last decades several approaches were proposed to numerically predict and simulate the average engine knock in the framework of RANS simulations.Their development was driven to meet the increasing efficiency targets requested by legislation.Thermal efficiency and specific output power are raised to unprecedented levels, in order both to reduce fuel consumption and pollutant emissions and to simultaneously preserve the desired target performance levels.
Under high thermal loads abnormal combustion events are promoted, the most harmful of which is engine knock.Knock is the consequence of the self-ignition of a portion of unburnt mixture ahead of the main propagating flame front, and its occurrence is enhanced by the mentioned strategies, as outlined in [1].Given the need to operate as close as possible to the theoretical optimum of the regular combustion range [2] and the simultaneous random nature of the turbulent combustion typical of internal combustion engines, the occurrence of engine knock is a possibility that is always to be accounted for when the operating condition is experimentally calibrated.To this aim, in-cylinder pressure is monitored in order to observe the random presence of knocking events.This is assessed by the definition of knock indices, such as the commonly adopted MAPO (Maximum Amplitude of Pressure Oscillations) or IMPO (Integral Modulus of Pressure Oscillations) as widely surveyed by [3,4,5], as well as other indicators like the DKI (Dimensionless Knock Indicator) [6] or the time-derivative of the in-cylinder pressure trace.Irrespectively of the chosen indicator, an arbitrary threshold value is always present to discern between a soft and acceptable knock intensity and a heavy and damaging knock level.Threshold values are part of engine manufacturer know-how and standardized limits are not defined.
All the mentioned aspects motivate why a significant research effort was paid to knock prediction in the recent years.Quasi-dimensional and three-dimensional CFD models for average knock prediction were developed and validated against experiments by several research groups.Vancoillie and co-authors used a fuel-specific Arrhenius formulation for the reaction rate of methanol and ethanol fuels in [7], and the average ignition delay was used to integrate a knock precursor species.A similar modelling approach was used by Forte et al. [8] and Corti et al. in [9] for gasoline fuels.
However, the ensemble average approach to knock modelling through the use of RANS simulations suffers of the inability to reproduce the intrinsically stochastic nature of knock; this is a strong limitation for this type of models.The dramatic impact of cycle-to-cycle variability (CCV) on all the in-cylinder physical processes, such as fuel-air mixing, combustion initiation and turbulent burn rate, motivates the adoption of more refined approaches.In fact, since engine knock depends on all the preceding processes, it is itself a typically stochastic and cycle-dependent phenomenon whose accurate prediction is therefore extremely complex.A rigorous analysis of CCV can only be carried out through the use of Large-Eddy Simulation (LES), where the largest flow structures are resolved allowing the simulation of flow unsteadiness deriving from large-scale turbulence.Despite the still demanding cost of this type of analyses, several promising studies of this kind were presented in the recent years, such as the works by Robert et al. [10,11].They showed that the simulated combustion CCV was able to replicate the degree on instability measured at the test-bench for a premixed isooctane-air engine at several spark timings.Large-Eddy Simulation was used to predict knock occurrence in a turbocharged GDI unit by the authors in previous studies [12,13,14], and the cycle-dependent knock-signature well correlated with the outcomes from the experimental test-bench for both Knock-Limited Spark Advance (KLSA) and for a knocking regime with an advanced spark-timing.
These examples showed the investigation insight made possible by LES and the possibility to explain individual misfiring cycles or cycle-specific knocking events, thus allowing a direct comparison between simulation results and engine test-bench output.However, the application of LES on production engines still suffers from the severe computational cost, preventing a full application in the design process of current units.
In this context the definition of a new approach for knock modelling emerges as a necessary bridge between the poorly representative RANS mean knock prediction and the relevant CPU effort of a multiple cycle LES study.This is based on the RANS formalism for average quantities, combined with the use of transport equations for variances of physical conditions, allowing to estimate a knock probability or a fraction of knocking cycles.The statistical RANS knock model proposed in this paper relies on transport equations for mixture fraction and enthalpy variances; detailed chemistry is used to calculate an accurate ignition delay of a gasoline surrogate model for the unburnt mixture.The variances of the variable are used as a basis of a multi-variate Gaussian model of the unburnt fluid cell, from which information on both the average reaction rate and its deviation are used to infer a presumed distribution of knocking events around the mean knock onset.An innovative definition for a probability of knocking cycles is proposed based on the same statistical basis deriving from the model equations, differently from what previously proposed by Linse et al. [15].An initial application of the PDF-knock model was presented by the authors in [16], and the application on a knock-limited turbocharged GDI engine successfully predicted 6% of knocking cycles while the ensemble average realization was knock-safe.If a traditional RANS knock model was used, a knock-safe condition would have emerged, with no further information available on knock probability and lack of correlation with the experimental acquisitions.Conversely, the use of the presented PDF-knock model gave a quantitative information regarding the presumed fraction of knocking cycles affecting the mean simulation for a given operating condition, thus enhancing the meaning of a RANS simulation of knock with a typical test-bench acquisition dataset.This initial application motivated the development of the statistical knock model and the application on a research single-cylinder engine presented in this paper.
In the next section the details of the knock model are presented, which is based on detailed chemistry tabulation to accurately reproduce the local reaction rate.A statistically-based treatment accounting for presumed turbulence-chemistry interaction is presented and transport equations for the local perturbation of the thermal and mixing state are introduced.Finally, the derivation of the mean knock precursor and of its variance are presented.The model is applied to an optically accessible GDI engine, and the knock prediction given by the presented model is compared to the experimental outcomes in terms of frequency of knocking cycles.A criterion is also proposed to correlate the results from the Probability Density Function (PDF) based knock model with the number of knocking cycles, and the potentiality and the limitations of the presented model are critically presented and discussed.

Cell-Average Reaction Rate
The first step of the presented model is based on the calculation of a cell-wise average reaction rate.A procedure for the calculation of the autoignition (AI) delay is presented in [17,18] and it is deputed to the interpolation of a cell-specific delay time from a pre-calculated database of calculated delay times .Multidimensional interpolation is carried out considering the local Favre-averaged physical conditions, i.e. the input vector  for the delay time interpolation considers the density-average values for each of absolute pressure, unburnt temperature, equivalence ratio and residuals mass fraction (Eq.1).

𝜑 = 𝜑(𝑝 ̃, 𝑇 ̃𝑢, Φ ̃, 𝑌 ̃𝐸𝐺𝑅 )
(1) The  vector lists the independent variables that govern a multiple interpolation technique whose result is the AI delay time .In this study the Andrae et al. [19] Toluene Reference Fuel (TRF) mechanism is adopted to generate a detailed look-up table of AI delays reproducing the autoignition behavior of a commercial RON95 European gasoline, corresponding to the fuel quality used in the experiments.Once the local ̃ is known, it can be time-integrated in the Livengood and Wu [20] or in the Lafossas et al. [21] models.
The limitation of the crude use of the cell-wise average ̃ value to represent the local reaction rate is the concept that a perfectly uniform value of pressure, unburnt temperature, equivalence ratio and residuals mass fraction is assumed in the cell.As a consequence, every fluid cell is considered as a laminar well-stirred reactor.Even in a RANS framework, ensemble average turbulence may affect the local thermal and mixing states around the mean value which could, in turn, induce variations in the reactivity of the unburnt charge.
Such an effect would be completely neglected by using a unique ̃ value.Therefore the model is extended to provide information regarding not only the average knock onset but also its dispersion around the mean value; details are presented in the next sections.

PDF-Knock Model Equations
Two additional transport equations are introduced to account for the statistical reconstruction of the physical states which may simultaneously be present in each fluid cell.
The   variable links the local flow turbulence to the variance dissipation rate.To this aim, it is expressed as a function of the local turbulent Reynolds number   as proposed in [23] and illustrated in Figure 1.A moderate dependence on the local   is visible in Figure 1, and a first order approximation of this formulation would be to consider   = 2.2.For low turbulent conditions (i.e.low   ), a small value of the   (  ) parameter is calculated, resulting in a long turbulent relaxation time-scale   (  ) from Eq. 4. As a consequence, the variance destruction operated by turbulence-operated mixing in Eq. 2 and 3 is slow and a high probability to find in-cell far from average states is accounted for.The opposite is verified for highly turbulent conditions.The perturbation of the local fluid state is not arbitrarily imposed, but it is derived from transport equations originating from turbulence intensity itself, hence no artificial or user-imposed variation of flow variables is introduced.
Equations 2 and 3 constitute the fundaments of the statistical treatment of the presented knock model and their application will be described in the next sections.
Finally, despite the absence of explicit spray-related terms in Eq. 2 and 3, the fuel spray affects variance fields by promoting gradients in the mean field of both  ̃ and ℎ  ̃ due to fuel evaporation and evaporative cooling, respectively.This leads to variance production, as accounted for by the first term on the RHS of Eq.
2 and 3. Also, the spray-induced turbulence acts as a local thermo-mechanical mixer, which is considered by the dissipation terms on the RHS of Eq. 2 and 3.

Statistical Description of In-Cell Reaction Rate
Both mean values and their variances are considered for enthalpy ℎ and mixture fraction .As for cellaveraged values, these are calculated from standard Favre-averaged Navier-Stokes transport equations, whereas local variance values are derived from the presented two additional transport equations (Eq. 2 and 3).
As a first step, the mixture fraction  is examined.Given the mean  ̃ and the variance " 2 ̃, a normal Gaussian distribution around the mean value in the mixture fraction space is assumed, whose spreading is represented by " 2 ̃.The probability () to find the generic  mixture quality in the cell volume is represented by Eq. 5.
If a single variable fluctuation was of interest, an analogous treatment could be carried out for enthalpy ℎ.
However, in this model both mixture fraction  and unburnt enthalpy ℎ  are affected by turbulence, as expressed by the variance values ℎ"  2 ̃ and " 2 ̃.Their joint effect has to be accounted for in each single cell and it is evaluated using a multi-variate Gaussian distribution, reported in its general form in Eq. 6 and considering a correlation coefficient  ℎ whose meaning will be described later.
A simplification of Eq. 6 could be introduced by assuming mixture fraction Z and unburnt enthalpy ℎ  as uncorrelated variables, implying  ℎ = 0 and leading to the simplified form in Eq. 7.This approximation is only used at this initial stage to describe the multi-variate model for the cell as a function of the local turbulence level.
As a final step, the mixture fraction space is converted into an equivalence ratio one and the unburnt temperature   is substituted to enthalpy through Eq. 8 and 9.
In Eq. 8   is the isobaric mixture specific heat, while in Eq. 9   is the stoichiometric air-to-fuel ratio of the fuel-air mixture and  ̃2 and  ̃2 are the mass fractions of oxygen and nitrogen respectively.From a mathematical point of view, this treatment stands as a presumed statistical reconstruction of all the possible combinations of (Φ,   ) which may exist in the cell and whose dispersion around the mean value is given by the local turbulent time scale.The expression reported in Eq. 7 represents the probability of a given (Φ,   ) state to be present in the cell volume.The probability is maximum for the mean value pair (Φ ̃,  ̃), while it is progressively reduced for far-from-average states, although it is not null for these; this anticipates the limitation of a single AI delay to describe the whole cell reactor.The statistical two-dimensional model for the fluid cell is graphically resumed in Figure 2 for different levels of turbulence, which is reflected by the amplitude of the multi-variate Gaussian model.The use of normal distribution for mixture fraction is motivated to keep the model assumptions to a minimum and to allow the use of a well-established multivariate Gaussian distribution on  and ℎ  .However, beta-distribution is another candidate choice for mixture fraction statistical representation, although it relevantly complicates the definition of a bi-variate distribution function such as the one in Eq. 6.
Figure 2. Multi-variate Gaussian-PDF distribution of physical states for the fluid cell in the presumed-PDF knock model.For the same most probable condition (Φ ̃, T ̃u), increasing turbulence intensity levels (from left to right) lead to a more effective mixing and to a probability reduction for far from average states to exist.
For the sake of numerical implementation, the multi-variate Gaussian distribution of physical conditions (Φ,   ) is discretized in an arbitrary number of physical states.In the present study, a cell-wise discretization step of half a standard deviation is adopted for both equivalence ratio and temperature, i.e. it is based on the local values of √ "  2 ̃2 ⁄ and √ Φ" 2 ̃2 ⁄ .Finally, a clipping distance from the mean value must be chosen in order to bound a finite in-cell physical space: in the present analysis a clipping is adopted at two standard deviations of each variable, i.e. ±2 • √ "  2 ̃ and ±2 • √ Φ" 2 ̃.A consequence of this is that approx.95% of the overall probability of mixture states is accounted for.Each of the two independent variables is therefore divided into 9 discrete values and the distribution of the in-cell states counts 9 2 conditions.Sensitivity analyses showed that this is a balanced compromise between the resolution of the discretized distribution and the amplitude of the simulated states.The result of the outlined procedure is a discrete clipped multi-variate Gaussian distribution.Due to the discretization operation and the boundary truncation, a final renormalization is carried out to re-normalize the sum of all the represented discrete states to unity.
The presence of the correlation coefficient  ℎ in Eq. 6 is used here to account the degree of relationship between unburnt temperature variation and mixture fraction.This is needed since this is a likely scenario in a modern GDI unit, where intense fuel stratification is observed.Charge non-homogeneity may persist until the end of the compression stroke, and it causes temperature inhomogeneity due to the dependency of the specific heat on mixture quality.The correlation coefficient  ℎ is calculated in the model at each iteration through the analysis of the in-cylinder instantaneous  ̃ and ℎ  ̃ fields and it is modelled using a Pearson-like formulation as in Eq. 10: 2  (10) In Eq. 10, the   ̃ and  ℎ ̃ terms are the standard deviations of the mean Favre-average in-cylinder  ̃ and ℎ  ̃ fields.Therefore, the   coefficient is dynamically calculated at each iteration based on the instantaneous modeled mean  ̃ and ℎ  ̃ fields, from which the spatial average 〈 ̃〉 and 〈ℎ  ̃〉 values are calculated.The  ℎ term is found to be always negative: high  ̃ cell-values (i.e.rich-in-fuel regions) are more likely associated to low ℎ  ̃.This is a consequence of the relationship between  ̃ and ℎ  ̃ deriving from the mixture isobaric specific heat.The instantaneous local  ℎ correlation coefficient modifies the bi-variate Gaussian model for the in-cell conditions as illustrated in Figure 2 for several  ℎ parameter values.In Figure 3 the equivalence ratio is used instead of  ̃ to represent fuel concentration.Hotter states are associated with leaner mixtures (bottom-right side in Figure 3) and the same is for cooler and richer conditions (upper-left side in Figure 3), while hot and rich (upper-right) or cool and lean (lower-left) combinations are less probable.

Knock Precursor Variance
From a general point of view, once a global reaction rate  ̃= ̃− 1 is known a knock precursor growth rate can be calculated in the same way as the Livengood and Wu knock integral function  ̃.This can be transported as a generic scalar with an equation as Eq.11.
The difference introduced by the presumed-PDF treatment lies in the definition of the  ̃ term.Since a variety of physical states is statistically possible in the single fluid cell, a distribution of reaction rates is also to be considered in the cell volume.The theoretical global reaction rate follows a more complex definition and it is expressed as in Eq. 12.
The calculation of the integral can be numerically challenging since an analytical function of  ̃(Φ, ) is not known a priori, hence an assumption is made considering the dispersion of the in-cell reaction rates as the sum of a mean term and a contribution due to fluctuations.This last can be either positive or negative, i.e. it can accelerate or slow down the global reaction rate depending on the considered local physical state.Since the focus of this study is on the earliest probability that a portion of a fluid cell experiences autoignition, just the faster than average part of the reaction rate distribution is of interest.
These arguments lead to the representation of the knock inceptor reaction rate  , by means of an average value  , and an accelerating contribution given by its root mean square (rms)   value (Eq.13).
, =  , +   (13) The first term on the right side is obtained by a PDF-weighted averaging operation of the population of reaction rates calculated based on the discrete multi-variate PDF space, hence the  , name (Eq.14).
The second term is calculated as the difference between the  , term and the faster than average reacting state considered in the discrete cell representation.It represents the root mean square of the faster than average reaction rates within the fluid cell.Following an analogous ±2 • √ " 2 ̃ clipping for AI delays, this reads as Eq. 15 and corresponds to the reaction rate of the 2.5% of the fastest-reacting portion.
Finally, the   term is calculated as in Eq. 16 and it represents the net increase in reaction rate due to the accelerating contribution of the statistically present faster than average states.
Once all the needed terms are available, the time integration of a knock integral function can be calculated considering Eq. 13, and it is manipulated as in Eq. 17 to split the integration of the mean reaction rate and its statistically faster reacting portion: The decomposition in Eq. 17 immediately leads to the independent calculation of two knock precursors,  ̃, and  ̃ respectively (Eq.18).The former expresses the average chemical reaction rate towards autoignition, while the latter is the precursor variance contribution given by the turbulence intensity.
Since the model aim is to track an autoignition probability for the fluid cell, the heat release due to potential knock is purposely not simulated.If autoignition heat was simulated, when AI is met for a portion of the cell due to faster than average states it would affect all the other realizations (e.g. the average knock onset) by varying the local thermo-physical conditions, while it is more interesting to transport both average and maximum probability states within the Gaussian-based model and independently track their time-history.

Experimental Apparatus and Engine Knock Characterization
Measurements were performed on a single-cylinder optically accessible DISI engine; whose main specifications are carefully detailed in [24] and here briefly resumed for the sake of completeness in Table 1 along with the operating conditions.The crank angle reference is made to the TDC at the end of compression.The engine is equipped with the cylinder head of a 1.4 litre currently made SI turbocharged power unit.The wall-guided fuel injection system features a side-mounted injector with a six-hole configuration and the spark plug is centrally located, as reported in Figure 4. Optical access is ensured through an 18 mm-thick fused silica window fixed on the piston crown featuring a Bowditch design [25] with a 45 degree UVenhanced mirror.Self-lubricating piston rings ensured oil-free operation, thus avoiding contamination of the visible field of view.More details into the application of optical techniques on this engine and other specific challenges are available in [26,27].Coolant and lubricant temperature were monitored and maintained at 330 K using a thermal conditioning unit.Engine speed was set at 2000 rpm, while start of injection was triggered at 300 CA bTDC with a single-pulse strategy at a pressure equal to 100 bar.The overall air-to-fuel ratio was set close to stoichiometry (λ≈1.05) and monitored using an oxygen sensor on the exhaust line, with an accuracy of ±1%.A turbocharged operating condition was examined in this study, with 0.5 bar boost pressure and intake manifold temperature around 315 K. Spark timing was set at 15 CA bTDC, and this also constituted the trigger for recording the optical measurements.The adoption of an instrumented GDI turbocharged unit represents an optimal testing for the presented knock model, allowing to include in the analysis the effects of stratified mixture distribution and high end-gas thermal loading common to most of the modern SI production units.A dataset of in-cylinder pressure measured (with a resolution of 0.2 CA) during 173 consecutive firing cycles were recorded through a piezo-electric transducer (that featured an accuracy of ±1%) flush-mounted on one side of the combustion chamber between an intake and an exhaust valve.The experiments were carried out using a commercial RON95 gasoline, and a knock-affected condition is observed from the cycle-resolved in-cylinder pressure derivative reported in Figure 5 (left).An arbitrary limit of pressure derivative was chosen equal to 5 bar/CA to discern between knocking and non-knocking cycles.The result of this filtering is that 109 out of 173 cycles exceed this criterion, depicting a 63% fraction of experimentally knocking cycles for this operation.The CA of knock onset for this subset of cycles is reported in Figure 5 (right), and the mean CA for knock onset for this portion of cycles is +11.9CA aTDC, with a standard deviation equal to 3.2 CA.This illustrates the 'run-away' behavior of knock, with the specific pressure oscillations being more evident for cycles at the end of the recorded dataset.the combustion visualizations even for images with low signal to noise ratio.As sketched in Fig. 6 (a-e), after the application of an appropriate circular mask (a->b) to cut the spurious light from reflections at the boundaries of the optical window of the piston crown, a look-up table (LUT) transformation [28] was used to adjust the brightness at 174, contrast at 73 and the gamma value at 0.76 (b->c).Then, thresholding was performed by fixing the minimum image intensity at 21 on 256 grey-scales (Fig. c->d).Thus, images were segmented into two regions, foreground and background respectively, obtaining a binarized image [29].

Displacement
Finally, the contours of flames were defined (d->e) and the related coordinates in pixels were stored.In order to perform a frequency map of autoignition events, only the border coordinates of the flames in the end-gas detected at fixed delay from the spark timing were considered.It should be noted that the difference between the optical window and the engine bore determined a 5.25 mm thick blind circular crown; therefore, only autoignition flames sufficiently large to cover the distance between the wall and the optical crown limit can be considered for the evaluation of autoignition maps.The experimentally measured knock probability is reported in Figure 7 through a spatial map of the normalized autoignition occurrence.The probability for the end-gases to undergo autoignition is assessed by means of the optical analysis carried out through the transparent piston window which allowed to identify the distribution of the end-gas regions subjected to knock; the peak probability was found on the exhaust side (top part in the picture), although knock is also probable on the intake side (low side in the picture).

3D CFD simulations
The 3D CFD analyses presented in this paper are carried out by means of a customized version of STAR-CD v4.22.Time varying pressure boundary conditions derive from the experiments and they are used to validate a 1D model of the engine, from which the corresponding temperature trace is extracted.Turbulence is modelled though the k-ε RNG turbulence model for compressible flows.The grid adopted for the simulations is reported in Figure 8 and it reproduces the whole combustion chamber and both the intake and exhaust ports.A close-up of the spark plug geometry and of the injector region are presented in Figure 8 as well.The total number of cells is approx.1.48 M and 430000 at BDC and TDC respectively, while the average cellsize is about 0.55 mm throughout the simulation.Combustion is modelled with the Flame Surface Density (FSD) ECFM-3Z model [30], coupled with a relatively simple algebraic ignition model based on a flame profile deposition to account for flame kernel formation [31].The fuel injector is a 6-hole full-cone GDI one whose nominal data are reported in Table 2 and whose nozzle orientation is sketched in Figure 9.The multi-hole liquid spray is modelled using a Lagrangian approach, where the fuel atomization is replaced by a Rosin-Rammler droplet distribution function.Nozzle-specific mass flow rate is prescribed as to reproduce experimental flow unbalance between the nozzles.The effective nozzle diameter is evaluated using the Kuensberg 1D model [32].The secondary break-up is modelled by the Reitz and Diwakar approach [33].Finally, spray is validated against experiments carried out in a spray bomb at an injection pressure of 100 bar; the spray morphology and the penetration curve are reported in Figure 10 at two instants after Start Of Injection (SOI), 400 ms and 800 ms respectively.In Figure 10 the penetration curve of the simulated spray is also reported and it confirms the satisfactory agreement with the experiments.

Knock Prediction Simulation
The presented knock model is applied to the simulated mean realization of combustion.In Figure 14 the unburnt temperature and the equivalence ratio field are represented at +10 CA aTDC.A detailed analysis of the end-gas condition is carried out for this CA.End-gases are identified through a conditioning on the Favre-average reaction progress variable , by filtering cells whose ̃ value is below 0.8. Figure 15 shows a scatter plot of the mixture composition against unburnt temperature, highlighting the promoted heating of the lean portion of peripheral mixture visible in Figure 14.The peak reaction rate and its average counterpart are illustrated in the scatter plots in Figure 17 for the endgas region at TDC, +10 CA aTDC and +20 CA aTDC.It is interesting to observe that the promotion of autoignition tendency given by the faster than average reacting states is present at all combustion stages, as stated by the population of fluid cells cleary lying above the 45-degree line.This is verified for relatively low reaction rates (e.g. at TDC) as well as for highly reacting conditions (e.g.+10 CA aTDC).
The field of the average knock precursor  ̃, and of its root mean square  ̃ , as calculated by Eq. 19 and 20, are represented in Figure 18.The average knock precursor field shows an evenly distributed population at both the exhaust and the intake sides, although the exhaust side of the combustion chamber (left side in figures) appears slightly more prone to knock.This is confirmed by the experiments, where knock onset locations were measured on both sides of the optical access.Furthermore, a  ̃, value below unity allows to consider the whole chamber as knock safe from an average point of view at such crank angle.Nonetheless, the magnitude of the precursor  ̃ field is relevant compared to the mean value.In particular, its additional contribution on the exhaust side points out that turbulence-induced perturbations of the reaction rate lead to more probable knock in that region of the combustion chamber, despite the average realization is not knocking at +10 CA aTDC.The degree of statistical knock tendency of the end-gases is further analyzed through the ratio of peak over average reaction rates as a function of the Favre-averaged fields of unburnt temperature and of the equivalence ratio, illustrated in Figure 19.Such analysis clearly shows that the leanest and hottest portion of the end gases, located on the exhaust side of the cylinder, are also those more subject to turbulence-induced reaction rate increase.This is accounted for by the  ℎ correlation coefficient calculated as in Eq. 10.This confirms the experimental evidence indicating the exhaust-side region as a knock critical area.Finally, the volume of the end-gas region where the autoignition criterion is met for both the average knock precursor  ̃, and the peak precursor  ̃, is illustrated in Figure 10 together with the simulated mass fraction burnt (MFB) curve.

Definition of Knock Probability in RANS
The knock prediction given by the presented model is analysed to reconstruct a fraction of knocking cycles to be compared with the experimental evidence.Knock occurrence is observed through the mass fraction of fuel which is burnt at knock onset (hereafter   ).If the cell-average ̃ value alone was used, the only information available from the analysed engine would be that the ensemble average cycle is knocking after +15CA aTDC, but no conclusions could be inferred regarding the dispersion around this value.
The additional information given by the presented presumed-PDF model is the estimation of a probability function around the mean value, which is obtained from the knock occurrence of the peak knock precursor  ̃, .This is used to define a second knock phasing indicator, i.e.  , .Given the stochastic occurrence of knock, it is reasonable to assume a normal Gaussian distribution for knock occurrence, which is centred in  , with a standard deviation    derived from  , .In the considered case,  , is equal to 83.8% while    is 39.8%.These values are based on a minimum non-null value of autoigniting volumes in the end-gas imposed as 1 mm 3 needed to avoid the spurious autoignition of individual cells to be considered.Based on these indicators a probability function of   is reconstructed following Eq.21 and it is illustrated in Figure 21.Given the above reconstruction for knock occurrence distribution, it is possible to calculate the cumulative probability of cycles exhibiting knock before the completion of regular combustion, i.e. 100% of burnt fuel.This is calculated as the Cumulative Distribution Function (CDF) of the above distribution, which is reported calculated in Eq. 22 and illustrated in Figure 22.
Figure 22.Cumulative Distribution Function of the percentage of burnt fuel at knock onset   .
Since the region of interest is limited to the probability to have a knocking cycle before the regular combustion finalization, the cumulative probability to verify this condition is given by Eq. 22 calculated at the limit value, i.e. (  )|   =100% .In this case the knock frequency value is 65.8%, stating that approximately 66% of the possible realizations reach the conditions for autoignition in a portion of their volume, while the remaining 34% is knock-safe.This result is in very good agreement with the experimental evidence showing that a fraction of 63% of the measured cycles is knocking.This information is inferred from a single RANS simulation, while the average knock prediction alone would be unable to estimate any knock dispersion around the mean value.
It is important to underline that the aim of the presented model is neither to substitute a multiple cycle LES simulation, which remains the only way to properly simulate most of the CCV-promoting processes, nor to identify the exact fraction of knocking cycles measured during the experiments.The presented presumed-PDF model is a numerical tool to identify the regions of the combustion chamber which are statistically more prone to autoignition, and this information is inferred by a combination of an average knock precursor and a statistical description of the reaction rate deviation produced by local gradients and dissipated by local turbulent intensity.Since RANS simulations are the most appropriate tool to investigate the mean behaviour of a fluid system (e.g. a new engine concept), the indications from the presumed-PDF model may provide useful statistical information on the probability to trigger potentially damaging knocking events.

Conclusions
In this paper a statistics-based knock model is presented in the context of RANS combustion simulation, which couples RANS traditional equations with the transport of variances for the physical conditions affecting local reaction rate.The model is based on separate transport equations for both the mixing and the thermal variance originated by local mean gradients and turbulent scales.They are combined to reconstruct a statistical model for the in-cell reaction rate, which is represented though a clipped multi-variate Gaussian distribution of probability.Two independent knock precursors are transported, in order to account for both the ensemble average knock proximity and its variance around the mean value.The combined use of the precursors is able to outline a statistically-based description of the in-cell reaction rate distribution.
When applied to a laboratory GDI engine with optical access, the results from the presumed-PDF knock model assesses knock onset on the exhaust side of the combustion chamber.This is the area where flame visualization indicated the highest knock probability, and numerical simulations showed that this region suffers of lean and hot end-gases, promoting knock onset.Moreover, turbulence-induced dispersion promoted reaction rate increase of the unburnt mixture in this region, making it also the highest probability knock location.An overall fraction of 66% of knocking cycles is calculated, which is in close agreement with the experimental 63%.equations allows to quantify the probability of knocking events given by turbulence-originated fluctuation of end-gas pockets, and it allows RANS simulations to be directly correlated with engine test-bench acquisition.

Figure 1 .
Figure 1.Variance dissipation rate   parameter as a function of the local turbulent Reynolds number   .

Figure 4 .
Figure 4. Experimental apparatus for the optically accessible SI engine.

Figure 5 .
Figure 5. Cycle-resolved in-cylinder pressure derivative (left) and CA of knock onset for the subset of knocking cycles (right).Cycle resolved visualization was performed for 100 consecutive cycles out of the 173 via a CMOS camera (Optronis CamRecord 5000 -512 x 512 pixel, 8-bit pixel digitization and 5000 frame per second at full chip) equipped with a 50 mm focal Nikon lens.The exposure time was fixed at 166.7 μs and the dwell time between two consecutives images was 200 μs (2.4 CAD at 2000 rev/min).A custom designed image processing using Vision Assistant of National Instruments allowed to retrieve quantitative information from

Figure 6 .
Figure 6.Sketch of the image processing steps.

Figure 7 .
Figure 7. Map of normalized autoignition location (left) and CA of knock onset for the subset of knocking cycles (right).

Figure 8 .
Figure 8. Computational grid (top) and spark plug and injector region detail (bottom).

Table 2 .
is modelled by means of the presented presumed-PDF knock model, which is coupled with the STAR-CD solver through in-house developed user-coding.Nominal injector data.

Figure 9 .
Figure 9. Sketch of the nozzle holes position and orientation.

Figure 10 .
Figure 10.Simulated spray shape against experiments (left side) at 400 ms (top row) and 800 ms (bottom row) after SOI; comparison between experimental and simulated penetration curve (right side).

Figure 15 .
Figure 15.Scatter plot of equivalence ratio against unburnt temperature for the end-gas region at +10 CA aTDC.The field of average AI delay times and of its root mean square are reported in Figure16.While the former field (left side), related to the gas phase average reaction rate (expressed as the inverse of the local AI delay time), shows local minima at the periphery of the intake valve and around the injector cavity, the latter (right side) visually suggests that local turbulence induces a more relevant degree of reactivity fluctuation on the exhaust side of the combustion chamber.

Figure 16 .
Figure 16.Average AI delay time field (left side) and delay root mean square field (right side) at +10 CA aTDC.

Figure 19 .
Figure 19.Scatter plot of the ratio of peak over average reaction rate for the unburnt mixture at +10 CA aTDC as a function of the unburnt temperature (left) and of the equivalence ratio (right).

Figure 20 .
Figure 20.Autoignition volume as predicted by the average  ̃, and the peak  ̃, knock precursors (black lines), together with the ensemble average MFB curve (red line).

Figure 21 .
Figure 21.Assumed Gaussian-PDF of the percentage of burnt fuel at knock onset   .
The use of a PDF-based model allows to give a quantitative estimation of the knock probability events associated with an average knock occurrence condition.The statistical model is based on transport equations for the variance of physical conditions (enthalpy and mixture fraction) and tracks their variability in terms of knock onset variance.The presented PDF-based knock model combines both the accuracy of a detailed chemical mechanism used to calculate the autoignition delay times and the applicability of a RANS-based model.The presented model aims at filling the existing gap between the scarcely representative average knock prediction given by traditional RANS models and the cycle-resolved knock simulation possible with CPU-intensive multi-cycle LES simulations.The combination of RANS models with variance transport

Table 1 .
Single cylinder optically accessible SI engine characteristics.