Abstract
Considering the strict regulations on the transport sector emissions, predictive models for engine emissions are essential tools to optimize high-efficient low-emission internal combustion engines (ICE) for vehicles. This aspect is of major importance, especially for developing new combustion concepts (e.g. lean, pre-chamber) or using alternative fuels. Among the gaseous emissions from spark-ignition (SI) engines, unburned hydrocarbons (uHC) are the most challenging species to model due to the complexity of the formation mechanisms. Phenomenological models are successfully used in these cases to predict emissions with a reduced computational effort. In this work, uHC phenomenological model approaches by the authors are further developed to improve the model predictivity for multiple variations including engine design, engine operating parameters, as well as different fuels and ignition methods. The model accounts for uHC contributions from piston top-land crevice, wall flame quenching, oil film fuel adsorption/desorption and features a tabulated-chemistry approach to describe uHC post-oxidation. With the support of 3D-CFD simulations, multiple novel modelling assumptions are developed and verified. The model is validated against an extensive measurement database obtained with two small-bore single-cylinder engines (SCE) fuelled with gasoline-like fuel, one with SI and one with pre-chamber, as well as against data from two different ultra-lean large-bore engines fuelled with natural gas (one equipped with a pre-chamber and one dual-fuel with a diesel pilot). The model correctly predicts the trends and absolute values of uHC emissions for all the operating conditions and the engines with an accuracy on average of 11.4%. The results demonstrate the general applicability of the model to different engine designs, the correct description of the main mechanisms contributing to fuel partial oxidation, and the potential to be extended to predict unburned fuel emissions with other fuels.
Keywords
Introduction
The internal combustion engine (ICE) has and will have a central role in the transportation sector in the coming decades,1–5 especially in combination with non-fossil fuels 2 and in the so-called ‘hard-to-electrify’ or ‘hard-to-abate’ sectors.6,7 Virtual engine development tools have become essential to effectively design new engine concepts and tackle efficiently the associated engineering targets. 8 To contain the environmental impacts of ICE, the research is strongly focused on the reduction of exhaust emissions. This applies to non-polluting species contributing to global warming like carbon dioxide (CO2), and to pollutant emissions like mainly soot, carbon monoxide, nitrogen oxides and unburned fuel, that is, unburned hydrocarbons (uHC) in the case of hydrocarbon-based fuel. Regarding CO2 reduction, research is focused on the design of new combustion concepts that can achieve high efficiencies possibly in combination with future-relevant renewable non-fossil fuels (e.g. green methanol, hydrogen). To tackle the pollutant issue, the new developed combustion concepts are optimized to generate low engine-out emissions in order to achieve ‘zero-impact emissions’ 9 in combination with a highly efficient after-treatment system. In this context, reliable exhaust emission models allow to optimize engine designs and minimize pollutant emissions for new engine concepts.
Among the gaseous engine emission species, unburned hydrocarbons are the most challenging species to simulate because of the heterogeneity and complexity of the multi-physical mechanisms contributing to their formation. Indeed, uHCs result from multiple three-dimensional (3D) phenomena taking place in the combustion chamber like the fuel injection, spray-wall interactions, oil and fuel films, complex geometry, flame-wall interaction and chemical kinetics. These phenomena are impossible to be fully modelled all together with the complete physical description even in physics-based and accurate 3D-CFD simulation. 10 Phenomenological models offer a valuable possibility to estimate uHC contribution from physical phenomena at different engine operating conditions. In particular, the possibility to implement this kind of models within engine 1D/0D simulation environments offers the chance to refine the engine design and verify it on the whole operating map with reduced computational effort.11–13
This work builds-up on the authors previous studies.11–13 In the work from Esposito et al.,11,12 extensive validation of an uHC model over a large experimental database from an SI direct injection (DI) small-bore SCE has been presented. A good overall capability of the uHC model has been demonstrated at stoichiometric and warm-engine operation. However, underestimation of uHC engine-out emissions has been shown at lean and cold-engine conditions. On the other hand, the models from Malfi et al. 13 have been tested only on a reduced dataset from ultra-lean large-bore engines fuelled with natural gas.
In this work, previous approaches, knowledge and data gathered from the authors11–13 have been merged and further developed to introduce new modelling features and enhance the accuracy of uHC phenomenological models as well as enlarge the prediction capability over relevant engine operating parameters and engine architectures.
Sources and modelling of uHC
Unburned hydrocarbons are a result of hydrocarbon fuel incomplete oxidation. The mechanisms leading to unburned hydrocarbon production in premixed-charge engines (e.g. SI, dual fuel, or pre-chamber concepts) are fundamentally alike and originate from the partial escape of a part of the fuel from the main flame propagation. According to a vast literature focused on the investigation and quantification of the most relevant sources of uHC,14–18 the main mechanisms contributing to uHC emissions are 12 : combustion chamber crevices, wall flame quenching and oil film. The process leading to the uHC engine-out emission values can be subdivided into ‘formation’ (escape of unburned fuel from the main combustion phase according to the above-mentioned sources) and ‘post-oxidation’ (inside the cylinder, the uHC released in the combustion chamber, once mixed with the hot bulk gases, can be partially oxidized). 12
Methodology
Experimental data
In this study, data from four different engines (A, B, C, D) are considered for the validation of the developed model. The engines have different characteristics, sizes, combustion systems and were operated with different fuels. Engines A and B are both small-bore single-cylinder engines (SCE) with very similar bore-to-stroke ratios, but different compression ratios, fuelling and ignition systems. Indeed, engine A is a conventional DI SI engine designed to operate stoichiometric with reduced enleanment capabilities in this setup, while engine B has an active pre-chamber (PC) fuelled with gasoline designed to reach ultra-lean conditions. Engines C and D are large-bore engines for marine or power generation applications, designed to operate at ultra-lean conditions and with natural gas (NG) port-fuel injection. Both engines have been tested with two different piston designs (different piston top-land heights), labelled as ‘low top land’ (LTL) and ‘high top land’ (HTL). 13 While engine C ignites the charge with a pilot diesel injection, engine D is equipped with an active PC with an auxiliary NG fuel injection. A summary of the considered engines is reported in Table 1. This includes also the previously published works where the experimental methodology for each engine has been explained in detail.
Main characteristics and data of the four engines considered in this study.
RON95E5 & a corresponding surrogate (toluene, iso-octane, n-heptane, ethanol). 48
Values refer to LTL and HTL piston designs.
The main development of the uHC model has been performed on the base of engine A because of the availability of a very extensive engine database, and 3D-CFD simulations.
All the considered measurements have been performed at stationary conditions, by setting a certain level of engine load in terms of indicated mean effective pressure (IMEP), engine speed (n E ), relative air-fuel ratio (λ), coolant and oil temperature (T cool = T oil ). For engine A, the extensive emission database consists also of emission values at different sampling locations in the exhaust runner of the SCE. In total, for each selected operating point, 1000 consecutive cycles were recorded. Figure 1 reports a schematic of the Engine A SCE with the emission measuring positions relevant for the evaluation of uHC emissions in this study. Both values will be used for comparison against the simulated results.

