*Motivation*

The present work aims to (1) introduce and validate an improved version of the tabulated chemistry solver initially presented in [33]; (2) assess the solvers' combustion and emissions predictions under Diesel and gasoline engine conditions within the 0-D SRM framework; (3) demonstrate the computational e fficiency of the method, which allows large design optimization studies while taking detailed chemistry e ffects into account. With respect to the analyses discussed in [33], the present work includes numerous performance and accuracy refinements on both the tabulation and the engine simulation codes, phenomenological formulations of the turbulence models for CI and SI engines, a novel approach for thermal NO source terms tabulation and a method validation under SI engine conditions in addition to a broader CI engine simulation campaign. To the authors' knowledge, a comprehensive SRM-based tool chain featuring tabulated chemistry, physical turbulence models and computational e fficiency comparable to multi-zone models (i.e., [7]) has not been demonstrated before.

The paper is structured as follows: In Section 2, a concise description of the SRM framework is given together with a description of the turbulence models and progress variable-based solver. In Section 3, an accuracy assessment investigation of the refined storage/retrieval strategy under zero-dimensional constant pressure reactors in shown. The validation of the development is demonstrated in Sections 4–6. Online and tabulated chemistry-based solvers are compared using the SRM under Diesel (Section 4) and spark ignited engine (Section 5) conditions. Finally, the accuracy of the introduced solvers is discussed with focus on combustion and emissions predictions along with remarks on computational performances in Section 6. For validation, the following criteria were defined: (1) Is the solver accuracy consistent with the general model accuracy? (2) Is the CPU time of the developed software acceptable for engine development, engine optimization and driving cycle analysis?

## **2. Simulation Methods**

The SRM has been coupled with two di fferent chemistry solvers: (1) online, where chemical source terms are calculated during run-time; (2) Combustion Progress Variable (CPV), which uses a pre-calculated look-up table to retrieve the source terms for combustion as well as thermal NO formation. In the first sub-section, a brief description of the SRM engine modelling framework is presented for better readability. In the following sub-sections, aspects of the di fferent phenomenological turbulence models adopted for the SI and CI simulation campaigns are briefly presented as well as a description of the tabulated chemistry solver.

#### *2.1. The Stochastic Reactor Model*

In the SRM computational domain, the in-cylinder trapped mass is discretized as an ensemble of notional particles that can mix with each other and exchange heat with the cylinder walls. A schematic visualization of the concept, along with an exemplary distribution of the particles' property (i.e., enthalpy or gas composition) is shown in Figure 1. Depending on the initial conditions assigned at Intake Valve Closure (IVC), a given chemical composition, temperature, and mass are assigned to each particle, while pressure is assumed to be the same across all particles.

**Figure 1.** Schematic visualization of the Stochastic Reactor Model (SRM) concept; (**a**) engine physical space; (**b**) 0-D particles in the SRM computational space and an exemplary distribution of particle properties at a given time-step t.

A probability density function is used to describe the in-cylinder content and each particle contributes to the realization of the given PDF at each time-step. Since all stochastic particles in the SRM represent a portion of the in-cylinder mass, a Mass Density Function (MDF) is used to solve the main transport equation of the 0-D SRM. The joint scalar vector of the MDF can be expressed as reported in Equation (1).

