Numerical Study of Natural Gas/Diesel Reactivity Controlled Compression Ignition Combustion with Large Eddy Simulation and Reynolds-Averaged Navier–Stokes Model

In the current study, a comparative study is performed using Large Eddy Simulation (LES) and Reynolds-averaged Navier–Stokes (RANS) turbulence models on a natural gas/diesel Reactivity Controlled Compression Ignition (RCCI) engine. The numerical results are validated against the available research work in the literature. The RNG (Re-Normalization Group) k− ε and dynamic structure models are employed to model turbulent flow for RANS and LES simulations, respectively. Parameters like the premixed natural gas mass fraction, the second start of injection timing (SOI2) of diesel and the engine speed are studied to compare performance of RANS and LES models on combustion and pollutant emissions prediction. The results obtained showed that the LES and RANS model give almost similar predictions of cylinder pressure and heat release rate at lower natural gas mass fractions and late SOI2 timings. However, the LES showed improved capability to predict the natural gas auto-ignition and pollutant emissions prediction compared to RANS model especially at higher natural gas mass fractions.


Introduction
In last decade, the use of new and advanced combustion strategies along with alternative fuels in internal combustion engines (ICEs) are of global interest to achieve higher fuel economy and lower pollutant emissions. It has been reported that simultaneous reduction in particulate matter (PM) and nitrogen oxides (NOx) emissions can be obtained with low temperature combustion strategies like reactivity controlled compression ignition (RCCI) combustion [1][2][3]. In RCCI combustion, in-cylinder reactivity gradient is generated by port-injection of a fuel with low reactivity, like natural gas or gasoline and direct-injection of a high reactivity fuel, like diesel which enables combustion and heat release rate control of the engine. Generally, RCCI combustion results in lower NOx and PM emissions and higher gross indicated efficiency (GIE) compared to conventional diesel combustion (CDC) and homogenous charge compression ignition (HCCI) combustion [4][5][6][7][8][9][10]. There are comprehensive and thorough review papers in the previous studies about RCCI combustion, that can be referred for more detailed information regarding various aspects of RCCI combustion [11,12].

Computational Model
In the current research RCCI engine simulations carried out utilizing the CONVERGE CFD tool [32]. Port fuel injected hydrogen and natural gas proposed to be homogeneously combined at IVC and the blob injection model was used to simulate the diesel injection procedure [33]. The primary and secondary breakups of the injected fuel were modelled employing the hybrid KH-RT model [34]. The No Time Counter (NTC) collision approach [35] was utilized for collision of fuel droplets and droplet drag was described using the dynamic drag model [36]. The RNG (Re-Normalization Group) k − ε model with wall-functions was employed to model turbulent flow [37,38]. A well-defined chemistry method utilized via the SAGE solver for combustion modelling. Tetradecane (C14H30) was used as a fuel surrogate to model physical characteristics of diesel fuel.
The simulations conducted utilizing the multi-zone chemistry solver to decrease the run time [39]. A 76 species and 464 reaction mechanism of GRI-Mech 3.0 was utilized [40]. The NOx formation was modelled utilizing an extended Zeldovich procedure [41]. Hiroyasu-NSC soot model was employed to predict soot according to the method of Hiroyasu [42] which uses acetylene as an inception species, allowing coupling of the soot model and the detailed chemical kinetics solver. A crevice model based on the model of Namazian and Heywood [43] was used which is a substitute to analysing the crevice regions in the CFD grid directly. The model includes two rings and five crevice regions and calculates the relevant parameters needed to determine the amount of mass entering and leaving the cylinder. A flowchart representing the overview of the solvers and models was depicted in Figure 1.