uHC measurement positions: the distances are approximately d1 = 120 mm and d2 = 400 mm.
Simulation model
3D-CFD RANS simulations
The 3D-CFD simulations considered in this study are the same as presented in Ref., 12 performed according to the methodology explained in Ref. 10 These 3D-CFD RANS investigations have been carried out with the software CONVERGE CFD v2.4 mainly to investigate the uHC formation mechanism from the piston top-land crevice. The data has been then post-processed to provide information to support and validate the modelling hypotheses and approach for the phenomenological uHC model. In particular, a geometrical post-processing of the 3D results from CFD has been coded in MATLAB and used to extract the trends of variables of interest in specific zones within the cylinder domain. A schematic and overview of the different zones evaluated in the post-processing is provided in Figure 2.

Illustration of the zone post-processing of the 3D-CFD results with details of the piston top-land crevice zones (dimensions are given in mm). Reworked from Ref. 12 with permission of SAE International.
In Figure 2, z c0 represents the base of the piston top-land crevice while z c1 and z c2 are the lowest and highest piston height respectively. The code for the zone post-processing uses as input the 3D-CFD cell-based results at a specific crank angle, the zones are identified based on cell coordinates, and the calculation of average quantities of interest is done in three separate zones around the TL crevices as depicted in Figure 2. All of them have ring-like forms. The zone ‘Crevice1’ consists of the lowest and narrowest zone of the TL crevice. The zone ‘Crevice2’ is directly above Crevice1 and reaches the piston shoulder. Since the piston is not flat, Crevice2 is a real crevice zone only for approximately 50% of the piston circumferential length. On the yz plane of Figure 2, it can be observed how Crevice2 really represents the gap between the piston and the liner while on xz plane, Crevice2 is a crevice outlet zone. The zone ‘Crevice-out’ is above Crevice2 and as wide as its outlet and it extends for 3 mm above the highest piston surface and allows the evaluation of the crevice outflow and oxidation thermodynamic state. The operating point simulated and considered for the development of the uHC phenomenological models are summarized in Table 2 and assume warm engine operation (T cool = T oil = 90°C).
Engine A operating conditions simulated in 3D-CFD.
Base 1D engine models
The GT-SUITE models used for the implementation of the uHC modelling include the geometry of the combustion systems between the low-pressure indication transducers, respectively in the intake and the exhaust runners as explained in Refs.11,12 Measured instantaneous pressures and average temperatures are used as boundary conditions, as well as the rotational speed. The injected fuel is controlled by a PID controller to match the experimental air-to-fuel ratio. More details about the base model settings for engine A can be found in Refs.,11,12 for engine B in Ref., 49 and for engines C and D in Ref. 13
The combustion model together with the uHC sub-models are both implemented in the GT-SUITE framework under the form of ‘user coding’. The combustion chamber is described by a two-zone approach, where the burned and unburned gases are subdivided into two zones at different temperatures and at the same pressure. To focus only on the predictivity of the uHC models, a combustion profile obtained from the analysis of the measurements through a three-pressure analysis (TPA) has been used to describe the combustion for the selected operating points. The uHC model is coupled to the combustion model by a burn rate manipulation process dependent on the release and post-oxidation uHC phenomena. 13 The heat transfer calculation is managed using a modified Woschni correlation 19 identified in GT-SUITE as WoschniGT with assumed wall temperatures for engine A as described in Ref., 11 3D-estimated temperatures for engine B and experimentally-derived temperatures for engines C and D.
Unburned hydrocarbon model
Among the main uHC mechanisms, 12 this work focuses on the piston top-land filling and emptying mechanism,50–55 the wall flame quenching,34,56 the adsorption and desorption of fuel from the oil layer42,57 with the addition of a novel approach to estimate the fuel direct entry in the oil film for engine A.
In the following, both uHC formation and post-oxidation processes will be described for the four different formation sub-models considered in this study.
where E is the total energy in the crevice volume, δQW/dt is the wall heat transfer, dmcrev/dt is the mass flow rate from the cylinder, and h flux represents the enthalpy of the volume where the flow flux originates. The mass and species balances are:
where Ycrev,k represents the mass fraction corresponding to the k-th species transferred between the cylinder and the crevice volume. Regarding the total internal energy, it can be calculated by multiplying the mass and the respective specific internal energy:
Therefore, if the mixture is considered as an ideal gas, the internal energy can be written as:
In which c v represents the specific heat at a constant volume which generally relies on both the gas temperature, T, and its composition, x k . Once assigned the crevice volume, V tl , the gas density can be also evaluated. The thermodynamic state of the gas trapped in the crevices is univocally determined by the knowledge of instantaneous density and specific internal energy. Using the gas state equation, the pressure in the top land crevice can also be estimated.
The heat transfer coefficient, σ crev , in the crevices is calculated through the definition of the Nusselt number, considering that the TL geometry predominantly extends in a specific direction:
with k crev representing the thermal conductivity of the trapped gas in the crevices, and D crev indicating the crevice’s hydraulic diameter, defined as:
Here, A crev represents the cross-section area of the crevices, taken perpendicular to its primary dimension, and P crev is the cross-section wetted perimeter. The Nusselt number is calculated with the empirical correlation proposed in Ref. 58 The solution of mass and energy balances within crevices allows to estimate the instantaneous profiles of temperature and pressure, reported in Figure 3. The pressure trend is almost equal to that in the cylinder, while the temperature evolution strictly depends on the size and geometry of the crevices volume. The gas trapped in the TL crevice follows the unburned gas temperature evolution, with some delay and dimming related to the filling and emptying process and to the wall heat transfer.

