**Validation of a RANS 3D-CFD Gaseous Emission Model with Space-, Species-, and Cycle-Resolved Measurements from an SI DI Engine**

## **Stefania Esposito 1,\*,†, Max Mally 1,‡, Liming Cai 2, Heinz Pitsch <sup>2</sup> and Stefan Pischinger <sup>1</sup>**


Received: 19 July 2020; Accepted: 8 August 2020; Published: 19 August 2020

**Abstract:** Reynolds-averaged Navier–Stokes (RANS) three-dimensional (3D) computational fluid dynamics (CFD) simulations of gaseous emissions from combustion engines are very demanding due to the complex geometry, the emissions formation mechanisms, and the transient processes inside the cylinders. The validation of emission simulation is challenging because of modeling simplifications, fundamental differences from reality (e.g., fuel surrogates), and difficulty in the comparison with measured emission values, which depend on the measuring position. In this study, detailed gaseous emission data were acquired for a spark ignition (SI) direct-injection (DI) single-cylinder engine (SCE) fueled with a toluene reference fuel (TRF) surrogate to allow precise comparison with simulations. Multiple devices in different sampling locations were used for the measurement of average emission concentration, as well as hydrocarbon (HC) cycle- and species-resolved values. A RANS 3D-CFD methodology to predict gaseous pollutants was developed and validated with this experimental database. For precise validation, the emission comparison was performed in the exact same locations as the pollutants were measured. Additionally, the same surrogate fuel used in the measurements was defined in the simulation. To focus on the emission prediction, the pressure and heat release traces were reproduced by calibrating a *G*-equation flame propagation model. The differences of simulation results with measurements were within 4% for CO2, while for O2 and NO, the deviations were within 26%. CO emissions were generally overestimated probably because of inaccuracies in mixture formation. For HC emissions, deviations up to 50% were observed possibly due to inexact estimation of the influence of the piston-ring crevice geometry. The reasonable prediction accuracy in the RANS context makes the method a useful framework for the analysis of emissions from SI engines, as well as for mechanism validation under engine relevant conditions.

**Keywords:** internal combustion engine; combustion; emission; RANS simulation

#### **1. Introduction**

The modeling of pollutant emissions from internal combustion engines is of great importance for the development of future low-emission powertrains. The emission analysis within 3D-CFD simulation allows the optimization of the combustion system [1] or the operating strategy [2]. Additionally, 3D-CFD simulation can provide insights into the emission formation mechanisms as a function of in-cylinder phenomena such as fuel distribution, residual gas fraction, etc. [3–11]. The simulation of gaseous pollutants for spark ignition (SI) direct-injection (DI) engines in 3D-CFD require the correct

prediction of many physical and chemical aspects, above all mixture formation, flame propagation, and chemical kinetics. Various studies have shown the usage of 3D-CFD to analyze gaseous emission formation with 3D-CFD simulation in SI DI engines [6,12–17].

The validation of emission simulation models with experiments is very challenging due to some fundamental differences. One aspect regards the surrogate fuel used in CFD that, due to its simpler composition, fails in the prediction of all the chemical and physical characteristics of market real fuels (e.g., distillation curve) [18]. For this reason, a few simulation studies used surrogate fuels with a large number of components to reproduce the evaporation curve of real fuels [8]. Additionally, validation requires a comparison in the location where the emission measurements are performed. It was shown in previous works [18,19] that the emission measuring position has a strong impact on the average emission level and cycle-resolved trends. However, validation studies reported in the literature rarely provide details on how the computed emissions were evaluated when compared to experiments. Often, the simulation domain does not cover the emission measuring position, and little information is given on what kind of computed values (e.g., average, instantaneous at a specific time) are evaluated where (e.g., in-cylinder, in the whole exhaust region, in a specific location) for the comparison with measurement data.

