© The Author(s), published by EDP Sciences, 2022

Nomenclature

AHRR: Apparent Heat Release Rate

AKTIM: Arc and Kernel Tracking Ignition Model

COV: Coefficient Of Variance

CCV: Cycle-to-Cycle Variability

CFD: Computational Fluid Dynamics

: Laminar flame thickness

Δ: Grid or filter size

DI: Direct Injection

DPIK: Discrete Particle Ignition Kernel

E: Wrinkling factor (or efficiency function)

ECFM: Extended Coherent Flame Model

EGR: Exhaust Gas Re-circulation

F: Thickening factor

F_max: Maximum thickening factor

HRR: Heat Release Rate

ISSIM: Imposed Stretch Spark-Ignition Model

LDV: Light-Duty Vehicle

LES: Large-Eddy Simulation

LESI: Lagrangian–Eulerian Spark-Ignition

n_res: Number of grid points across the flame front

N.Cyc.: Total Number of Cycles

PFI: Port Fuel Injection

PIV: Particle Image Velocimetry

RANS: Reynolds-Averaged Navier–Stokes

SI: Spark-Ignition

SparkCIMM: Spark Channel Ignition Monitoring Model

S_l: Laminar Flame Speed

S_t: Turbulent Flame Speed

TFM: Thickened Flame Model

UDF: User-Defined Function

WSR: Well-Stirred Reactor

1 Introduction

The transportation sector in the United States is the largest producer of CO₂ emissions, followed closely by the electrical sector [1]. Within the transportation sector, Light-Duty Vehicles (LDVs) consume the largest share of energy, which mostly originates from gasoline, as shown in Figure 1 [1]. It naturally follows that LDVs are responsible for the most CO₂ gas emissions within the sector [1]. Most LDVs are powered by four-stroke gasoline Spark-Ignition (SI) Internal Combustion Engines (ICEs) due to a combination of low production and maintenance costs and availability of gasoline distribution infrastructure. As a result, efficiency improvements and CO₂ emission reductions in SI engines can potentially lead to national and global scale emission reductions, which are favorable for global warming and climate change.

Fig. 1 History and projections of the United States transportation sector energy consumption: (a) by transportation mode and (b) by fuel source [1].

On the other hand, electrification of LDVs by Original Equipment Manufacturers (OEMs) and recent government-driven policy changes cast uncertainty on the future of SI engines in LDVs. Despite predictions by the U.S. Energy Information Administration (EIA) that SI engines will remain a large part of the market by 2050 [1], recent changes predict an unfavorable future for SI engines. However, even in the most unfavorable scenarios, spark ignition systems can still have applications in hybrid vehicles, heavy duty, and even off-road engine applications such as rail, marine, and aviation.

In the meantime, to reduce fleet CO₂ emissions, engine automakers have turned to advanced and unconventional engine operation [2]. Highly dilute, stratified charge, and boosted operation, among others, have the potential to increase SI engine efficiency [2]. At these operating conditions, the engine experiences high levels of dilution through Exhaust Gas Re-circulation (EGR) or overall lean mixture, which leads to higher efficiency at the cost of elevated cyclic variability and combustion instability. To compensate, ignition systems must provide higher energy for a longer duration compared to traditional engine operation to ensure stable combustion. Reliable and accurate spark ignition models are necessary to design ignition systems that reduce cyclic variability [3], making ignition a relevant research item for the foreseeable future.

Computational spark ignition models have become synonymous with spark channel modeling and early flame formation. In addition, while flame transition from laminar kernel to fully developed turbulent propagation is in principle an ignition problem, it is currently being handled by combustion models [4, 5]. In general, ignition begins with plasma energy deposition, followed by laminar kernel formation, and finally laminar-to-turbulent kernel to flame transition. In modern engine simulations, there is overlap between the ignition model and combustion (or turbulent flame propagation) model, and thus the coupling between the two is important.

In general, ignition models are classified as either Eulerian [6–8] or Lagrangian [9–11]. In Eulerian models, usually (but not exclusively), a heat or energy source is deposited on a refined finite volume Eulerian grid that is equivalent to the spark (or plasma) energy delivered to the gas mixture [12]. Then, pressure and temperature are allowed to expand and convect with the flow field, making use of the grid resolved velocity. Initially, Eulerian approaches were limited to simple electrode geometries with one energy source [6]. As computational power advanced, Eulerian models became more viable with smaller grid sizes and can now handle full engine geometries with multiple sources corresponding to the different stages of ignition [8]. The energy source can have any pre-determined shape represented by a group of computational cells: a sphere (industry standard), cylinder, and even a line or cube in some Computational Fluid Dynamics (CFD) solvers. The source can be stationary or allowed to move with the flow and can account for gas ionization. Energy deposition is a Eulerian process even though the ignition model can have Lagrangian aspects for arc tracking [13, 14].

Lagrangian models, on the other hand, were the more common type of ignition models 5–20 years ago. The general approach is to track the spark channel or arc through a series of points or kernels in space and time, using flame speed correlations. In addition, transport properties are solved to keep track of the size of the kernel. When the flame kernel reaches a critical radius or volume, a Eulerian combustion model, usually a flamelet model, handles flame propagation. In general, Lagrangian models are designed for a specific turbulent combustion model [4, 15]. Most well-known Lagrangian ignition models are Discrete Particle Ignition Kernel (DPIK) [16], Arc and Kernel Tracking Ignition Model (AKTIM) [10], and Spark Channel Ignition Monitoring Model (SparkCIMM) [15], which generally rely on a Lagrangian methodology to track the plasma or arc channel using flame speed expressions and initiate a flame kernel when certain criteria are met (such as Karlovitz number).

