Numerical Prediction of Cyclic Variability in a Spark Ignition Engine Using a Parallel Large Eddy Simulation Approach

Cycle-to-cycle variability (CCV) is a phenomenon, primarily in spark ignition (SI) engines, in which significant variations in the in-cylinder pressure are observed from cycle to cycle. CCV is important for two major reasons. First, the optimum engine settings are designed for the average cycle, and hence the faster and slower engine cycles will lead to losses in power and efficiency. Second, when CCV is extremely high, this can lead to partial burn, misfire, or engine knock and an overall reduction in the reliability of the engine. With the advent of advanced modes of combustion including low temperature combustion and dilute operation, the effect of CCV is observed to be higher.

Prior studies [1] have shown that CCV is caused by a combination of factors including variations in in-cylinder flow pattern, mixture inhomogeneity, turbulence intensity and spark discharge characteristics. The relative importance of these factors will depend on the engine geometry and operating conditions. There have been several experimental studies, which investigated the causes of CCV. Recently, many experimental studies have employed particle image velocimetry (PIV) in optically accessible engines to characterize the flow structures present inside the combustion chamber and their cyclic variabilities. Reuss [2] utilized two-dimensional (2D) PIV to measure the cyclic variability in the flow structures at the top dead center (TDC) of a motored transparent combustion chamber (TCC) engine. It was shown that depending on the valve design, the cyclic variability in the flow structures could be minimal or extremely large. Funk et al. [3] extended this study and showed that at low swirl ratios, the flow is dominated by turbulence at the large scales, while the turbulent Reynolds numbers are much lower at high swirl ratios. The TCC engine geometry has undergone several revisions and has resulted in the TCC-III engine benchmark [4] with PIV measurements available at several crank angles, measurement planes, and engine speeds, and has been used by several modelers to validate their simulations.

Several experimental measurements have been made in other optical engine set-ups as well. Baum et al. [5] measured the three-dimensional (3D) flow field in a motored optical engine using a tomographic PIV technique. These measurements provided insights into the formation and orientation of vertical structures inside the engine. Voisine et al. [6] performed high-speed PIV and analyzed the spatial structure of the flow and its temporal evolution during series of consecutive cycles in a research engine of moderate tumble ratio. They showed that approximately 30% of the fluctuating kinetic energy was generated due to cyclic fluctuations in the chamber near TDC. There have been a few other experimental studies [7,8], which focused on the large-scale structures and their cyclic variability in SI engines. However, in experimental studies, it is difficult to isolate the different coupled processes and hence cannot be used to explain the causes of CCV.

Numerical prediction of CCV is extremely challenging for two key reasons: (i) high-fidelity methods such as large eddy simulation (LES) are required to accurately capture the in-cylinder turbulent flow field, and (ii) CCV is experienced over long timescales and hence the simulations need to be performed for hundreds of consecutive cycles. In spite of these challenges, there have been several numerical studies in recent times where multicycle LES has been used to predict CCV in engine-like geometries and understand the causes for CCV. These include simplified geometries resembling piston and/or valve motions [9–13], realistic engine geometries under motored [14,15] and fired conditions [16–18]. Moureau et al. [9] performed 2D and 3D LES of multiple cycles of a simplified square piston engine and compared the simulation results with experimental PIV measurements. It was shown that even though the mean flow is essentially two-dimensional during the intake phase, only the 3D simulation was able to predict the breakdown of the tumble during the compression phase. Another simplified engine geometry that has been widely studied is the “Imperial College” engine, which consists of a single, centrally located valve, a circular cylinder and a flat piston. Direct numerical simulation (DNS) [11] and LES [12,13] studies have been performed for this engine geometry. These studies showed that the flow-field variability is dominated by the dynamics of the intake jet and the vortex ring it creates. Yang et al. [14] compared the predictions of cycle-averaged velocity and turbulence predictions by Reynolds-averaged Navier–Stokes (RANS) and LES for the motored TCC engine. They showed that 2–3 cycles of RANS were able to capture the overall qualitative flow trends; however, LES was required to obtain accurate estimates of velocity magnitudes, flow structures, turbulence magnitudes, and its distribution.

Recently, improved computational capabilities and improved solvers have enabled the use of multicycle LES to capture CCV in fired engines. Schmitt et al. [17] performed multicycle LES of the flow-field evolution and combustion process in a direct injected methane engine. It was shown that the simulations reproduced the average pressure trace and the cyclic fluctuations from the experiment. They also showed that the global tumble ratio strongly correlated with the turbulent kinetic energy and this variability was a major cause for the CCV.