$$\mathcal{F}\_{\Phi}(\boldsymbol{\upmu},\mathbf{t}) = \mathcal{F}\_{\Phi}(\boldsymbol{\upmu}\_{1'}, \dots, \boldsymbol{\upmu}\_{N\_{\rm Sc}}; \mathbf{t}) \tag{1}$$

The vector ψ is the realization of the vector of local scalar variables noted herein as **Φ**, while NSc is the number of scalars transported in the computational domain. Depending on the chemistry solver used in the simulation framework, the number and type of transported scalars varies significantly.

$$\text{Online chemistsy solver } \Phi(\mathbf{t}) = \begin{pmatrix} \mathbf{Y}\_1 \ \mathbf{Y}\_2 \ \dots \ \mathbf{Y}\_{\text{N}\_{\text{Sp}}} \text{ h} \ \mathbf{t} \end{pmatrix} \tag{2}$$

$$\text{Tabulated chemistry solver (CPV) } \Phi(\mathbf{t}) = \left(\phi\_{\prime} \text{ h}\_{\prime} \text{ w}\_{\text{q}} \text{ } \text{T}\_{\text{p} \prime} \text{ yEGR, h}\_{298} \text{:} \mathbf{t}\right) \tag{3}$$

In the online chemistry solver formulation (see Equation (2)), the specific enthalpy (h) as well as the full vector of species mass fractions **Y** must be transported in order to correctly compute the chemical reactions during run time. The size of the mass fraction vector is therefore dependent on the number of species (NSp) defined in the chemical mechanism adopted. With respect to the tabulated chemistry solver, only a reduced set of quantities need to be transported in the computational domain independently on the size of the mechanism used to generate the look-up table. These are the equivalence ratio (φ), the specific enthalpy (h), the mean molecular weight (wq), the vector (**Tp**) containing the thermodynamic polynomial coefficients, the EGR mass fraction (yEGR) and the latent enthalpy (h298) of the in-cylinder mixture. Once initial conditions are assigned, the transport equation of the joint scalar vector of the MDF is numerically solved using a Monte Carlo particle method (e.g., Pope [34]) with the operator splitting technique as previously presented by Maigaard et al. [35]. A series of sub-models are employed to sequentially solve the different physical and chemical processes occurring in the combustion chamber, as summarized in Equation (4):

$$\frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} = \left. \frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} \right|\_{\Delta V} + \left. \frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} \right|\_{\text{lir} \mathbf{j}} + \left. \frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} \right|\_{\text{FP}} + \left. \frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} \right|\_{\text{mix}} \left. \frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} \right|\_{\text{chem}} + \left. \frac{\partial \mathbf{F}\_{\Phi}}{\partial \mathbf{t}} \right|\_{\text{ht}} \tag{4}$$

In the equation above, the subscript ΔV represents the piston movement, inj fuel injection, FP flame propagation (considered only for SI engine simulations), mix the turbulence and particle interaction process, chem the chemical reactions and ht heat transfer to the walls. Detailed descriptions of the

sub-models used for piston movement, pressure correction, online chemistry solver and fuel injection can be found in the work from Tuner [36]. The treatment of the flame propagation has been previously introduced by Bjerkborn et al. [37] and broadly validated by Pasternak et al. [38] and Netzer [39]. The Woschni model [40] is used to evaluate total wall heat transfer, while the distribution of the heat transfer over the SRM particles is calculated using a stochastic approach, explained by Tuner [36] and further validated by Franken et al. [41] using Direct Numerical Simulation (DNS) results from Schmitt et al. [42]. A short overview of the turbulence models adopted in the present work is given in the following sub-section.

#### *2.2. Phenomenological Turbulence Models and Particle Interaction*

The expression of the mixing term in Equation (4) is presented below:

$$\frac{\partial \mathbf{F}\_{\Phi}}{\partial t}\Big|\_{\text{mix}} = \frac{\mathbf{C}\_{\Phi}\beta}{\pi} \left[ \int\_{\Delta\psi} \mathbf{F}\_{\Phi}(\psi - \Delta\psi, \mathbf{t}) \mathbf{F}\_{\Phi}(\psi + \Delta\psi, \mathbf{t}) \mathbf{d}(\Delta\psi) - \mathbf{F}\_{\Phi}(\psi, \mathbf{t}) \right] \tag{5}$$

C**Φ** and β are two model parameters that in the present study have been set to 2.0 and 1.0, respectively. The mixing time history, τ in Equation (5), is the main input parameter of the SRM. The mixing time is needed to model the turbulent mixing occurring in the combustion chamber due to various phenomena such as: fuel injection, swirl motion, tumble motion, etc. The shape and magnitude of the mixing time profile controls how intense the SRM particle mixing process is. Since τ influences mixture inhomogeneities in both composition and enthalpy spaces, a strong influence on the auto-ignition process, the local rates of heat release and pollutant formation are seen when the mixing time is changed.

A simple approach to construct the mixing time history is by using a set of empirical constants as done by Pasternak et al. [43,44] for Diesel combustion studies. A more comprehensive approach is to calculate τ during run-time by employing a k-ε or K-k based turbulence model. Depending on whether Diesel or gasoline combustion are considered, different approaches have been implemented in the present work. For Diesel combustion, the approach initially proposed by Kožuch [45] was adopted and validated by Franken et al. [46] across a large set of operating points. For gasoline engine conditions, the K-k model proposed by Dulbecco et al. [47] was implemented and validated in the SRM framework for different gasoline fuel surrogate mixtures in [48,49]. Examples of the mixing time histories that are computed with the mentioned turbulence models are shown in Figure 2 as function of crank angle degrees (CAD) assuming 0 as firing top dead center (TDCf).

