Predicting Cycle-to-Cycle Variation With Concurrent Cycles in a Gasoline Direct Injected Engine With Large Eddy Simulations

Cycle-to-cycle variation (CCV) in gasoline engines has undesirable effects, including noise and vibration, engine damage, and poor drivability. Moderating CCV is important to the viability of many advanced combustion concepts, such as low-temperature combustion strategies, which offer the potential of low emissions and high efficiency [1–4].

Research into CCV has been undertaken both experimentally and computationally. In computational work, computational fluid dynamics (CFD) has employed both unsteady Reynolds-averaged Navier–Stokes (uRANS) and large eddy simulation (LES) turbulence models. LESs are typically more computationally expensive than uRANS due to smaller computational cell sizes needed [5].

Simulations using uRANS have been able to qualitatively predict CCV [6,7]. Scarcelli et al. [7] simulated 21 consecutive cycles for a gasoline direct injection (GDI) engine at two operating modes. The engine and operating points were the same as used in the current study. It was shown that the cycle-to-cycle variations generated by different initial conditions at the start of every cycle did not damp out after a large number of cycles. At both stoichiometric nondilute and dilute operating modes, the uRANS predicted CCV of the two modes qualitatively. The simulations correctly predicted the higher CCV at the dilute operating point versus the nondilute. Furthermore, reduction in numerical accuracy of the uRANS simulation reduced the predicted CCV. Finally, it was shown that the uRANS result could qualitatively predict the change in combustion stability for differences in test cases.

LESs have been able to quantitatively predict CCV; however, cell sizes required for LESs are smaller than uRANS and typically the simulations are more computationally expensive. Vermorel et al. [8] acquired nine consecutive simulation cycles that showed the same level of CCV as experiment. Cell sizes of 0.2 mm–1.5 mm were used. Variations that are damped out in uRANS from the ensemble averaging, persist in LES, since large eddies, on the order of the size of the grid, are directly resolved. These eddies are influenced by the flow field behavior. Smaller eddies are modeled by subgrid models.

In several studies, where 30–100 consecutive cycles were obtained, LESs are compared favorably with experimental CCV [9–11].

Granet et al. [12] found that LES was able to distinguish stable and unstable operating conditions, where 25 cycles were required to capture the stable operating condition and 50 cycles were required to capture the unstable operating condition. Cell sizes of 0.2 mm–0.8 mm were used in the cylinder, with cells of 0.04 mm at the valve seats. Fontanesi et al. [13,14] showed that the CCV captured by LES was closer to the experimental data than uRANS, with the uRANS underpredicting the experimental CCV more. Cell sizes of 0.55 mm were used during the combustion process with 0.08 mm cells around the spark plug. Larger cells, increasing to a maximum of 0.7 mm, were used during induction and expansion. The experimental operating mode was stable combustion for a GDI engine where ten consecutive simulation cycles were acquired.

Coarse LES, which uses grid sizes similar to uRANS, has also been able to predict CCV while having the advantage of affordable computational expense, compared with uRANS [15]. In the study by Kodavasal et al. [15], the G-equation combustion model was used to model combustion with the dynamic structure LES subgrid turbulence model. The G-equation model constants were tuned to match the mean of the experimental pressure. No tuning was performed to match the CCV. For a nondilute condition, the CCV was predicted well, while the CCV for the dilute condition was underpredicted. The reason for the underprediction at the dilute condition was unknown but expected to be either from the limited size of the CFD domain or the simplified inclusion of the dilute fraction in the G-equation flame speed correlation.

The sources of CCV have been observed to primarily be bulk flow field variation occurring within the simulation domain from cycle to cycle. In some cases, particularly when there is an external exhaust gas recirculation (EGR) loop that is not included in the computational domain, CCV will be underpredicted. This will be particularly evident when there is a CCV correlation between subsequent cycles: a poor combustion cycle will feedback through the EGR loop to influence the subsequent cycle or cycles. The underresolved flow field behavior will also typically cause predicted CCV to be lower than seen in the experiment. It is expected that the quantitative CCV predicted by simulation will increase as the fidelity is increased, moving from uRANS, coarse LES, to fine LES, and ultimately to direct numerical simulation.

CCV is typically measured in statistical terms using the coefficient of variance (COV) of indicated mean effective pressure (IMEP) or peak pressure. The standard deviation can also be used to quantify CCV. The acquired number of cycles used to calculate the variation is important to accurately represent the CCV.