The implementation of chemical kinetics is fundamental for emission predictions. RANS cannot fully describe the turbulence-flame interaction, and the calculation of flame propagation only by solving species transport equations with detailed chemistry is often unsuccessful, as shown by Yang et al. [20]. Thus, a flame propagation model that takes into account the turbulence enhancement of the burning velocity is necessary. Among the different approaches [21], the level-set method (*G*-equation) formulated by Peters [22] for turbulent flames proved to be reliable for combustion simulation of SI engines [7,8,17,20,23,24]. To allow the interaction with chemical kinetics, Yang and Reitz [20] presented a criterion based on the ratio between turbulent and chemical time-scales to select whether fuel oxidation should follow the *G*-equation or the kinetics.

This study aims to provide a rigorous validation of a RANS CFD methodology for gaseous emission predictions. The simulations reproduce a specific experimental setup at a single-cylinder engine test bench. A toluene reference fuel (TRF) surrogate with ethanol (TRF+E) was considered in both experiments and simulation to minimize the impact of fuel on the validation results. The CFD geometry covers the emission sampling positions, allowing the comparison of computed emission values with experimental data at the same locations. The extensive comparison involves not only average gaseous emissions, but also species-, space-, and cycle-resolved hydrocarbon (HC) measurements. In this study, the chemical kinetics is used in combination with the *G*-equation combustion model, and a novel approach for their interaction is presented. Additional information on the described methodology can be found in [25].

#### **2. Methodology**

#### *2.1. Engine*

In this study, the measurement data came from a specific experimental setup, which was already presented by the authors in [18,19]. The main specifications of the single-cylinder engine are reported in Table 1.

Both the DI injector and the spark plug are in a central position in the combustion chamber. The spark plug is located between the exhaust valves, while the fuel injector is between the intake valves. The intake ports have a symmetrical high-tumble design.


**Table 1.** Single-cylinder engine technical data and operating parameter.

CA = crank angle; aTDCF after firing top dead center.

#### *2.2. Test-Bench Instrumentation*

#### 2.2.1. Thermodynamic Measurements

The cylinder pressure was measured with two pressure transducers (Kistler 6045B), flush-mounted in the combustion chamber roof, between the intake and the exhaust valve seat rings. The cylinder pressure signal was sampled via Kistler 5064 charge amplifiers and acquired with an FEV indication system (FEVIS) with a 0.1◦ crank angle (CA) resolution. The cycle-resolved intake and exhaust pressures were measured with Kistler 4045 A5 pressure transducers, sampled via Kistler 4665 and Kistler 4603 charge amplifiers with 0.1◦ CA resolution as well. For each operating point, one-thousand consecutive cycles were acquired. The static pressures and temperatures were measured with standard pressure transducers and thermocouples. The values were averaged over an interval of 30 s. The oil and coolant conditioning systems allowed steady-state operation. The intake air was conditioned to 25 ◦C upstream of the throttle flap. The pressure upstream of the throttle flap and in the exhaust manifold was controlled to 1.013 bar during throttled operation. During boosted operation, to simulate turbocharging, the pressure in the exhaust manifold was imposed equal to the pressure upstream of the throttle flap. The engine coupling with an eddy-current brake and an electric dynamometer allowed maintaining the desired engine speed with an accuracy of ±1 1/min, independently of the operating point. An ultrasonic air mass meter and a Coriolis-type mass flow sensor were used to measure the intake air mass flow and the fuel mass flow, respectively.

#### 2.2.2. Emission Measurements

In these experiments, detailed emission measurements were conducted in two sampling positions in the engine exhaust manifold and with multiple devices. One position was in the exhaust port (Pos. 1) and the second position (Pos. 2) further downstream. As measurement devices, a fast flame ionization detector (FFID) was used for cycle-resolved measurements of total-HC (THC) emissions (THC are referenced to C3H8), while an ion molecule reaction mass spectrometer (IMR-MS) was applied to measure average concentration of selected HC species. The FFID device used was a Cambustion HFR500, which was equipped with two probes of 210 mm, with an internal diameter of 1.07 mm (0.042 inches). The two sampling lines were heated to 200 ◦C and had a length of 45 cm to the FID measuring head and had an estimated response time (*t*10−90%) of 1.8–2 ms. The FFID signal output (0–10 V) was connected to the FEV indication system (FEVIS) system and acquired with the resolution of 0.1◦ CA. The raw FFID signals need to be reconstructed for delay in the probe and the filtering of the dynamics inside the probe [26]. More details about the FFID signal reconstruction can be found in [19]. The FFID accuracy was expected to be lower and the drift higher than conventional FID devices due to the selected larger scales (due to high instantaneous values especially in Pos. 1) and pressure sensitivity.

