3.1. Turbulence Models
In computational fluid dynamics (CFD), turbulence modeling can generally be divided into direct numerical simulation (DNS), large eddy simulation (LES), and the Reynolds-averaged Navier-Stokes (RANS) model. Due to the challenges associated with applying DNS to real-world engineering problems, especially in the field of automotive engineering, LES and the RANS model are more frequently utilized.
Turbulence models, especially simpler ones like RANS models, often struggle to accurately predict not only flow separation and reattachment, which are critical in engine intake and exhaust processes but also turbulence length scale which is crucial for determining not only the rate of energy dissipation and mixing efficiency in engine flows but also combustion processes.
Moreover, the flow fields encountered in PCSI engines are highly complex, involving phenomena such as impinging flows, jet-like flow, flow separation, strong swirl, highly variation of turbulent length scale, and vortex shedding (refer to
Section 2). Accurately and universally modeling the all these turbulent flow characteristics remains a formidable challenge.
As shown in
Figure 4, in the analysis of 52 publications reviewed, , the majority, 52.8% studies, relied solely on the RNG k-ε turbulence model, while 11.91% studies opted for LES alone. Moreover, 7.14% studies integrated both LES and RANS models in their research. Common RANS turbulence models cited were standard k-ε, renormalization group (RNG) k-ε, shear stress transport (SST) k-ω. Sub-grid scale models for LES included Smagorinsky only.
3.1.1. RANS Turbulence Models
The RANS-type turbulence models are simple to use, computationally inexpensive, and economical [
113,
114].
Hence, given its simplicity and shorter computation time, the RANS-based model, which represents the average behavior of important parameters, is typically chosen over the LES (Large Eddy Simulation) model. However, the RANS-type k- ε turbulence models are based on Boussinesq’s isotropic eddy viscosity assumptions, and it is well known that have several problems with deteriorated prediction performance in cases of swirling flow, separation and reattached flow and flows with large rapid extra strains. Therefore, numerous variants of the k- ε model have been researched, reported, and are currently in use to overcome these challenges [
115,
116,
117]. The existence of numerous k- ε model variants is due to the ε equation with 3-4 model coefficients, which are empirically derived or determined through ad-hoc methods based on various turbulent flow patterns [*]. Additionally, turbulent flows can vary significantly in their characteristics depending on the geometry, flow velocity, pressure gradients, and other factors. To accurately capture these diverse turbulent flow patterns, modifications to the standard k-epsilon model, particularly in the ε equation, are necessary [
115,
116,
117]. Each variant attempts to address specific deficiencies of the standard model in certain flow regimes, such as free shear flows, boundary layers, or rotating flows. Consequently, a universally accurate RANS-type turbulence model has not yet been achieved. Therefore, engineers face the challenge of selecting a turbulence model that is appropriate for the specific turbulent flow characteristics of the geometry of interest.
The RNG (Re-Normalization Group) k-ε model [
118], which is the most frequently adopted in previous studies for PCSI system for various engine types [
115,
116,
117], is one of the variations of the standard k-ε model used for turbulence modeling in various kinds of CFD simulations. It incorporates additional terms that accounts for the interaction between turbulence and mean strain rate, which is not present in the standard k-ε model and modifications to the standard k-ε model to improve accuracy for certain flow conditions, particularly for swirling and highly strained flows. To achieve this goal, this model contains a strain-dependent correction term in the constant C
1ε of the production term in the RNG k-ε model’s dissipation rate(ε) equation [
115,
116,
117,
119]. The RNG k-ε model is particularly useful for simulations involving complex flow features such as swirling flows, recirculating flows, and flows with high strain rates, making it suitable for jet-like flows and complex industrial applications in which the velocity gradients are significant causing intense mixing and variations in velocity. [
54,
115,
116,
117]. In PCSI internal combustion engines, the scavenging processes in pre-chamber, along with turbulent flame jets and combustion, create regions of high strain rates. Particularly, turbulent jets issuing from nozzles of pre-chamber are typical examples of flows with high strain rates. Therefore, for these reasons, applying RNG k-ε model to PCSI engines extensively is considered practical for engine simulation, as it provides a good balance between accuracy and computational efficiency. However, noteworthy demerit of RNG k- ε model is the fact that near-wall treatment of this model can struggle with accurately predicting flows close to walls, particularly in cases involving adverse pressure gradients or separation [
54,
115,
116,
121,
122,
123]. Hence it should be careful to use RNG k- ε model if the turbulent jet interacts significantly with cylinder walls or piston head surface. Additionally, it is clear that RNG k- ε model is less accurate for detailed Structures: and does not capture the detailed eddy structures as well as LES(Large Eddy Simulation) [
115,
116,
117,
119,
120]
The k- ζ -f turbulence model [
121,
122,
123] is the second most commonly used turbulence model in PCSI engine CFD simulation after the RNG model. The k-ζ-f model is an extension of the eddy-viscosity concept and includes three transport equations for turbulence quantities, namely turbulent kinetic energy(k), Turbulence Frequency (ζ) and Dissipation Rate (f). This is selected as the turbulence model in previous studies due to its high accuracy and convergence stability. This model, optimized from Durbin’s near-wall turbulence closure model [
124,
125], is a variant of RSM(Reynolds Stress Model) turbulence models [
124,
125,
128]. The k-ζ-f model is particularly useful in complex flow simulations, including those with significant near-wall effects and flow separation [
117,
121]. It enhances the standard k-ε turbulence model by introducing the wall-normal velocity fluctuation v
2 and its source term f, which incorporate near-wall turbulence anisotropy and non-local pressure-strain effects. The careful introduction of these relaxation terms eliminates the need for damping functions. Additionally, this model improves numerical stability over the original v
2-f model by solving a transport equation for the velocity scale ratio ζ=v
2/k instead of the velocity scale v
2. Moreover, This model demonstrates superior predictive accuracy in heat transfer, surface friction at low Reynolds numbers, adverse pressure gradients, and recirculation regions compared to traditional k-ε turbulence models [
126,
127]. While it shares characteristics with low-Reynolds-number models, it uniquely eliminates the need for wall functions by being applicable near the wall. Instead, it introduces a transport mechanism for turbulent energy from the wall, effectively representing near-wall viscous damping effects through an elliptic relaxation equation [
121,
122,
123]. Hence, this model may offer superior performance in capturing the detailed flow features and interactions in not only the intake and exhaust processes but also scavenging process inside pre-chamber, crucial for predicting flame propagation, heat transfer, and emissions. Its improved near-wall treatment is advantageous for accurately predicting heat transfer inside orifices of pre-chamber and flow separation around valves and ports.
3.1.2. LES Turbulence Model
Due to the substantial turbulent energy and significant influence on momentum transfer and turbulent mixing carried by large eddies, LES methodologies offer superior accuracy compared to RANS turbulence models. [
129]
LES captures flow structures from the domain scale down to the filter scale, necessitating significant resolution of high-frequency turbulent fluctuations. This requires the use of either high-order numerical methods or fine grid resolution when employing lower-order numerical techniques. Therefore, the implementation of Large Eddy Simulation (LES) methods in automotive and mechanical engineering necessitates finely resolved grids with grid points positioned in close proximity to the boundary layers. This results in significantly heightened computational costs compared to Reynolds-Averaged Navier-Stokes (RANS) methods [ over serveral applications [
130]. Consequently, previous studies have reported only a limited number of simulations utilizing LES, with the majority being conducted on Rapid Compression Machines (RCEMs) [
102,
105] rather than on full-scale metal engines [
67].
Recently, 3D CFD analysis using LES turbulence models was performed to investigate the engine characteristics of active and passive PCSI engines fueled by natural gas for large ships under a single operating condition, both at stoichiometric and lean conditions (λ=2). Additionally, the heat release rate curves obtained from the analysis were compared with experimental results [
67].
However, to perform accurate CFD simulations of PCSI engines using the LES turbulence model, it is necessary to understand the limitations and characteristics of various sub-models.
The most important noteworthy point is that the accuracy of LES heavily relies on the Sub-grid Scale (SGS) models used to represent the unresolved scales. This is because the computational grid limits the size of eddies that can be physically represented. Despite their presence in the flow field, these eddies cannot be resolved because the CFD mesh lacks the resolution to accurately capture and depict them in CFD simulations. In LES simulation, we are particularly concerned with the eddies that are just larger than the mesh size. These eddies are too large to be broken down by molecular viscosity. Therefore, we need to find a way to model and remove these eddies from the grid. If we do not, the turbulent kinetic energy predicted in our large eddy simulation will be too high. These eddies are removed by applying an additional stress term to the Navier-Stokes equations, known as the sub-grid scale (SGS) stress. While various SGS models such as Smagorinsky model [
113], dynamic Smagorinsky model [
55,
56], WALE model [
137], and others are available, none are universally applicable to all types of flows, leading to potential inaccuracies in specific scenarios. Different SGS models all they do is provide different methods for calculating the subgrid kinematic viscosity,
.
One of the shortcomings of the Smagorinsky Subgrid Scale model which is the most frequently adopted in CFD simulation of PCSI combustion is that it contains a model constant Cs called Smagorinsky coefficient that is not universal and depends on the local flow conditions and the fraction of the cell size that gives the sub-grid length scale. Cs for homogeneous isotropic turbulence is around 0.17. This value is very crucial for accurate simulation and has come from an analytical mathematical derivation based on a homogeneous isotropic turbulence. It is noted that modern CFD codes use different values of Cs. In STAR-CD [
55], Fire [
56], Fluent [
131] Cs=0.1 and in PHOENICS [
132], Cs=0.17. However, no matter what Cs values are, this is not true in case of rotation or near wall because of too much dissipation near the wall. Additionally, in Smagorinsky model, sub-grid stress is not damped close to wall. Hence, some kind of modification to the model should be needed. Therefore, various types of sub-grid scale models have been developed. The most frequently used near-wall treatment of Smagorinsky model is the Van Driest damping function which is a damping function for proper results in wall-bounded flows [
55,
56,
131].
Figure 9 illustrates the necessary adjustments to the Smagorinsky model for accurate near-wall treatment. To address this, the sub-grid scale stress should ideally approach zero as the wall is approached. This can be achieved through various methods. One option is to implement an entirely different model, such as one based on sub-grid scale kinetic energy, replacing the traditional Smagorinsky model. Alternatively, we can modify the length scale near the wall, reducing it to zero to account for sub-grid scale eddies. Another approach is to adjust the velocity scale, making it less dependent on strain rate. The primary goal of these adjustments is to decrease the sub-grid kinematic viscosity near the wall, which will, in turn, reduce the sub-grid scale stress, effectively simulating the damping effects on the eddies.
Instead of a single user-defined constant Cs, modern CFD codes like Fire, STAR-CD implemented the Dynamic Smagorinsky Subgrid Scale model (134, 135) which computes a local time varying Cs value by test-filtering the flow field at a length scale greater than the grid length scale, which allows it to compute the correct result for wall-bounded flows without the use of damping functions.
The WALE (Wall-Adapting Local-Eddy Viscosity) Subgrid Scale model (136) is a more modern subgrid scale model that uses a novel form of the velocity gradient tensor in its formulation and widely adopted by CFD codes like STAR-CCM+ and OpenForm []. Similar to the Smagorinsky Subgrid Scale model, it suffers from the limitation that the model coefficient Cw is not universal. It is known that the WALE model is seemingly less sensitive to the value of this coefficient than the Smagorinsky model. Another advantage of the WALE model is that it does not require any form of near-wall damping—it automatically gives accurate scaling at walls [
136].
Figure 5 schematically shows various options for correcting sub-grid scale kinematic viscosity near the wall. Here, S
ij is the strain rate of the resolved eddies on the CFD mesh.
Additionally, a Coherent structure model(CSM) as a subgrid-scale model is applied [
64,
67] was adopted in Fire code [
56]. It is reported that the CSM gives good predictions and is almost the same performance the dynamic Smagorinsky model for various complex geometries [
137].
Additionally, applying appropriate boundary conditions for LES is complex. Inflow boundary conditions, in particular, need to accurately represent turbulent fluctuations, which is challenging to achieve in practice. Due to the aforementioned factors, achieving dependable LES simulations necessitates a higher level of expertise and experience compared to the requirements for utilizing the RANS model [
138].
Bolla et al. [
100,
102] executed numerical studies of RANS-LES comparison using an RCEM to analyze an automotive-sized scavenged pre-chamber, aiming to compare the two turbulence models’ ability to predict turbulence and fuel-air mixing. In this study, a Smagorinsky-type sub-grid scale model [
139] was used to compute the unresolved turbulent scales for the LES turbulence model, while the RANS turbulence model employed the time scale Bounded k−ε Turbulence Model in VECTIS, which is an enhanced version of the standard k−ε model for high strain rates or strong adverse pressure gradients. [
57]. The outcomes indicated that the RANS-based model could effectively reproduce the key ensemble-averaged flow patterns seen in LES for two pre-chamber setups. However, in RANS often falls short in accurately capturing the radial fuel-air mixing compared to LES.
3.2. Physical Phenomena and Combustion Models of PCSI Engines
The most commonly used combustion models for engine simulations are based on the flamelet combustion regime, particularly for spark ignition (S.I.) engines. In this regime, turbulence can distort and increase the surface area of the flame front while maintaining its inner laminar structure and flame speed. This condition is characterized by a Damköhler number (Da) greater than 1 and a Karlovitz number (Ka) less than 1.
Numerous previous studies [
17,
25,
29,
45,
140,
141,
142,
143,
144,
145,
146] have investigated the operating mechanisms of TJI in PCSI engines, which are schematically illustrated in the above figure. The characteristics of TJI can be categorized into enrichment, thermal effects, and chemical effects. Key heat transfer phenomena during flame propagation include thermal quenching, which occurs due to rapid heat transfer to solid surfaces as the flame passes through the nozzle, and hydrodynamic quenching, which happens when the flame mixes with the cool, lean mixture as it enters the main chamber. For a detailed explanation of TJI, please refer to the referenced literature [
45]. Recent studies [
95,
141,
142,
143,
144,
145,
146] have validated that the two-stage combustion process in the main chamber comprises a jet-dominant phase and mixing controlled, which is influenced by the combustion intensity within the prechamber, and a flame propagation phase, which depends on the reactivity of the mixture in the main chamber driven by the in-cylinder bulk flow and the associated turbulence.
As described above, capturing the flame dynamics across the pre-chamber nozzles is a significant challenge due to the complexity of flame propagation through the orifices [
147]. Most previous premixed combustion models were based on the corrugated flame zone; however, recent studies have demonstrated that PCSI engines operate at highly diluted condition, which approach the thickened flame regime. Consequently, flamelet-based combustion models, conventionally used for SI engine combustion analysis, are not suitable for application in PCSI engines operating under these ultra-lean conditions. These models fail to capture the combustion behavior of flame quenching and stretching through the orifices, potentially shortening ignition delays in the main chamber [
116].
In the studies conducted thus far on 3D CFD combustion analysis of PCSI engines, the combustion models applied have been flamelet-based premixed combustion models originally used for premixed combustion in homogeneously operated spark ignition engines. However, the PCSI engines operate under very lean conditions (λ > 2) and in high Karlovitz number (Ka > 1) regimes due to the presence of turbulent jets and turbocharging, which result in extremely high turbulence intensities. Therefore, to accurately predict the combustion characteristics of PCSI engines, a combustion model must effectively capture multiple and distributed ignition points within the main chamber, covering both premixed and partially premixed combustion regimes. Furthermore, under the typical operating conditions of TJI systems, characterized by lean mixtures and highly turbulent flow fields, the foundational assumptions of flamelet-based models may no longer be applicable [
93]. Therefore, recent studies have raised questions about the predictive performance of flamelet-based combustion models and have made significant efforts to find countermeasures. In this section, we discuss the application cases, comparative predictive performance, and limitations of combustion models that have been applied to the combustion analysis of PCSI engines. Additionally, we review various methods to enhance predictive performance and overcome the limitations of flamelet-based models.
3.2.1. Flamelet Assumption
As mentioned above, combustion models for PCSI engines are most often based on the so-called flamelet assumption [
54,
55,
56,
57,
147,
148], which has frequently been adopted for premixed combustion in homogeneously operated spark ignition engines [
148]. The combustion models based on flamelet assumption simplify the complex interactions between turbulence and chemical reactions by assuming that the flame can be represented as an ensemble of thin, locally laminar flame structures, or “flamelets,” embedded within the turbulent flow. They assume a clear separation of scales between the turbulent eddies and the flame thickness. This allows the detailed chemical reactions to be precomputed and stored in a database, which can be accessed during the simulation. The chemical reactions are solved in a laminar flame configuration under varying conditions of temperature, pressure, and mixture composition. The results are stored in flamelet libraries, which provide information about species concentrations, temperature, and reaction rates as functions of mixture fraction and scalar dissipation rate.
In summary, flamelet models strike a balance between accuracy and computational cost, making them suitable for capturing the essential features of premixed combustion in S.I. engines. However, it’s essential to recognize their limitations, especially when dealing with non-premixed or partially premixed combustion regimes. Researchers continue to refine these models and explore more detailed approaches to improve engine combustion simulations.
3.2.2. G-Equation Model
A flamelet-based combustion model, the G-equation, is one of the widely adopted combustion models for simulating the combustion processes of PCSI engines in CFD simulations within the engine modeling community. This approach utilizes a level-set method, which represents moving interfaces or boundaries on a fixed computational mesh. It is particularly useful for problems where the computational domain is divided into two regions separated by an interface. The Level Set modeling technique allows the fluid-fluid interface to move within any given velocity field [
1,
148,
149]. Detailed information on G-equation model can be found in the literature [
61,
71,
81,
86,
90,
91,
94,
100,
101,
102,
103,
104], and only the brief descriptions are provided here. In order to obtain a formulation that is consistent with the well-established use of Favre averages in premixed turbulent combustion, we split G and the velocity vector v into Favre means and fluctuations. Using a number of additional closure assumptions described in Peters [
4], one finally obtains governing equations for
and its variance
are defined [
1,
2,
3,
4]. In order to solve
equation, a model for the turbulent flame velocity must be provided. The turbulent flame speed is a key parameter and is typically modeled as a function of the laminar flame speed and turbulence characteristics such as turbulent intensity and length scales. This allows the model to incorporate the effects of turbulence on flame propagation without directly solving the detailed turbulence-chemistry interactions.
Within the G-equation context, several correlations for turbulent burning velocity from literatures [
4,
56] are presented and evaluated in the previous studies [
152]. Among the most popular correlation is Peters’s correlation [
150] which is valid for both large-scale and small-scale turbulence as:
where S
L is laminar flame speed, u
’ the fluctuating turbulent component, δ
L laminarflame thickness,
the integral length scale, b
1 and b
3 are model constants corresponding to large and small-scale turbulence enhancement, respectively, and Da Damkohler number which is a ratio between the flow time scale(
over the chemical time scale(
The laminar flame speed, S
L depends upon the local pressure, the fresh gas temperature, the local unburned fuel-air equivalence ratio using the Metghalchi and Keck correlation [
103,
107] and the chemical time arising due to the flame stretching [
56,
151,
153,
154,
155]. These common correlations for S
L are equations derived from fitting forms based on combustion experiments conducted over various temperature and pressure ranges. Therefore, outside the range of these correlations, the S
L is calculated using extrapolation methods, which inherently introduce relevant input errors into any combustion model.
Another way to get S
L is using tabulated values which was created based on the 30-species skeletal mechanism developed by Lu and Law [
91].
The laminar flame thickness is calculated from the temperature profile along the normal direction of the flame front and also from the chemical time. The chemical time is calculated from the characteristic time of the laminar flame using Zeldovich Number which depends on the activation temperature of the fuel oxidation.
3.2.3. The Extended Coherent Flame Model(ECFM)
The Extended Coherent Flame Model (ECFM) [
56,
151] builds upon the basic principles of the Coherent Flame Model (CFM) by incorporating additional features to handle more complex combustion scenarios which focused on turbulence and flame interaction and includes detailed modeling of how turbulent eddies affect flame stretch and flamelet behavior. This involves correcting for the effects of turbulence on flame stretch, considering the turbulence intensity, and adjusting for the curvature and thermal expansion of flamelets [
156]. Thus, the flame stretch is influenced by turbulence, as well as the ratios of turbulent to laminar flame velocities and lengths. It is adjusted for curvature and thermal expansion effects caused by laminar combustion in flamelets, based on the assumption of isotropic flame distribution [
56].
The model calculates the rate of fuel consumption based on the flame surface density (FSD) and the reaction rate per unit flame surface area.
In the case of the coherent flame model, flame surface area per unit volume, defined as
Using flamelet assumption, the mean turbulent reaction rate is computed as the product of the flame surface density Σ and the laminar burning velocity S
L via:
In this ωL as the mean laminar fuel consumption rate per unit surface along the flame front, ρfu,fr the partial fuel density of the fresh gas, ρ the density of the fresh gas and yfu,fr is the fuel mass fraction in the fresh gas.
When combustion starts, several new terms have to be computed. Amongst them are source terms and two quantities in order to use equation (Eq. 1): Σ and S
L.
The first term of left hand side is the time dependent component and second term is the convective transport of the FSD. The first term of right hand side is source term which represents the production of flame surface density comes essentially from the turbulent net flame stretch, the second term is sink term which represents the quenching effect referring to the local extinction or reduction of the flame surface density due to unfavorable conditions, such as excessive strain, heat loss, or insufficient reactants. Hence, the FSD transport equation incorporates these effects into a combined source term , S
∑ which includes both the creation and destruction mechanisms:
Here, Sother represents additional source terms that might be relevant depending on the specific combustion scenario.
This approach allows for detailed tracking of how turbulence affects the flame surface and, consequently, the combustion process.
Recently, the ECFM-3Z model has been extensively adopted for 3D CFD analysis of PCSI engines using AVL’s Fire CFD code [
56]. The ECFM-3Z model is is an extension of the ECFM combustion model based on FSD transport equation and mixing model that can describe inhomogeneous turbulence premixed and diffustion combustion and operates within three distinct zones or regions: fuel, air, and the air-fuel mixture. In this model, the fuel can be represented as a mixture of various components. Both the burnt and unburnt gases are categorized into these three zones.
The extent of mixing among these zones is determined using a characteristic time scale, which is computed based on the k-zeta-f turbulence model [
61,
65,
68,
69]. The ECFM-3Z model assumes that the composition of unburnt gases, including air and Exhaust Gas Recirculation (EGR), remains consistent across both mixed and unmixed zones. The properties of the burnt gases are calculated based on the reaction progress variable.
In conclusion, this combustion model focuses on flame propagation and the interaction between turbulent flow and flame, while utilizing simplified global kinetics for chemical kinetics.
3.2.4. The Multizone Well-Stirred Reactor (MZ-WSR) Model
A homogeneous reactor-type combustion model, MZ-WSR, operates on the premise that sharp gradients in temperature and density are unlikely to occur within a cell and models combustion as an ignition-based phenomenon. This model divides the reactor into several well-mixed zones, each of which is assumed to be perfectly mixed with uniform composition. Additionally, chemical reactions are assumed to occur instantaneously within each zone. The governing equation of the MZ-WSR (Multi-Zone Well-Stirred Reactor) combustion model can be represented as follows:
where, γ is the multiplier, Y
i the mass fraction of species i,
the reaction rate of species i,
the mass flow rate, and
the mass fraction of species i in the inflow.
The MZ-WSR model, which is coupled with detailed kinetic calculations, is particularly suitable for modeling the jet ignition process. Therefore, most of the studies [
83,
84,
87,
88,
89,
90,
91,
93,
94,
106] investigating the PCSI system using WSR models have considered detailed chemical mechanisms derived from GRI Mech 3.0 by Lu and Law [
158] by utilizing the integrated chemistry solver known as SAGE [
78,
88,
157], which is included in the commercial code CONVERGE [
54]. However, the use of the SAGE model for combustion involves a considerable simplification. In the simulations, turbulence models are not applied to the mean chemical production terms in the governing equations. This means the potential impacts of turbulent fluctuations on these terms are not accounted for in the simulations. It should be noted that the influence of turbulence, as modeled using the RANS k-
ɛ turbulence approach discussed earlier, is applied exclusively to the transport equations of mass, energy, and momentum in the averaged equations. Prior numerical investigation into the jet ignition process [
157] have demonstrated that the WSR assumption can lead to cooler temperatures in cells where thin flamelets are present. This issue can become more pronounced in ultra-lean mixtures. Recently, several papers [
90,
93] have evaluated the prediction performance of this model by comparing it with the G-equation model for the combustion processes in an both passive [
93] and active PCSI engine [
90] fueled by natural gas. The results of these literatures showed that MZ-WSR model predicts faster combustion rates in the main chamber than the G-equation model fails to match the pre-chamber combustion phase in the right place. However, the prediction of the combustion rate in the main chamber by this model matched well with the experimental results. The MZ-WSR model requires that each cell be treated as an individual well-stirred reactor. For this to be accurate, the characteristic time of the turbulence in the cell must be significantly smaller than the characteristic time of the combustion chemistry. In other words, the Damköhler number must be much less than 1 to assume a well-stirred reaction. Therefore, in the case of PCSI combustion, where lean and highly heterogeneous turbulent flows exist, this model can result in significant errors.
Despite its limitations, this model has been extensively used for analyzing the combustion processes of PCSI engines within the RANS framework up to the present [
26,
27]. This is because the MZ-WSR model excels in providing a detailed and chemically accurate representation of combustion processes, and allows for more flexibility in adjusting the model to account for different fuels such as natural gas [
27,
28] and gasoline [
26] and dual-fuel combustion conditions [] due to its detailed chemical kinetics. This makes it particularly useful for research and development where precise emission predictions and understanding of combustion chemistry are crucial. The experimental validation of this model was conducted using in-cylinder pressure and heat release profiles, and it demonstrated good predictive accuracy for both pressure and the combustion process within range of 1.6 < λ < 2. [
78,
88,
157].
In recent, When the G-equation is utilized, the MZ-WSR model is integrated before, during, and after the flame front to calculate the intermediate and post-reaction species instead of using simplified global kinetics [
91,
93,
104]. However, its effects have not been quantitatively proven.
3.3. Turbulence-Chemistry Interaction
The interaction between turbulence and chemistry-related quantities plays a fundamental role in determining combustion characteristics of PCSI applications, including ignition, flame propagation speed in both the prechamber and main chamber, and the burn rate trend. This is because, compared to traditional SI engines, the turbulent flow field in PCSI systems is highly inhomogeneous, has large spatial gradients at the jet boundaries, and exhibits rapid temporal evolution. Additionally, a high level of turbulent fluctuation near the spark plug, due to jet-to-jet interaction produced by the throttling effect during the scavenging process, is one of the key factors distinguishing the flame evolution characteristics of the PCSI system from conventional SI engines. To simulate these complex turbulence-chemistry interactions, previous literatures utilized the source terms of combustion models that were developed for pre-mixed or non-premixed conventional gasoline S.I engines.
As previously explained, the ECFM-3Z and G-equation combustion models are the most widely and frequently used for combustion analysis of I.C engines including PCSI system. These two models have different approaches to representing turbulence-chemistry interactions. Therefore, they calculate different laminar flame speeds, turbulence intensities, and spatiotemporal scales, resulting in different turbulent flame speeds. Turbulent flame speed significantly affects engine output, fuel consumption rate, and pollutant formation. Thus, the predictive accuracy of these two combustion models ultimately depends on the accurate calculation of turbulent flame speed 159,160]. The G-equation combustion model calculates the turbulent flame speed using explicit correlations, as shown in Equation (1), whereas the ECFM-3Z model derives it from the FSD equation. By examining the changes in velocity scales ratio(
and length scale ratio(
during the combustion process for an optically accessible GDI engine under stoichiometric conditions using the two combustion models(G-equation with three different turbulent flame speed correlations), as depicted in the Borghi-Peters diagram in
Figure 8 below, clear differences can be observed as combustion progresses [
161]. The differing results from the two combustion models can be attributed to their distinct approaches in handling the flame brush, particularly in terms of the spatial distribution of the reaction volume predicted by each model [
161]. Specifically, the differing flame brush volumes and notable deviations between two models are due to different turbulence intensities and scales, namely mean turbulence-flame interaction.
The flamelet-based combustion models undergo significant differences in the combustion process in terms of turbulence/flame interaction. However, validation results compared to experimental ensemble-averaged in-cylinder and burn rate traces do not show substantial differences. The figure below compares pressure traces on the left and heat release rate curves on the right from experiments and various combustion models under the same engine and operating conditions [
161]. From an engineering perspective, the results show overall good agreement with the experimental ensemble average traces for all models, exhibiting only minor deviations. Moreover, a noteworthy point from these results is that the results of the ECFM-3Z and the G-equation with Bradley correlation are almost identical. Consequently, it is concluded that comparing the predictive performance of various combustion models using ensemble-averaged pressure traces and burning rate curves is not appropriate. Thus, the overall burn rate only partially reflects the accuracy of the simulation framework concerning turbulent flame speed (S
T). This is because combustion models must consider a broad range of interactions between flow and chemical reactions [
159,
160,
161].
The FSD transport equation indeed includes source and sink terms that account for both the stretching and quenching effects. These terms are crucial for accurately modeling the dynamic behavior of the flame front in both laminar and turbulent combustion conditions. By incorporating these effects, the FSD equation can provide a comprehensive description of the flame surface evolution, which is essential for predicting combustion performance and stability. As shown in equation (7), the stretching effect in the FSD transport equation accounts for the deformation of the flame surface due to the flow field. This can be caused by both laminar and turbulent strain rates. In turbulent flows, this term can be more complex and might include contributions from both large-scale and small-scale turbulent eddies that stretch the flame front. To simulate stretching and quenching conditions in term S
∑, of equation(7), the intermittent turbulent net flame stretches (ITNFS) model [
159] was used. The ITNFS model concept involves using an advanced model to characterize the interaction between a single vortex and the flame front, and then extrapolating this to represent the complete turbulent flow. This model postulates that each turbulence scale influences the flame independently, implying no interaction between different turbulence scales. It is based on the idea that the overall effects of turbulent fluctuations can be predicted from the behavior of individual scales in the unburned gas mixture as shown in
Figure 9.
By applying this model to a complete turbulent flow, it is assumed that the cumulative effect of all turbulent fluctuations can be inferred from the behavior of each individual scale. The limitation of this model is that it cannot account for vortex interactions. It is clear that this model has limitations when applied to PCSI engines, where strong turbulence intensity is distributed near the flame front of the hot turbulent jet and highly inhomogeneous turbulence exists. The production of flame surface density primarily results from the net flame stretch due to turbulence. This flame stretch is expressed as the large-scale characteristic strain (ε/k), adjusted by a function C
t, which considers the size of turbulence scales, and viscous and transient effects [
160]. C
t depends on turbulence parameters and the properties of the laminar flame. Hence, the right-hand side of Equation (3) can be expressed as follows.
Kt is a very important property since it influences the source term for the flame surface and therefore the mean turbulent reaction rate. The coefficients of α known as stretching factor and β in equation(7) are arbitrary tuning constants used in ECFM.
However, it is important to better understand the contexts in which the ITNFS model can be effectively applied and where additional considerations or alternative models may be necessary, particularly due to the lack of consideration for nonlinear effects in the interaction of turbulence scales with the flame front, flame-generated turbulence, and reignition of fresh gases crossing a locally quenched flame front [*].
*159 C. Meneveau and T. Poinsot, “ Stretching and quenching of flamelets in Premixed turbulent combustion”, Combustion and Flame, Vol.86, pp.311-332, 1991.
[**] 160 Meneveau, C. and Sreenivasan, K.R. “The Multifractal Nature of Turbulent Energy Dissipation.” Journal of Fluid Mechanics 224 (1991): 429-484.
Several previous studies [
11,
12,
150,
162,
163] have reported the average turbulence/flame interaction experienced by engine flames in PCSI engines and recently developed highly downsized engines on the Borghi-Peters diagram, confirming the occurrence of multiple regimes during S.I. combustion. The results of these studies emphasize the necessity for combustion models to be predictive across all potential combustion regimes [
150,
160,
163]. Namely, the high levels of turbulence generated by the pre-chamber (PC), combined with the reduction in laminar flame speed due to dilution, significantly impact the combustion regime of PC-initiated combustion systems. These interactions are illustrated in the schematic Borghi diagram in
Figure 10, where the Karlovitz (Ka) and Damköhler (Da) numbers are used to compare relevant timescales of turbulent combustion to determine the combustion regime. Consequently, the PC-initiated jet ignition combustion regime shifts into the thin reaction zone, bringing it much closer to the stability limits.
The two figures below are based on the same passive PCSI engine with a port fuel injection system in a gasoline engine.
Figure 10(a) shows the progression of the combustion region at 2000 rpm on the Borghi-Peters diagram for λ = 1 and λ = 2 cases [
12].
Figure 10 (b) illustrates the change in flame structure during the combustion period at 4500 rpm by increasing the air-fuel ratio using air and EGR [
11]. From these results, it is clear that lean combustion shifts the combustion characteristics to the thickened flame regime. As is well known, this occurs because the eddies become smaller than the flame thickness, allowing some eddies to penetrate the pre-heat zone of the flame. Another important piece of information evident from
Figure 10 (b) is that EGR dilution exhibits higher sensitivity compared to air dilution [
11,
12].
In very recent, gasoline-fueled passive pre-chamber engine was numerically modelled using RNG k-ε turbulence model and MZ-WSR with GAGE combustion model and investigated combustion charateristics for λ=1.0 and 1.2 at 4000rpm [
83].
Figure 11 shows the evolution of the turbulent regimes in pre-chamber and main chamber for two operating conditions on the Borghi-Peters diagram. Noteworthy feature of this figure is that almost entire pre-chamber combustion, even in lean case, evolve in the thin reaction regime due to high level of turbulent turbulence intensity produced by strong jet to jet interaction during compression stroke. As shown in
Figure 11(a), during the initial stages of combustion, the velocity scale ratio decreases due to the weakening of the initial turbulence intensity. Subsequently, it recovers as the residual gas decreases. The length scale ratio shows that the integral length scale(l
l) remains constant, and the laminar flame thickness(δ
L) gradually decreases and then stabilizes. Under lean conditions, the laminar flame thickness increases and the laminar flame speed decreases, causing the curve to shift upward. From these results, it can be deduced that the MZ-WSR model is suitable for pre-chamber combustion.
Figure 11(b) shows that compared to pre-chamber combustion, main chamber combustion occupies a wider area on the Borghi–Peters diagram, starting from the border of the broken reaction zone, passing through the thin reaction zone, and moving into the wrinkled flame zone. This occurs because the high turbulence intensity is distributed across the flame front when the hot turbulent jet is ejected. As the intensity of the jet and turbulence decreases, the combustion quickly transitions through the thin reaction zone and linearly moves into the corrugated reaction zone. Therefore, combustion in the main chamber experiences three combustion regimes due to the rapid changes in velocity scale ratio and length scale ratio, indicating the presence of complex turbulence-flame interactions. This figure suggests that for PCSI engines, a multi-mode combustion model that can cover a wide range of the Borghi–Peters diagram needs to be developed. It is important to note that
Figure 10 and
Figure 11 are based on the RANS type k-epsilon turbulence model using Boussinesq’s isotropic eddy viscosity assumptions and flamelet-based combustion models, which may introduce some level of error.
3.2.5. The well-Tuned Versions of Combustion Models
As previously explained, to accurately predict the combustion in a PCSI engine, it is essential to develop a model capable of calculating multi-mode combustion. Currently, utilizing a flamelet-based combustion model is a practical alternative. Therefore, recent research has focused on tuning the model coefficients b
1, b
3 in equation (1), which represent the large and small-scale turbulence enhancements in the G-equation model described by equation (1) [
91,
93,
104]. This tuning process is performed ad hoc.
In the case of combustion models using Flame Surface Density (FSD) or MZ-WSR, efforts have been made to adjust the flame stretch factor,
α in equation (7) for FSD and multiplier, γ in equation (6) ad hoc to align with experimentally obtained in-cylinder pressure traces and heat release rate profiles [
164]. However, these treatments could not be adequate to achieve a precise correlation between simulations and experimental results in both the pre-chamber and the main combustion chamber [
93,
104,
164]. Such studies, as described above, attempt to address the limitations of describing the turbulence-flame interaction using RANS-type turbulence models and flamelet-based combustion models by calibrating the model coefficients included in the equations that represent combustion speed or burn rate.
Kim et al. [
93] tuned the multiplier of the MZ-WSR model within the range of 1.0 to 1.4. Additionally, for the G-equation model, they adjusted the coefficients b
1 and b
3 between 1.0 and 2.5. for two air-fuel ratios at 1200rpm. The engine used in this study is a 1.86-liter passive PCSI PFI engine fueled by natural gas. They validated the improvement in prediction accuracy by comparing the in-cylinder pressure traces and heat release rate profiles with experimental data.
Figure 5 presents a comparison between the cylinder pressure and heat release rate obtained from experimental data and the simulation results using the tuned models. The experimental data includes results from 300 cycles (represented by light gray lines) and the averaged pressure (depicted by a black bold line).
Figure 12 and 13 compares the simulation results without tuning the model coefficients to the corresponding experimental results. From the results, it is evident that without tuning the model coefficients, it is impossible to accurately predict the combustion characteristics of the PCSI engine. Additionally, the discrepancies are more pronounced in the lean region, likely due to the inability to account for the transition to the thin reaction zone as the flame thickness increases and the laminar burning velocity decreases.
The tuning of model coefficients significantly improves the predictive accuracy of the G-equation model. However, the MZ-WSR model fails to improve predictive accuracy with coefficient tuning, as it does not account for the effects of small-scale turbulence on the reaction rates.
Silver M et al. [
104] conducted a detailed analysis on a 2.1-liter active pre-chamber engine using natural gas as fuel, focusing on the impact of turbulent jets ejected from the pre-chamber on the burn rate in the main chamber. This study utilized the G-equation combustion model and RANS-type turbulence model for two different orifice diameters. In this work, the turbulent flame speed equation was tuned with coefficients b
1 and b
3 set to 0.78 and 2.0, respectively. The figure below compares the pressure traces in the pre-chamber and main chamber at 1200 rpm for the two orifice diameters between the calculated and experimental results. As observed, there is a excellent degree of agreement between the experimental and predicted data.
In this computation, a highly sophisticated 0D/1D wave model [**] was employed to obtain initial conditions and temperatures at the walls, as well as temperature and pressure values at the inflow and outflow boundaries for the 3D-CFD analysis. However, this model requires input conditions such as lift curves of intake and exhaust valves, pressure pulsations in the intake and exhaust pipes, and fuel lines. Additionally, experimentally obtained pressure traces are needed to analyze combustion phenomena in the pre-chamber and main chamber.
Therefore, without preceding experiments, it is challenging to perform precise 3D CFD analysis of the PCSI engine. Moreover, the tuning of model constants in the combustion model based on experimental values is necessary, which varies depending on the engine displacement, fuel used, pre-chamber geometry, and fuel supply system.
To the best of the authors’ knowledge, there have been only a few studies utilizing well-calibrated combustion models [
46,
62]. Although the physical phenomena of TJI combustion are not yet fully understood, model constants in the correlation or source term for various combustion models adjusted to correct the laminar flame speed and align with experimental pressure traces [
41,
48,
49]. However, the ability to accurately predict all phenomena involved in TJI combustion remains uncertain [
46].
Figure 15 shows the recent share of (a)combustion modelsused and fuels for PCSI simulations. The majority, 45.9% studies, relied on the MZ-WSR model, while 32.4% studies opted for G-equation model. Subsequently, the frequency of usage for ECFM(3Z) model accounted for 16.2%, while less frequently used combustion models, such as the Weller model, were grouped as ‘etc.’ Additionally, the fuels used in the CFD simulation for the PCSI engine were categorized and shown in
Figure 15(b). As illustrated in the figure, 50% of PCSI engine studies utilized natural gas fuel, followed by gasoline fuel. Recently, research has emerged on using carbon-free fuels such as ammonia and hydrogen in PCSI engines. These studies employ a dual fueling strategy, utilizing natural gas and diesel fuel for initial stage of combustion, and have been grouped under dual fuel.