Validation of Model
To improve computational efficiency, closed-cycle calculations on sector grids via periodic boundaries are conducted. The liner temperature of 450 K, piston temperature of 550 K and head temperature of 525 K is taken in the simulation for heat transfer model. The law-of-wall boundary condition is applied to the boundaries like piston, liner and head. The grid represents one of the six holes of the direct injector, so a sixty-degree sector geometry is created. The CFD code uses a structured Cartesian grid with base cell size of 1.0 mm for RANS and LES calculations as reported in previous studies [25,31]. Boundary, nozzle and cylinder embedding along with adaptive mesh refinement (AMR) are employed to provide mesh resolution only when and where it is needed to minimize the cell count. Injector embedding was used since the amount of injected diesel is low and the injector angle is also towards the squish region, it is not so distinguishable. Fixed embedding scale of four at the piston and cylinder head boundaries and near the injector nozzle and an AMR scaling of four according to fixed thresholds for the spatial gradients in velocity and temperature are used (see Figure 2). The numerical findings are validated against the Nieman et al. [17] results. The specifications of Caterpillar Single-Cylinder Oil Test Engine (SCOTE) 3401E are presented in Table 1. This heavy-duty diesel engine is a 2.4 L single-cylinder version of the Caterpillar 3406 six-cylinder, on-road engine and representative of Caterpillar's engine technology for the 1994-1997 model years. Validations are performed at low load 4 bar indicated mean effective pressure (IMEP) and medium load 9 bar IMEP as demonstrated in Table 2.  Figure 3 shows cylinder pressure traces and heat release rate at two engine loads compared to Reference [17] using RANS turbulence model. As observed, the prediction of combustion phasing and pressure traces are good and the model also captures the low temperature heat release region at −24 and −23 for the two loads which is the result of the energy release of the decomposition of early injected pilot fuel (low temperature reactions). The absolute errors in the peak in-cylinder pressure are 0.3956 and 0.0842 MPa and the absolute errors in the relevant crank angle at the peak in-cylinder pressure are 1.2448 and 1.15 CA for the IMEP 4 and 9 bar, respectively. The errors for other outputs are included in Table 3. As presented in Table 3, the model has successfully predicted emissions and performance parameters trends at medium and low loads of RCCI engine. It should be noted that the difference in the emission results is maybe due to the difference in the used chemical mechanism.

LES Model Description
Traditionally, in RANS models an effective turbulent viscosity is used to model the Reynolds stress term. Thus, the turbulent convective mixing is modelled by additional turbulent diffusion (i.e., diffusive mixing). The common RANS model which has been extensively used for turbulence modelling in internal combustion engines is Rapid Distortion RNG k − ε [37]. A key difference between LES and RANS models is how the fields are decomposed for modelling. For a RANS approach, the field is decomposed into an ensemble mean and a fluctuating component. In the LES approach, the field is decomposed into a resolved field and a sub-grid field as follows: where the over-bar indicates the resolved field and the prime indicates the sub-grid field. The resolved velocity field is defined as a spatial average of the actual velocity field. This differs from the RANS approach, where the mean velocity field is an ensemble average. Because of these decomposition differences, the LES filter has properties, unlike RANS, given by, If the LES decomposition is applied to the conservation of momentum equation, the following LES equation can be derived: where, Most LES models focus on modelling the expression for the sub-grid stress tensor, ij, given above. However, it is also possible to allow upwinding to be used as the LES model. In this case, no additional term is added for the sub-grid stress tensor [32].

Dynamic Structure Model (One-Equation)
In CONVERGE, there are two classes of LES models; zero-equation and one-equation models. For zero-equation models, CONVERGE does not solve any additional transport equations. For one-equation models, an additional transport equation is added for the sub-grid kinetic energy [32].
In the dynamic structure model, a transport equation is added for the sub-grid kinetic energy to enforce a budget on the energy flow between the resolved and the sub-grid scales. The dynamic structure model does not use turbulent viscosity to model the sub-grid stress tensor [44]. The sub-grid turbulent kinetic energy needs the sub-grid stress tensor models, which are given by: where the test level kinetic energy is defined by: The test and grid level kinetic energies are related by the trace of the Leonard term, thus another transport equation for K is not required: Substituting these models for the two stress tensors into the Germano identity yields the following: Alternatively, in an algebraic model, in which the tensor coefficient is removed from the integral and then solved for ij, the model for the sub-grid tensor becomes: Sufficient grid resolution would be used to resolve the flow into the viscous sub-layer, when running an LES simulation. In this case, walls are modelled directly with no-slip boundary conditions and Dirichlet temperature conditions. However, in many cases it is not feasible to add enough resolution to resolve the viscous sub-layer and thus wall models must be employed. The Werner and Wengle [45] wall model is designed to work with LES models. This model is more complicated than the standard law-of-the-wall but it does not require iterations. Like the other wall models there are two regimes. If the magnitude of the tangential velocity satisfies the inequality given as, Then the shear speed is given as, where A is 8.3 and B is 1/7. If the inequality is not satisfied, then the shear speed is given by, Once the shear speed is calculated, the stress tensor can be determined as, In addition, the Amsden [46] temperature law-of-the-wall is used for energy law-of-the-wall. Regarding the run time for LES simulations, it took longer time compared to RANS model.