Exemplary model results for Engine A at operating condition 2500@3 - λ1 (cfr. Table 2): (a) average cylinder temperature with unburned and burned gas as well as crevices’ gas temperature and (b) cylinder pressure and pressure in the crevices.
Usually, in literature, the crevice temperature is fixed at a constant value equal to that of the piston 30 or liner, or an average of the two temperatures. In a previous work, 12 the crevice temperature was calculated over the crank angle as a function of the temperatures of piston, liner and unburned gas with a specific formula that has been validated against 3D data. In this work, an alternative and more general approach is implemented in which the crevice temperature is calculated during the entire engine cycle from the energy balance equation. The calculated temperatures have been compared to the ones obtained through the 3D simulations to prove the approach reliability. In particular, both the Crevice1 and Crevice2 zones evaluated from CFD with geometric post-processing described depicted in Figure 2 are taken as reference. In Figure 4, the temperature comparisons have been depicted for all the available data. The crevice temperature from Esposito et al. 12 (black dotted line) reproduced a temperature profile similar to Crevice1 (cfr. Figure 2), which is the deeper and colder part of the TL crevice. The approach adopted in this study (black continuous line) overall shows similar temperature profiles to Ref. 12 with slightly higher levels in correspondence to the pressure peak, which can be considered more realistic considering also the temperature of the Crevice2 zone in CFD.

