*2.2. Numerical Model*

The open software OpenFOAM was employed in the present work. The mesh is indicated in Figure 2. In order to implement the movement of the piston and valves, a deforming mesh was used. In particular, Figure 2a represents the tridimensional mesh, Figure 2b a cross-section at BDC (bottom dead center), i.e., 180◦ or 540◦ CA ATDC, and Figure 2c a cross-section at TDC, i.e., 0◦ or 360◦ CA ATDC. Several meshes with di fferent elements were tested in order to verify that the results are independent of the mesh size. Table 2 indicates the error obtained between experimental and numerical results of pressure and fractions using a mesh with 501,769 elements at BDC (mesh 1), as well as 802,527 elements (mesh 2) and 1,264,873 (mesh 3). As can be seen, there is no di fference between the meshes 2 and 3. For this reason, the mesh 2 was chosen for the present work.

**Figure 2.** (**a**) Tridimensional mesh at BDC; (**b**) Cross-section mesh at BDC; (**c**) Cross-section mesh at TDC.


**Table 2.** Error (%) at 100% load obtained using di fferent mesh sizes.

A new in-house solver was programmed using C++. Briefly, this solver is based on the RANS (Reynolds-averaged Navier Stokes) equations of conservation of mass, momentum, and energy. The k-ε turbulence model was chosen.

The standard Kelvin-Helmholtz and Rayleigh-Taylor breakup (KH-RT) model [37] was employed for fuel droplet breakup, and the Dukowicz model [38] for the heat-up and evaporation. A comprehensible analysis about the adequacy of breakup models can be found in the literature. Compared to other breakup models such as WAVE, TAB (Taylor Analogy Breakup), etc., the KH-RT model is more suitable for the high injection pressures that take place in diesel engines [39]. As indicated previously, ammonia and water injection were modeled through an injector equipped with a dual nozzle with separate needles for water/ammonia and fuel.

In order to solve the chemical kinetics, a reaction mechanism was programmed by adding the three kinetic schemes described in Sections 2.1–2.3 for combustion (131 reactions and 41 species), NOx formation (43 reactions and 20 species) and NOx reduction (131 reactions and 41 species), respectively. Several additional equations must be added to model chemical kinetics. Given a set of *N* species and *m* reactions, Equation (1), the local mass fraction of each species, *fk*, can be expressed by using Equation (2).

$$\sum\_{k=1}^{N} v\_{kj}^{\prime} \mathcal{M}\_k \xleftarrow{k\_j} \sum\_{k=1}^{N} v\_{kj}^{\prime\prime} \mathcal{M}\_k \ j = 1, 2, \dots, \; m \tag{1}$$

$$\frac{\partial}{\partial t}(\rho f\_k) + \frac{\partial}{\partial \mathbf{x}\_l}(\rho u\_l f\_k) = \frac{\partial}{\partial \mathbf{x}\_l} \left(\frac{\mu\_l}{\mathbf{S}c\_l} \frac{\partial f\_k}{\partial \mathbf{x}\_l}\right) + S\_k \tag{2}$$

In the equations above, ρ is the density, *u* the velocity, *vkj* the stoichiometric coefficients of the reactant species *Mk* in the reaction *j*, *vkj* the stoichiometric coefficients of the product species *Mk* in the reaction *j, Sct* the turbulent Smidth number, and *Sk* the net rate of production of the species *Mk* by chemical reaction, given by the molecular weight multiplied by the production rate of the species, Equation (3).

$$S\_k = M\mathcal{W}\_k \frac{d[\mathcal{M}\_k]}{dt} \tag{3}$$

where *MWk* is the molecular weight of the species *Mk* and [*Mk*] its concentration. The net progress rate is given by the production of the species *Mk* minus the destruction of the species *Mk* along the *m* reactions:

$$\frac{d[\mathcal{M}\_k]}{dt} = \sum\_{j=1}^{m} \left\{ (\upsilon'\_{kj} - \upsilon'\_{kj}) \left[ k\_{fj} \prod\_{k=1}^{N} \left[ \mathcal{M}\_k \right]^{\upsilon'\_{kj}} - k\_{bj} \prod\_{k=1}^{N} \left[ \mathcal{M}\_k \right]^{\upsilon'\_{kj}} \right] \right\} \tag{4}$$

where *kfj* and *kfb* are the forward and backward reaction rate constants for each reaction *j*.

The simulation started at 360◦ CA ATDC and the whole cycle was analyzed. The initial pressure was 1.31 bar, obtained from experimental measurements. Concerning boundary conditions, the intake valve pressure and temperature after the turbocharger were 2.72 bar and 514 K, respectively. The heat transfer from the cylinder to the cooling water was modeled as a combined convection-radiation type, Equation (5). Previous investigations [40] demonstrated the accuracy of this type of boundary condition in comparison with adiabatic or constant temperature:

$$q = h(T\_{\text{gas}} - T\_{\text{water}}) \tag{5}$$

where *q* is the heat transferred, *Tgas* the in-cylinder temperature, *Twater* the cooling water temperature (78 ◦C), and *h* the heat transfer coefficient, given by the following expression [41]:

$$h = 10.4 \, k b^{-1/4} \, (u\_{piston}/\nu)^{3/4} \tag{6}$$

where *b* is the cylinder bore, *k* the thermal conductivity of the gas, *upiston* the mean piston speed, and ν the kinematic viscosity of the gas. Substituting values into the above equation yields *h* = 4151 <sup>W</sup>/m2K.