These numerical studies have shown the importance of the in-cylinder flow field and its cyclic variability on CCV. However, the numerical simulation of CCV using multicycle LES is still extremely expensive and performing these calculations on computing clusters would require a few months of computational time and thus would not be realistically used in the engine design process. To solve this problem, a few different ideas have been proposed to reduce the turnaround time by decomposing the long simulations into several shorter simulations.

Goryntsev [19] introduced a parallelization technique based on variation of the boundary conditions in each of the individual cycles. The perturbations were introduced at the inflow boundary during the intake stroke. Further work by Goryntsev et al. [20,21] has applied this technique in simulating the CCV in realistic direct-injection spark ignition engines. By applying an initial perturbation to the inflow velocity profile, ten consecutive engine cycles were simulated on five individual simulations resulting altogether in 50 individual engine cycles. It was shown that this technique resulted in an estimated statistical error for mean velocity and its standard deviation of 5% and 10%, respectively. It was also shown that this technique led to a 5× reduction in the computational expense.

Finney et al. [22] proposed the idea of using concurrent single-cycle simulations to replicate the long timescale dynamics of engine combustion by employing a simple one-dimensional model. This approach was later extended to replace the one-dimensional model with concurrent RANS simulations [23,24] by perturbing global parameters like exhaust pressure, exhaust temperature, and spray mass. They showed that perturbing these variables led to variations in the cylinder pressure. However, the range of variation observed from this approach has not been validated against experimental data. It is also unclear whether this approach of branching out the simulations leads to similar CCV as that from the traditional approach of simulating consecutive cycles.

Ameen et al. [25] recently developed the parallel perturbation model (PPM), a methodology to reduce the turnaround time in simulating CCV without sacrificing the accuracy. The strategy is to perform multiple parallel simulations, each of which encompasses 2–3 cycles, by effectively perturbing the simulation parameters such as the initial and boundary conditions. The methodology is schematically shown in Fig. 1. More details of this methodology are explained in Ameen et al. [25]. This methodology was validated for the prediction of CCV due to gas exchange in a motored TCC engine by comparing the flow-field statistics with PIV measurements. It was shown that by perturbing the initial velocity field effectively based on the intensity of the in-cylinder turbulence, the mean and variance of the in-cylinder flow field were captured reasonably well.

Recently, Zhao et al. [18] performed LES of 50 consecutive cycles for a fired port fuel injected (PFI) SI engine. They employed the G-equation approach to model the flame propagation and showed that the simulations were able to accurately predict the cyclic variability from the experiments. They also showed that, by decoupling the effects of the velocity field and the equivalence ratio field, the velocity field and not the equivalence ratio field is what dominated the CCV for this engine configuration. This motivates the use of PPM to simulate the CCV for this engine configuration, since PPM was shown to be effective in capturing the cyclic variability of the in-cylinder flow field [25].

In this paper, the PPM approach is extended to simulate the CCV in a fired PFI SI engine previously studied by Zhao et al. [18]. Two operating conditions are considered—a medium CCV operating case corresponding to 2500 rpm and 16 bar brake mean effective pressure (BMEP), which was studied by Zhao et al. [18] and a low CCV case corresponding to 4000 rpm and 12 bar BMEP. In Experimental and Numerical Setup section, the details about the experimental and numerical setups employed will be discussed briefly. This will be followed by a comparison of the simulated results with the experiments as well as with the consecutive cycle LES results from Zhao et al. [18]. The goal is to demonstrate that the PPM can predict CCV for different operating conditions and then document the reduction in computational time due to this approach. The successful demonstration of this approach can result in further use of simulations in the engine design process and can result in significant cost-saving for engine manufacturers.

The experiments were performed on a four-cylinder SI PFI engine. Two intake and two exhaust valves were present per cylinder. Table 1 summarizes the major details about the engine, and Table 2 describes the details of the two operating conditions under which the engine was run. The fuel employed in the experiments was the unleaded 95 RON (research octane number) European gasoline. Crank-angle resolved in-cylinder pressure measurements were made for 1000 consecutive engine cycles for one of the four cylinders. Further details about the experimental setup and the measurement techniques are described in Mirzaeian et al. [26]. Both cases A and B operate at stoichiometric conditions and correspond to stable operating points. Case A corresponds to a late spark timing, medium CCV condition, and case B corresponds to a low CCV condition.