Typically, 25 cycles has been considered the minimum to predict the mean, while 50 cycles has been considered the minimum to predict CCV [12]. However, in an engine that is displaying moderate to high CCV (of greater than 2% COV of IMEP), more engine cycles may be necessary to predict CCV without unacceptably high uncertainty. The uncertainty in a statistical measurement of mean or variation is directly related to the measured variation in the sample size.

In experimental datasets, often 100–500 engine cycles will be acquired in order to estimate CCV accurately [9,10].

Running multiple cycles in order to predict CCV can be computationally expensive and leads to long runtimes, on the order of months for a sufficient number of engine cycles. This is unacceptably long for many industry research timelines. With the availability of large high performance computing clusters, running many cases at once has become possible, and an approach to concurrently obtain valid realizations of engine cycles has been explored.

A method to run cycles concurrently was proposed by Ameen et al. [16], where each cycle is initialized from a single 3D flow field result that is perturbed slightly. The perturbation method is as follows: The initialization is obtained from running a case for two engine cycles initially to wash out initial assumptions and then to save a point cloud map file with cell velocity, thermodynamic, and species information. The map file is perturbed for each case to be run concurrently in order to avoid the exact same solution. The concurrent cases are run for two cycles, where the second cycle of the concurrent cases is analyzed as a distinct cycle representative result. The perturbation was calculated using root-mean-square (RMS) velocity applied with an energy spectra methodology to create a physically realistic perturbation noise on the velocity field.

The use of a map file will “lose” the grid information in a case where the grid has adaptive mesh refinement (AMR), meaning that any AMR would be discarded upon the initialization of the concurrent case. This will wash out some fine scale flow structures and potentially bias the next cycle result.

Alternatively, a concurrent case could be initialized with a restart file, which is a binary file saving the detailed flow structure as well as the grid (everything necessary to restart a case). Thus, a concurrent case begun with a perturbed restart file may not be biased as much as a case started from a map file.

This parallel perturbation method was found to yield valid realizations of engine cycles, such that the statistics from the concurrent cycles yielded the same result as consecutively obtained engine cycles, when flow field snapshots were examined. The engine case used was an spark-ignited (SI) port fuel injection (PFI) engine with low CCV.

The engine used in this work is a single-cylinder GDI engine located at Argonne National Laboratories. It has been used as the basis for the previous studies by Scarcelli et al. [6,7] and Kodavasal et al. [15]. The engine specifications are listed in Table 1.

Two operating conditions were examined: a stoichiometric nondilute case with 0% EGR and a stoichiometric dilute case with 18% EGR. The IMEP was maintained the same for the two operating conditions by adjusting the intake manifold pressure. Table 2 summarizes the two operating conditions.

For each operating condition, 500 engine cycles were acquired experimentally.

The commercial CFD software converge, version 2.3.25, was used for engine simulations [17]. The code uses a modified cut-cell Cartesian mesh at run time that eliminates user meshing. The mesh is defined by user-defined settings, including a base mesh size, which is the largest cell size in the domain, and local embedding, which can define smaller cell sizes in certain regions or along boundaries, such as the spark gap or at the fuel injector nozzles.

AMR allows local areas of mesh refinement where the subgrid field is largest in variables such as velocity and temperature. The subgrid scale, defined by the user, designates when the AMR will be activated and deactivated [18].

The LES dynamic structure model was used to simulate the turbulent subgrid scale stresses. This is a one-equation, nonviscosity dynamic LES model [19]. The energy flow from the resolved scales to the subgrid scales is modeled by applying the subgrid scale kinetic energy transport equation.

The fuel spray was modeled with a conventional Euler-Lagrangian discrete droplet method. The liquid phase was modeled as discrete parcels (Lagrangian), and the gas phase was modeled as continuous (Eulerian). The “blob” injection method was used where the initial droplet diameter is initialized as the size of the effective nozzle diameter. The Kelvin–Helmholtz Rayleigh–Taylor model was employed for the spray breakup. The Rayleigh–Taylor model was used without a breakup length [20]. The Frossling model was used for modeling evaporation [21]. The no time counter parcel collision method was used [18] with post-collision outcomes [22] and the dynamic droplet drag model [23] with the O'Rourke turbulent dispersion model [21]. Wall film was modeled with the O'Rourke film splash model based on Weber number, film thickness, and viscosity [24].

