2.4.1. Geometry and Meshing
The simulations were performed with the commercial software CONVERGE CFD (v2.4). The computational mesh domain was divided into three regions: the intake region, which included the intake runners and ports, the combustion chamber, and the exhaust region, which included exhaust ports and runners. The combustion chamber geometry included also the geometry of the piston top-land crevice up to the first compression-ring. The cold-geometry diameter (
D) gap on the piston crown was 1 mm, with
of 75 mm and
of 74 mm. The gap in warm operation can be reduced due to the thermal expansion. However, the amount of thermal expansion is difficult to estimate and depends on the operating points. For this reason, in
Section 3.5.2, a sensitivity analysis of the computed emissions on the piston-liner gap is reported. The CFD geometry included also the tips of all the sampling probes applied in the two emission measurement positions, in order to take into account possible flow disturbances.
Figure 1 shows the CFD geometry including the emission sampling locations for the different devices.
The software CONVERGE has an automatic meshing at simulation time. The base size of the structured Cartesian mesh uses 2 mm cells in the intake and exhaust region and 1 mm in the cylinder region. Mesh refinements are implemented in the valve gaps during gas exchange and around the spark plug electrode gap at spark timing. Additionally, an automatic mesh refinement (AMR) reduces the mesh size up to 0.5 mm in the intake and exhaust regions if temperature and/or velocity gradients overcome a selected maximum value.
More details on the meshing refinement settings, boundary and initial conditions, and standard models adopted in the simulation are provided in
Appendix A.2.
2.4.3. Injection and Wall-Film Model
For the spray break-up, the Kelvin–Helmholtz/Rayleigh–Taylor (KH-RT) hybrid model [
32] was adopted. The model was calibrated in order to achieve the same spray penetration measured with static injector optical measurements. The calibration was performed by changing the time constant parameter [
33] of the KH-RT model [
32]. Typical values of this constant are in the range 5–100 [
33]. In this work, it was calibrated to a value of 30. For the film splash, the O’Rourke model [
34] was selected.
A strong accumulation of the fuel wall film in the piston top-land area was observed, especially in the case of the simulation of multiple consecutive cycles. Results of the film accumulation in the piston top-land area are shown in
Figure 2a, while a schematization of the problem is reported in
Figure 2b.
This effect had a strong impact on the HC emissions, increasing them by almost a factor of two. This accumulation effect was likely overestimated because the fuel-film in the piston-ring crevice would be mainly incorporated in the oil-film, transported in the ring pack, or reduced by the effect of the blow-by in the real engine operation. This was considered here by removing the residual wall film mass at −20 CA before top dead center (TDC).
2.4.4. Chemical Kinetics and Combustion Model
The kinetic solver available in CONVERGE CFD (SAGE) [
35] is active when a cell temperature exceeds 500 K and when the molar concentration of HC plus CO is higher than 0.1 ppm. The gasoline surrogate mechanism from [
36] was taken into account in this study. Since the focus is on gasoline combustion in SI engines and the prediction of its knocking behavior is out of scope, the mechanism was reduced for improved computational cost by removing the low temperature chemistry of long-chain components,
i-C
8H
18 and
n-C
7H
16. The reduced mechanism consisted of 239 species and 1068 reactions, including the detailed NO
formation mechanism from Lamoureux et al. [
37]. It was validated against the experimental laminar burning velocities available in the literature for gasoline fuels and their surrogates. An example of the mechanism validation is reported in
Appendix A (
Figure A1), and the kinetic mechanism itself is provided in the
Supplementary Materials.
In order to focus on the prediction of the emissions, an accurate reproduction of the heat release and pressure trace is necessary, and it was achieved with the calibration of a
G-equation flame propagation model. The model from Gülder [
38] was selected for the calculation of the laminar flame speed. In this work, a novel methodology for the interaction of
G-equation and kinetics with the scope of correct emission simulation was developed. As presented from Yang and Reitz [
20], it is important to detect if the combustion is turbulence- or chemistry-dominated, especially in the late combustion phases. Indeed, the usage of the
G-equation when the combustion is no longer turbulence-controlled does not allow for a correct evaluation of HC and CO emissions, since the flame propagation continues in the expansion stroke as long as
G burnt (
) and unburnt (
) cells exist. An alternative to the approach of Yang and Reitz [
20] was developed, in which the
G-equation model was deactivated in the last phase of the combustion, when the turbulence level was low and further oxidation was kinetically controlled. The deactivation was achieved by means of a re-initialization of the
G-scalar to negative values in the cylinder, and the later burn-out phase was then calculated solely by the SAGE solver. The timing at which the
G-equation was deactivated was selected on the basis of the considerations of the heat release profile. The point at which the heat release rate slowed down was when the flame front touched the liner walls. This was when the flame propagation stopped being essentially turbulence-dominated and the combustion rate was expected to be mainly chemistry-dominated. This point was analytically determined on the basis of the measured heat release profiles. A sensitivity analysis on the effect of the
G-equation deactivation point on the emission results will be presented in
Section 3.5.2. The crank angle taken as deactivation angle
is where the highest maximum of the second derivative (
) of the heat release rate (
) is reached in the second half of the combustion duration.
Figure 3 reports graphical examples of the determination of
for LS2 and LS3.
In
Figure 3, it can be observed how this analytical procedure identifies the point at which the heat release slows down for the two operating points.
Due to the selected methodology, the SAGE kinetic solver is used not only in the burnt zone, but also in the unburnt zone for the entire simulation time. Even if the low-temperature chemistry was not included, this was done for two reasons. First, after the G-equation deactivation and re-initialization of G to negative values, the whole combustion chamber is considered as the unburnt zone, and the SAGE solver is necessary to continue the calculation of the heat release. Additionally, the simulation of multiple consecutive cycles results in pollutants being present in the unburnt zone. The use of the SAGE solver allows the kinetic determination of species change in the unburnt zone during the time when the G-equation flame propagation model is active.
At the flame front (cells with
), the species concentration is changed to burned conditions represented by chemical equilibrium to allow feasible computational times. This chemical equilibrium chemistry is only active as long as the
G-equation model is active, which is for a relatively short time of about 30
CA over the entire cycle. The assumption of equilibrium chemistry during the main part of combustion, in which the highest temperature is reached, is a good approximation for almost all the species relevant for this study. The only exception is NO because it has longer chemical time scales, also at high temperatures. It has been verified that this approximation is acceptable for stoichiometric operation because the kinetically determined NO concentration come close to equilibrium levels during combustion, especially at higher loads, as also mentioned by Heywood [
39]. This aspect was investigated for the selected operating points with 0D chemistry calculation in the burnt zone with a two-zone combustion chamber model in GT-POWER (v2019). These results are reported in
Figure A2 in the
Appendix A. For OPs LS2, LS3, and HP, the overestimation of NO concentration with equilibrium chemistry in comparison to chemical kinetics was approximately 15% at the time of
G-equation deactivation. For OP LS1, a stronger overestimation of NO, over 50%, was observed due to lower combustion temperature. However, the magnitude of these deviations observed with 0D-chemistry cannot be directly carried to 3D-CFD. In conclusion, the adopted combustion methodology can show an overestimation tendency regarding NO that is expected to be quite limited in stoichiometric operation at mid-high loads, higher at stoichiometric operation at low loads, and very strong at lean operation.
To start the combustion with the
G-equation model, the variable
G was set to positive values at spark timing in a sphere of 0.6 mm in radius centered in the spark gap. This simplified ignition model resulted in a shorter early kernel development phase. Thus, the spark timing was slightly retarded to compensate for this aspect. The calibration of the combustion speed was achieved by modification of the
parameter of the
G-equation model (
according to [
22]). In
Table 3, the parameters resulting from the combustion calibration are reported.
It is possible to note how the calibrated
is on a similar level for all the operating points. The burnt fuel fraction for
G-equation
was around 0.9 for all the operating points, with slightly higher values for the higher loads. In
Table 3, the difference between the simulation and measurement ST is reported in degrees CA (
) and in ms (
). The values are consistent with the hypothesis that the ST retardation needs to compensate for the burn delay, which is artificially shortened from the ignition methodology. Indeed,
follows a physical trend, as it increases at lower loads and at higher speeds.
2.4.5. Emission Post-Processing
Regarding the emission post-processing, monitoring points were introduced at the tips of all the modeled measurement probes. The average emission measurements (FEVER, IMR-MS, and FFID average values) were compared with the time averaged results from simulations in the same points as measured. In order to compare the emission values from CFD to the measurements, it is necessary to properly post-process them. For emission species like O
, CO
, and CO that are measured as dry (without water vapor), it is necessary to correct the CFD emissions (wet) according to:
in which
is the molar fraction of a generic emission species
i and
is the molar fraction of water vapor.
Regarding HC emissions, an online evaluation of the total-HC concentration was done during the simulation. Indeed, the recording of all the species at each time step and a successive post-processing would have been unfeasible. In particular, C
-equivalent HC values were calculated as user-defined passive quantities, in order to allow the direct comparison with the test bench results. The theoretical FID factor
for the generic HC species
i can be defined as:
where
represents the carbon atoms of the species
i. The total-HC C
-equivalent concentration
including all the HC species
i can be calculated according to:
The real FID factor of each species is different from the theoretical one
, considered in the post-processing, especially for oxygenated species. A comparison between real and theoretical factors was experimentally verified by the author for selected species and is reported in
Appendix A in
Table A4. However, even if a few real factors were known, the majority of them could not be estimated. For this reason, it was preferred to use the theoretical factors for all species. This was considered a possible source of deviations between CFD and measurements in terms of total-HC emissions.