The AVL optical engine was modeled and the numerical simulations were executed by the density-based finite volume CFD solver CONVERGE T M 3.0 utilizing its multi-physics capabilities Richards et al. (2021). The computational domain is illustrated in Figure 2. The Reynolds-averaged Navier-Stokes (RANS) formulation was used with the RNG κ − ϵ closure model Han and Reitz (1995). Pressure-implicit with splitting of operators (PISO) algorithm was employed to couple pressure and velocity dynamics Issa et al. (1986). Wall heat transfer was modeled with the GruMo-UniMORE model developed by Berni et al. (2017). The Redlich-Kwong equation of state was used along with real gas properties and mixture-averaged diffusion coefficients. Turbulent Prandtl and Schmidt numbers were set as 0.9 and 0.78, respectively. The multi-zone well-stirred reactor (MZ-WSR), referred to as SAGE, combustion model was selected as a baseline combustion submodel, with the kinetic mechanism by Zhang et al. (2021). Furthermore, G-equation combustion model is assessed. The initial field variables were mapped from a prior simulation to replicate the effects of field initialization. The simulations began at −175 CAD and continued until 900 CAD, thereby covering a full cycle. Variable time-step algorithm was implemented with the minimum time step governed by the critical Courant-Friedrichs-Lewy number.

Boundary conditions were imposed based on the experimental measurements Uddeen et al. (2023a); Uddeen et al. (2023b); Uddeen et al. (2023c). The mixture composition was stoichiometric and assumed as a perfectly mixed inflow boundary condition. Constant Dirichlet pressure inlet and outlet boundary conditions were set, simplifying the exchange gases process, yet capturing the correct delivered mass per cycle. Piston, liner and head walls temperatures were kept at 450 K, while the intake and exhaust systems temperatures were maintained at 333 K and 500 K, respectively. To model ignition, a spherical energy source was placed at each spark plug electrode gap. The ST was set to − 30 CAD. For G-equation, ignition was achieved using passive source modeling with a G = 0 surface initialized at the spark location as a spherical surface of 0.5   m m . Alternatively, the multi-zone well-stirred reactor model used a high- and low-energy source to capture ignition, depositing 0.04   J and 0.02   J in a spherical region of 0.5   m m radius over 0.5 CAD and 10 CAD, respectively. The five-spark combustion configuration was selected for evaluation because it operates at a low COV, making the experimental data a reliable baseline for modeling. A base mesh size of 2.5   m m was employed throughout the entire domain. Furthermore, adaptive mesh refinement (AMR) was applied in regions with strong gradients in temperature and velocity, and embedded mesh refinement was used to capture near-wall flow characteristics. Spherical embeddings were implemented around the spark plug regions. Mesh features are shown in Figure 2, and details are provided in Table 2 for the five different configurations that were evaluated.

The first step in the calibration of the engine model was related to the proper representation of the compression and expansion of the in-cylinder gases. Optical engines experience a significant level of leakage due to hardware modifications and the setup configuration. Therefore, to assure accurate reproduction of the pressure history, a proper crevice model is crucial. Following the studies by Irimescu et al. (2013); Irimescu et al. (2018); De Renzis et al. (2021), 7.5% of the total inducted mass per cycle was allowed to leak to the crankcase. The model was further calibrated while maintaining the squish height as in the engine, leading to a geometrical compression ratio of 12. The results are presented in Figure 3, revealing close agreement with the experimental data, particularly during the expansion stroke. This level of accuracy is especially noteworthy given the optical configuration of the engine under investigation.

Subsequently, combustion cases are calibrated, starting from mesh sensitivity; different mesh strategies were tested and as summarized in Table 2. The baseline case setup is established based on literature Yang et al. (2023), whereas the additional mesh features were considered to improve the description of the initial kernel growth with the refinement of the mesh in the spark plug region. Additionally, a cylindrical embedding of the combustion chamber was explored in meshes three and four to observe the keys aspects of the main chamber turbulence dissipation, along with the AMR to refine the flame zones. Finally, the geometry of all the spark plugs was introduced in Mesh four and its effects on the combustion behavior will be explained in detail below.