**Figure 2.** Exemplary mixing time histories computed using turbulence models for (**a**) Diesel engine conditions [47] and (**b**) gasoline engine conditions [49].

Once the mixing time history has been computed, an additional sub-model is employed to decide which and how many particles are selected in each mixing step. In this work the Curl model [50] has been used to describe the particle interaction for both the Diesel and the gasoline engine simulations.

#### *2.3. Combustion Progress Variable Model*

The CPV model relies on the pre-tabulation of the auto-ignition and emission formation processes for a broad range of initial conditions in terms of exhaust gas recirculation rates (yEGR), pressure (p), unburned temperature (Tu) and equivalence ratio ( φ). The reconstruction of the thermo-chemical state during run-time is then performed by means of an appropriately chosen progress variable (herein noted as C). As discussed in detail by Niu et al. [51] and Ihme et al. [52], any progress variable-based method needs to fulfill the following conditions to be mathematically sound: (1) strict monotonicity with time; (2) non-negativity and normalization (C = 0 for unreacted mixture and C = 1 for fully burned mixture); (3) if reactive scalars (i.e., chemical species) are used in the definition of the progress variable, said scalars should evolve on comparable time scales [52].

Several progress variable definitions have been proposed by di fferent groups for both homogeneous and flamelet-based tabulation frameworks [10,13,21]. The most common approach consists in formulating the progress variable via a combination of chemical species (typically including the surrogate fuel molecule, O2, CO, CO2, H2O, OH, CH2O and others). The choice of the species to use and their weighting factors within the progress variable definition is usually optimized so that both low and high temperature combustion regimes are captured [51,52]. In the present work, latent enthalpy (enthalpy of formation at standard state, h298) is used to define the reaction progress variable C as reported below in a normalized fashion.

$$\mathbf{C}(\mathbf{t}) = \frac{\mathbf{h}\_{298}(\mathbf{t}) - \mathbf{h}\_{298,\mathbf{u}}}{\mathbf{h}\_{298,\text{max}} - \mathbf{h}\_{298,\mathbf{u}}} \tag{6}$$

In Equation (6), h298 is the latent enthalpy calculated at 298 K, and summed over all species, subscript u denotes unburned state, and max denotes the most reacted state, which is assumed to be where the maximum of the accumulated chemical heat is released. As discussed in [15,33], this progress variable can be used to track both low and high temperature reactions, which is important when tabulation methods are developed for fuels exhibiting cool flames. Figure 3 shows an exemplary h298 profile together with temperature for a constant pressure reactor calculation.

**Figure 3.** Temperature (black) and latent enthalpy at 298 K (red) as function of time for a constant pressure reactor calculation at 10 bar and 750 K for stochiometric *n*-heptane/air mixture.

It can be discussed that in a constant pressure reactor, progress variables based on temperature (i.e., as proposed by Knop et al. [24]) or on latent enthalpy will be equal. However, unlike temperature the latent enthalpy will not be influenced by pressure changes, fuel vaporization or mixing. As h298 is a conserved scalar, a transport equation, as well as a spray source term, can be easily formulated. For the majority of the tabulated conditions, the maximum of the accumulated chemical heat release coincides with the chemical equilibrium state. However, for fuel-rich conditions as well as for states that lead to pronounced endothermic reactions (i.e., high pressure or high EGR conditions), the mixture relaxation towards equilibrium is not taken into account by Equation (6). Nevertheless, this definition was considered acceptable for engine applications in the SRM since it a ffects a small fraction of the particles for a limited portion of the cycle in direct injected-engine simulations. Over the whole range of computed conditions, this definition was found to be convenient to correctly capture both low and high temperature combustion regimes when generating the table through adiabatic homogeneous constant pressure reactor calculations. To keep the table size within a feasible range for engine applications, storage of the progress variable source terms was done over 25 points for each reactor simulation.

As introduced in [33], the coupling of the CPV solver with the SRM required substantial changes to the pre-existing online chemistry-based code. As discussed in the previous section, the tabulated chemistry solver requires only a reduced set of scalars as opposed to the full vector of species mass fraction needed for the online solver. It is important to notice that the normalized reaction progress, as reported in (6), is not transport directly. The latent enthalpy ( h298) is transported instead and the progress variable is calculated, whenever a table look-up is needed. This practice facilitates the formulation of the fuel injection model and the resulting evaporation source term and the particle–particle interaction sub-model. Beside the treatment of the chemistry step, the same set of sub-models noted in Equation (4) and their relative constants were applied for both the tabulated chemistry runs the online chemistry simulations. Hence potential di fferences between the two solver solutions are to be interpreted as due to the chemistry step only.