The IMR-MS used in these investigations was a V&F TwinSense, which is a soft-ionization mass spectrometer. This ionizes the sample with a lower ionization energy (IE) in comparison to traditional electron-ionization mass spectrometry [27–30]. The lower IE results in lower fragmentation of HC molecules and makes this kind of mass spectrometry more suitable for application to engine exhaust gas. The IMR-MS has two channels available and can measure in two different positions at once, as the FFID. The sampling lines were insulated metallic capillaries of 0.5 mm of diameter, heated to 100–120 ◦C, approximatively 1 m long. The device measured single species concentration, and the results were averaged over the total measuring time (≥1000 cycles).

The standard gaseous emission measurements were performed by means of an FEV emission rate (FEVER) measurement system. This device contains analyzers to measure cycle-averaged concentrations of main emission species in the exhaust gases. The emission measurements were performed on a partial mass flow of exhaust gases, which was sampled downstream of the cylinder head flange. The sampling line was heated to 193 ◦C. The technical details of the FEVER system are reported in Appendix A in Table A1.

#### *2.3. Fuel*

Kinetic mechanisms and fuel surrogates are used in 3D-CFD simulation to emulate the combustion behavior of real fuels. The surrogate formulations are usually calculated in order to match relevant fuel chemical indicators (e.g., knocking tendency). However, emission formation mechanisms can be affected from the surrogate formulation, especially regarding HC emissions. Thus, the comparison of 3D-CFD simulated emissions with emission measurement data with real fuel (e.g., gasoline) can lead to deviations. The authors verified [18] that the average global emissions (CO, CO2, NO, total-HC) of a TRF+E surrogate as defined in Cai et al. [31] are almost identical to that of its corresponding RON95E5 gasoline. However, high discrepancies were observed in terms of specific HC species, especially for the species that were also fuel components (toluene and n-heptane). In this study, in order to allow the comparison of CFD simulation with the species-resolved HC emissions, measurements at the single-cylinder engine with a TRF+E surrogate were taken into account. The surrogate composition in mass fraction was the following [18]: 51.8% *i*-C8H18, 13.3% *n*-C7H16, 28.8% C7H8, and 6.1% C2H5OH. The same fuel was used in 3D-CFD.

#### *2.4. RANS 3D-CFD Simulation*

#### 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 *D*liner of 75 mm and *D*piston 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.

**Figure 1.** Top view of the CFD geometry with an overview of the emission sampling positions (Pos.) and length of the intake (Int.) and exhaust (Exh.) regions. FFID, fast flame ionization detector; IMR-MS, ion molecule reaction mass spectrometer; FEVER, FEV emission rate.

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.2. Operating Points

Four operating points (OPs) were selected for simulation. A load sweep (LS) at *n*<sup>E</sup> = 2500 1/min was simulated, as well as a high power (HP) operating point at *n*<sup>E</sup> = 4000 1/min. Table 2 reports the operating parameters of the selected operating points.


**Table 2.** Simulated operating points (OPs). LS, load sweep; HP, high power.

#### 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.

**Figure 2.** Fuel-film distribution in the combustion for the OP HP at 20◦ CA before TDCF. (**a**) CFD rusults of the fuel-film formation on different planes, (**b**) schematization of the film positioning in the piston top-land crevice [25]. .

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*-C8H18 and *n*-C7H16. The reduced mechanism consisted of 239 species and 1068 reactions, including the detailed NOx 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 (*G* > 0) and unburnt (*G* < 0) 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 *α*G0 is where the highest maximum of the second derivative (*x*˙ <sup>b</sup> ) of the heat release rate (*x*˙b) is reached in the second half of the combustion duration. Figure 3 reports graphical examples of the determination of *α*G0 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.