Ignition models, regardless of type, can include additional sub-models to improve the characterization of the physics of ignition and therefore increase the computational accuracy. An on-line secondary circuit predicts the secondary circuit voltage, current, and power based on the spark channel elongation. A short-circuit sub-model monitors the voltage across the spark channel and predicts the occurrence of a short-circuit that reduces the spark channel length and creates a drop in the voltage. A blowout sub-model monitors the current supplied to the spark channel and predicts a blowout if the current drops below a level deemed sufficient to maintain the spark. If the remaining energy in the system is more than the breakdown energy, then a re-strike is possible [17–19].

As computational power continues to grow, Eulerian models are becoming more appealing due to their simplicity and increasing accuracy with finer grid sizes. A shift from Lagrangian to Eulerian ignition models has occurred. Even historically Lagrangian ignition models have adopted a Eulerian formulation (for example AKTIM–Euler).

Earlier ignition research at Argonne concluded that detailed energy deposition can predict ignition success and misfires in quiescent conditions, assuming the correct input parameters are used with accurate source geometry, detailed chemistry, and conjugate heat transfer [20]. An advanced energy deposition hybrid Lagrangian–Eulerian Spark-Ignition (LESI) model was then developed to model the glow phase of ignition in cross-flow conditions [21]. LESI tracks the arc using Lagrangian methods and deposits energy using a Eulerian approach. The pre-glow (or breakdown) phase is modeled using offline equilibrium calculations to initialize pressure and temperature and account for species ionization.

This paper presents results from on-going work using the LESI model making it usable in the community of the engine modelers. LESI is coupled with finite-rate chemistry-based combustion models (WSR and TFM) and flamelet based combustion models (g-equation) to model spark channel elongation and kernel and flame development in a four-valve Direct-Injection Spark-Ignition (DISI) engine, using CONVERGE CFD solver [22]. First, in Section 2, the LESI model is described in detail. Second, in Section 3, an overview of the engine domain is provided, including the experimental and computational setups. Then, in Section 4, the simulation results are compared to experimental results, including global engine metrics and spark channel elongation images. Finally, in Section 5, the paper concludes, and main findings and future direction are discussed.

2 Method – LESI model

LESI is a line-source spark-ignition model for the glow phase of ignition, developed at Argonne within the CONVERGE CFD framework by using User Defined Functions (UDFs) for the glow phase of ignition. It has been discussed and evaluated in several previous publications [14, 21, 23]. In this section, the main features of LESI are listed and the coupling procedure with the combustion models is described.

2.1 Arc tracking and energy deposition

The LESI model is initialized based on user input as a line connected to both electrodes in the spark gap. The model creates a line of Lagrangian particles, with the general recommendation being one particle per computational cell. The model then tracks the motion of every particle in a Lagrangian fashion using an algorithm that relies on the resolved Eulerian velocity field. At the end of each time-step, energy is deposited in the Eulerian framework at the location of each Lagrangian particle. The total source energy, less any losses (for example resistance Ohmic losses), is distributed in the computational cells based on the length of segments that belong to that cell. Then, the model iterates to the next time-step and, again, starts by tracking the motion of the particles. The energy is deposited as a volumetric power source term in the energy equation of the computational cell that exists at the location of the corresponding Lagrangian point according to equation (1):(1)where ρ is the density, h is the enthalpy, u is the velocity, P is the pressure, σ_ij is the stress tensor, λ is the thermal conductivity, D is the molecular diffusivity, Y_m is the mass fraction of individual species, and S_ign is the energy source term introduced in specific cells, by the LESI model (J/m³ s).

The main features of LESI include:

• Lagrangian particles movement derived from the local resolved flow, where the level of elongation can be controlled by the user to mimic the behavior of electrical current. In addition, the algorithm utilizes the velocity of adjacent or neighboring Lagrangian particles to determine the motion of every individual particle, to smooth local instabilities out. Once two particles move too far apart, the model inserts new particles in between.

• The end particles are always connected to the electrode surfaces. The particles are allowed to move along the surface itself but the distance between the particle and the surface cannot exceed a minimum distance set by the user, as shown in Figure 2 where the red arrow denotes the movement of the top electrode end particle. In addition, if a center Lagrangian particle gets closer than the end particle to the electrode surface, the spark channel is truncated and the center particle becomes the new end particle, as can be seen in Figure 3, where the red and yellow arrows denote the location of the old and new end particles, respectively. All the particles in between are then removed.

Fig. 2 Attachment and movement of upper electrode end particle, denoted by the red arrow.

Fig. 3 Truncation of spark channel between frames 1 and 2, when a center Lagrangian particle gets too close to the lower electrode surface. Red and yellow arrows denote the location of the old and new end particles, respectively.

The LESI model leverages energy deposition on fine Eulerian grids and, when coupled with finite-rate chemistry, does not need explicit criteria for ignition, transition to turbulent flame, or quenching.

2.2 Coupling with combustion models