#### **3. Look-Up Table Generation and Testing**

Prior to the engine simulation campaigns, a grid density investigation together with an assessment of the interpolation error was performed so that the CPV model could be verified. In the following sub-sections, the details of the table construction process as well as the outcome of the interpolation error analysis are discussed.

#### *3.1. Look-Up Table Generation*

Detailed chemistry simulations under adiabatic constant pressure reactor conditions were performed using the software LOGEtable [53]. Initial mixture unburned temperature ( Tu), pressure (p), equivalence ratio ( φ) and mass fraction EGR (yEGR) are the input variables of the look-up table generator. Table grid points were chosen so that the typical range of thermodynamic conditions encountered in conventional direct injected engines that can be expected to be found in conventional direct injected engines. The table grid adopted in this work relies consists of 205,200 points which are distributed across tabulation parameters as outlined in Table 1.


**Table 1.** Tabulation grid adopted for all tabulated chemistry runs.

Except for the EGR mass fraction space, the points are distributed in a non-equidistant manner so that typical engine relevant conditions (i.e., around stoichiometry in equivalence ratio space) are better resolved. The EGR stream is assumed to have a fixed composition comprising combustion products (CO2, H2O and N2) of the given fuel mixture at stoichiometry.

For each of the table grid points, the progress variable C introduced in Equation (6) is used to trace the reaction trajectory, and its source terms (dC/dt) are stored in the table at 25 di fferent points between unburned C = 0 and fully burned C = 1.0 state. In addition, mean molecular weight, the thermodynamic polynomial coe fficients and any chemical species that the user decides to monitor are also stored across the mentioned progress variable grid. In the present work, the following species have been included in the tabulation process so that major engine-out emissions and standard engine performance analyses could be done: fuel molecule, i.e., n-C7H16, O2, N2, CO2, H2O, CO, H2, C2H2, C2H4 and unburned hydrocarbons (uHC) (defined as the sum of all species containing an H and C or H, C and O atoms). The resulting size on disk was approximately 1 GB.

#### *3.2. Table Density*/*Accuracy Investigation*

Before applying the CPV solver to engine simulations in the SRM, the interpolation strategy for the progress variable source term retrieval was verified. For each of the data entry points included in the tabulation grid (see Table 1), the auto-ignition solution from the constant pressure reactor calculation using the online chemistry solver has been compared to the respective solution retrieved from the table using linear interpolation. The discrepancy between the two solution has been quantified based on the mean relative di fference of the vector composed by the errors at 5, 10, 50 and 90 percent of energy released location (noted as E5%, E10%, E50% and E90% in Figure 4) at all given reactor points.

**Figure 4.** Online and tabulated chemistry solutions for the progress variable (based on Equation (6)) evolution of a stoichiometric *n*-heptane/air mixture at 10 bar and 750 K.

The global interpolation error has been then defined based on the mean square di fference (MSD) of the location of E5%, E10%, E50% and E90% points. In Equations (7), where i = 5, 10, 50, 90, the exact definition of the interpolation error is reported.

$$
\varepsilon\_{\rm i} = \frac{\mathbf{E}\_{\rm i, \rm Tabulated} - \mathbf{E}\_{\rm i, \rm ordine}}{\mathbf{E}\_{\rm i, \rm ordine}} \times 100 \,\mathbf{\tilde{z}} = \sum\_{i=1}^{\rm N\_{\rm crru}} |\varepsilon\_{\rm i}| / \mathbf{N\_{\rm crrus}} \,\varepsilon\_{\rm MSD} = \sqrt{\sum\_{\rm i} \left( \overline{\varepsilon} - |\varepsilon\_{\rm i}| \right)^{2}} \tag{7}
$$

The interpolation errors have been assessed for a large set of operating conditions covering lean to rich ( φ = 0.6−3.0) conditions as well as from low to high unburned temperatures (T = 600–1400 K) under constant pressure reactor conditions. *n*-Heptane was used as fuel molecule, and the mechanism developed by Zeuch et al. [54] was employed. The main goal of this campaign was to quantify the accuracy of the interpolation routine under the simplest reactor conditions, so to have a higher confidence in results assessment during the subsequent engine simulation campaign. In the SRM simulation, numerous sub-models are coming into play, and hence, error compensation e ffects may arise. In Figure 5, exemplary results of the error quantification campaign are displayed in pixel plot format for conditions at 1 and 35 bar pressure.