**Figure 3.** Example of *G*-equation deactivation angle calculation on the basis of the measured (TPA) burning functions (*x*b) for two operating points [25]. (**a**) LS2; (**b**) LS3.

At the flame front (cells with *G* = 0), 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 *b*<sup>1</sup> parameter of the *G*-equation model (*b*<sup>1</sup> = 2.0 according to [22]). In Table 3, the parameters resulting from the combustion calibration are reported.

It is possible to note how the calibrated *b*<sup>1</sup> is on a similar level for all the operating points. The burnt fuel fraction for *G*-equation *x*b(*αG*<sup>0</sup> ) 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 (*Δα*ST) and in ms (*tΔα*ST ). 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, *tΔα*ST follows a physical trend, as it increases at lower loads and at higher speeds.


**Table 3.** Calibrated *b*<sup>1</sup> values, *G*-equation deactivation angle *α*G0 , and fuel burnt mass fraction at *G*-equation deactivation *x*b(*αG*<sup>0</sup> ); spark timing (ST) in simulation versus the one in the experiments for all operating points, together with the difference in CA (*Δα*ST) degree as well as in ms (*tΔα*ST ).

#### 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 O2, CO2, and CO that are measured as dry (without water vapor), it is necessary to correct the CFD emissions (wet) according to:

$$Y\_{i, \text{dry}} = \frac{Y\_{i, \text{wet}}}{1 - Y\_{\text{H}\_2\text{O}}} \tag{1}$$

in which *Yi* is the molar fraction of a generic emission species *i* and *Y*H2O 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, C3-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 *f*thFID,*<sup>i</sup>* for the generic HC species *i* can be defined as:

$$f\_{\text{thFID},i} = \frac{c\_i}{3} \tag{2}$$

where *ci* represents the carbon atoms of the species *i*. The total-HC C3-equivalent concentration *Y*THC including all the HC species *i* can be calculated according to:

$$Y\_{\rm THC} = \sum\_{i} Y\_{i} \cdot f\_{\rm thFID,i} \tag{3}$$

The real FID factor of each species is different from the theoretical one *f*thFID,*i*, 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.

#### **3. Results and Discussion**

#### *3.1. Mixture Formation*

Since the mixture formation affects the emissions strongly, the simulation results of the mixture formation are analyzed first. In Figure 4, the relative air-to-fuel ratio *λ* distribution at −20◦ CA before TDC for all the simulated operating point is depicted. Rich tails can be observed for all the operating points. In the piston-ring crevice zone, a richer mixture is usually observed since the injected fuel droplets enter that region.

**Figure 4.** Lambda distribution for all the operating points at 20◦ CA before TDCF: mass fraction distribution over *λ*-bins and central section of the combustion chamber. Operating points: (**a**) LS1, (**b**) LS2, (**c**) LS3, (**d**) HP.

As the starting value for the fuel mass, the injected mass is calculated on the basis of the measured fuel mass flow. However, since the cylinder filling and the trapped air mass in simulation can slightly differ from the trapped air of the real engine, the injected fuel mass was adjusted accordingly to achieve stoichiometric operation. Additionally, the simulation of multiple consecutive cycles is affected by CFD model inaccuracies (wall-film, injection, etc.), and small deviations in the simulation average *λ* up to 1.2% are observable. However, the accuracy of the *λ* determination in the measurements is typically around 1%. Figure 5 shows the deviation in terms of average *λ* for the different operating points.

**Figure 5.** Relative deviation between simulation average *λ* and measured *λ* (according to Spindt [40]).

#### *3.2. Combustion and Heat Release*

Figure 6 shows the comparison of the measured and simulated pressure curves (*p*cyl) and the heat release rate (*Q*˙ cyl) for all simulated operating points. In all the cases, a good agreement in the pressure traces and heat release is achieved. The above-described *G*-equation deactivation analytical procedure results in a deactivation time at a point where the burn rate slows down. The heat release rate is accurately predicted even after the deactivation of the *G*-equation model, which is an indication of the validity of the proposed methodology.