Comparison among the 3D-CFD average temperatures in the zones Crevice1 (dashed grey line) and Crevice2 (dashed-dotted grey line), the simulated crevice temperature Tcrev in this study (black continuous line) and calculated temperature Tcrev from Esposito et al. 12 (dotted black line) for all the operating points simulated in 3D-CFD.

Combustion chamber areas schematization for wall flame quenching.
In Figure 5, representative flame configurations are depicted during four different instants of the flame progression. In the initial stage, the flame remains entirely detached from the cylinder walls (1); as it expands, a section of the flame front touches the cylinder head (2). Subsequently, the flame extends to reach the piston crown (3) and finally, the cylinder liner (4). The quenched volume variation can be estimated by multiplying the i-th combustion chamber wall areas (head, liner and piston) reached by the flame front in a time step by the quenching distances, δ q,i = Peq,i· δ fl , which is evaluated through the Péclet number (Peq,i) definition. Among all the terms cited above contributing to the evaluation of the quenching contribution, detailed in previous author work, 13 a focus on the Péclet number estimation is posed. This term can be related to the normalized wall heat flux, φ q , according to the correlation presented by Boust et al.56,60:
In the model presented, the normalized wall heat flux φ q is held constant at 0.3, a value recommended in the situation of a laminar head-on quenching at ambient conditions. 34 The wall temperature Tw,i varies according to the specific region reached by the flame front, whether it is the cylinder head, liner, or piston. To correct the inaccuracies of the presented correlation, a tuning constant, c q , is introduced in equation (8). By combining the quenching distances δ q,i , and the area subject to the flame quenching, dAq,i, with the gas density, it is possible to calculate the finite quenched mass in a time step. As a result, this mass, dm q , is integrated over the engine cycle to compute the uHC contribution from wall quenching.
To quantify this additional uHC source, 3D-CFD calculations available for some representative operating points10,12 for engine A are post-processed and the maximum percentage of fuel impinging the liner walls is evaluated and assumed as the amount of fuel that enters the oil film at cold engine conditions. Subsequently, based on 3D-CFD data, a correlation to estimate the injected fuel percentage forming the liquid film is developed depending on engine IMEP and λ values:
Figure 6 reports the agreement between the fraction of injected fuel impinging the liner as evaluated in 3D-CFD and the estimation through equation (9). It can be observed that the amount is in the order of magnitude of 1%–4% and the developed correlation reproduces the 3D-CFD trend quite accurately for the different operating points.

Estimation of fuel direct entry in the oil layer: comparison between 3D-CFD estimation and developed correlation from equation (9).
In this work, the boundary layer (bl) state is individually computed for various uHC-forming regions, namely crevices, wall quenching and oil layer. The calculations in those regions adopt a common mathematical approach. Regarding the uHC emissions from the TL crevice, the boundary layer represents a volume near the cylinder wall accumulating the mass exiting from crevices during the expansion phase. During this phase, the mass released from the TL crevice mixes with a part of the burned gas originating from the bulk burned gas region of the cylinder.
The temperature and volume of the boundary layer are evaluated on a crank angle basis by solving specific mass and energy balances expressed as:
Here, m bl represents the mass contained in the boundary layer, ṁ crev the mass flow exiting from the top land crevice, ṁ b is the entrained burned gas mass flow rate, e bl is the specific internal energy of the boundary layer and V bl is its volume, p cyl represents the in-cylinder pressure, and h crev , h b and h bl are the specific enthalpies of top land crevices, burned gas and boundary layer, respectively. ṁ b characterizes the burned gas flux coming from the bulk region of the cylinder. This is defined to be proportional to the mass flux coming from the crevices through the tuning constant k oxid (koxid_crev for the TL crevice):
The energy balance for the boundary layer allows to estimate the temperature evolution for this cylinder region. In Figure 7, the calculated temperature trends for the crevices boundary layer are compared to the post-oxidation temperature modelled in Esposito et al. 12 and the temperature evaluated in the Crevice-out zone from 3D-CFD simulations (see Figure 2). In the proposed model, the post-oxidation process initialization corresponds to the rise of the boundary layer temperature. Hence, this happens after the top land crevice begins to release gases, after the cylinder peak pressure. This explains the absence of the BL temperature trend in the early stage of the expansion stroke, represented by the grey bands in Figure 7. The rest of the diagram represents the time at which the post-oxidation process occurs and hence the boundary layer region in the phenomenological model aims to reproduce a similar temperature. Indeed, it can be observed that the BL temperature quickly rises reaching a level that is in line with that in the Crevice-out zone from 3D-CFD simulations and similar to the post-oxidation temperature calculated in Ref. 12

