Development of a 0D multi-zone model for fast and accurate prediction of homogeneous charge compression ignition (HCCI) engine

Homogeneous Charge Compression Ignition (HCCI) is a promising advanced combustion mode, featured by both high thermal efficiency and low emissions. In this context, a 0D multi-zone model has been developed, where the thermal stratification in the combustion chamber has been taken into account. The model is based on a control mass Lagrangian multi-zone approach. In addition, a procedure based on a tabulated approach (Tabulated Kinetic of Ignition TKI) has been developed, to perform an accurate and fast prediction of the air/fuel mixture auto-ignition. This methodology allows combining the accuracy of detailed chemistry with a negligible computational effort. The tabulated procedure has been preliminarily verified through the comparison with the results of a commercial software (GT-PowerTM). In this assessment, single zone simulations have been performed comparing the TKI strategy to a conventional chemical kinetics one, in four different cases at varying the intake temperature and the equivalence ratio. Then, the proposed 0D multizone model has been validated against experimental data available in the literature. The analyses are carried out with reference to an HCCI engine fuelled with pure hydrogen and working in a single operating point, namely 1500 rpm, 2.2 bar IMEP and with a fuel/air equivalence ratio of 0.24. Three different temperatures, i.e., 373, 383, and 393 K, have been considered for the intake air. The experimental/numerical comparisons of pressure cycles and burn rates proved that the proposed numerical approach can reproduce the experiments with good accuracy, without the need for case-by-case tuning.


Introduction
In the present-day scenario, the automotive industry is facing a challenging path consisting in the design of low emission and low fuel-consumption vehicles, through the development of innovative Internal Combustion Engine (ICE) architectures. To this aim, an increasing interest is addressed to innovative combustion modes. Among these options, the Homogeneous Charge Compression Ignition (HCCI), as a low-temperature combustion (LTC) concept proposed by Onishi et al. [1], represents a reliable technique, able to produce ultra-low NOx levels and to provide greater fuel conversion efficiencies, compared to those of conventional ICEs [2][3][4].
This combustion mode merges the benefits of the two conventional engine combustion processes. Indeed, the presence of homogeneous air/fuel charge preparation, like SI combustion, removes the fuel-rich region, resulting in low smoke emission. The mixture burns in a controlled auto-ignition, like DI combustion, allowing the adoption of a lean air/fuel mixture and higher compression ratios with relevant improvements for fuel consumption. However, HCCI combustion presents some drawbacks which should be overcome [4][5][6][7][8]. They include the lack of any direct control method for combustion phasing; the high production of HC and CO emissions; limitations in achieving high engine loads (knock issue); misfiring; and cold start problems. To overcome the above-listed issues, the HCCI concept is still under investigation and development with the aim to become a feasible and reliable alternative to conventional combustion engines. In this context, both 1D and 3D-CFD numerical approaches represent valuable tools to support the study and the improvement of HCCI. The engine modelling and simulation studies are applied to enhance the engine design, improve engine calibration, developing the optimal engine control strategies. The combustion process of conventional SI engines in 1D codes is successfully modelled employing a phenomenological model based on a two-zone schematization of the combustion chamber, namely the burned and unburned zone separated by a turbulent propagating flame [9,10]. Instead, since the HCCI engines are characterized by the absence of the flame front propagation, the combustion phenomenon is described through the evolution of the mixture composition in the cylinder via a detailed chemical kinetic mechanism [11][12][13]. For this reason, the most common HCCI modelling is based on a single zone approach [14][15][16]. However, this schematization is not able to accurately predict the combustion duration, the heat release rate, and the pressure rise rate since a unique auto-ignition event of all mixture instantaneously occurs. In this case, an over-prediction of the pressure and temperature peaks is obtained [17][18][19][20][21][22][23]. To overcome these limits, phenomenological multi-zone models have been developed, where the combustion chamber is subdivided into smaller volumes. The number of zones and their configurations are determined through a balance between the accuracy and the computational effort. Either two zone, the core and the boundary layer [24]) or many zones (up to thirty [25]) can be specified. The multi-zone schematization improves the prediction capabilities since it allows to compute a temperature gradient between the zones, ensuring a progressive ignition of the charge [26][27][28][29][30][31][32][33][34][35][36][37][38][39]49]. It is known that the thermal stratification before ignition is primarily developed by heat loss to the cylinder walls, with little contribution from mixing with internal residuals [35]. However, the multi-zone approach has got some issues: one of them is the computational time, due to the resolution of chemical kinetics in each zone. In order to reduce the computational effort in the HCCI combustion modelling, two approaches can be followed. The first one consists of simplifying the kinetic mechanism and finally use a scheme with few species. Because of the complex nature of hydrocarbon fuels, those simplifications typically produce misleading results. The second uncommon approach is the so-called Tabulated Kinetic of Ignition (TKI). It is based on the tabulation and interpolation methods, which record the essential parameters of an adiabatic reactor simulation in a database to be employed at runtime for the engine simulation [36]. In this framework, this study takes place focusing on the development of a multi-zone model, based on the TKI approach, to reproduce the combustion in an HCCI engine.
The paper is arranged as follows. First, the balance equations for the multi-zone model have been described, to derive the temporal evolution of the thermodynamic properties during the engine cycle. Then, the adopted tabulated procedure is verified in a single-zone schematization, through the comparison with results obtained by a commercial software employing a chemical kinetic simulator, under various operating conditions. Once validated, the model is utilized in its multi-zone version, and a sensitivity analysis is performed on the tuning parameters of the temperature stratification in the combustion chamber. Finally, the whole procedure is applied to reproduce the experimental data related to an HCCI engine fuelled with hydrogen.