The simulations were performed using the commercial CFD code CONVERGE v2.3 [27], which employs a Cartesian cut-cell approach to mesh generation. The modeling was performed for the one cylinder only, and Fig. 2 shows the computational domain and identifies the location of the injector, ports, valves, and spark plug. Figure 2 also shows the location of the Y = 0 tumble plane, which will be used for further analysis in Figs. 3–7. Table 3 summarizes the important details of the numerical setup employed in the present study. A second-order central difference scheme was used for spatial discretization and first-order implicit scheme for temporal discretization. A maximum grid size of 0.7 mm was employed inside the cylinder. Boundary embedding was added to refine the grid to 0.35 mm in the intake and exhaust valve openings to accurately model the flow near the valves. Additional refinement of 0.175 mm was enforced near the spark plug region. Further, adaptive mesh refinement (AMR) was activated based on velocity and temperature gradients to generate cell sizes of 0.35 mm in regions having AMR.

In this study, the liquid spray (from the PFI injector in the intake port) is modeled as a discrete phase in the Lagrangian framework. The details about the spray submodels are summarized in Table 3 and discussed by Senecal et al. [28]. The flame kernel growth and flame propagation were modeled using the G-equation combustion model [29,30], which is a level set-based model, which ignores detailed chemical kinetics. The position of the flame front is at G = 0, and the burned region has G > 0, where G is the level-set variable. Ignition is modeled by adding a volumetric source for G in a sphere of radius 0.45 mm around the spark plug. The turbulent flame speed closure employed and other details of the G-equation model employed in this study are explained in Zhao et al. [18].

In the present study, two different operating conditions are employed, the details of which are summarized in Table 2. For both these conditions, 100 parallel cycles were simulated using the PPM approach as shown in Fig. 1. 50 consecutive LES cycles were simulated for case A by Zhao et al. [18], and these results will also be used for further validation of the PPM approach.

Finney et al. [31] provided a comprehensive review of the recent developments in the understanding of CCV in SI engines. They showed that depending on the operating conditions, CCV can exhibit either stochastic or deterministic behavior. When the engine is operating at a relatively stable operating point with low cyclic variability, the features of CCV are inherently random with no short-term predictability. Under these stochastic conditions, there are minimal correlations between individual engine cycles, and CCV can be modeled as a stochastic process by applying appropriate perturbations to the model parameters. However, when the engine reaches limits of stability, CCV can exhibit a deterministic behavior with appearance of low-dimensional deterministic structures among the individual engine cycles. To determine the stochastic and deterministic aspects of CCV, Kaul et al. [32] determined the relationship between the indicated works of consecutive cycles. They showed that when the indicated works of consecutive cycles exhibited minimal correlation with each other, the engine operation can be expected to exhibit stochastic CCV features. When the engine operation reaches limits of stability at higher exhaust gas recirculation levels, the consecutive cycles exhibit strong correlation with each other (for instance, a low cycle could be followed by a high cycle), and the engine CCV is expected to be deterministic in nature. Figure 8 compares the indicated work of consecutive cycles for both cases A and B. It can be seen that there are no correlations between consecutive cycles and thus both these conditions exhibit purely stochastic features of CCV. Thus, the extension of the PPM approach, which models CCV as stochastic process, to the fired PFI engine conditions described in Table 2 is justified.