Figure 4a presents the pressure and apparent heat release rate (AHRR) evolutions for the different evaluated cases, with 100 experimental combustion cycles shown in gray for reference. At first glance, the Baseline, Mesh 1, and Mesh two simulations underpredict the pressure rise despite accurately capturing the ignition delay, indicating a slower heat release rate compared to the experiments. This discrepancy can be attributed to the coarse mesh resolution in the spark plug region and insufficient refinement during the early stages of combustion. In particular, coarse meshes promote a rapid energy dissipation, which reduces the reactivity of the cell-based ignition in the SAGE combustion model.

Additional evidence appears in Figure 4b, which shows how the temperature in the SPC vicinity rapidly decays for the coarse meshes. In the Baseline mesh, the peak temperature is even lower due to this rapid energy dissipation. A second temperature peak occurs once combustion initiates; notably, this peak appears earlier for Mesh three and takes longer in the coarser meshes, indicating a slower heat release process.

In contrast, Mesh three overpredicts the heat release rate and the peak pressure reaches higher values compared to experiments. Although the mesh in the spark plug surroundings is refined, as well as the combustion chamber, the early flame development is shorter and a rapid rise of the in-cylinder pressure is substantial. To understand how the local flow field and early flame development are affected by the presence of spark plugs geometry (J-gap type), their geometrical details are included in Mesh 4. Remarkably, the prediction is much improved, suggesting the importance of capturing the near-spark flow field during the early combustion development. To obtain additional insights about the importance of the detailed geometry, additional test simulations were done with the spark plugs rotated (results not shown in the document), but marginal changes were observed in the numerical results. It is thus concluded that the most important aspect is the local flow field near the wall.

These findings suggest that typical engine CFD meshing practices for hydrocarbon fuels may not be fully applicable to NH 3 -fueled engines, primarily due to NH 3 ’s lower reactivity and narrower flammability range. Under coarse mesh resolutions, numerical dissipation can quickly dilute the ignition kernel and reduce reactivity, making it more difficult to accurately capture the onset and early stages of NH 3 combustion.

Different heat transfer models may be selected for the law-of-the-wall temperature boundary condition treatments available in CONVERGE, such as GruMo-UniMORE by Berni et al. (2017), Amsden, 1997, and Han and Reitz (1997). A proper heat transfer description is necessary to correctly describe the flow physics in the combustion chamber, especially at the end of the compression stroke and during combustion period when the walls are experiencing higher heat transfer rates. As the viscous sublayer is not resolved, it is necessary to model the temperature profile near the wall. In addition to the heat transfer equation, the transition of the logarithmic layer, y + , is also important. In general, the GruMo-UniMORE model is recommended for highly-charged/highly-downsized spark-ignited engines and account for different transient elements during the power cycle; Han and Reitz incorporates compressible effects and offers a transient description of the heat transfer processes; O’Rourke and Amsden rely on the quasi-steady assumption, where the heat transfer process has reached a local equilibrium.

Figure 5 shows that the implementation of the three different models lead to a similar behavior in the pressure and AHRR evolution, with the O’Rourke model predicting a slightly higher peak pressure and AHRR. This suggests that the heat losses are minimal when using this model. On the other hand, the AHRR peak in GruMO-UniMORE and Han and Reitz is similar, with the latter predicting a longer combustion duration. In addition, the pressure peak is lower for the Han and Reitz results, indicating higher heat transfer losses, observed as well in the lower pressure captured towards the end of the compression stroke. Considering the small variations, for this specific engine configuration running with NH 3 as fuel, any of the three heat transfer models are found to be adequate in predicting the engine combustion behavior.

The reference ST was set at −30 CAD aTDC. Given the lack of data in optical NH 3 engines, however, the uncertainties in arc breakdown timing were assessed by varying the ST by ± 2.5 CAD, as shown in Figure 6. Although the AHRR rise begins at a similar time for earlier STs, the retarded ST exhibits a significant ignition delay. Furthermore, the three cases demonstrate distinct pressure evolutions. A monotonic reduction in the pressure peak and an increased ignition delay are observed with retarded STs. These results indicate that while the HRR for earlier STs remains similar, the retarded case shows a significantly lower peak. Due to the slow flame speed of NH 3 , delaying the ST results in a substantial reduction in both pressure peak and AHRR, along with a further delayed ignition. However, more data is needed to draw a conclusive understanding of the sensitivity to ST.