Naturally, Eulerian energy deposition models are well suited for finite-rate (or detailed) chemistry models such as the Well-Stirred Reactor (WSR) and the Thickened Flame Model (TFM). The model generates a local increase in temperature at the site of ignition (usually the spark plug gap) through an energy source equivalent to the total ignition energy in the experiment, less any Ohmic losses. Generally, in engines, the ignition source is on the order of mJ and is deposited in an L-shape profile to model the different phases of ignition. The temperature rise generated by the energy deposition increases the Arrhenius source term in individual species transport equations , which then leads to reactant consumption. After that, the flame is self-sustaining through heat and species diffusion and convection. The flame–flow interaction that controls kernel formation and subsequent flame propagation necessitates a fine grid for Eulerian models to be viable, especially when coupled with combustion and turbulence models. Since the LESI model relies on an Eulerian energy deposition scheme, the coupling with finite rate chemistry models follows the same scheme. From a practical standpoint, a volumetric power source term (J/m³ s) is introduced in the energy equation of the computational cells that overlap with the Lagrangian particles. The source term is equivalent to the total source energy divided by the duration of ignition and distributed among the Lagrangian particles based on segment length, as described in Section 2.1.

While Eulerian models couple seamlessly with finite-rate flame propagation models, they can still be coupled with flamelet concepts such as the g-equation model, which is perhaps the SI combustion model most widely used in industry (together with the Extended Coherent Flamelet Model [ECFM]). The hand-off between ignition and flame propagation occurs through temperature criteria or g-sourcing algorithms. In the former, the g = 0 iso-surface, which identifies the flame front, is initialized in computational cells with temperatures higher than a user set value (typically 3000 K). Then, the flame front is allowed to propagate at the turbulent flame speed, which is evaluated through a predefined expression (Eq. (2)). In the latter, an empirically determined G source term is introduced at the spark plug gap for the duration of ignition. The resulting g = 0 initialization is similar to that of the temperature criteria, even though the G source criteria is less physics-based. Using either approach, the early flame growth can be overestimated since immediate turbulent kernel behavior is assumed and the laminar-to-turbulent flame transition is overlooked. Hence, faster early flame propagation is expected with this approach. Laminar-to-turbulent transition functions for the g-equation model exist in literature and rely on an algebraic smoothing function, between laminar and turbulent behavior, rather than a physics-based function [5]. In this study, the coupling between the g-equation model and LESI is done through the temperature criteria with T_cut = 3000 K. A side effect of this treatment could be a reduced combustion sensitivity to ignition variability.

3 Method – DISI engine

In this paper, the LESI model is coupled with WSR, TFM, and g-equation models and used to carry multi-cycle Large-Eddy Simulations (LES) in a DISI engine.

The DISI engine is a pent-roof four-stroke four-valve engine with a Bowditch type piston and is located at Sandia National Laboratories [24]. The engine is equipped with a variable intake flow configuration to control in-cylinder large scale flow conditions, such as tumble and swirl. The engine can be seen in Figure 4 and the engine dimensions are listed in Table 1. The engine can be run in an all-metal configuration for continuous fired operation or in an optical configuration with a quartz window in the piston bowl and pent-roof gables.

Fig. 4 Sandia DISI engine: cross-sectional view of the combustion chamber. The valves, pent roof, spark plug and injector location and orientation are visible.

3.1 Experimental setup

Multiple experiments have been carried out in the DISI engine with final data sets including Particle Image Velocimetry (PIV) and direct spark-plasma and flame imaging. For this work, the engine runs at 1000 rpm with one intake valve deactivated to increase swirl in the combustion chamber. The eight-hole injector is slightly offset from the cylinder axis and is 5° from the vertical axis. The 60° injector produces eight spray plumes in a circular pattern. Fuel is delivered into the combustion chamber in three injection pulses, each delivering around 6 mg of E30 fuel (gasoline with 30% ethanol). The first injection has a SOIe of −298° CA after Top-Dead Center (aTDC) with a dwell of 15° CA between the three injections. The secondary ignition coil delivers around 106 mJ of ignition energy starting at −12.6° CA aTDC. A summary of the engine operating conditions is provided in Table 2.

3.2 Numerical setup

LES of the operating condition described in Section 3.1 were carried for ten consecutive engine cycles (for each setup) out using CONVERGE CFD solver v3.0 [22]. The computational domain includes the engine without the intake and exhaust plenums and can be seen in Figure 5. CONVERGE is a general-purpose CFD solver that can handle compressible and incompressible reacting flows in complex three-dimensional geometries. An orthogonal grid with a 4 mm base size is used in the computational domain. Several embedded refinements are included to achieve a minimum grid size of 0.125 mm at the spark plug gap. The grid size at the injection site is 0.25 mm. A dynamic structure turbulence model is used within the LES framework. A snapshot of the grid during ignition can be seen in Figure 6.

Fig. 5 Engine computational domain as visualized in CONVERGE STUDIO and includes the piston, cylinder walls, pent roof, valves, spark plug, valve ports and runners.

Fig. 6 Engine computational grid during ignition with the refined regions around the spark plug, piston, and pent roof visible. The grid size away from the spark plug gap and the cylinder and cylinder head walls is 1 mm.