**Figure 5.** Error maps for the interpolation error (<sup>ε</sup>MSD as defined in Equation (7)) between online and tabulated chemistry solver under constant pressure reactor conditions (**a**) at 1 bar and (**b**) 35 bar, both at 0% exhaust gas recirculation using *n*-heptane as fuel.

For a more effective error visualization as well as to better highlight the exact areas where interpolation errors may arise, pixel contours have been conditionally formatted so that every case having an error below 0.5% is displayed in white. All cases presenting an error above 0.5% will have contours consistent with the RGB color bar instead. For the whole range of computed conditions, the solution retrieved from the table did not exceed a 6% discrepancy on a single point basis, while the overall average error is approximately 1.5%. As expected, interpolation errors are more visible in the Negative Temperature Coefficient (NTC) regimes, which is the most challenging area to parametrize when progress variable methods are concerned [24,55]. While a tighter tabulation grid in temperature space could have delivered even lower errors, the authors considered this level of accuracy to be acceptable. This decision was based on the globally low errors but also on method usability considerations.

A tighter table grid results in a larger file size on disk, which in turn affects the random-access memory (RAM) requirement for the engine simulations during run-time. Hence, the present configuration was considered to be a good trade-off.

#### **4. Compression Ignition Engine Simulations**

In this section, results of a heavy-duty Diesel engine simulation campaign are presented and discussed. In the first sub-section, experimental data and simulation setups are presented. Secondly, a result comparison between experimental data and the two chemistry solvers is shown and discussed with respect to engine performance parameters and engine-out emissions. All the simulation results presented in this section were obtained using the commercial software LOGEengine version 3.2.1 [53].

#### *4.1. Engine Data and Simulation Setup*

The simulation setups were constructed based on experimental data from a heavy-duty (HD) 13.0 L Diesel engine that features a direct injection system capable of injection pressures up to 2000bar. Although the engine has an in-built external EGR system, the analyzed engine conditions did not include external EGR. A total of 10 operating points from 1000 rpm to 1700 rpm and 6 bar to 22 bar indicated mean effective pressure (IMEP) are used for the present simulation campaign. Details of the operating conditions are outlined in Table 2.


**Table 2.** Heavy-duty engine operating conditions.

As presented in Figure 6, operating point 10 features a double injection rate profile (pilot + main) while the rest have a single injection event. The crank-angle resolved pressure profiles were measured for one cylinder and used to calibrate the SRM engine model.

**Figure 6.** Fuel injection rates as a function of crank angle degree (assuming 0 as firing TDC) for operating conditions (**a**) from 1–5 and (**b**) 6–10, as listed in Table 2.

Commercial Diesel fuel was used during experiments, whereas in the simulations a blend of *n*-decane and α-methylnaphthalene (75–25% mass fractions respectively, also noted as IDEA surrogate fuel) was employed as surrogate fuel model. A liquid fuel properties comparison is summarized in Table 3.


**Table 3.** Liquid properties of the experimental and surrogate fuel mixture used in the heavy-duty Diesel engine simulation campaign.

The chemical kinetic scheme employed in this study has been taken from the LOGEfuel database [53]. The mechanism is an improved version of the detailed model from Wang [56]. It features oxidation models for *n*-decane, α-methyl-naphthalene and methyl decanoate as main fuel species as well as a detailed PAH growth mechanism [57] and thermal NOx model. The detailed reaction scheme was reduced to a size of 189 species using the method described in [12].

The measured pressure history was analyzed using the thermodynamic analysis of LOGEengine [53]. This procedure provided chemical kinetics-based estimations of wall temperature, in-cylinder temperature at IVC, internal EGR fraction and the apparent rate of heat release (RoHR). In Tables 4 and 5, the main SRM model parameters and the calibrated k-ε model constants are presented, respectively [45].


**Table 4.** SRM main model settings for the heavy-duty Diesel engine simulation campaign.



The SRM model calibration for the presented operating conditions was carried out using the procedure described in [58], and the Curl [50] particle interaction sub-model was used. To ensure consistency during chemistry solver comparisons, the same set of model parameters and constants were applied to both the online and tabulated chemistry solver runs without any re-calibration.