Ameen et al. [25] explored different strategies to parallelize the multicycle LES by applying perturbations to initial flow field, initial temperature, initial pressure, and boundary pressures. They showed that applying perturbations to the initial flow field by superimposing multiple instances of a synthetic flow field was sufficient to obtain similar flow field statistics as that of the consecutive cycle LES. Also, Zhao et al. [18] showed that the cyclic variability in the flow field was the primary contributor to the CCV under the operating conditions described in Table 2. Hence, in this study, the PPM LES approach was applied to the conditions described in Table 2 by applying perturbations to the flow field. No perturbations were applied for other initial or boundary conditions. If a more advanced combustion model is employed, it is possible that neglecting equivalence ratio variations may introduce errors. In such a case, PPM approach should be enhanced by introducing perturbations in the equivalence ratio fields as well. Initially, LES is performed from the TDC of compression for three full cycles to remove the effect of the initial conditions. Multiple parallel simulations are then performed by starting from the flow field at the end of cycle three by superimposing synthetic flow fields as shown in Fig. 1. The synthetic flow fields were generated by assuming a Batchelor spectrum with prescribed length and velocity scales. Our rule of thumb is to use the mean clearance height at top dead center as the length scale and the mean piston speed as the velocity scale. These estimates were based on the observations of Heywood [1] that the length and velocity scales of turbulence at TDC were approximately equal to the mean clearance height and the mean piston speed, respectively. The length and velocity scales of this synthetic turbulence were chosen to be 11.8 mm and 7.0 m/s (11.2 m/s for case B), respectively. 100 cycles were simulated using PPM LES for both cases A and B. The simulation time was approximately 3 days on 96 cores for each parallel cycle. The turnaround time for simulating 100 PPM LES cycles is limited only by the availability of computing resources, and can be as short as 3 days if sufficient resources are available. By comparison, the 50 consecutive LES cycles took more than 3 months of simulation time on 96 cores. The turnaround time for the parallel cycles is about ten times smaller than the consecutive cycle LES. In the rest of this paper, these results will be compared with experimentally measured pressure traces as well as the consecutive cycle LES presented by Zhao et al. [18].

The primary objective of applying the PPM LES approach is to drastically cut down the turnaround time to numerically predict CCV without sacrificing the computational accuracy. The validity of the PPM LES for fired engines is first verified by comparing the simulation results from PPM LES with those from the consecutive LES cycles for case A.

Figure 9 compares the trapped mass from 50 consecutive LES cycles with those from 50 PPM LES cycles for case A. It can be seen that both the approaches predict less than 0.3% cyclic variability in the trapped mass. The mean trapped mass predicted by both the approaches are very similar to each other, i.e., the consecutive LES predicted 590.28 mg and PPM LES predicted 589.28 mg. The standard deviation of trapped mass predicted by the consecutive LES was 0.537 mg, whereas it was 0.387 mg by the PPM LES. Although the PPM LES predicts a slightly less cyclic variability in trapped mass than the consecutive LES, these differences are expected to play a minimal role in predicting the CCV of combustion.

Ameen et al. [25] showed that the PPM LES approach was able to accurately replicate the flow-field statistics predicted by the consecutive LES cycles. The validity of this conclusion for fired engines is verified in Figs. 3 and 4 for case A. Although the objective of the PPM approach is to capture the cyclic variability of the flow field and combustion, it is important to ensure that the mean fields are predicted correctly as well. Since the perturbations are added to the flow field at the end of cycle 3, and not the mean flow field, there could be differences in the mean fields predicted by the PPM and the consecutive cycles approach. Figure 3 compares the mean velocity magnitudes and directions for consecutive and PPM LES cycles at 1 deg before spark timing along the Y = 0 plane that passes through the spark plug. The location of this plane with respect to the engine configuration is shown in Fig. 2. 50 cycles were used to compute the statistics for both the consecutive LES and PPM LES. It can be seen that there is an excellent agreement in the large-scale flow structures predicted by the two approaches. The simulations show the presence of a dominant tumble flow pattern at this crank angle position. Figure 4 compares the RMS velocity magnitudes predicted by the two approaches for the same measurement plane. The peak values of the RMS magnitudes predicted by both the approaches are seen to be very similar to each other. Both the approaches also predict peak RMS velocities close to the tumble center. Although there are minor differences in the RMS velocity distributions, these differences are not expected to appreciably affect the CCV of combustion.

Figures 5 and 6 compare the statistics of the equivalence ratio distributions for case A at the same measurement plane and crank angle location as those considered in Figs. 3 and 4. Figure 5 shows that both the consecutive and PPM LES approaches predict very similar distribution of equivalence ratio at this measurement plane. The simulations also show stratification in the mean equivalence ratio distribution with a slightly richer mixture near the exhaust valve side (left side of the measurement plane). Figure 6 compares the corresponding distributions of the RMS of the equivalence ratio. Both the approaches predict that the cyclic variability in the equivalence ratio distribution is less than 20%. It is also observed that there are noticeable differences in the distributions of the RMS of equivalence ratio between the two approaches. However, Zhao et al. [18] had shown that the cyclic variability of the composition field had only a minor influence in the cyclic variability of combustion rates. Hence, the differences in the standard deviation of the equivalence ratio field between the two approaches are not expected to affect the CCV of combustion significantly.