Ignition is split into two sources. The first deposits 20 mJ in 15 μs to mimic the breakdown phase of ignition. The second source deposits 86 mJ in 2.3 ms, represents the glow phase of ignition, and is handled by the LESI model. The line source is initialized with nine (9) Lagrangian particles (to keep with the general recommendation of one particle per computational cell) and moves with the flow while depositing energy at every time-step (refer to Sect. 2.1). LESI requires that the spark plug surface file be split into at least two surfaces, to allow for proper identification of the lower and upper electrodes, as shown in Figure 7. During ignition, time-stepping is derived from a dynamic scheme that relies on user defined CFL numbers (velocity, viscosity, and speed of sound). For the results shown in this work the time step during ignition varied between 0.1 and 2 μs depending on the combustion model used.

Fig. 7 Spark plug surface split into more than one surface, as required by the LESI model.

WSR, TFM, and g-equation models are used for combustion modeling. WSR is a detailed (or finite-rate) chemistry model that solves the chemistry at every computational cell and the model output is highly dependent on grid size. TFM locates the flame front using a sensor function. The flame is then resolved in the computational cells identified by the sensor function using an artificial thickening scheme. Here, the sensor is based on the Heat Release Rate (HRR) and not the fuel species. The sensor uses Jaravel’s correction to improve the accuracy of the flame front location. The artificial thickening, however, reduces the flame wrinkling due to turbulence. Here, Charlette’s model [25] of the wrinkling factor is used to account for Sub-Grid Scale (SGS) flame-turbulence interaction and reduced flame wrinkling. Here, adaptive mesh refinement is not coupled with TFM, even though CONVERGE provides that capability [26]. Finally, the g-equation model solves the transport equation of the g variable which locates the flame front at g = 0. The Pitsch expression [27] is used for LES and can be seen in equation (2):(2)where S_l is the laminar flame speed, u′ is the sub-grid scale velocity, μ_t is the turbulent viscosity, μ is the viscosity, and b₁ and b₃ are modeling constants.

For both TFM and g-equation, tabulated laminar flame speed values are used, which were generated offline using CONVERGE 1D laminar flame speed chemistry solver, for the pressure, temperature, equivalence ratio, and dilution levels present in the combustion chamber at this operating condition. The tabulated flame speed values were generated using the same chemical mechanism used in TFM and WSR models. The numerical setup of the simulations is summarized in Table 3, where E30 fuel is a gasoline surrogate composed of C₂H₅OH, C₆H₅CH₃, IC₈, and NC₇H₁₆. The fuel mechanism and liquid properties were taken from previous studies on the same engine platform with the same fuel [28, 29].

4 Results

The results are split into two sections. The first section compares in cylinder metrics, such as pressure and AHRR, to experimental results and additional three-dimensional results for the different combustion models, for validation purposes. The second section compares arc tracking among the different combustion models with experimental images to assess the success of the LESI model in tracking the spark channel. Kernel formation among the different combustion models is also compared and discussed.

4.1 Well-Stirred Reactor model

The WSR model serves as the baseline for this study. Generally speaking, the results of such an approach are dependent on the grid size. In this study, a reaction multiplier was tuned to match peak experimental pressure traces. Within the LES framework, the coupling between LESI and the WSR model is successful, and the peak in-cylinder pressure is shown in Figure 8 while the AHRR is shown in Figure 9 for LES. Both combustion metrics are compared to experimental results, which are plotted in green. Prior simulations with Reynolds-Averaged Navier–Stokes (RANS) (not included here) showed lower levels of Cycle-to-Cycle Variability (CCV) than what is observed here, as is expected.

Fig. 8 WSR and LESI ignition model: peak in-cylinder pressure for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.

Fig. 9 WSR and LESI ignition model: in-cylinder AHRR for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.

The average experimental peak pressure in Figure 8 is approximately 23.3 bar, while that predicted by simulation is approximately 23.5 bar, a difference of 0.2 bar or 0.86%. While the tuning process matches the average peak in-cylinder pressure with its experimental counterpart, the WSR model struggles to match the pressure rise. Even when the peak pressure has the same magnitude, the WSR model predicts a delayed pressure rise and a retarded peak pressure. The shortcomings of the WSR model are likely due to inadequate grid resolution to fully resolve the flame front.

The mismatch in pressure rise is mirrored in the AHRR profile shown in Figure 9. HRR is underestimated prior to 10° CA aTDC and overestimated afterwards. The experimental average reaches a maximum HRR of approximately 24.8 J/CA at around 9° CA aTDC, while the simulation average reaches a maximum HRR of approximately 35 J/CA at around 18° CA aTDC. The pressure rise slope in Figure 8 exhibits the same behavior, where the slope is flatter than what is observed in the experiment prior to 10° CA aTDC. The early flame development and heat release are likely underestimated due to insufficient grid resolution.