Comparison among the calculated boundary layer temperature for the TL crevice developed in this work, the average 3D-CFD temperature in the zone Crevice-out, and the previous approach to compute the TL crevice oxidation temperature from Esposito et al. 12 for all the investigated points with CFD. Grey bands highlight the instants before post-oxidation (and top-land crevice boundary layer model) activation.
As mentioned previously, the determination of the boundary layer state for the uHCs generated by wall flame quenching and absorption-desorption from oil follows a comparable approach to that used for the TL crevice. Concerning the wall flame quenching, one difference lies in the replacement of the production term ṁ crev , attributed to the TL crevice, with the quenched mass rate, denoted as ṁ quench , derived from the model detailed in the section above. Another variation consists in the quenching boundary layer which is assumed to not mix with the burned gas (ṁ b = 0 in equations (10) and (11)), and the tuning constant k oxid , in this case (koxid_quench), represents a weighting factor used to define the temperature of the gas in the boundary layer, whose level is between the temperatures of combustion chamber walls and burned gas, according to:
The approach presented in equation (13) is different from the previous boundary layer model, 13 but similar to that presented in Ref. 12 with the difference that the weighting factor there was between 0.2 and 0.3. The choice to have a separate boundary layer computation for the wall flame quenching arises from the necessity to account for a different release timing and thermodynamic state of the quenched mixture compared to those of the mixture released from the TL crevice. Indeed, the wall flame quenching happens in a phase of the expansion stroke with high and rising temperatures. Conversely, the TL crevice mixture is released after the peak pressure and encounters comparatively lower and decreasing temperatures. Adopting the same formulation for the boundary layer temperature would result in an unrealistic complete oxidation of the quenched mixture due to the high temperatures. Physically, to undergo post-oxidation, the quenched mixture needs to diffuse into the burned zone and this process has been proven to be a limiting factor for the post-oxidation of quenched layers. 64 To account for this slow diffusion, Esposito et al. 12 introduced a factor to strongly limit the post-oxidation of the wall flame quenching. In this work, koxid_quench accounts for this diffusion delay of the quenched layers into the hot burned gases. As far as the oxidation of uHC released from the oil layer is concerned, an approach similar to that adopted for the uHC oxidation from wall flame quenching is applied. This is controlled by the tuning constant koxid_oil, which acts as expressed by equation (13). A summary of the of the boundary layer calculations is schematized in Figure 8.

Sketch of the boundary layer areas and calculations considered in this work.
As an advance of the existing models, a tabulated kinetic approach is here adopted for the oxidation process. This allows considering the exact conditions present in the boundary layer, both the thermodynamic state and the composition, improving the evaluation of the oxidation rate and allowing the application of the model within the thermochemical conditions used for the tabulation.
Similarly to previous works of the authors, 65 the chemistry tabulation is synthesized in a progress variable (and its derivative). This is defined as the formation enthalpy variation of the evolving mixture, computed with respect to a reference temperature of 298 K. The formation enthalpy (h f ) evolution is normalized by its variation between the final and beginning stages of the reaction process:
The progress variable time derivative, ċ i , is tabulated in offline chemistry simulations in an adiabatic homogeneous reactor at 99 different and prefixed values of the progress variable, c i , for various initial conditions of temperature, pressure, equivalence ratio and residual gas content. During engine simulations, the C-atom-based fuel equivalent oxidation fraction (r ox ), as already presented in previous work, 12 is calculated according to:
where Δt is the simulation time step and values of
Results
Model calibration and sensitivity analysis
Before presenting the results of the uHC model, this section describes how the model can be tuned and what are the effects of its constants. Based on the description provided in the previous paragraphs, the constants to be selected for each engine are: c q from equation (8) for the Péclet number estimation, and the three oxidation constants for TL crevice, wall flame quenching and oil layer, respectively (koxid_crev, koxid_quench, koxid_oil). To show the effect of the calibration constants on the model results, a sensitivity analysis is presented in Figure 9 for Engine A. In particular, the trend of the released uHCs (before post-oxidation) and post-oxidized concentration values are shown over crank angles for all the sub-models, (a) for the TL crevice, (b) for the wall flame quenching and (c) for the oil layer. Additionally, also the r ox variable, which is evaluated in the boundary layer of each sub-model, is reported (right y-axis).