Effect of Natural Gas Mass Fraction
The medium load case (i.e., 9 bar IMEP) in Table 2 was selected as the base case condition in the current paper. Figure 4 depicts the cylinder heat release rate (HRR) and pressure with different natural gas percentages for LES and RANS models. As it is observed, when mass fraction of natural gas increases, combustion delays, heat release rate decreases and ignition delay (ID) increases. Increasing the diesel fraction increases the local fuel reactivity and so, results in advanced start of combustion (SOC) and decreased ignition delay, which improves combustion efficiency but it has a detrimental effect on soot (see Figure 5). It can be seen that in low natural gas mass fractions, both models have insignificant difference but this difference between two models is much more distinctive with increasing of natural gas percentage. Figure 5 shows ignition delay, gross indicated efficiency (GIE), ringing intensity (RI) and CA50 variations for RANS and LES models. It is apparent that wherever natural gas mass fraction increases, the ignition delay rises [47] and the LES and RANS give similar trend. The detailed discussion about this behaviour can be found in Ref. [14]. In addition, LES model predicts misfiring at 89% natural gas mass fraction and descending RI trend compared to RANS, thus that LES is much more reliable in combustion simulation of natural gas/diesel RCCI engine.  In-cylinder temperature contours are shown in Figure 6 in order to discuss more details of the in-cylinder combustion processes in RCCI engine. As known, the early-injected diesel targets the squish region whereas the late-injected diesel targets the piston bowl and plays as an ignition source, which results in the high temperature combustion happening in both bowl and squish regions. It is apparent from Figure 6 that the ignition zone of 89% natural gas case was smaller than that of 85% natural gas case, indicating that with by increasing natural gas mass fraction, ignition delay increases. According to Figure 6, more precise temperature contours are captured by the LES model due to the point that LES can specify more detailed turbulence structures, fuel/air combination and temperature fluctuations. This is much more distinctive in the chemically dominated premixed auto-ignition step. The main drawback of the RANS model is its over-estimation of the turbulent kinetic energy, so, the RANS results in a larger Damköhler number compared to the LES and leads to a lower reaction ratio [31].
Comparison of engine-out emissions for LES and RANS models is depicted in Figure 7. As shown, by increasing natural gas mass fraction ratio, the CO emission decreases. Low CO is achieved when ignition happens near TDC (i.e., low CA50). UHC emission slightly increases at higher natural gas mass fraction in the case of LES model due to the part of fuel not burning resulting of the high combustion delay time. In addition, previous studies have showed that the UHC depend more on crevices and squish volumes in RCCI combustion. As shown in Figure 7, NOx emissions decrease by increasing natural gas fraction due to the in-cylinder temperature drop along with reduction of the residence time of the mixture. Soot emission also decreases by increasing of the natural gas, since the amount of soot produced is highly related to the mass of direct injected diesel fuel and mainly depends on the local equivalence ratio. The detailed discussion about this behaviour can be found in reference [14]. As can be seen, both LES and RANS models successfully captured the emission trends but LES predicts them more precisely than RANS particularly at high natural gas mass fractions.