The objective of this work is to couple the LESI model with commonly used combustion models for SI engine simulations. It follows that, despite the shortcomings of the WSR model in predicting in-cylinder metrics with the presented numerical setup, the coupling with the LESI model is successful. LESI models the ignition event in the DISI engine without incurring significant computational costs, if any, and is compatible with grid changes originating from Adaptive Mesh Refinement (AMR) or fixed embedding zones. The main advantage of LESI over state-of-the-art models, is its ability to track the arc or spark channel as it moves and stretches with the flow. To that end, Figure 10 shows the arc movement as predicted by LESI for two representative cycles at different stages of ignition. Recall that ignition starts at −12.6° CA aTDC and ends around 1.1° CA aTDC. For cycle 1 in Figure 10, the arc stretches inwards towards the lower electrode (or the J-gap). The Lagrangian particles close to the lower end of the arc experience a stronger velocity field, as denoted by the color of the spheres used to visualize the arc, and thus stretch creating a unique arc shape. For cycle 2, on the other hand, the arc moves in a completely different manner. The elongation experienced by the arc is smaller compared to that in cycle 1, and the arc moves inwards then to the left. The local flow field that convects the arc can have many sources, which are not mutually exclusive: leftover turbulence from large-scale flow structures originating from the intake stroke, directional velocity field originating from the injection event, or instantaneous velocity changes due to flame kernel formation and growth. Regardless of the origin of these changes to the local velocity field, LESI captures those variations and models their effect on the spark channel. While the operating condition evaluated in this paper is characterized by high swirl and moderate tumble, LESI captures cyclic variability in arc elongation and stretch and transfers this variability to kernel formation and subsequently combustion, thus achieving its intended goal.

Fig. 10 WSR and LESI ignition model: arc movement as predicted by LESI for two representative cycles at different stages of ignition. Individual Lagrangian points are colored by the flow velocity at the particle location.

4.2 Thickened Flame Model

Given that the shortcomings in predicting engine performance using the WSR model in Section 4.1 are most likely due to inadequate flame resolution using the current grid, TFM could improve the results. Remaining within detailed chemistry models, TFM builds on top of WSR to resolve the flame front through manipulation of mixing and chemistry terms.

Since TFM alters the heat release properties of the flame, it is not recommended to be turned on immediately with energy deposition as it might lead to slow flame kernel growth and even premature quenching. As a result, TFM is turned on at −10° CA aTDC, roughly 2° CA after spark timing. In addition, for these LES results, TFM relies on a laminar flame speed table generated for E30, a flame sensor with Jaravel’s adjustment, and Charlette’s model for the wrinkling factor (sometimes called efficiency function), with an exponential power factor of 0.6. The flame speed table also contains flame thickness and reference reaction rates to be used in the flame sensor and efficiency function formulations. Ten grid points are utilized to cover the flame front, or in other words, the minimum grid size required to fully resolve the flame is 10% of the laminar flame thickness, as shown in equation (3). TFM does not use an explicit turbulent flame speed expression:(3)

TFM is expected to sufficiently resolve the flame front, and that can be seen in Figure 11, where the average simulation peak pressure is approximately 24.1 bar, a difference of 0.8 bar or 3.4% versus the experiment. In addition, on average, TFM overestimates the pressure rise at times earlier than 10° CA aTDC, which is a reversal of what was seen with the WSR in the previous section. With TFM, the flame front is better resolved, and the effect of turbulence is accounted for, which leads to higher initial pressure rise.

Fig. 11 TFM and LESI ignition model: peak in-cylinder pressure for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.

The changes in in-cylinder pressure prediction by TFM over WSR is mirrored in the AHRR results shown in Figure 12. TFM predicts a peak HRR of approximately 27.6 J/CA at around 8° CA aTDC, while the experimental average reaches a maximum HRR of approximately 24.8 J/CA at around 9° CA aTDC. TFM improves the overall shape of the HRR plots, better predicting the location and magnitude of the peak as well as the profile shape. While the initial heat release can be reduced by delaying the onset of the TFM activation, fine tuning the models to provide an exact match with experiment does not significantly contribute to the scope of this work. The purpose of this research is to present successful coupling of LESI with common turbulent flame propagation models.

Fig. 12 TFM and LESI ignition model: in-cylinder AHRR for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.

In addition to improving the simulation results over the WSR model, TFM is compatible and couples seamlessly with the LESI model. Iso-contours of the thickening factor F and the wrinkling factor at −5° CA aTDC shown in Figure 13 indicate that TFM is working as intended alongside LESI. First, the yellow shading in Figure 13, which denotes the flame front as identified by the flame sensor, overlaps the high temperature region. Second, the magnitude of F changes abruptly between two zones (light blue color versus orange color), which indicates that the grid size has changed due to fixed embedding and TFM is responding adequately. In addition, the magnitude of F ranges between 50 and 100 which is typical and consistent with the grid size and n_res of 10. Finally, is concentrated at the flame front and has a larger magnitude at the leading edge of the flame where F has a larger magnitude and turbulence is expected to be higher, which is consistent with the expectations and is another indication that TFM is working as intended. In addition, at this time, the shape of the arc is different than what has been seen with the WSR model, another indication that LESI is capturing the flow variability on the spark channel.

Fig. 13 TFM and LESI ignition model, from left to right: iso-contours of temperature (in K), iso-contours of the thickening factor F, iso-contours of the wrinkling factor , and arc shape as predicted by LESI model at −5° CA aTDC. The yellow shading denotes the flame front as identified by the TFM sensor.

4.3 g-equation model

The g-eq model is a popular flame propagation model and widely used by automakers. Hence, coupling the g-eq model with an advanced ignition model is inherently of interest. The model relies on solving transport equations for G, where g = 0 is defined as the flame front. The propagation of g is mainly controlled by the turbulent flame speed S_t. In the LES framework, S_t is based on the expression by Pitsch [27]. Here, the expression is tuned with b₁ = 30 and b₃ = 1.65. The g-eq model also relies on the chemistry equilibrium solver in the burnt and flame regions, which allows for a fast engine simulation.