To further demonstrate the fidelity of the model, the instantaneous flame development and topology is analyzed in Figure 7. In the first row (a) optical experiments are presented, whereas the second row (b) shows the computational volumetric rendering of the temperature field. Upon initial inspection, the CFD results demonstrate a qualitatively accurate prediction of flame development from the five different spark locations. However, the numerically predicted flame growth appears slightly faster than in the experiments, which may relate to the higher AHRR observed in Mesh 4 (see Figure 4). Furthermore, it should be noted that the experimental optical setup does not capture the entire outer region of the combustion chamber, and the single-shot image from the experiment may differ from the ensemble-averaged RANS solution. In the initial stages of combustion at −10 CAD, the flame surface areas are small, and then grow towards top dead center (TDC), increasing the pressure and temperature of the combustion chamber. At 10 CAD, the different flame fronts have merged towards the late combustion, consuming nearly the entire unburned premixed charge. As reported in previous studies by Huang and Liu (2024), the multi-spark methodology speeds up the overall charge consumption due to the simultaneous creation of multiple flame fronts. One important observation is the flame resultant from the central spark (SPC) growing with a faster rate. Furthermore, it is observed that the topology for the different flames differ among them, as further discussed later.

The volumetric-averaged turbulent kinetic energy (TKE) at each spark location is presented in Figure 8. There are two important peaks during the cycle. The first one corresponds to the intake valve opening event, promoting the increase in the TKE as a consequence of the entering turbulent inflow mixture charging the main chamber. The second peak occurs at the ST, where the sudden rise in the pressure and the expansion of the flames disturb the local flow field. It is seen that TKE is higher at SPC before the start of ignition, which promotes the central flame to grow faster. Additionally, higher fluctuations are observed in this region. There is irregular behaviors in the different regions during the intake stroke. However, the level of turbulence decays once the intake valves are closed. Towards the compression stroke, the fluctuations in the flow field remain at similar low values.

It is evident that each flame front behaves distinctly as a consequence the different local flow field conditions at ST and afterwards; this implies the need for appropriate mesh selections at different regions to capture the detailed turbulent flow characteristics in addition to the spark geometry. Furthermore, due to the turbulent flows and the swirl-type engine, the in-cylinder fluid dynamics varies spatially, translating in a complex flame-flame interaction. The present results show significant differences from those obtained by Huang and Liu (2024), where a simplified initial setup resulted in multiple flame fronts propagating in comparable shapes and sizes. The present study suggests that capturing local variations in the flow and flame behavior may lead to a significant difference in the prediction of the overall engine combustion characteristics.

More details of the flow field are shown in Figure 9 along with the TKE map. The higher TKE found in the center of the combustion chamber is found to be caused by the swirl motion. The rotational motion promotes an increase in the pressure close to the walls, where the velocity is zero due to the non-slip condition; the boundary layers of the piston and head induce the flow to move inwards to the center of the combustion chamber. As a consequence of the different streams collapsing, velocity fluctuations are concentrated in that region. The iso-surface for T = 1700 K represents the flames in Figure 9. Due to the growth of the multiple flames and the heat release during combustion, the gas expansions at different locations affects the flow patterns, where multiple locus of flow stagnation are observed as a result of the flames propagating onto the fresh charge.

Turbulence varies significantly both spatially and temporally throughout the cycle. These variations have a crucial influence on flame dynamics and need to be properly captured in the model. The SAGE combustion model treats each computational cell as a homogeneous (well-stirred) reactor. Accordingly, within the RANS framework, turbulence-chemistry interactions (TCI) are effectively neglected, as local reaction rates are computed solely from the mean flow field and temperature without accounting for sub-grid fluctuations Mohan and Haworth (2015). To properly account for TCI in the context of flame propagation, the G-Equation combustion model is being employed, where a transport equation for the scalar field variable G is solved and the iso-scalar surface G ( x , t ) = G 0 represents the flame surface at each instantaneous solution field, separating the burned and unburned gases. The proposed transport equation per Peters (2001) is written as: ∂ ρ G ̃ ∂ t + ∂ ρ u ̃ i G ̃ ∂ x i = ρ u S t ∂ G ̃ ∂ x i − ρ D t ′ κ ∂ G ̃ ∂ x i (1)