Methodology
In the study of HCCI engines, 0D multi-zone modelling presents many advantages over the 0D single-zone approach and the 3D CFD one. It must be noted that the 3D CFD model for HCCI combustion treats any cell as a homogeneous reactor and therefore it does not account for turbulence-chemistry interactions (TCI) [37][38][39]. Since the HCCI combustion is primarily controlled by chemistry, the 0D multi-zone model is an excellent option when a broad variety of parameters have to be modified, thanks to the low computational cost. There are two macro-categories of quasi-dimensional multi-zone models: Eulerian controlvolume and Lagrangian control-mass approaches. Both assume that the pressure in the combustion chamber is uniform, and both require a correlation for the overall heat transfer with the combustion chamber wall, such as the Hohenberg correlation. In the Eulerian control-volume model [29][30][31][32], the knowledge of the relative zonal size and position is assumed. Since the pressure is uniform, the mass that crosses the zonal boundaries is calculated using the ideal gas law. Instead, the Lagrangian model [33][34][35][36][37][38][39] is developed on a control mass basis (each zone has a specified constant mass) and no information is given on zone's position inside the combustion chamber. Here, the ideal gas law allows the computation of the zonal volume. The model proposed in this paper is based on the latter framework. It allows allocating a mass fraction quantity in each zones, based, for example, on a normal distribution, as reported in Fig. 1. This distribution assumes an almost constant level for very large values of the standard deviation, . Further, in the Eulerian model the heat transfer to the walls only occurs in the thermal boundary layers (outer zone) and the heat is transferred between the inner zones through interzonal heat transfers. Since there are multiple mechanisms involved in this model, such as the diffusion and convective exchanges, several empirical correlations must be applied. This tuning condition severely limits the implementation and predictability of Eulerian models, although this model appears to be closer to reality. The Lagrangian approach has been conceived to avoid the specification of interzonal mass and heat transfer. However, it still requires tuning for the assigning a predefined zone-byzone heat transfer fraction. Indeed, the thermal stratification is achieved by setting different fractions of the total cylinder's heat transfer to each zone, following the rule that the charge closer to the boundaries (typically outer regions and boundary layer) loses more heat than the charge further away (typically inner regions). A fully adiabatic inner zone could be also specified. The above heat transfer fraction distribution can be derived by experiments or 3D-CFD outcomes and was verified to be generalized to a broad variety of operating conditions [35,36]. Fig. 2a depicts such a distribution of heat transfer fraction per unit mass versus the cumulative mass fraction (CMF). With this curve, each infinitesimal mass layer, starting from the core of the combustion chamber (CMF = 0) to its walls (CMF = 1), exchanges a certain amount of the overall heat to create the desired thermal stratification in the combustion chamber. For consistency, the integral of this curve must be equal to 1. The continuous curve versus the CMF in Fig. 2a cannot be directly used during the simulation, since the model needs a discrete distribution based on the number of zones and zonal mass distribution. Once the number of zones and the zonal mass fraction distribution are known, the heat transfer fraction for a specific zone, called , , can be easily derived by integrating the curve in Fig. 2a in the appropriate interval. Fig. 2b depicts the integration results for each zone. Although no spatial information is required, the first zone is presumably located in the core of the combustion chamber, while the last one is adjacent to the cylinder wall, as the schematization of the boundary layer. Table 1 summarizes a comparison between the described approaches.