The peak in-cylinder pressure as predicted by the g-eq model combined with the LESI ignition model is shown in Figure 14. The average LES peak pressure is approximately 21.9 bar, a difference of 1.4 bar or 6% versus experiment. In addition, the g-eq model improves the slope prediction at early times (earlier than −10° CA aTDC), an improvement over both TFM and WSR models. The individual pressure cycles, however, are less spread which is an indication of lower CCV being predicted (more on this in Sect. 4.4).

Fig. 14 g-equation model and LESI ignition model: peak in-cylinder pressure for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.

The AHRR as predicted by the g-eq model combined with the LESI ignition model is shown in Figure 15. The g-eq predicts a peak HRR of approximately 24.7 J/CA at around 9° CA aTDC, while the experimental average reaches a maximum HRR of approximately 24.8 J/CA also at around 9° CA aTDC. Qualitatively, the g-eq model reproduces the heat release profile observed in experiment, a reflection of the improvement on the pressure rise slope seen in Figure 14, when compared to TFM and WSR models.

Fig. 15 g-equation model and LESI ignition model: in-cylinder AHRR for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.

Similarly to the WSR and TFM models, here the coupling between the g-eq and LESI models is discussed to ensure the models work properly alongside each other. Iso-contours of the g progress variable and turbulent flame speed at −5° CA aTDC shown in Figure 16 indicate that g-eq is working as intended alongside LESI. The flame front as identified by the g = 0 surface is thin when compared with TFM, as the g-eq model is designed to track the flame front. The g flame front coincides exactly with the high temperature region. In addition, the turbulent flame speed is lower in the burnt region when compared to the flame front and unburnt regions. Finally, at this time, the shape of the arc is different than what has been seen with the TFM and WSR models (moving away from the spark plug gap), an indication that LESI is capturing the local flow variability on the spark channel between cycles.

Fig. 16 g-eq model and LESI ignition model, from left to right: iso-contours of temperature, iso-contours of the g progress variable, iso-contours of the turbulent flame speed St, and arc shape as predicted by LESI model at −5° CA aTDC.

4.4 Arc Tracking and Flame Initialization

In Sections 4.1–4.3, engine simulation results with the WSR, TFM, and g-eq models were validated against experimental pressure and AHRR traces. In addition, the coupling between LESI and the combustion models was shown to be successful, where the models are behaving as expected alongside each other. Here, the arc tracking as predicted by LESI is validated against optical engine data. In addition, kernel formation and CCV under the different combustion models are inspected.

Figure 17 shows the spark channel (or arc) and flame propagation for the three combustion models and LESI compared to experimental images for two representative cycles at three different times during which ignition is active. In the experimental images, the flame front is represented by the blue shading. The arc as predicted by the LESI model is represented by the yellow line in the vicinity of the spark gap while the flame surface is represented by the 1500 K iso-surface. The engine operating condition is characterized by high levels of swirl and moderate levels of tumble. The spark channel is not expected to greatly elongate outwards of the spark plug gap but rather to convect sideways with minimal levels of elongation depending on the direction of the swirl flow in the individual cycle. The experimental images in Figure 17 show that this is indeed the case, where for cycle 1 the arc convects to the left of the spark plug while for cycle 2 the arc convects to the right and back to the left of the spark plug, with minimal elongation. In addition, the simulation images for all combustion models show a similar behavior being predicted by LESI. In cycles WSR-1 and g-eq-1 the arc moves to the left of the spark plug, while in cycles WSR-2, TFM-1 and g-eq-2 the arc moves to the right of the spark plug. In cycle TFM-2 the arc shows minimal movement to the left then returns to the center with spark channel elongation being parallel to the plane of view. In summary, the LESI model is capable of fully capturing the spark channel (or arc) shape, location, and elongation for this operating condition.

Fig. 17 Spark channel (or arc) and flame propagation for the three combustion models (WSR, TFM, and g-eq) and LESI compared to experimental images for two representative cycles at three different times during which ignition is active. The flame surface is represented by the T = 1500 K iso-surface and the arc is represented by a yellow line as predicted by the LESI model.

Given LESI’s improved arc tracking capabilities over traditional spherical energy deposition, it is expected to provide more predictive results in lean and dilute operating conditions that exhibit higher CCV levels than the operating condition presented here. Nevertheless, the flame kernel formation can be inspected, and conclusions can be made about the validity of each model combination for CCV prediction. Figure 18 shows temperature iso-contours for the three models at different times during ignition in three representative cycles, as an analogy for flame kernel growth. The g-eq model shows the fastest flame kernel formation and growth when compared to TFM and WSR, and, as previously hinted at, is due to the assumption of immediate turbulent behavior after flame front initialization using the temperature cut-off method. The flame kernel at its early stages grows according to the laminar flame speed S_l. Once the flame kernel reaches a large enough size to interact with turbulent eddies that wrinkle the flame front, it propagates at S_t. Hence, with the g-eq model, early flame development is overestimated using the temperature cut-off method and a laminar to turbulent transition scheme, as proposed in [5], is potentially needed.