Combustion was modeled using the G-equation level-set method, which is a nondetailed chemistry approach [25]. The G-equation combustion model defines the flame front as the interface between the burned and unburned gases with the scalar passive G = 0, where G indicates the distance from the flame front. Positive values of G are within the burned region, while negative values of G are in the unburned region. The flame speed is calculated based on a correlation with user-defined tuning parameters [15]. A key distinction between detailed chemistry and the G-equation combustion model is that the G-equation method does not require the flame front to be highly resolved.

A 4-mm base mesh was used in this study, with embedding in the combustion cylinder and around the valves to 0.5 mm. A minimum cell size of 0.125 mm was used at the spark plug gap. Velocity AMR was used to refine the mesh to 0.5 mm with a subgrid scale of 1.0 m/s, where this value was chosen based on being 1/100 of the expected peak quantities of flow velocities (100 m/s). AMR based on the temperature was not used, since the G-equation combustion solver was used. The peak cell count was 2.5 million cells.

The geometry of the GDI engine is shown in Fig. 1. The external EGR loop was not included in the CFD domain; for the dilute case, the EGR fraction was imposed with species fractions for the inflow boundary condition.

The simulations were run in parallel using the message passing interface. Wall clock time for each cycle was approximately 18 h on 64 cores or 51 h on 16 cores.

For the consecutive cycles, the first cycle was discarded to ensure assumed initial conditions are washed out.

Perturbations were applied to initial conditions using either a velocity perturbation or a numerical perturbation, following the methodology of Cao [26].

A numerical perturbation was generated by changing the random number seed to converge in the standard input files. This random number seed changes the random numbers used by the spray models in converge. The effect of this is very subtle on the spray injection event, but does numerically perturb the simulation. It is worth pointing out that numerical perturbations that are miniscule are sufficient to send a simulation onto different valid realizations. There are several ways to numerically perturb a CFD simulation including running the simulation on a different number of cores, changing convergence tolerances slightly, or slightly altering load balance parameters.

The start of each concurrent case was initialized from a restart file. The restart file in converge contains a complete description of the case allowing it to start with the same physical characteristics (velocity, thermodynamic, spray parcels, and turbulence field) as well as the exact grid condition (including regions of refinement from AMR). In contrast to this, a map file that is a point cloud of the initial data resets the grid to base grid and fixed embedding, so any local refinement from AMR will be lost at the start of the concurrent case. By losing the fine scale cell refinement of the AMR, velocity fluctuations will be washed out at the start when using a map file.

The restart file used to initialize the concurrent cases was taken from the consecutive simulation after at least one engine cycle had been simulated to wash out the initial boundary conditions. The restart file after the sixth engine cycle was used. In addition, restart files from this same case were generated at exhaust valve open (EVO) (−600 degrees after top dead center (dATDC)), Intake valve open (IVO) (−390 dATDC), start of injection (SOI) (−300 dATDC), intake valve close (IVC) (−135 dATDC), and spark timing (−26 dATDC) to allow concurrent cases to be launched at these times. Note the fuel injection SOI occurred during the intake stroke, and this was not close to the spark timing of either case. The timing map of the engine is shown in Fig. 2.

The uncertainty analysis used in this paper applies to the data with a normal distribution. It is well accepted that the CCV engine data are normally distributed, and the 500 cycles of the experimental data examined in this work do demonstrate a normal distribution for the peak pressure.

The dilute and nondilute cases were each simulated for 100 engine cycles consecutively (omitting the first cycle to wash out initial conditions). The cases had several restarts, between 5 and 10, due to server downtime and other unavoidable reasons. It was found by Kodavasal et al. [15] that a case with frequent restarts produces the same predicted CCV as an uninterrupted case. Thus, this was considered inconsequential.

The consecutive cases ran in 18 h/cycle on 64 cores. The 100 consecutive cycles took approximately 2.5 months to complete, not including unavoidable server downtime.

Figure 3 shows the 100 cycles for the nondilute case (0% EGR) in comparison with the 500 experimental data cycles.

The mean peak pressure of the experiment was 40.927 ± 0.311 bar (95% CI), while the COV of IMEP was 1.42% with + 0.10%/−0.09% uncertainty. The simulation results had a mean peak pressure of 40.576 ± 0.5495 bar (95% CI), with a COV of IMEP of 1.04% with +0.17%/−0.13% uncertainty. The simulation mean peak pressure matched the experimental data within the uncertainty of the sampling. The simulation COV of IMEP underpredicted the experimental variability more than the uncertainty in the sampling. This is likely due to the coarse LES resolution.