Lagrangian Multi-zone model: balance equations
The proposed HCCI model is implemented in the commercial software GT-Power through user coding, and it is utilized during the closed valve phase for engine simulations. The multi-zone model switch to single-zone model during the open-valve period, coupled to the 1D resolution of the intake/exhaust pipes. The calculation of the initial conditions at Inlet Valve Closure (IVC) is provided by the 1D flow model at each engine cycle. Starting from the initial conditions at IVC, the average properties and the composition of each zone (pressure, temperature, and gas composition) can be initialized. Furthermore, an initial temperature stratification among the zones can be also defined, maintaining the mean energy content of the cylinder. For the considered model, several assumptions are used: 1) each specie in the combustion chamber is treated as an ideal gas, taking into account the variability of gas properties; 2) since no interzonal mass transfer is allowed in the Lagrangian approach, the mass fraction of each zone remains constant during the simulation; 3) blow-by and crevices effects are neglected. The zonal energy evolves during time due to liquid species evaporation, combustion, heat transfer to the wall, and interzonal work exchange (i.e., work done by one zone on another). The zonal composition, on the other hand, only varies due to the auto-ignition process. Eq. (1) depicts the first law of thermodynamics for each zone.
The term is the overall heat transfer and is calculated using the Hohenberg correlation [40], used with the average cylinder temperature. Then, it is subdivided into zones through the coefficient , , computed as specified in Fig. 2b. Moreover, if the engine is fuelled with gaseous fuel, the evaporation heat is zero. The last term in Eq. (2) takes into account the energy released by the variation in the composition of the species within the zone related to the auto-ignition event.
For any specified condition, the following equation must be solved for each chemical specie: being the molar mass of species i, being the density of the zone z, and ̇, being the molar consumption/production rate of specie i per unit volume [ 3 ] during combustion. Typically, this term is provided by the solution of a chemical kinetic scheme at each time step and for each zone. As known, the complex chemical processes occurring inside an HCCI engine would require a large-scale mechanism, able to take into account both low and high temperature reactions. For a typical hydrocarbon fuel, such as gasoline, more than hundred species are usually required and, in this case, the computational time becomes too long for a practical application. Thus, the multi-zone chemical kinetics models offer no substantial advantages over the CFD ones. Furthermore, the detailed knowledge of species concentrations in a single zone is not mandatory; indeed, only the development of the most important species (fuel, O2, N2, CO2, H2O) is required for the computation of the zonal heat release. For these reasons, a methodology, known as Tabulated Kinetic of Ignition -TKI [41][42][43][44][45][46][47], is designed to perform an accurate and fast prediction of the air/fuel mixture auto-ignition.
The TKI technique will be briefly described and tested in the next section.