**Figure 6.** Comparison of measured and simulated pressure curves (*p*cyl), as well as the heat release rate (*Q*˙ cyl) obtained from three-pressure analysis (TPA) [41] and simulation. Operating points: (**a**) LS1, (**b**) LS2, (**c**) LS3, (**d**) HP.

#### *3.3. Average Emissions*

Figure 7 illustrates the comparison of the simulated average emission concentration with the FEVER measurements in Pos. 2. Together with the absolute values, the relative deviation of the simulation in comparison to the measurements is plotted.

**Figure 7.** Comparison of average emissions in Pos. 2 between simulation and FEVER measurements (bars, left axis) with relative deviation (line plot, right axis): (**a**) CO2, (**b**) O2, (**c**) CO, (**d**) THC, and (**e**) NO.

Regarding CO2 results (Figure 7a), a very high accuracy in simulation is achieved, with a maximum deviation of 4%. The O2 (Figure 7b) content in the exhaust is in relatively good agreement as well, with a maximum deviation below 10%, with the exception of OP LS3, where the deviation reaches about 25%. The relative deviation trends of CO2, O2, CO are consistent with each other and suggest a mismatch of the calculated mixture formation with the real one. Since the injector was calibrated only with static measurements, it is possible that the spray-tumble interaction at different operating points is not correctly predicted and leads to inaccurate predictions of the *λ*-field. Above all, CO (Figure 7c) is the most sensitive species to mixture formation, especially around stoichiometric operation and hence shows overestimation up to 130% and incorrect trends over the operating points. However, the deviation in CO emissions can be partly explained with the trend of the average *λ* deviation of Figure 5.

The best overall agreement for CO2, O2, CO is observed for OP LS1, which shows a more homogeneous mixture in Figure 4a than the other OPs. In Section 3.5.1, a sensitivity analysis on global *λ* for OP LS3 will be presented in order to evaluate the impact on the CO emissions.

Regarding THC emissions (Figure 7d), for three out of four points, a relatively good agreement with deviations within 30% is observed. The LS2 value is overestimated by about 50%. In the CFD simulations, the main source of HC emissions is from the crevices, since the mesh is not fine enough to resolve the flame-wall quenching. Furthermore, fuel wall film is removed, and other HC mechanisms are not modeled. Thus, a possible reason for the overestimation of HC can result from a too large piston-ring crevice volume. In Section 3.5.2, a sensitivity analysis of the piston-ring crevice volume will be presented for OP LS2. A further reason for the possible overestimation could be a too rich mixture in the crevice, e.g., due to inaccuracies in the calculation of the mixture formation. Indeed, a richer mixture in the piston-ring crevice results in higher HC emissions, since the mass in the crevice expands back in the main combustion chamber during the expansion phase. Figure 4 shows that the OP LS2 (b) has the highest inhomogeneity among the simulated points and a very rich tail that corresponds partly to the mixture in the piston-ring crevice.

As far as NO is concerned, good predictions are achieved in all operating points with maximum deviations of 23%. The overestimation observed for LS1 can result from the combustion methodology and the expected overestimation tendency of the chemical equilibrium calculation at low loads. For the other operating points, specific deviations are difficult to be traced back to single sources. Since the *λ*-sensor used in the measurements has a certain inaccuracy (normally ≈ 1%) as well, a deviation in global *λ* could be a possible source. The sensitivity analysis on *λ* presented in Section 3.5.1 will give an evaluation of the variability of NO with *λ*.

Figure 8 shows the comparison of the non-standard emission measurement devices (FFID and IMR-MS) in Pos. 1 and Pos. 2 for two operating points (LS1 and LS2). The selected quantities are the average FFID THC values (a, e) and three selected HC species measured with the IMR-MS, C2H2 (b, f), C6H6 (c, g), and C7H16 (d, h).