Figure 4 shows 100 consecutively obtained cycles for the dilute (18% EGR) operating condition in comparison with the 500 experimental data cycles.

The dilute case was not tuned differently than the nondilute case. All modeling settings and parameters were maintained the same between modes. The mean peak pressure of the experiment was 43.558 ± 0.311 bar (95% CI), while the COV of IMEP was 7.79% with +0.57%/−0.51% uncertainty. The simulation results had a mean peak pressure of 40.982 ± 0.668 bar (95% CI), with a COV of IMEP of 2.20% with +0.37%/−0.28% uncertainty.

The simulated COV for the dilute case was significantly lower than the experiment, beyond the uncertainty in the measurement. It is not conclusively known why the dilute variability is underpredicted more significantly than the nondilute. The coarse resolution of the LES simulations may be a factor affecting both operating points.

The simulation results were considered accurate enough for the purposes of the remainder of the paper, namely comparing the prediction of CCV of consecutively obtained simulations with concurrently obtained simulations with perturbations to the initialization.

Several situations were evaluated for running concurrent cases with a perturbed initialization. For all cases, the first combustion cycle after initialization was analyzed. For each case, 100 concurrent cycles were completed to get a statistical measurement.

The concurrent cases were run in batches of 100 on 16 cores each. The cases that began closer to combustion ran in less time, due to the shorter simulation duration. The batches of cases took 51 h for EVO, 37 h for IVO, 30 h for SOI, 14.5 h for SOI, and 12 h for spark angle. In total, approximately 1 million core-hours were used for the study.

When to start the simulation was evaluated, i.e., how close to the combustion event could the perturbed cycles be initiated, in order to minimize the required runtime, while matching the statistics of the consecutive cycles. Starting times of EVO, IVO, SOI, IVC, and spark angle were considered. The fuel injection was during the intake event, and thus, it was not close to combustion as in other direct injection engines, such as conventional diesel engines. Figure 5 shows the predicted standard deviation of peak pressure for the cases started at EVO, IVO, SOI, IVC, and spark angle. These cases were perturbed by reading a restart file and perturbing the velocity field by adding a 0.355 m/s RMS velocity field.

The consecutive case result was matched within the statistical uncertainty by the concurrent cycles when started early enough before the combustion event. At timings of EVO, IVO, or SOI, the predicted standard deviation of peak pressure was close to the standard deviation predicted by the consecutive case. When the timing was at IVC or spark, the predicted standard deviation of peak pressure was lower than the consecutive case.

For comparison with the consecutive results, 100 concurrent cycles for the different timings are shown in Figs. 6–10.

First, the 100 concurrent cycles started at EVO are shown in Fig. 6. This result is statistically similar to the consecutive cycles result shown previously in Fig. 3.

Then, Fig. 7 shows the 100 concurrent cycles started at IVO.

Figure 8 shows the 100 concurrent cycles started at SOI.

Figure 9 shows the 100 concurrent cycles for the case started at IVC. The cyclic variation is less than the cyclic variation for the cases started at earlier timings, as well as the consecutively run case.

The 100 concurrent cycles for the case started at the spark angle are shown in Fig. 10. The variability in the peak pressure for these cases is very low, presumably because the differences in the flow field that developed in the other cases had insufficient time to develop in this case.

The COV of IMEP for the cases is shown in Fig. 11. The same outcome can be seen here as was evident with the standard deviation of peak pressure: a sufficient amount of time is needed after the velocity perturbation is applied to allow the concurrent cycles to develop into distinct valid realizations.

Finally, the mean peak pressure is shown in Fig. 12. The mean peak pressure was predicted within the statistical uncertainty for all timings except the spark angle. The higher mean was also evident in the plot of all 100 cycles, as shown in Fig. 10.

It is believed that for the case started at spark, the concurrent cases did not have sufficient time to develop into a valid distinct flow field after the initialization, and thus, the 100 cycles were biased toward high combustion pressure. If different initializations were used, the cases may have been biased differently, anywhere between the min and max cycles seen in the consecutive data.

When a sufficient amount of time is available before combustion (EVO, IVO, and SOI), the concurrent cases predicted the same mean and cyclic variation as the consecutively obtained cycles.