Instantaneous trends of uHC concentration before (released, dotted black line) and after the post-oxidation process (full black line) together with the rox variable that describes the post-oxidation (red dashed line). (a) Piston TL crevice, (b) wall flame quenching and (c) oil layer adsorption/desorption. Engine A: nE = 1500 min−1, IMEP = 3 bar, λ = 1.1.
Figure 9(a) shows that the TL crevice uHCs are released after the peak pressure when the combustion is almost complete, and oxidize intensely only for about 20° CA. This mechanism limits the production of uHC from this source until about 35° CA aTDCF. Moreover, it is visible that the uHCs from flame wall quenching (Figure 9(b)) start to accumulate first, with a quite weak post-oxidation rate. The release of uHC from the oil layer (Figure 9(c)) also starts when the combustion ends, having an almost constant production rate and quite low levels.
A sensitivity analysis is superposed to all the representative trends shown in Figure 9. All four model constants are varied by ±10%. The variation of c q is shown in the sub-plot (b) and highlighted through the light blue scatter band. For the k oxid constants, the grey scatter bands highlight the sensitivity related to the ±10% variation. It can be observed that only the wall flame quenching model is quite sensitive to the variation of the model constants (namely c q and koxid_quench) with strongly nonlinear behaviour due to the complex and high-temperature sensitivity of the wall quenching and its post-oxidation. The TL crevice sub-model shows approximately ±6% on the final uHC concentration with ±10% of koxid_crev, while the oil layer sub-model is almost insensitive to the variation of koxid_oil. This is due to the very late release of uHC when the temperature is too low for oxidation.
Regarding the tuning procedure, some steps are applied to identify the most appropriate set of constants, which are based on the following general observations. For the operating points at stoichiometric conditions and especially at high loads, the uHC production is expected to be dominated by the contribution of the crevices. The mass exchange between the top-land crevice and the cylinder is only determined by the instantaneous in-cylinder pressure and the geometric characteristics of the crevice volume, and no calibration parameters can control this phenomenon. Therefore, the final uHC contribution from crevices can be varied only by changing the intensity of the post-oxidation process through the constant koxid_crev, whose effect is depicted in Figure 9(a). This constant has been calibrated with a focus on the high-load points at stoichiometric conditions and kept constant for all other operating points. In the second step, the uHC production due to flame wall quenching is controlled by the c q constant. This parameter allows increasing the quenching layer thickness at the wall and, as a consequence, it raises the quenching uHC emission especially at low loads and at diluted conditions. Then, the tuning on the post-oxidation of this uHC contribution through the constant koxid_quench as highlighted in Figure 9(b) allowed to improve the trends of uHCs, especially over λ. Finally, concerning the production of uHC from the oil layer, no actual tuning is applied, and also its post-oxidation is generally marginal whatever the value of koxid_oil as demonstrated in Figure 9(c).
A summary of all the optimal values selected for the model constants for the considered engines is reported in Table 3. The model constants are kept fixed in all the operating conditions of each engine. It is worthwhile considering that the engines under study are very different in terms of dimensions, type of fuel, air/fuel proportions, mixture ignition mode and crevice size. Regarding the variability in c q , this seems correlated with the mixture air-fuel ratio because engine A, which runs stoichiometric or close to stoichiometric, has a lower value than the other ultra-lean engines. Further differences among engines B, C and D can be related to the flame propagation in combination with the different geometries. Regarding koxid_crev, this looks strongly correlated with the TL crevice volume reported in Table 1. It is interesting to note that engines C and D, which have similar designs and the same main fuel (natural gas) have the same values for koxid_crev and koxid_quench, while the difference in c q can be connected to the differences in the mixture ignition mode (pre-chamber fuelled with NG and diesel pilot, respectively) and consequent flame propagation. Both engines have no sub-model for oil layer adsorption/desorption because of the very high solubility of natural gas in oil, resulting in an early uHC release during expansion and an almost complete post-oxidation.
Model constants selected for the considered engines.
uHC prediction results
Figure 10 shows the comparison between predicted and measured cylinder-out uHC concentrations over the extensive database available for Engine A with variations in engine load, λ, valve timing and spark timing. The last load variation, λ variation and the spark timing variation are at cold engine conditions, where the coolant and oil temperature are maintained at 30°C instead of the standard level of 90°C, typical of warm engine operation. The spark timing variation at cold engine conditions represents the so-called ‘catalyst heating operation’ which has spark-timing retardation from optimal to very delayed combustion (until 10% coefficient of variation of IMEP is reached) at low load and low speed (nE = 1200 min− 1, IMEP = 3 bar, λ = 1.0). uHC measurement data (black lines) refers to two measurement positions along the exhaust line (see Figure 1). In particular, the overall model results (red continuous line) are shown together with the individual sub-models contributions (dashed blue line for TL crevice, dashed-dotted orange line for wall flame quenching, dashed-dotted green line for oil layer). For reference, the red dotted line shows the predicted values of the model from Esposito et al. 12 to show the improvements achieved in this study. It is quite visible how the prediction results follow the trends and the absolute levels of the measured uHC for all the operating conditions. As an example, the model predicts the reduction of uHC production at increasing load and speed (with stoichiometric air/fuel ratio), which is mainly due to the more intense post-oxidation process caused by higher in-cylinder temperatures. The model also detects the increase of uHC concentration in the case of lean air/fuel mixtures, which is related to the growth of wall flame quenching intensity (thicker quenching layer). The contribution of the oil layer is marginal in most cases but becomes essential to enhance the model predictivity at cold engine conditions. The simulation is also able to predict the trends and levels of uHC at catalyst heating operation, which is a very extreme operating condition. Compared to previous studies,11,12 the actual model is able to capture better the trends over IMEP and λ, and above all, reaches very high accuracy at cold engine operation. This aspect is of major importance since the uHC prediction especially at a cold start is of central interest. It can be noted that for the EVC variations (lower row of Figure 10), the model shows a tendency to overestimate uHCs and it does not seem to perceive the effect of EVC retardation (increase of valve overlap). This might be due to 3D gas exchange and mixing effects with residual gases, which are not described in the current model. Globally, the presented model demonstrates to be accurate with an average deviation of 11.7% considering all the operating conditions of Figure 10.