Effect of SOI2 Timing
In this section, SOI2 timing is varied in the range of −60 to −15 ATDC, meanwhile the SOI1 timing is kept constant at −80 ATDC. Figure 8 shows the predicted cylinder pressure and HRR for the SOI2 timing sweep for both RANS and LES models. It can be seen that, as SOI2 is advanced from 0 to −30 ATDC, the cylinder HRR and pressure are raised and combustion phasing is advanced, leading to higher in-cylinder temperatures (see Figure 9). But, as the SOI2 timing is further advanced to −60 ATDC, a contradictory behaviour is observed. With advanced SOI, the combustion is reactivity controlled; in contrast, the combustion with late SOI2 timing is close to mixing controlled. Generally, an early SOI2 timing is recommended to achieve low RI and emissions. The detailed discussion about this behaviour can be also found in Reference [14]. The interesting point is that the difference between LES and RANS models are much distinctive and obvious as the SOI2 timing is advanced. The comparison of RANS and LES models for SOI2 effect on ID, GIE, RI and CA50 are illustrated in Figure 10. Regarding the ignition delay, both models have predicted it precisely but LES is more robust and reliable in prediction of performance parameters.    Figure 11 shows the emissions trends as the SOI2 timing varies. As the SOI2 timing is delayed, CA50 is near to TDC, which increases the temperature of combustion and consequently increases NOx emissions. When ignition happens close to TDC, the lowest CO emission is achieved, since the in-cylinder temperature and pressure are high. Regarding UHC and CO emissions, both models have similar prediction with this difference that LES predicts the amount of these two pollutants more precisely than RANS. Figure 11 indicates that by using RANS model soot emission does not change with SOI2 timing sweep but LES can predict soot trend well. Generally, both LES and RANS models show the similar behaviour during SOI2 sweep but LES more exactly predicts the emissions trends.

Effect of Engine Speed
The engine speed effect on cylinder pressure and heat release rate is illustrated in Figure 12. The influence of engine speed on ID, GIE, RI and CA50 is also shown for LES and RANS models in Figure 13. It is known that, at lower engine speed, the needed time for initiation of combustion reactions is higher and combustion happens at lower CA50 and the higher temperature which results (see Figure 14) in higher NOx emissions and lower UHC and CO emissions. At higher engine speeds, the reduced available time for chemical reactions leads to more incomplete combustion and so UHC and CO soot are considerably increased (see Figure 15) [8]. This delay in combustion for LES model is more than RANS. Generally, LES and RANS model give similar predictions for heat release rate, temperature, pressure inside the cylinder and performance parameters; however, LES provides quite accurate results in prediction of engine emissions.

Conclusions
In the present research, the comparison of LES and RANS model in a natural gas/diesel RCCI combustion was studied by using a 3D-CFD modelling. Various parameters involving natural gas mass fraction, diesel SOI2 timing and engine speed were selected to explore effects of the two turbulence modelling approach on RCCI combustion prediction. The results obtained were validated by the available ones in the literature. It was found that as natural gas mass fraction increases, the SOC was delayed. Generally, a good level of agreement was observed between LES and RANS models for cylinder pressure and heat release ratios prediction for all cases but the LES gives much better prediction at higher natural gas mass fractions, while RANS model under-predicts emission levels in specific cases. Moreover, the RANS has difficulty predicting the cases containing natural gas auto-ignition. In conclusion, RANS simulations are good enough for capturing overall qualitative flow trends, however in order to get reasonably accurate estimates of velocity magnitudes, flow structures, turbulence magnitudes and its distribution, LES simulations must be used. For future research, use of open-cycle mesh with intake and exhaust valves and detailed study of flow structure considering cyclic variations may provide better predictions of in-cylinder RCCI combustion process.
Acknowledgments: The authors thank Mr. Pourya Rahnama for his support, and Convergent Science Inc. for providing a free version of their CONVERGE package for academic purposes.
Author Contributions: Amin Paykani planned the study, developed the methodology and modelling aspects and drafted the manuscript. Parvaneh Jafari performed the numerical simulations; Amir-Hasan Kakaee collaborated on editing/finalizing the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.