The dilute operating point was evaluated for the same timing points: EVO, IVO, SOI, IVC, and spark timing. Figure 13 shows the predicted standard deviation of peak pressure for the same cases as shown above in Fig. 5 for the nondilute case. Similarly, the cases were started by reading a restart file and perturbing the velocity field by adding a 0.355 m/s RMS velocity field.

The dilute operating point supports the same observation as the nondilute: the perturbed concurrent cases predict the same cyclic variability as the consecutively obtained cycles when a sufficient amount of time before combustion is allowed. When the concurrent cycles were initiated close to combustion, at IVC or spark timing, the cyclic variability was underpredicted in comparison with the consecutively obtained cycles. When concurrent cycles were initiated at timings of EVO, IVO, and SOI, the same cyclic variability was predicted within the uncertainty of the statistics.

Figure 14 shows the individual 100 cycles for the EVO timing case. This figure shows similarity with the dilute consecutive cases, as shown previously in Fig. 4.

The individual cycle results for the dilute IVO case are shown in Fig. 15. They appear similar to Figs. 14 and 4 and predicted the same variability within the uncertainty.

The 100 concurrent cycles started at SOI are shown in Fig. 16.

The 100 concurrent cycles started at IVC are shown in Fig. 17, which shows less variability than the consecutively run cycles.

Figure 18 shows the 100 concurrent cycles started at spark timing. The variability in the peak pressure for these cases is very low, presumably because differences in the flow field did not have time to develop. In addition, all the cycles underpredict the mean, presumably because the restart file they began with was correlated with a low-pressure combustion cycle.

The COV of IMEP is shown in Fig. 19 for the concurrently obtained cases started at EVO, IVO, SOI, IVC, and spark angle, in comparison with the consecutively obtained cycles. The COV is predicted within the statistical uncertainty for timings of SOI and earlier. When started at IVC or spark, the COV was underpredicted.

The mean peak pressure is shown in Fig. 20. The concurrent cases matched the mean of the consecutive cycles within the statistical uncertainty.

The results of the nondilute and dilute engine operating conditions supported the same conclusions on using concurrently obtained engine cycles to predict cyclic variation. When sufficient time is allowed after initialization, the flow field generates distinct and valid realizations of individual cycles. In this type of GDI engine, including the cylinder filling intake event is crucial to the flow field development for simulating distinct realizations of combustion cycles and predicting CCV.

By simulating CCV with perturbed concurrent simulations, researchers can shorten the calendar time required for studies. Simulating 100 cycles consecutively required 2–3 months to complete, while 100 concurrent cycles could be completed in 1–2 days. In addition, the concurrent cycles were restarted part-way through the cycle, so the required core-hours for each cycle was less. On 16 cores, a full engine cycle took 51 h, while a cycle simulation starting at SOI took 30 h. This reduces the total required core-hours by up to 41%.

How the perturbation should be applied was evaluated for several select cases and methods. These comparisons were conducted for the nondilute operating point.

A numerical perturbation was compared with a velocity perturbation at SOI. The numerical perturbation was a random number seed for the random number generator used by the spray models. The velocity perturbation for this comparison was a 0.355 RMS velocity, the same as used in section Concurrent Cases. Both types of perturbations were applied after reading the same restart file.

Figure 21 shows the standard deviation of peak pressure for the concurrent cases initiated at SOI compared with the consecutively obtained cycles. The predicted cyclic variability is not different when the concurrent cycles are perturbed with a numerical or velocity method. The results indicate that either type of perturbation will predict the same variability within the statistical uncertainty of the cases.

Figure 22 shows the 100 concurrent cycles perturbed with the random seed. The comparable plot for the 0.355 m/s RMS velocity perturbation at SOI was shown previously in Fig. 8.

A map file (a point cloud initialization without saving the grid) was compared with a restart file (grid retained) with an identical velocity perturbation at IVO and SOI. A 0.355 m/s RMS velocity perturbation was applied at IVO, while a numerical perturbation was applied at SOI via a change of the random number seed.

Figure 23 shows the standard deviation of peak pressure for the IVO and SOI cases, with the consecutively obtained cycles as a comparison. The map file case at IVO overpredicts the variability slightly, beyond the sampling uncertainty.

The 100 concurrent cycles for the IVO map file case are shown in Fig. 24. The variability appears larger, particularly due to several low cycles, when compared with the consecutively obtained cycles in Fig. 3 and the concurrent restart case in Fig. 7.