Experimental/simulation comparison of cylinder-out uHC concentration of Engine A. Upper row: load and λ variations; lower row: EVC and spark advance (SA) variation at catalyst heating conditions. Dotted red lines are the total predicted uHC from Esposito et al. 12 reported to highlight the improvement in the results achieved in this work.
Figure 11 presents the results of the uHC model for engines B, C and D. In this case, only one measurement dataset is used for comparison. As previously, single contributions of the uHC sub-models are reported in the figure. As far as engine B is concerned, a λ sweep at fixed engine speed and load is used for validation. Higher deviations are observed only at the first two richer points at λ=1 and 1.1, while the rest of the points show a very high accuracy of the predictions. On average, the deviation is 12.7% while it drops to 4.4% excluding the first two points. It is interesting to note how the trend of increasing λ is reproduced as a result of the quenching contribution.

Experimental/simulation comparison of cylinder-out uHC concentration of Engine B (λ sweep), C and D (load sweep).
As far as engines C and D are concerned, the results consist of load variations at ultra-lean conditions for both ‘high top-land’ (HTL) and ‘low top-land’ (LTL) piston ring-pack configurations. The model predicts the increase of uHC production when the load decreases. This is mainly caused by an increase in the flame quenching layer thickness at the cylinder walls depending on the higher residual content, the mixture air-fuel ratio and on the lower cylinder pressure and temperature levels. On the other hand, the uHC production from the top-land slightly decreases with the load reduction, since a smaller amount of fuel is trapped in the crevice region when the in-cylinder peak pressure is lower. The simulation detects the effects of the ring pack dimension up to a certain extent since the piston design with the smallest top land height (LTL) leads to slightly reduced uHC emissions as shown in the experiments. The accuracy of the predictions is very high, and the average deviation observed is about 10% for both engines C and D and all the considered operating conditions.
As a final remark, it is worthwhile observing that the uHC model performs quite well for all studied engines, adopting similar calibration constants (see Table 3). This is a very important point in light of the substantial differences among the engines in terms of size, supplied fuel, rotational speed, air/fuel ratio and ignition system. This also allows to foresee the possibility of successfully applying the uHC model to further engines and fuels to evaluate their unburned fuel emissions in a pure predictive mode.
Assessment of model upgrades
In this final section, the advantages of the presented model improvements are highlighted. These can be summarized in three points:
1) Fuel inclusion in the oil layer at cold conditions.
2) Energy equation in the crevices replacing the isothermal hypothesis.
3) Post-oxidation computation by a tabulated approach instead of one-step chemistry.
The advantages of the first improvement can be inferred from the observation of Figure 10, where results for the Engine A conditions at 3 and 8 bar of IMEP and 1500 min−1 are compared for two different coolant temperatures. This comparison highlights that uHC contributions from crevices and flame wall quenching have similar magnitude at cold and warm engine operation as to be expected based on thermodynamics (same IMEP means similar pressure and gas temperatures). However, as also well-known from literature,14,66 the oil layer contribution becomes crucial at cold engine conditions because of the increased oil film thickness and increased fuel solubility in oil. Without the introduced enhancing mechanism of fuel entrainment within the oil layer, a poor uHC prediction would have been found under cold conditions, resulting in uHC levels similar to those estimated under warm operations as visible also in the results from Esposito et al. 12
Figure 12 shows the impact of the model improvements 2) and 3) for three load variations and four λ variations for the Engine A. Concerning the introduction of the energy equation to estimate the temperature into the crevices, the advantages can be observed in Figure 12 comparing the actual final results of this work (red lines) with results obtained by fixing the temperature of the gas into the crevices equal to that of the piston (dashed teal lines). The comparison is performed by superimposing the crevice boundary layer temperatures since these temperatures have been verified against CFD results as shown in Figure 7. The average deviation on the selected points varies from 7.2% to 29.7% demonstrating the improvement associated to the crevice temperature calculation.