**Figure 8.** Comparison of average HC-emissions in Pos. 1 and in Pos. 2 between simulation (hatched bars) and measurements (FFID and IMR-MS) for all the operating points. FFID: (**a**) Pos. 1, (**e**) Pos. 2. IMR-MS: C2H2 (**b**) Pos. 1 and (**f**) Pos. 2, C6H6 (**c**) Pos. 1 and (**g**) Pos. 2, C7H16 (**d**) Pos. 1 and (**h**) Pos. 2.

In comparison to the average FFID values, the simulated THC are underestimated in Pos. 1 and overestimated in Pos. 2 (similarly to what is observed with the FEVER THC comparison in Figure 7c). The THC concentration in Pos. 1 is found to be strongly dependent on the fuel-film in the piston-ring crevice, and elimination of this film resulted in a reduction of the THC concentration in Pos. 1.

Regarding the HC species measured with the IMR-MS, overall similar levels in simulations are observed. This aspect by itself is a very positive indication that from one side, the experimental methodology to measure HC species with the IMR-MS is reliable. Indeed, the correct determination of HC species concentration in the engine exhaust is challenging [18,19]. From the other side, similar levels of some specific HC species (fuel components like C7H16 and partly oxidized species) are a further confirmation that the CFD methodology (combustion model, kinetics) is reproducing realistic results. Higher deviations are observed in Pos. 1 due to the high sensitivity to the piston-ring crevice mechanism. In Pos. 2, the deviation is reduced, and in several cases, a very high accuracy is achieved.

#### *3.4. Crank Angle-Resolved HC Curves*

Another detailed validation step is provided by the comparison of simulation results with the crank angle-resolved FFID measurements. Indeed, this can verify if the dynamics of the THC emissions in the CFD simulation are similar to what was measured on the test-bench. In Figure 9, the simulated instantaneous THC concentration is compared with the measured FFID signals in Pos. 1 and Pos. 2 for the OP LS2. The raw FFID signals (all 500 measured cycles and the average) are shown together with the signal reconstructed as explained in [19]. The exhaust valve opening and closing angles (EVO and EVC) are indicated as well.

Regarding the comparison in Pos. 1, a quite good agreement in trend and absolute values is observed. The simulated value during the time when the exhaust valve are closed has a high influence on the time-averaged value. The value appears to be lower than the average FFID signal (even if still in the scatter-band), while with the exhaust valve open, similar or slightly higher values are observed. In Pos. 2, a similar trend in simulation can be qualitatively observed even if the absolute values are higher overall. The timing of simulated and measured THC peak value are close to each other. This peak results from the passage of the HC stored in the exhaust port from the previous cycle, shortly after EVO, as visible in the 3D picture of Figure 9 at 170◦ CA aTDC. However, the peak shape and intensity differ between simulation and measurement, possibly because of the peak smoothing effect of the FFID probe.

**Figure 9.** Top (diagram): comparison between FFID cycle-resolved measurements and simulated THC concentration in Pos. 1 and Pos. 2; FFID raw signals (all 500 measured cycles and the average) and the reconstructed curve, compared to the CFD prediction. Bottom: CFD pictures of the total-HC (THC) distribution in the combustion chamber and in the exhaust system (see the pictogram) at different crank angles. Operating point LS2.

#### *3.5. Sensitivity Analyses*

In this section, sensitivity analyses of specific parameters identified as relevant for the gaseous emission predictions are reported.

#### 3.5.1. Average Relative Air-to-Fuel Ratio *λ*

O2, CO, and NO emissions are supposed to be very sensitive to the average *λ* value. This is verified in OP LS3, where the highest deviation in CO is observed. Three simulations were calculated starting from 20◦ CA before TDCF from the same 3D solution, in which the average *λ* was adjusted by means of a scaling of the fuel components in order to achieve *λ* = 1.00 (base value), *λ* = 0.98, and *λ* = 1.02. In Figure 10, the effect of this variation on the emission concentration in Pos. 2 is evaluated during the exhaust stroke.

**Figure 10.** Effect of ±2% *λ* on simulated emissions concentrations in Pos. 2 compared with measurements (FEVER): (**a**) CO2, (**b**) CO, (**c**) O2, (**d**) NO, and (**e**) THC. Operating point LS3.