Fig. 18 Kernel growth for the three combustion models and LESI at three times during ignition, for representative cycles.

When comparing Figures 17 and 18 with the in-cylinder pressure and heat release rate results in Sections 4.1–4.3, a discrepancy arises because the g-equation model provides the best agreement of global results with experimental results among the three combustion models (Figs. 14 and 15). However, the underlying formulation of the combustion models is fundamentally different. The g-eq model employs a turbulent flame speed expression that controls the flame kernel growth. A chemical equilibrium solver with only the global species is used in the burnt and flame regions to estimate the temperature rise and subsequent pressure rise. Hence, the tuning of the expression compensates for deficiencies in the chemistry solver to match the experimental in-cylinder pressure and heat release. As a result, while the global engine results may show a high level of agreement with experimental results, the kernel size and perhaps the species mass fractions would show a lower level of agreement, which is likely the case here. On the other hand, with WSR and TFM the kernel growth and temperature rise are controlled by a reduced chemical mechanism, rather than a turbulent flame speed expression. As such, while the results will be sensitive to the grid resolution, these models, at least in theory, should provide better agreement of kernel growth against experimental images or engine optical data, since laminar-to-turbulent flame transition is modeled to a certain extent, at least compared to the g-eq model. In other words, in the g-eq model a tuned turbulent flame speed expression dictates the flame kernel growth to predict in-cylinder pressure and heat release, regardless of the chemistry. The flame kernel growth predicted by CFD in the g-eq case does not qualitatively look like the optical images in Figure 17, hinting that locally the model lacks the necessary level of detail. For the WSR and TFM, the detailed chemistry solver predicts the kernel growth and temperature rise (refer to Fig. 17), which then dictate the in-cylinder pressure and heat release and given adequate grid resolution, provides accurate flame kernel growth.

Overestimation of flame growth with the g-eq model would lead to a reduction in the effects early flame growth has on overall cycle combustion, and that includes cyclic variability. In other words, the higher fidelity level brought by advanced ignition models such as LESI would not transfer to higher variability levels observed at full cycle combustion. Table 4 lists the range of peak in-cylinder pressures predicted by each simulation setup presented in this work. Even though the WSR model did not have the best agreement in heat release (Fig. 9 in Sect. 4.1), here, it shows the best agreement (8.16 bar) compared to the experiment (7.1 bar). The g-equation model shows a smaller range (less variability: 3.63 bar) while TFM exhibits a larger range (more variability: 13.25 bar, likely due to the slow burn cycle predicted by TFM. The pressure range would be around 8 bar if the slow burn cycle is ignored, which is in line with WSR and experiment). Furthermore, the results discussed so far are mirrored in the combustion phasing behavior of every case. Figure 19 shows seven combustion phasing metrics (CA1, CA2, CA5, CA10, CA25, CA50, and CA90) for every case, compared to the experiment. Every scatter point on the plots represents an individual cycle. For the WSR model, early combustion is retarded while late combustion is advanced, which reinforces the conclusions of Section 4.1. On the other hand, for TFM and g-eq combustion phasing is more aligned with experiment, as discussed in Sections 4.2 and 4.3. The spread of the scatter plots at each CA in Figure 19 also confirms that the WSR model shows a level of variability comparable to experiment, while g-eq model underestimates variability and TFM slightly overestimates it, as discussed earlier in Table 4.

Fig. 19 Different combustion phasing metrics for every simulation case (10 consecutive cycles each) compared to the same metrics from 35 experimental cycles. Every scatter point represents one cycle.

In summary, while the g-eq model accurately predicts global pressure and heat release profiles, when coupling it with advanced ignition models, a transition function is necessary to ensure early flame growth is well captured to better predict cyclic variability, and by proxy, ignition effects on early flame kernel formation and growth. For detailed chemistry models a transition function might not be necessary but sufficiently fine grids or model setup are required to resolve kernel formation and growth and subsequently heat release profiles. The importance of capturing ignition effects on early flame growth is even greater in operating conditions where the spark channel is expected to stretch, elongate, short circuit, and re-strike. In other words, unconventional operating conditions that promise efficiency improvements and CO₂ emissions reductions will likely lead to strenuous ignition behavior (such as re-strikes). Providing models that not only capture this behavior but carry it downstream in the simulation is important for the design of such engines and for realizing CO₂ reductions.

5 Conclusion

In this paper, the LESI advanced energy deposition ignition model was coupled with turbulent flame propagation models, more specifically the WSR, TFM, and g-equation model, and used to carry out LES of the Sandia DISI engine.

First, the main features of the LESI model were highlighted, which include: a Lagrangian particle arc tracking algorithm which is based on the local velocity field and an Eulerian energy deposition scheme which deposits fractional energy based on the segment length of the respective Lagrangian particle. Second, the coupling of the LESI model with other combustion models was discussed, where LESI couples with the detailed chemistry models seamlessly and uses a temperature cut-off approach to couple to the g-eq model.

Then the DISI engine is discussed in detail with the experimental operating condition and the numerical setup being described. The main engine parameters are an intake-induced swirl motion and stoichiometric operating condition with a three-pulse injection of E30 fuel at 1000 rpm. The simulations rely on a dynamic structure turbulence model within a LES framework with AMR and fixed embedding and combination of LESI with the different combustion models.