Comparison of uHC results from this work (red lines) against results in isothermal hypothesis for crevice gas (dashed teal lines) and results obtained with one-step chemistry for uHC post-oxidation (dot-dashed light blue lines).
Finally, the model enhancement deriving from the new post-oxidation approach can be observed again in Figure 12 by comparing the actual final results of this work (red lines) with results obtained through the application of the well-assessed one-step chemistry approach26,27,30,31,43–45 (dot-dashed light blue lines). Once again, these results have been obtained maintaining the same crevice boundary layer temperature. In this case, the average deviation of the predicted uHC for the considered points increases from 7.2% up to 83.9% demonstrating the potential need for an arbitrary recalibration of the one-step approach and the advantage of the here proposed post-oxidation method based solely on tabulated chemistry.
Concerning the same type of comparisons for engines B, C and D, regarding model improvement 2), the average deviations of the predicted uHC increase from 12.7% up to 15.0%, from 10.7% up to 11.7% and from 9.5% up to 29.7%. Regarding model improvement 3), the average deviations for engines B, C and D, increase from 12.7% up to 63.2%, from 10.7% up to 26.0% and from 9.5% up to 20.5%. Also, for these engines, the results have been obtained maintaining the same crevice boundary layer temperature to evaluate model improvements 2) and 3). The obtained results confirm increased predictivity of the here proposed model.
Summary and conclusions
This work has presented further developments to a model for the prediction of unburned hydrocarbons emitted by internal combustion engines. Starting from the authors’ previous work,11–13 model approaches have been refined in order to increase the level of physical description of the adopted phenomenological models. New modelling aspects like the direct fuel inclusion in the oil film have been considered. The model is developed with the help of 3D-CFD simulations and validated against four different engines including two small-bore single-cylinder engines fuelled with gasoline, and two large-bore naval engines fuelled with natural gas. The model capabilities to predict uHC are very high at a very wide range of operating conditions including cold start, very different engine designs, fuels and ignition sources. In total, 83 different engine operating conditions from the four engines have been considered for the validation and the average deviation is 11.4%.
The model proves to correctly describe the major effects contributing to uHC and incomplete fuel combustion in a physical manner. Therefore, the model is suitable for the future virtual development of new combustion concepts. With the right update of model constants, the model could be extended to predict unburned fuel emissions from engines fuelled with alternative potentially non-fossil fuels like methanol, hydrogen, or ammonia. Future work will focus on the extension and validation of the models in combination with those fuels.
Footnotes
Appendix
Acknowledgements
The authors are thankful to Convergent Science Inc. for providing licences for the CONVERGE software.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was partially funded by ‘Deutsche Forschungsgemeinschaft’ (DFG – German Research Foundation) – GRK 1856. The 3D-CFD simulations were performed with computing resources granted by RWTH Aachen University under project ‘rwth0239’.