Description and validation of TKI
As mentioned, the TKI approach aims to describe the heat release due to the auto-ignition process based on the results obtained from a priori simpler chemistry calculations. The burn rate prediction in engine simulation can be achieved through this approach by exploiting a phenomenon called Burn Rate Equality [36]. According to this property, the burn rate of an auto-ignition combustion engine, at each time step, equals the burn rate of a homogeneous adiabatic constant-volume (CV) or constant-pressure (CP) reactor with the appropriate initial pressure, temperature, composition, and the same combustion progress variable [36].
One of the challenges in the TKI approach is the definition of the combustion progress variable, which defines the location of the engine chemical state between the initial pre-reaction state and the final equilibrium one. It is used to find a similar and unique state between the engine simulation and the adiabatic reactor one. Indeed, the progress variable definition has been widely explored in literature works [41][42][43][44][45][46][47]. Typically, it is based on the time evolution of species concentration (CO2, CO, O2, etc.), formation enthalpy or temperature.
In this work, the definition of the progress variable proposed by Lehtiniemi et al. [45] has been adopted. It defines the c variable as a non-dimensional heat released by combustion, calculated as: ℎ being the reactor enthalpy of formation at reference temperature (298 K), ℎ 298 . The ℎ (0) and ℎ ( ) represent the initial and the final values of reactor enthalpy of formation, respectively. Each point in the thermochemical state space is uniquely characterize by the definition of ℎ : Preliminary, a database for engine simulations is built from the solution of the chemistry in the constant-volume or constant-pressure reactors with specified initial conditions. At the end of each reactor calculation, progress variable reaction rates () are stored as a function of the discrete values of the normalized progress variable c. The initial conditions in the reactor simulation are selected to replicate all the possible incylinder conditions. They are defined based on two parameters describing the thermodynamic state (initial pressure, 0 , and fresh gas temperature, 0 ) and two variables defining the mixture composition (fuel/air equivalence ratio, , and EGR fraction, ). For each parameter, a variation range for the initial conditions is selected to properly match the local sensitivity of chemistry to that parameter. Table 2 shows the number of tabulated values, and the maximum and minimum levels for each parameter in the look-up tables. At each time step of the simulation, the solver, based on thermodynamic engine state and the table used (CP or CV), identifies the appropriate initial properties ( 0 , 0 , and ), and inquires a specific row of the table. Here, the value of tabulated c closest to the engine simulation one is identified and the corresponding ̇value is retrieved from the table, to integrate the progress variable in each zone. This datum is also used in the energy balance equation to calculate the amount of heat released by the on-going auto-ignition (AI) process, as a consequence of the species composition variation. During a time step, a fraction of fuel in the specified zone, calculated by , ∝̇, ∝ ̇, , is assumed to react with the air, leading to the formation of chemical equilibrium species in the in-cylinder mixture.
To assess the reliability of the c integration and the TKI method, Fig. 3 reports the pressure trend of a single-zone model, related to four different cases at varying operating conditions (i.e., intake temperature and fuel/air equivalence ratio). These single zone simulations are performed to compare the online chemical kinetics, realized using a conventional chemical kinetics software, and the tabulated kinetics (TKI approach), either performed in a homogeneous adiabatic reactor at constant pressure or constant volume. The simulations are carried out considering an engine fuelled with hydrogen, modelled with Hong's kinetic scheme [48], with 10 species and 20 reactions. Fig. 3: Comparison of the pressure results between the online chemical kinetics approach (the red line, "CHEMKIN"), the TKI approach (the blue and green lines, based, respectively, on the CPreactor and CV-reactor), on different intake temperatures and relative fuel/air equivalence ratios.
In Fig. 3 a perfect agreement between the results of online chemical kinetics and TKI can be observed. These results highlight that the TKI-based HCCI model enables fast computations and accuracy comparable to one of the conventional chemical kinetics models. Furthermore, the adopted tabulated approach, although characterized by a scheme of just 10 species and 20 reactions, is three times faster than the conventional chemical kinetics. This disparity would be much greater if a more detailed kinetic scheme is utilized, as, for example, occurring for a hydrocarbon fuel. Although both CV and CP approaches can be successfully used with an appropriate interpolation, CP simulations are chosen for the remainder of the paper. Conversely to traditional chemical kinetics methods, TKI model is characterized by three different benefits. First, the numerical method is non-stiff, which makes the simulation much faster. Second, the simulation time for the pre-generated data is much faster than solving the chemical kinetics in the engine simulation. Third, a more accurate and detailed chemical kinetics mechanisms can be used to generate the relative tables, without impacting the computational effort, even if the table build-up will be more time demanding.
In this paragraph, a sensitivity analysis is performed to evaluate and quantify the impact of tuning parameters on the multi-zone HCCI model. These results will be employed during the model calibration to reproduce the experimental data, as deeply treated in the following paragraph.
The tuning parameters of the proposed multi-zone model are the number of zones; the zonal mass fraction distribution; the profile for the heat transfer fraction distribution (Fig. 2a); and initial temperature distribution between the various zones. The analyses are performed at a fixed engine speed of 1500 rpm, an intake pressure of 1 bar, an intake temperature of 383 K, and a fuel/air equivalence ratio of 0.3. Except otherwise specified, the set of tuning parameters is selected as listed in the second column of Table 3, while the third column describes the values used in the sensitivity analyses. Table 3. The set of tuning parameters in the sensitivity analyses.