The results showed that while the WSR model was tuned to agree with the peak in-cylinder pressure, it was not able to match the AHRR and slope of the pressure rise, likely due to insufficient grid resolution. TFM, being an add-on model to WSR, improved the agreement with experiment in peak in-cylinder pressure magnitude and location, as well as the slope of the pressure rise and AHRR. Finally, g-eq model was tuned to provide adequate agreement in pressure and HRR but exhibited lower levels of variability, than the WSR and TFM. All three models coupled well with LESI and operated within the respective expected behavior. While all models require improvement in flame kernel growth, WSR and TFM showed better results than the g-eq model. On the other side, the g-eq model had the best pressure and heat release profiles.

The spark channel elongation was then inspected in all the simulations and compared to experimental images, where it was seen that LESI successfully predicted arc movement in the spark gap for this engine operating condition. Finally, it was shown quantitatively that the g-eq model predicted less CCV than the WSR, TFM, and experiment, which was attributed to capturing less of the early flame growth due to the assumption of immediate turbulent behavior. It was concluded when coupling the g-eq model with advanced ignition models, a transition function would be necessary to ensure early flame growth is well captured, and by proxy, ignition effects on early flame kernel formation and growth.

In conclusion, the LESI model has been shown to be a successful advanced ignition model that couples well with the most commonly used combustion models for SI engine simulations. LESI accounts for the impact of the flow field variability on the ignition event, and by proxy, can account for any spatial variability (e.g., turbulence, equivalence ratio stratification, dilution, etc.) Future work for LESI includes similar engine studies in high-tumble, lean, and/or dilute operating conditions.

Acknowledgments

This manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). 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 is funded by DOE’s Vehicle Technologies Program, Office of Energy Efficiency and Renewable Energy, under the FY19 Technology Commercialization Fund (TCF) program. The LESI model has been developed and continues to be further improved as part of the Light-Duty National Laboratory Consortium PACE (Partnership to Advance Combustion Engines). The authors would like to express their gratitude to Gurpreet Singh, Mike Weismiller, and Kevin Stork, program managers at DOE, for their support. Numerical simulations were run on the Blues and Bebop Clusters at the LCRC, Argonne National Laboratory. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.

References

All Tables

All Figures

	Fig. 1 History and projections of the United States transportation sector energy consumption: (a) by transportation mode and (b) by fuel source [1].
In the text

	Fig. 2 Attachment and movement of upper electrode end particle, denoted by the red arrow.
In the text

	Fig. 3 Truncation of spark channel between frames 1 and 2, when a center Lagrangian particle gets too close to the lower electrode surface. Red and yellow arrows denote the location of the old and new end particles, respectively.
In the text

	Fig. 4 Sandia DISI engine: cross-sectional view of the combustion chamber. The valves, pent roof, spark plug and injector location and orientation are visible.
In the text

	Fig. 5 Engine computational domain as visualized in CONVERGE STUDIO and includes the piston, cylinder walls, pent roof, valves, spark plug, valve ports and runners.
In the text

	Fig. 6 Engine computational grid during ignition with the refined regions around the spark plug, piston, and pent roof visible. The grid size away from the spark plug gap and the cylinder and cylinder head walls is 1 mm.
In the text

	Fig. 7 Spark plug surface split into more than one surface, as required by the LESI model.
In the text

	Fig. 8 WSR and LESI ignition model: peak in-cylinder pressure for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.
In the text

	Fig. 9 WSR and LESI ignition model: in-cylinder AHRR for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.
In the text

	Fig. 10 WSR and LESI ignition model: arc movement as predicted by LESI for two representative cycles at different stages of ignition. Individual Lagrangian points are colored by the flow velocity at the particle location.
In the text

	Fig. 11 TFM and LESI ignition model: peak in-cylinder pressure for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.
In the text

	Fig. 12 TFM and LESI ignition model: in-cylinder AHRR for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.
In the text

	Fig. 13 TFM and LESI ignition model, from left to right: iso-contours of temperature (in K), iso-contours of the thickening factor F, iso-contours of the wrinkling factor , and arc shape as predicted by LESI model at −5° CA aTDC. The yellow shading denotes the flame front as identified by the TFM sensor.
In the text

	Fig. 14 g-equation model and LESI ignition model: peak in-cylinder pressure for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.
In the text

	Fig. 15 g-equation model and LESI ignition model: in-cylinder AHRR for LES (blue) and experiments (green). Individual cycles are denoted by the dashed lines and the ensemble average is denoted by bold lines.
In the text

	Fig. 16 g-eq model and LESI ignition model, from left to right: iso-contours of temperature, iso-contours of the g progress variable, iso-contours of the turbulent flame speed St, and arc shape as predicted by LESI model at −5° CA aTDC.
In the text

	Fig. 17 Spark channel (or arc) and flame propagation for the three combustion models (WSR, TFM, and g-eq) and LESI compared to experimental images for two representative cycles at three different times during which ignition is active. The flame surface is represented by the T = 1500 K iso-surface and the arc is represented by a yellow line as predicted by the LESI model.
In the text

	Fig. 18 Kernel growth for the three combustion models and LESI at three times during ignition, for representative cycles.
In the text

	Fig. 19 Different combustion phasing metrics for every simulation case (10 consecutive cycles each) compared to the same metrics from 35 experimental cycles. Every scatter point represents one cycle.
In the text