It can be observed that while CO2 and THC are weakly influenced by the average *λ*, the effect on the other emission species is strongly non-linear, and a scatter-band up to ±30% is shown. This is a confirmation that matching the average *λ* in a very narrow band is of great importance for emission simulations and a great challenge due to the interaction of the measurement inaccuracies with all the CFD models (spray, wall film, evaporation).

#### 3.5.2. Crevice Volume and *G*-Equation Deactivation

The sensitivity to the *G*-equation deactivation angle and to the crevice volume were investigated for OP LS2. To verify that the developed combustion methodology does not affect the emissions, a simulation in which the *G*-equation was deactivated 4.5◦ CA later (20 instead of 15.5◦ CA aTDC) was performed. Additionally, the effect of a 20% reduction of the crevice volume (by reducing the piston-liner gap from 1 mm to 0.8 mm) was investigated. The effect of these two changes was evaluated by calculating the relative change in the mass of a certain pollutant species released through the exhaust valves during the exhaust stroke. Figure 11 shows the results of these investigations.

The effect of later *G*-equation deactivation is negligible. This is a confirmation that the adopted combustion methodology is robust.

The reduction of 20% of the piston-liner gap results in a decrease of ≈15% for CO, O2 and HC. This result correlates with the less than linear dependency of HC emissions with piston-liner gap verified by Min and Cheng [42].

**Figure 11.** Relative change in species mass flows through the exhaust valves in the exhaust stroke in case of: later *G*-equation deactivation (*G*-equ. deactiv.) (hatched bars) and -20% piston-liner gap (grey bars). Operating point LS2.

#### **4. Conclusions**

In this study, a RANS 3D-CFD simulation model was used to reproduce a specific experimental setup and to enable precise validation with space-, species-, and cycle-resolved emission measurements. A newly developed approach for the interaction of the calibrated *G*-equation flame propagation model with kinetic mechanism was shown to be robust and to successfully reproduce measured in-cylinder pressure and heat release. The validation of emission predictions was based on four operating points with varying load and engine speed. The deviation of CO2 was below 4%, O2 below 26%, and NO below 23% for the investigated points. For CO, a high deviation at *n*<sup>E</sup> = 2500 1/min and *p*mi = 16 bar of 130% was calculated, together with a high deviation of 25% in O2. This effect suggests a mismatch in the simulated mixture field in comparison to the test-bench. Regarding THC, two operating points showed deviations below 3%, while for mid-high loads and lower speed deviation, up to 50% was observed. This behavior can derive from additional HC formation mechanisms not modeled in the simulation (e.g., fuel-film, quenching), which affect by a different magnitude the operating points. Additionally, an overestimation of the piston-ring crevice contribution to the THC formation is not excluded.

The comparison with single HC-species measurements showed that the simulated values were in the same order of magnitude, with better agreement in the measurement position further downstream the exhaust valve. Although deviations were present, the combined experimental-simulative methodology shows potential for further validation of kinetic mechanisms under engine conditions. Lastly, the transient trace of THC emissions was validated by the measurement results, which showed a similar shape and timings in two sampling points. Sensitivity analyses showed that the effect of the timing at which the *G*-equation model was deactivated was negligible, while the global *λ* led to changes in CO, O2, and NO of up to ±30% when varied by ±2%. The reduction of the crevice volume resulted in a less than linear reduction of THC, CO, and O2 emissions.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1996-1073/13/17/4287/s1. The kinetic mechanism in CHEMKIN format is made available as Supplementary Material.

**Author Contributions:** Conceptualization, S.E., M.M., S.P.; Methodology, S.E., M.M., L.C.; Software, S.E. and M.M.; Validation, S.E. and L.C.; Formal Analysis, S.E.; Investigation, S.E.; Writing—Original Draft Preparation, S.E. and M.M.; Writing—Review & Editing, S.E., M.M., L.C., H.P. and S.P.; Visualization, S.E. and M.M.; Supervision, L.C., H.P. and S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by "Deutsche Forschungsgemeinschaft" (DFG)—GRK 1856.