where H_y is the angular momentum of the flow about the y-axis, M_y is the moment of inertia about the y-axis, and ω is the engine speed in revolutions per second. Figure 7 compares the mean and standard deviations of the Y-tumble ratio between the two approaches. It is seen that the PPM LES approach is able to accurately mimic the variation of the mean tumble ratio for all crank angle locations. Both the simulation approaches are able to predict the tumble breakdown process. The magnitude of the tumble ratio initially increases and is followed by a strong reduction during the compression stroke (−180 deg to 0 deg). In a numerical study conducted in similar engine geometry, Vermorel et al. [16] had shown that the magnitude of this tumble breakdown is correlated to the combustion phasing, i.e., cycles with large tumble breakdown had shorter combustion durations. Figure 7 shows that the PPM LES is able to accurately predict the qualitative as well as quantitative variabilities in the global flow patterns.

Figure 10 compares the mean pressures from the consecutive and PPM LES cycles with the experimentally measured pressure traces. It can be seen that both the simulation approaches predict similar spread in the pressures as the experimental measurements. The intensity of CCV is typically characterized in terms of the coefficient of variation (COV). COV of any measured or simulated quantity is defined as the ratio of its standard deviation to its mean. Figure 11 compares the COV of maximum pressure, indicated mean effective pressure (IMEP), and three combustion phasing quantities, i.e., CAign-10, CAign-50, and CA10-75. Here, CAign-10 and CAign-50 are the crank angle durations from ignition timing for 10% of the fuel mass to burn and 50% of the fuel mass to burn, respectively. CA10-75 is the crank angle duration from the 10% fuel mass burn time to 75% fuel mass burn time. These are the same definitions employed by Zhao et al. [18]. A few important observations can be made from Fig. 11. First, both consecutive and PPM LES cycles overpredict the COV in almost all the quantities as compared to the experimental measurements. This could imply that the simulations are adding certain unphysical variabilities into the CCV. One possible source of additional variability could be due to the use of numerical restarts in the simulations. The computing clusters used for the current study enforce certain limits in the maximum wall-clock time for any simulation leading to the necessity to restart the simulations multiple times from checkpoints. These restarted simulations can introduce additional disturbances to the simulations (due to round-off errors in the restart files for instance) that could grow and manifest themselves as increased COV. Further studies are necessary to quantify the effect of restarts on CCV. The second major observation from Fig. 11 is that the COV of early stage of combustion, i.e., CAign-10, is underpredicted by both the LES approaches as compared to the measurements. This could be an artifact of using the G-equation turbulent combustion model. The G-equation model initializes the spark kernel as a spherical volume of fixed radius at the same location for all the cycles. This possibly enforces an unphysical constraint to the early stage flame development, thus leading to lower variability in CAign-10. The validity of this hypothesis can be verified by applying more accurate spark energy deposition and turbulent combustion models. Finally, it is also observed that the PPM LES consistently underpredicts the COV as compared to the consecutive LES cycles by 5–10%. This is also along expected lines as the PPM LES presented in this study only included the perturbations in the initial flow field. Overall, over- or underprediction of CCV is within ±20% from experiments, which we believe is reasonable for simulations.

It was shown in the Parallel Perturbation Model Large Eddy Simulation of Case A section that the PPM LES approach was able to mimic the flow field statistics as well as the CCV of combustion as predicted by the consecutive LES cycles reasonably well. Although the quality of any numerical model is determined by how close it is able to quantitatively match the experimentally measured quantities, it is equally important for the model to be able to accurately predict the qualitative trends in these quantities with changing operating conditions without the need to tune any model parameters. To determine whether the PPM LES is able to predict the effect of operating conditions on CCV, 100 PPM LES cycles were also performed for case B (refer to Table 2). Case B corresponds to a low-CCV operating condition at a higher RPM and earlier spark timing. Figure 12 compares the mean pressure predicted by PPM LES with the measured pressure traces. Comparing Figs. 10 and 12, it can be seen that the PPM LES is able to predict the reduction in CCV with changing operating conditions reasonably well. However, the discrepancy between the experimental and simulated pressures is slightly larger for case B. This could be because of the fact that the G-equation model employed in this study was optimized for case A [18]. Also, the temperature at the intake boundary was not measured for case B, and hence had to be approximated based on the measurements for case A. In spite of these uncertainties, the simulation approach is able to correctly predict the trend in CCV with changing operating conditions.