Equation 1 comprises several variables and constants, including the turbulent flame speed represented by S t , the unburned mixture density ρ u , the local density ρ , the Favre mean value of the scalar G represented by G ̃ , the curvature κ , and the turbulent diffusivity D t ′ . The first term on the right hand side of Equation 1 represents the mean turbulent burning rate, while the second term represents the impact of the flame front curvature. To solve Equation 1, it is necessary to close the system by first determining the turbulent flame speed, S t , which is is calculated using Peters’s correlation, defined as: S t = S l + u ′ − a 4 b 3 2 2 b 1 D a + a 4 b 3 2 2 b 1 D a 2 + a 4 b 3 2 D a , (2) where the turbulent Damköhler number is defined as D a = ℓ   S l   /   ( ℓ F   u ′ ) , which characterizes the ratio of flow timescales to chemical timescales; u ′ is the turbulence intensity (derived from the local turbulent kinetic energy), and ℓ is the integral length scale, obtained from the RANS formulation as ℓ = C μ 0.75   k 1.5 / ϵ . The laminar flame speed, S l , is tabulated by the built-in chemistry tool. The flame thickness, ℓ F , is estimated as ℓ F = λ m / c p 0   /   ρ   S l u , where λ m is the molecular thermal conductivity, c p is the specific heat, and ρ is the fluid density. Subscripts 0 and u represent in cell value, and unburned region value, respectively. Any additional constants that appear in (Equation 2) are specific to the turbulent flame speed model.

Unlike the SAGE model, which initiates combustion via an energy source that mimics the arc discharge to trigger reactions, standard ignition practice for the G-equation model is to initialize the corresponding G -field (i.e., the flame) at the spark plug electrode gap. In order to capture additional reactivities in both the unburned and burned regions, the SAGE model is coupled to the G-equation on both sides of the flame front. Due to the system’s low reactivity, an uncommon practice is adopted: both an energy source and a G -field passive source are initialized in the gap vicinity. This approach increases the local reactivity and ensures successful ignition at each spark plug location, consistent with observations in the optical experiments.

Both sources are spherical with a radius of 0.5 mm, activated simultaneously at the ST. Figure 10 shows the pressure and AHRR traces obtained using the G-equation combustion model, with experimental pressure traces for comparison. The combustion process is accurately captured, with small deviations observed in the numerical results, particularly at the early stages of combustion. Consequently, the model effectively characterizes the flame-flow interaction, providing a solid foundation for evaluating key aspects of the turbulence-chemistry interaction.

Figure 11 shows the visual flame development in the combustion chamber based on the iso-surface for the flame ( G = G 0 ) colored by the local D a at different stages of the combustion process. At the initial stage of combustion (CA < 1), the majority of the flame front surfaces propagate with a D a < 1 , indicating chemistry is slow and rate-determining. At this stage, the flames experience low levels of turbulence due to wall confinement, resulting in similar growth rates for all flames, as observed at CA1, where the different flame surfaces grow at a comparable rate. Once the flames surpass the spark obstacles, however, they are subjected to a wider range of turbulence level, causing D a to increase and vary in magnitude, thus affecting the flame propagation rate more significantly. At CA5, the wrinkling effect of the flame, induced by flame-turbulence interaction, accelerates flame growth, particularly for the central flame, which experiences higher turbulence and grows faster. Silva et al. (2023) reported that due to the low flame speed of NH 3 , most of the combustion occurs in the corrugated flamelets regime, where the eddy length scales are larger than the flame thickness, impacting flame growth only in a kinematic manner. The increase in D a , along with the enlarged surface area of the central flame, further enhances the propagation rate. Finally, after CA10, the different flame fronts collapse and complete the consumption of the reactant charge.

Huang and Liu (2024) suggested that the inclusion of multiple sparks does not enhance turbulence levels in the chamber; instead, the advantage of this approach lies in the ignition multi-locus adjoined propagation. This behavior is akin to advanced combustion systems like the pre-chamber, where multiple ignition locations are achieved. Unlike the pre-chamber method, however, the multi-spark approach does not enhance turbulence and mixing.