It is inconclusive if a map initialization can be used in place of a restart file when estimating cyclic variability from concurrently perturbed cycles.

A large velocity RMS perturbation (0.355 m/s) was compared with a small velocity RMS perturbation (0.003 m/s) at EVO and IVO. The size of the RMS velocity perturbation was set equal to 0.335 (large) or 0.003 m/s (small). To assess the magnitude of these velocity values, images of the velocity field were evaluated at IVO.

Figure 25 shows the velocity magnitude on a cutplane through the intake and exhaust runners at the IVO crank angle.

The velocity in the exhaust port is high, while the velocity in the intake port is low. Figure 26 shows the intake port velocities in more detail. Near the valve, the intake port has velocity magnitudes on the order of 1 m/s. Since the magnitudes of the mean velocity field are quite low near the entry to the intake valve, the 0.335 m/s perturbation is considered to be large. Reducing this by a factor of 100–0.003 m/s was considered to be small.

After the large velocity perturbation (0.355 m/s) was applied, a modification of the velocity field can be observed in Fig. 27. This cutplane was generated from the start of a concurrent case initiated at IVO.

A similar cutplane was generated with the small velocity perturbation (0.003 m/s), as shown in Fig. 28. Due to the small velocity perturbation, there was no observable difference from the nonperturbed velocity field.

The 100 concurrent cycles are shown in Fig. 29 compared with the consecutively obtained cycles for the large and small velocity perturbations at EVO and IVO.

The size of the velocity perturbation did not affect the prediction of variability within the statistical uncertainty when initiated at either EVO or IVO.

It was found that concurrent cycles can be used to evaluate CCV in coarse LES simulations for a GDI engine when a sufficient amount of time is allowed before the combustion event after a perturbation is applied. This was shown by acquiring 100 consecutively obtained cycles as a benchmark and comparing 100 concurrently obtained cycles for several different cases. Two operating points supported this finding: a nondilute and dilute.

The predicted CCV was compared with both the standard deviation of peak pressure and the COV of IMEP. The statistical analysis compared the results by calculating the uncertainty with a 95% confidence interval.

Several different ways to concurrently start simulations were evaluated. It was found that

• Concurrently obtained cycles could predict the CCV when started early enough before the combustion event. It was found that starting the concurrent simulations at timings of EVO, IVO, and SOI were early enough, while starting timings of IVC or Spark were not. This was demonstrated at both dilute and nondilute operating points using a velocity perturbation applied to a restart file.

• The type of perturbation was not found to be important, as long as there was time for the flow field to develop into a valid distinct realization. A large velocity perturbation (0.355 m/s) and a small velocity perturbation (0.003 m/s) were compared at two start times, IVO and SOI. The size of the velocity perturbation did not produce statistically different results. A numerical perturbation (changing the random number seed) was compared with a velocity perturbation at SOI and was also found to produce the same variability within the statistical uncertainty.

• The method of initialization was compared between a restart file and a map file. These two methods are the most common approaches to initialize concurrent cycles. A restart file restores a simulation with all variable values and modeling parameters, including grid refinement due to AMR, while a map file simply initializes the variables, such as temperature, velocity, and species, but does not save the grid refinement or other simulation characteristics. The results were inconclusive as to whether a map file was suitable to start concurrent simulations as close to the combustion event as a restart file, since a case initiated at IVO overpredicted the variability in comparison with the restart case.

By concurrently running cases, 100 engine cycles can be completed in significantly less wall clock time than consecutive simulations. Simulating 100 concurrent cycles took less than 2 days, while the consecutive simulations took 2–3 months. This increase in turnaround is valuable for researchers as it allows CCV to be studied more rapidly. In addition, the required core-hours is reduced since only a portion of the engine cycle needs to be solved with the concurrent method (for instance, starting at SOI). This can yield up to a 41% reduction in computational expense.

UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”), a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-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. This research was partially funded by DOE's Office of Vehicle Technologies, Office of Energy Efficiency and Renewable Energy (Contract No. DE-AC02-06CH11357). The authors wish to thank Gurpreet Singh and Michael Weismiller, program manager at DOE, for their support. We gratefully acknowledge the computing resources provided on the Bebop and Blues clusters, which are operated by the Laboratory Computing Resource Center at Argonne National Laboratory.