To quantitatively evaluate the effect of operating conditions on the simulation results, Fig. 13 compares the COV of the maximum pressure between cases A and B for both the measured and the PPM LES cycles versus number of cycles. It can be observed that 200 experimental cycles are required to attain a converged value for COV of the maximum pressure. However, the COV of maximum pressure predicted by PPM LES attains convergence at around 45 cycles. Similar number of cycles was sufficient for the COV of the burn rate related quantities as well. From the experimental pressure traces, it is observed that the COV reduces from 7.64% for case A to 4.14% for case B. The PPM LES predicts a reduction in COV of maximum pressure from 9.17% for case A to 5.99% for case B.

Figure 14 compares the COV of all pressure and burn rate related quantities between the experimental and PPM LES cycles. It can be seen that generally the conclusions obtained for case A in Fig. 11 holds for case B as well. The largest discrepancy is observed for the COV of the early stage of combustion. This could again imply the deficiencies of the G-equation combustion model, and the nonoptimal model parameters for this operating condition. Use of more sophisticated spark energy deposition and turbulent combustion models is again expected to improve the simulations. In spite of all these limitations, it is observed that the PPM LES is able to reasonably predict the measured COV as well as the trends in COV with changing operating conditions. We believe that this is the first demonstration that a parallel perturbation approach can capture trends in CCV accurately with changing operating conditions.

Ameen et al. [25] applied PPM LES for a motored TCC engine. In their study, they examined the convergence of the flow-field statistics with respect to the number of cycles performed for each condition with the PPM LES calculations. They had shown that simulating two consecutive cycles for each of the PPM LES calculations was sufficient to obtain convergence of flow-field statistics. In this section, the validity of this observation for fired engines is examined.

For all the results shown in the previous sections, two consecutive cycles were simulated for each of the PPM LES simulations and the COV was computed based on the second cycle. To examine the convergence of the PPM LES statistics, 50 additional PPM LES calculations were performed for case B in which three consecutive cycles were simulated for each simulation. Figure 15 compares the COV of P_max and CA10-75 predicted from cycles 2 and 3 of 50 PPM LES calculations. It can be seen that although there are differences in the initial transients of COV with increasing number of cycles, both the cycles predicts similar values of COV at 20 cycles. Thus, it can be concluded that two consecutive cycles are sufficient to attain convergence for the PPM LES approach under fired conditions as well.

In this work, the PPM approach introduced by Ameen et al. [25] for speeding up CCV simulations was extended to a fired PFI engine. The idea is to launch multiple parallel simulations simultaneously by effectively perturbing the initial flow field of each of these simulations. It is shown that the PPM approach is able to accurately predict the cyclic variability in the in-cylinder flow field, composition field, pressure, and burn rates. We demonstrated that the parallel LES approach is able to correctly predict the trends in the CCV with changing operating conditions. It is also shown that the results from the PPM LES match the results from consecutive cycle approach reasonably well. The turnaround time for the parallel LES approach is about ten times shorter than the consecutive cycle approach. This implies that these CCV calculations can be completed in a few days and can thus be realistically employed in the engine design process.

Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under contract No. DEAC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The authors wish to thank Gurpreet Singh, Kevin Stork, and Leo Breton, program managers at DOE, for their support. This research was conducted as part of the Co-Optimization of Fuels & Engines (Co-Optima) project sponsored by the U.S. Department of Energy (DOE) Office of Energy Efficiency and Renewable Energy (EERE), Bioenergy Technologies and Vehicle Technologies Offices. Co-Optima is a collaborative project of multiple National Laboratories initiated to simultaneously accelerate the introduction of affordable, scalable, and sustainable biofuels and high-efficiency, low-emission vehicle engines.

We would also like to acknowledge Ms. Emma Zhao from Michigan Technological University and Dr. Janardhan Kodavasal and Dr. Abdul Moiz from Argonne National Laboratory for help with the simulation setup.

The computing resources were provided by the Laboratory Computing Resource Center at Argonne National Laboratory. The authors would also like to acknowledge Convergent Sciences Inc. for providing CONVERGE licenses which were used for this work.

This research was partially funded by DOEs Office of Vehicle Technologies, Office of Energy Efficiency and Renewable Energy under Contract No. DE-AC02-06CH11357.