**Acknowledgments:** Simulations were performed with computing resources granted by RWTH Aachen University under the project "rwth0239". The authors are thankful to Convergent Science Inc. for providing licenses for CONVERGE.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**

#### *Appendix A.1. Standard Emission Analyzer (FEVER)*


**Table A1.** Specifications of the standard gaseous exhaust analyzer (FEVER) [18].

<sup>1</sup> Chemiluminescence detector; <sup>2</sup> non-dispersive infrared detector; <sup>3</sup> para-magnetic detector; <sup>4</sup> flame-ionization detector; <sup>5</sup> total-HC, referred to propane (C3); <sup>6</sup> time 10 to 90% (*t*10−90%).

#### *Appendix A.2. Additional Information on the Simulation Methodology*

**Mesh refinements** The mesh refinements in the CFD model, both user-defined and AMR (automatic mesh refinement), consist of a halving of base cell size *δ* of *n*ref times, with *n*ref the refinement degree (*δ*ref = *δ*/2*<sup>n</sup>* ref) [33]. Table A2 summarizes the meshing parameters of the simulation.

**Table A2.** Geometry meshing parameters. Base mesh size *δ* = 2 mm. AMR, automatic mesh refinement; EVO, exhaust valve opening; EVC, exhaust valve closing.


*<sup>a</sup>* Embedding sphere radius centered in the spark gap; *<sup>b</sup>* maximum difference between adjacent cells allowed from AMR.

**Boundary and initial conditions, heat transfer, and wall temperatures** The pressure and temperature boundary conditions were generated on the basis of the post-processing of the measured indicated pressures through a three-pressure-analysis (TPA) [41] with the software GT-POWER (v2019). Regarding the gas composition at the boundaries, on the intake side, fresh air 76.8% in mass of N2 and 23.2% of O2 were set, while on the exhaust side, the gas composition was calculated based on emission measurements. Since the real HC composition was unknown, the measured THC concentration was set at the exhaust boundary as C3H8.

The simulations were started at the second half of the exhaust stroke. The initial composition in the intake region was fresh air. In the combustion chamber and in the exhaust region, the gas composition of measured engine exhaust gas was used. However, for each simulation point, multiple consecutive cycles (≥3) were simulated in order to minimize the possible influence of initial conditions on simulated emission species.

For the heat transfer, the O'Rourke and Amsden model [43] was adopted. For the description of the temperature and velocity boundary layer, the law of the wall [44] was applied. The imposed wall temperatures on the simulation surfaces are reported in Table A3.


**Table A3.** Assumed wall temperatures in K for the different simulated operating points.

**Flow and evaporation modeling** The RANS equations were solved by using the standard *k* − [45] within CONVERGE [33]. The thermodynamic quantities of the gas were calculated as a function of temperature. The liquid fuel was considered as incompressible. The droplet evaporation was assumed as ideal and described by the Frossling correlation [46].

*Appendix A.3. Kinetic Mechanism Validation*

**Figure A1.** Burning velocities of a commercial gasoline TAE7000and its surrogate (13.7% n-heptane, 42.9% iso-octane, 43.4% toluene) with air mixtures with 15% ethanol addition. Symbols denote the experimental measurements by Dirrenberger et al. [47] at 1 atm and 358 K. Solid lines show the numerical results for the present kinetic model. Comparison analogous to that presented by the authors in [31].

*Appendix A.4. Comparison of NO Production between Kinetics and Equilibrium Chemistry*

**Figure A2.** Comparison of in-cylinder NO-production prediction with a 0D two-zone GT-POWER model (v2019) using chemical kinetics (same mechanism used in 3D-CFD and provided in the Supplementary Materials, red continuous line) and chemical equilibrium calculations (from [48]).

#### *Appendix A.5. Experimentally Determined FID Factors*


**Table A4.** FID response factors measured for selected species for the FEVER device and the FFID device (two lines). The theoretical response factor *f*thFID,*<sup>i</sup>* (s.Equation (2)) is compared with the measured one, and the relative deviation is reported [25].

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