Tuning Parameters Value Sensitivity value
Number of zones The outcomes of the sensitivity study on the number of zones and the mass fraction distribution within the zones are reported in Fig. 4. In this figure, the pressure has been plotted shifted of 10 bar on the y-axis for each simulation to ensure a better representation.  Here, the single-zone model (red line) burns with a great delay compared to the multi-zone model (even with 5 zones), since the single-zone model does not take into account the presence of the adiabatic core of the combustion chamber. Furthermore, in the multi-zone model, the pressure rises generated by the instantaneous auto-ignition of the first zone leads to the auto-ignition of the other zones. This mechanism allows to take into account the combustion duration. A smoother pressure rise can be obtained by increasing the number of zones, even though a higher computational effort is required. The pressure traces between the 25-and 50-zones simulations are essentially indistinguishable. The impact of the zonal mass fraction distribution are plotted in Fig. 4b. In these simulations, the standard deviation of the reference Gaussian distribution for the mass fraction between zones has been varied as listed in Table 3, following the trends reported in Fig. 1.
An increased mass fraction in the intermediate zones (small value) implies a corresponding higher ∆ . Similar behaviour can be noted in the evolution of the progress variable, although the overall results are very close to each other. These simulations are carried out with just 5 zones, since the differences are even lower by increasing the number of zones. A further sensitivity analysis is performed on the distribution of heat transfer between the different zones and on the initial temperature stratification at IVC. Concerning the heat transfer parameters, six different profiles have been tested, some of which derived from the literature papers (namely, the profile 1 and 3 [35,36] in Fig. 5). The functional form of the heat transfer distribution mainly depends on engine geometry than on operating conditions. For the engine analysed in this work, no experimental or numerical data are available to derive a thermal stratification profile. Basing on the above consideration, the profiles 1 and 3 were parameterized to create additional trends plotted in Fig. 5, which will be used in the calibration with the experimental data.   Fig. 6a reports the results of the application of the above defined profiles. With the yellow trend (#1), most of the heat transfer is located at the cylinder periphery and a quasiadiabatic inner zone arises with an early auto-ignition. On the contrary, the magenta trend (#6) corresponds, as expected, to a single-zone model, since all zones evolve with the same temperature. In Fig. 6b, a linear temperature distribution is imposed at IVC between the first and the last zone. The results of this analysis are similar to the ones of the heat transfer distribution. Once the initial temperature stratification is increased, the combustion start is advanced and its duration increases, as a consequence of the wider thermal stratification existing in the combustion chamber. After understanding the effects of these parameters and selecting the actual degrees of freedom of the model, a comparison of the outcomes of the multi-zone model and the experimental data can be attempted, as reported in the following section.

Model Validation
Here, the validation of the proposed model through the comparison with literature-derived experimental outcomes is presented. In particular, the experimental results reported by Ibrahim et al. [50] are considered. The adopted engine is a single-cylinder, four-stroke, water-cooled, naturally aspirated unit fuelled with hydrogen through a direct injector (Engine manufacturer: Kirloskar, Model: AV1). Main engine characteristics are reported in Table 4. Model validation is proposed in terms of the numerical/experimental assessment for the incylinder pressure traces and the burn rate profiles. The latter are derived by the authors' reverse procedure processing the measured pressure cycle. The analyses are carried out considering a single operating point, namely 1500 rpm, 2.2 bar IMEP, and with a fuel/air equivalence ratio of 0.24. Three different temperature levels at intake valve closing (IVC) were considered: 373 K, 383 K, and 393 K. The validation results are reported in Fig. 7, showing the numerical/experimental comparisons of in-cylinder pressure traces and burn rate profiles as a function of the engine crank angle by varying the intake temperature.
To reproduce the experimental data, the model was set with 40 zones, a constant mass fraction distribution, the yellow curve in Fig. 5 for the heat transfer distribution, and an initial temperature stratification at IVC of 2 K. Fig. 7 shows a quite satisfactory agreement between the model outcomes and experimental data. The proposed multi-zone approach, considering temperature stratification in the combustion chamber, accurately predicts the engine characteristics in terms of combustion phasing and burn rate. Moreover, it provides reasonable estimations of the maximum pressure rise rate; this parameter allows the calculation of the Ringing Intensity, which is used to properly forecast the load limit of HCCI engines. A greater error occurs for the lower intake temperature of 373 K (green line in Fig. 7). This case is close to the auto-ignition threshold and, in this condition, the model is extremely sensitive to the mixture temperature. Indeed, a slight temperature difference causes a relevant variation in the auto-ignition timing. The Fig. 8 shows the evolution of the overall combustion progress variable for the three considered cases. In each case, some unburned mixture is still present at EVO. None of the three trends reaches the unity, and, particularly, for the case with the intake temperature at 373 K around 20% of gases remains unburned. Thus, this model is able to predict the unburned hydrocarbons from unreacted unburned zones, although a more refined HC simulation also requires the characterization of the crevices flows and the post-oxidation. The observed incomplete combustion can be better clarified through the analysis of the temperature and zonal progress variable trends, reported in Fig. 9a-c.
Depending on the actual zone temperature, different auto-ignition events occur zone-byzone during time. However, in each operating condition, the very low temperature of the outer zones determines the absence or partial auto-ignition reactions, leaving some unburned hydrocarbons inside the combustion chamber.  Fig. 9: The zonal and bulk temperature (top) and the zonal and bulk combustion variable progress (bottom) for all considered cases in the experimental comparison From these figures, it can be noted that the pressure gradient, generated by the zone autoignition, is the only connection between the various zones. In particular, referring to the highest intake temperature of 393 K, Fig. 9a shows that the combustion occurs almost completely before the firing TDC, i.e., during the compression phase, and all zones burn close to each other. The close auto-ignition of the zones is due to the increase in pressure, and consequently in temperature, generated by the auto-ignitions of the previous zones. Instead, at lower intake temperature (third case of 373 K), Fig. 9c reports a much longer overall oxidation since it takes place entirely in the expansion phase. Here, the pressure increase produced by the early burning zones is mitigated by the expansion; thus, the colder near-wall zones are not able to reach the auto-ignition conditions.

Conclusion
This work proposes and validates a phenomenological 0D multi-zone model, implemented as an "user routine" in the commercial software GT-POWER, for a fast and accurate prediction of HCCI combustion. The model is based on the control-mass Lagrangian approach and on the TKI method. The TKI approach gives some significant benefits, such as the possibility to operate with a good accuracy and a larger number of zones, preserving the computational effort. Primarily, the TKI methodology has been validated by comparing the results of single-zone simulation to the ones of the conventional chemical kinetic model. Then, a parametric analysis on the tuning parameters of the multi-zone model has been performed. Here, the results show that it is required a larger number of zones for a smoother pressure trajectory, while the zonal mass distribution does not play a significant role. Furthermore, it is fundamental to describe the zonal heat transfer distribution and initial ∆ , since the thermal stratification is crucial to forecast the burn rate profile of HCCI combustion. Finally, the reliability of the HCCI multi-zone model has been proved through the comparisons with the experimental data in the case of an engine fuelled with hydrogen. It is worth underling that the comparison has been made by keeping constant the model tuning parameters at varying operating conditions. The model demonstrates an adequate capability to reproduce the behaviour of an HCCI engine, also perceiving the effect of the intake temperature. It can be usefully employed to help the development of new engines featured by HCCI combustion. As the next developments of this model, the research should be conducted in these four areas: 1) carrying out more experiments in various operating conditions to extend the model validation. This will allow showing the ability of the model in predicting operations in a wide range of fuel/air and EGR ratios; 2) developing a more accurate model for the heat transfer distribution. Currently, the tuning parameters are not correlated with any physical phenomenon in the combustion chamber. Thanks to 3D CFD calculations or experimental data, these values could be correlated to some variables, including turbulence values, walls' temperatures, the amount of EGR, etc; 3) calculating the NO emission and considering the blow-by and crevices effects for a better prediction of the HC and CO emissions; 4) extending the use of this model to other combustion strategies (SACI, RCCI), including the integration of this model with a conventional flame front propagation model.