Modeling of Hydrogen Combustion from a 0D/1D Analysis to Complete 3D-CFD Engine Simulations

Abstract

Hydrogen and its unique properties pose major challenges to the development of innovative combustion engines, while it represents a viable alternative when it is based on renewable energy sources. The present paper deals with the holistic approach of hydrogen combustion modeling from a 0D/1D reactor evaluation with Cantera up to complete engine simulations in the 3D-CFD tool QuickSim. The obtained results are referenced to the current literature and calibrated with experimental data. In particular, the engine simulations are validated against measurements of a single-cylinder research engine, which was specifically adapted for lean hydrogen operation and equipped with port fuel injection and a passive pre-chamber system. Special attention is hereby given to the influence of different engine loads and varying lambda operation. The focus of this work is the complementary numerical investigation of the hydrogen flame speed and its self-ignition resistance under the consideration of various reaction mechanisms. A detailed transfer from laminar propagation under laboratory conditions to turbulent flame development within the single-cylinder engine is hereby carried out. It is found that the relatively simple reaction kinetics of hydrogen can lead to acceptable results for all mechanisms, but there are particular effects with regard to the engine behavior. The laminar flame speed and induction time vary greatly with the inner cylinder conditions and significantly affect the entire engine’s operation. The 3D-CFD environment offers the opportunity to analyze the interactions between mixture formation and combustion progress, which are indispensable to evaluate advanced operating strategies and optimize the performance and efficiency, as well as the reliability, of the engine.

1. Introduction

A reduction in greenhouse gas requires a comprehensive approach that not only lowers the carbon footprints of new vehicles but also addresses emissions from the existing fleet if the objective of a 45% reduction in the global CO₂ emissions by 2030 (relative to 2010 levels) is to be achieved [1]. The scientific community underscores the importance of a technology-neutral approach to accelerate greenhouse gas (GHG) reduction, recognizing that relying solely on electric vehicles is not sufficient to achieve rapid progress [2]. A range of technologies must be adopted to ensure a faster and more effective transition, enabling significant emission cuts across all vehicle types and promoting the sustainable decarbonization of the transport sector.

Operating spark ignition (SI) engines with hydrogen is a key area of research in the development of sustainable combustion technologies. In order to adapt to the different physical and chemical properties of hydrogen compared to common gasoline, considerable opportunities and challenges must be addressed.

Table 1. Fuel properties of gasoline RON98 and hydrogen [3,4,5].

Property	Gasoline RON98	Hydrogen
Density [kg/m3] 1,2	750–770	0.09
Lower heating value [MJ/kg] 1	41.4	120
Energy density [MJ/m3] 1,2	31.7 × 103	10.7
Stoich. air-to-fuel ratio [-]	14.7	34.3
Caloric heating value (PFI) [MJ/m3] 1,3	3.76	3.19
Flammability limit [λ]	0.4–1.4	0.13–10
Minimum ignition energy [mJ]	0.24	0.02
Flame quenching distance [mm]	2	0.64

¹ at 273 K, ² at 1 bar, ³ at λ = 1.

offers a concise comparison of the fundamental properties of hydrogen and gasoline, especially considering parameters relevant to the performance of internal combustion engines (ICEs).

Hydrogen presents a much higher lower heating value (LHV) than gasoline considering its mass-based properties, but, due to its low density, it also has a 3000 times lower energy density in the standard state. As is generally known, liquid hydrocarbon fuels are nearly unmatched in terms of energy density at standard conditions, which is one of the main reasons for their wide range of applications in the transport sector. In addition to their low density, the higher stoichiometric air-to-fuel ratio leads to a low caloric heating value, which generates high volumetric efficiency losses once hydrogen is injected into the ports (PFI) of an engine [6,7,8]. For these reasons, a lack of power output can be expected compared to gasoline engines if no further actions are considered.

On the other hand, hydrogen presents favorable properties for a high-efficiency combustion process in ICEs. It offers wide flammability limits that enable an extended engine operation range with very lean combustion. Lean operation offers a number of advantages, such as reducing the throttle losses at lower loads and increasing the engine efficiency at high loads by minimizing wall heat losses. At the same time, due to high flame speeds, the combustion process occurs comparably quickly, leading to high knock resistance if only the chemical characteristics of hydrogen are considered. Indeed, because of its lower ignition energy (see Table 1), and if there exist any spots in the combustion chamber or even in the manifolds of a non-homogeneous mixture, with conditions prone to auto-ignition, hydrogen tends to generate abnormal combustion events, such as backfire and knock. With the purpose of reducing combustion anomalies, the engine can be run leanly and with the replacement of a PFI system with direct injection (DI), for example, to prevent backfire issues [3,9]. As a consequence, new challenges for the turbocharging system must be addressed to run the engine with sufficient air excess to avoid combustion anomalies. Packaging constraints may also prevent the implementation of higher excess air operation strategies, as the vehicle would need to be fitted with a suitably sized intercooler to regulate the air temperature before it enters the engine.

In hydrogen engines, a new balance between the key thermodynamic parameters must be established, including the injection strategy, air–fuel mixture, and combustion process. This is essential to harness the benefits of hydrogen while mitigating the drawbacks posed by this novel fuel.

Along with experiments, 3D-CFD simulations represent a powerful and more cost-effective solution for the optimization of the new ICE operating with hydrogen. Fuel modeling in 3D-CFD simulations for ICEs considers both the thermodynamic and chemical properties of the fuel. Key fuel-related parameters, such as the laminar flame speed (LFS) and ignition delay time (IDT), can be calculated using detailed chemical mechanisms, which provides the most precise way to evaluate these properties across the temperature and pressure ranges typical of ICE combustion [10,11,12]. For instance, while common fuels like gasoline or methane exhibit similar trends in LFS, hydrogen, at a stoichiometric air-to-fuel ratio (λ = 1), has a flame speed that is nearly ten times higher. Even under very lean conditions, hydrogen maintains a sufficiently high LFS. Accurately replicating hydrogen’s LFS is crucial for successful hydrogen combustion modeling in 3D-CFD codes. As mentioned, while hydrogen enables rapid combustion, its susceptibility to self-ignition also needs careful evaluation. The knock resistance, for example, can be assessed through detailed chemical mechanisms by calculating the IDT—the time required for the air–fuel mixture to self-ignite under a specific temperature, pressure, air-to-fuel ratio, and residual gas concentration. These calculations assume a perfect mixture of fuel and air. However, hydrogen’s low ignition energy means that even local ignition conditions can cause the fuel mixture to self-ignite, leading to backfire or pre-ignition events. Therefore, parameters like the LFS and IDT must be evaluated while considering the actual flow field within the engine. In this study, 3D-CFD simulations were conducted using the 3D-CFD tool QuickSim, which integrates both the chemical properties of hydrogen and engine thermodynamics into a single platform for comprehensive analysis [13].

The aim of this work is to present a methodology for the modeling of hydrogen’s thermodynamics and chemistry, specifically designed for the development of internal combustion engines. This approach strives to achieve a balance between accuracy—capturing key engine parameters like mixture formation and combustion—and computational efficiency, ensuring that the tool can be used effectively by the industry for the virtual development of high-efficiency engines within reasonable calculation times.

2. 0D/1D Combustion Analysis

In order to determine the detailed properties of hydrogen combustion, its reaction kinetics are analyzed first. Since the direct link between the detailed chemistry and the 3D-CFD simulation poses major challenges, investigations are initially carried out in a 0D/1D environment. Thanks to the simplified boundary conditions and reduced computational time, all important conditions can be calculated in advance and saved in look-up tables. With this practical approach, the required properties can then be obtained from these tables during the 3D-CFD simulation and do not have to be calculated parallelly.

2.1. Simulation Setup with Cantera

The freely available software Cantera (V3.0.0) was selected as the reaction kinetics solver for this work [14]. To cover all of the specific areas needed to describe the flame phenomena inside internal combustion engines, Cantera’s reactor models are used in combination with adapted coding in a Python environment. To determine all of the necessary combustion properties for H₂-ICE applications, three different reactor types are used in this work:

• Non-adiabatic 0D reactor—combustion products;
• Adiabatic 1D reactor—laminar flame speed (LFS);
• Adiabatic 0D reactor—ignition delay time (IDT).

Furthermore, Cantera makes it possible to use a large number of freely accessible reaction mechanisms, as the conversion of many file formats is supported. In

Table 2. Used reaction mechanisms and their properties.

Mechanism	# Species	# Reactions	Species
LLNL [12]	10	40	H, H2, O, O2, OH, H2O, N2, HO2, H2O2, AR
Polimi [15]	14	33	H, H2, O, O2, OH, H2O, N2, HO2, H2O2, AR, CO, CO2, HE, HCO
Keromnes [16]	15	48	H, H2, O, O2, OH, OH*, H2O, N2, HO2, H2O2, AR, CO, CO2, HE, HCO
Burke [17]	13	27	H, H2, O, O2, OH, H2O, N2, HO2, H2O2, AR, CO, CO2, HE
SanDiego [18]	9	21	H, H2, O, O2, OH, H2O, N2, HO2, H2O2

# Number, OH* chemiluminescence.

, the reaction mechanisms used in this work are listed. Each is specifically reduced for hydrogen reactions and excludes long hydrocarbon and nitrogen compounds. Only the mechanisms of Polimi, Keromnes, and Burke are suitable for CO oxidation, but this is not considered within this investigation (see Section 2.2). With less than 16 species, these compact mechanisms enable very short computational times compared to other fuel combustion calculations.

2.2. Mixture Composition Including Combustion Products

For the calculation of the LFS and IDT, the inputs of the reactors must first be determined. The mixture composition hereby plays a crucial role, because the wide flammability range of hydrogen–air mixtures causes great differences in the species composition. In the course of this work, the reactor inputs were reduced to the four main components of the mixture: H₂, O₂, H₂O and N₂ (see Figure 1). Reactions involving hydrocarbons and reactions of other nitrogen compounds were neglected, as they only play a minor role in the combustion.

For internal combustion engines, the residual gas concentration of the mixture must also be considered. For this reason, simulations were conducted with a 0D reactor to determine the relevant combustion products and their proportions. Within a homogeneous reactor volume, the state variables are a function of time and it is assumed that a thermodynamic (but not chemical) equilibrium is sought at all times while the hydrogen conversion is taking place. With regard to the residual gas composition in engines, the temperature is a very decisive factor. Therefore, the 0D reactor was not implemented in adiabatic conditions but with a defined wall heat loss to simulate a temperature gradient after combustion. The investigations showed that reactive combustion products such as OH, HO₂, and H₂O₂ were only stable at temperatures above 1500 K and progressively decayed for decreasing temperatures. By comparing the calculated temperature after combustion with that of a typical internal combustion engine during the expansion phase, it can be concluded that the burned mixture (becoming residual gas in the following engine cycle) most likely cools down to temperatures far below 1000 K before it takes part in the combustion in the following engine cycle. For this reason, it is valid to reduce the residual gas mixture to the four relevant main species, comparable to the fresh gas mixture composition, namely H₂, O₂, H₂O and N₂. Figure 1 shows an evaluation of the fresh mixture compositions at various air-to-fuel ratios and increasing residual gas concentrations.

2.3. Laminar Flame Speed

For the calculation of the LFS, an adapted 1D reactor in Cantera can be used to solve the governing conservation equations for the overall mass, species, and energy [14]. With regard to the validation of the approximated flame values, it is necessary to consider experimental data as a reference. Unfortunately, measurements of the LFS are only possible under limited boundary conditions, as it is necessary to keep the flame under laminar conditions without the onset of turbulent flame wrinkling for a correct evaluation. Optical measurements (e.g., schlieren) of spherical flame development are well suited for this purpose, as Beeckmann et al. [19] and Krejci et al. [20] have shown. Figure 2 (left) shows the results of the LFS measurements in comparison with simulated data with different reaction mechanisms. Considering the ambient pressure and temperature, all mechanisms show good agreement with the experimental data, as they have all been validated against measurements at such conditions. As no experimental data are available for ICE-relevant conditions (high pressures and temperatures), the properties of the mechanisms must be extrapolated for these boundaries. Figure 2 (right) shows the calculated LFS for relevant engine operations close to the ignition point during the engine compression phase and for different mechanisms. A variation in the input parameters was implemented at 60 bar, 800 K and a lambda of 2.0 for all considered reaction mechanisms. The differences among the mechanisms in the calculation of the LFS increase under these conditions. While the mechanism of Polimi predicts lower values, the results with Burke are increased by a considerable amount and are thus the highest in this comparison.

2.4. Ignition Delay Time

The mixture’s IDT was analyzed with a 0D constant-volume reactor. Within this volume, all state variables are homogeneously distributed and the states are a function of time, i.e., a generally non-stationary system. It is assumed that a thermodynamic equilibrium prevails in the entire reactor at all times. In this case, the reactor is also assumed to be isolated and therefore no heat or mass transfer with the environment occurs.

Although the correct measurement of the induction time under laboratory conditions is already possible in a wider area, it is not possible to evaluate it in real ICE conditions here either. The measurement procedure can be carried out either in a shock tube [20,21] or in a rapid compression machine (RCM), as also considered in the work of Villenave et al. [22] and Kéromnès et al. [16]. Figure 3 (left) shows the results of their measurements in comparison with simulated data of the different reaction mechanisms. As the induction times are very sensitive to the temperature, a logarithmic scale was chosen in the illustration. Considering the value of 32 bar and a lambda of 2.0, all mechanisms show good agreement with the experimental data, especially for higher temperatures and lower IDTs. For decreasing temperatures, a deviation is notable between the experimental and simulated data. The reason for this could be the increasing amount of wall heat loss during the experiments, which is not taken into account under the adiabatic conditions of the 0D reactor. Therefore, the timings for lower temperatures in the experiment become steadily higher, while a flatter gradient can be observed in the simulation. As no experimental data are available for more realistic ICE conditions, the properties of the mechanisms must be extrapolated for these areas as well. Figure 3 (right) shows the results for relevant engine operations close to the ignition point during an engine’s compression phase. A variation in the input parameters was implemented again at 60 bar, 800 K, and a lambda of 2.0 for all considered reaction mechanisms. A significant deviation between the results can be observed. While, with Keromnes, the IDT is the lowest and more critical in terms of self-ignition, the results with SanDiego are increased by a considerable amount and are thus the most conservative values in this comparison.

2.5. Residual Gas Sensitivity

The influence of various boundary conditions on the LFS and IDT was discussed in the previous sections. However, it is also necessary to consider the influence of the burned gas concentration, which is highly relevant for internal combustion engine applications. As described in Section 2.2, the residual gas that has already been combusted is added now to the fresh mixture. To better describe the reactivity of the residual gas, a fictitious air-to-fuel ratio λ_RG is introduced. This parameter is similar to the λ_fresh of the fresh hydrogen–air mixture and corresponds to its ratio before combustion. In order to model realistic engine conditions as far as possible, the λ_RG is always equal to that of the fresh mixture in this consideration. An increasing residual gas concentration leads to a dilution effect, which is clearly reflected in the characteristic values of the combustion process. Figure 4 shows the simulation results with increasing residual gas concentrations for all considered reaction mechanisms at 60 bar and 800 K over different lambdas. Although all five cases have different absolute values for the LFS (see also Figure 2), they all show a proportional reduction in the LFS for increasing residual gas. This reveals, for example, that, due to the dilution effect, a fresh mixture at lambda 2.5 has approximately the same LFS as a mixture at lambda 1.5 with 30% residual gas, as both exhibit a comparable proportion of H₂ molecules. The deviations in the IDT can be seen more clearly between the reaction mechanisms. Not only do each of the five cases show a different gradient over the air-to-fuel ratio, but increasing residual gas is also evident in different intensities. All mechanisms except one show a clear increase in the IDT with higher dilution; only the case with Keromnes is less sensitive to the residual gas concentration, and the IDT remains almost constant.

3. 3D-CFD Virtual Engine Environment

3.1. Simulation Methodology

The use of alternative fuels such as hydrogen, which have completely different properties compared to conventional liquid fuels, requires a fundamental new development for ICEs. Key challenges in the engine development process are the optimization of charge motion, mixture formation and combustion. Especially for fuels that have not been widely investigated yet, detailed analyses and a profound understanding of these complex processes are necessary. Complementarily to test bench investigations, 3D-CFD simulations offer almost unlimited potential for investigations. For this work, the 3D-CFD software QuickSim was selected, which was developed at the FKFS and IFS University Stuttgart. A more detailed description of the used model approaches can be found in [13,23,24,25] The following section provides a brief summary of the sub-models relevant to the present work.

Specifically designed for the virtual development of ICEs, QuickSim combines the commercial CFD solver of Star-CD with self-developed models for the injection, fuel, and combustion. The robust numerical framework is set up with a RANS approach in combination with a standard κ-ε turbulence model. Second-order spatial discretization (MARS) is selected with an implicit Euler temporal discretization scheme and for the pressure–velocity coupling, the Pressure-Implicit with Splitting of Operators (PISO) algorithm is used. Utilizing an innovative meshing and mesh motion methodology that does not necessitate remeshing during the complete engine cycle avoids mapping errors and leads to a significant reduction in computational time compared to other traditional 3D-CFD approaches. The combustion chamber is meshed with structured hexahedron cells, which are detailed enough to resolve the relevant flow effects so as to reduce numerical instabilities and enable larger time steps without compromising the accuracy [13,25]. By extending the simulation domain to the real test bench’s dimensions, and by calculating multiple consecutive operating cycles, the dependency on the initial and boundary conditions of the simulation is minimized. Figure 5 shows the whole simulation domain, which includes all of the peripheries from the real test bench to correctly reproduce its behavior in the simulation.

Since the mixture formation directly affects the investigation of the combustion process in ICEs, a reliable 3D-CFD injection model is essential. For gaseous fuels in general, there are mainly two different approaches applied to simulate injection: Either a detailed injector model, where the inner flow of the injector is fully resolved, or a more pragmatic and simplified Lagrangian approach can be considered [24]. Considering a PFI engine setup with multiple low-pressure injectors in the intake manifold, resolving the flow of the injectors would increase the computational effort significantly. For a more efficient development process, the macroscopic behavior, such as the propagation of the H₂ jet, the total penetration, and the displacement of air, must be prioritized. Therefore, a simplified Lagrangian injection approach is preferred for this simulation setup. Through this approach, the total number of cells can be reduced from many millions to only approximately 350,000. Due to these methodological simplifications, a good model calibration is required, for which optical measurements inside a spray chamber serve as a reference [24].

The accurate description of the combustion process ensures reliable simulations of multiple consecutive engine cycles. In the 3D-CFD tool QuickSim, the caloric properties, as well as the LFS and the IDT, are calculated in advance for all possible temperature, pressure, and mixture composition conditions by using a detailed reaction mechanism in Cantera and are stored in comprehensive look-up tables [23]. This approach avoids the necessity of performing a detailed chemical calculation for each cell in every simulated engine cycle. In reality, for most technical applications, there is not laminar but turbulent combustion, especially in ICEs. Turbulent eddies wrinkle the flame surface and accelerate the flame’s propagation due to the higher convection of reactants towards the flame [27]. The increase from the laminar flame speed S_L to the turbulent flame speed S_T is described by the wrinkling factor K_wrink, which accounts also for the increase in the effective surface of the wrinkled flame. In the context of this study, a semi-empirical formulation of K_wrink is used, which was originally developed by Herweg and Maly [28]:

$$K_{wrink}={\displaystyle \frac{S_T}{S_L}}=1+{\left({\displaystyle \frac{{\overline{v}}_T}{{\overline{v}}_T+S_L}}\right)}^{{\displaystyle \frac{1}{2}}}·{ \left(1-\mathit{exp}\left(-{\displaystyle \frac{r_K}{l_l}}\right)\right)}^{{\displaystyle \frac{1}{2}}}{· \left[1-\mathit{exp}\left(-{\displaystyle \frac{{\overline{v}}_T+S_L}{l_l}}\Delta t\right)\right]}^{{\displaystyle \frac{1}{2}}}· {\left({\displaystyle \frac{u^{\prime }}{S_L}}\right)}^{{\displaystyle \frac{5}{6}}}$$

(1)

This equation consists of multiple numerical inputs, where r_K is the kernel radius of the flame, l_l is the integral length scale of the flow field, and Δt is the elapsed time since the ignition point. The component ῡ_T represents the velocity relevant to the wrinkling factor, which, in this approach, is calculated via the combination of the main fluid’s kinetic energy and its turbulent intensity u’. The resulting term is hereby described as

$${\overline{v}}_T={ \left[{\left|\overline{v}\right|}^2+{\left(k_{turb}·u^{\prime }\right)}^2\right]}^{{\displaystyle \frac{1}{2}}} with u^{\prime }=\sqrt{{\displaystyle \frac{2 · k}{3}}}$$

(2)

By assuming u’ as an isotropic turbulent field, it is derived from the turbulence k. The factor k_turb is additionally introduced to consider numerical and fuel-specific influences on the conversion of turbulence effects to scale the turbulence intensity. While the chemical reactions for hydrogen are less complex than for fossil fuels, a wide range of air-to-fuel ratios has to be taken into account. Especially for increased lambdas, the laminar flame speed S_L rapidly declines, and, according to Equation (1), this results in the anti-proportional behavior of K_wrink with respect to S_L. This increased flame wrinkling phenomenon for leaner mixtures has also been observed in multiple studies, such as [19,20,29,30,31].

According to the model of one-step flame propagation proposed by Weller [32], and by introducing a flame progress variable, the position of the flame can be identified in the 3D-CFD mesh. This variable can take values between 0 and 1 and defines the extent to which a specific numerical cell has already been burnt. When using the approach of Weller in the context of a RANS simulation in combination with coarse meshes and relatively large time steps, an empirical formulation for flame wrinkling is required. The equation presented by Herweg and Maly [28] is a relatively simple semi-empirical formulation for flame wrinkling and attempts to account for the microscopic flame phenomena. A more comprehensive description of Equation (1), which explains in detail which term refers to each stage of flame progression, can be found in [13,28].

The 3D-CFD tool QuickSim provides additionally a 3D auto-ignition model, which was specifically developed for the reaction kinetics approach. With the help of the detailed look-up tables for the ignition delay times of the mixture, it is possible to locally predict abnormal combustion events such as knocking, pre-ignition, or backfiring. The most common approach for self-ignition models is based on the evaluation of the pre-reaction state of the mixture. The formulation used in this work was originally proposed by Livengood and Wu [33]:

$$1={\int }_{t=0}^{t=t_e}{\displaystyle \frac{1}{\tau }}dt$$

(3)

For this formulation, τ is the individual ignition delay timing of the mixture in its current conditions inside the 3D-CFD cell. The integral length t is the elapsed time, and t_e is the time at the end of integration. With an additional transport equation for this knock integral, the temporal changes in the conditions and the mixture properties are considered as well.

3.2. Engine Setup

The single-cylinder engine used to conduct these investigations was originally designed for lean gasoline combustion with a passive pre-chamber (PC) ignition system. In particular, the measurements were performed at Fraunhofer ICT and used for the validation of the simulations carried out by FKFS. The simulation model, shown in Figure 5, also illustrates the real test bench setup and the dimensions of the intake and exhaust system. For hydrogen injection, 8 PFI injectors are circularly mounted at the intake manifold with a 145 mm distance to the intake valves; see also Figure 5. Each injector can deliver approximately 0.7 g/s of hydrogen at 6 bar injection pressure. The used passive pre-chamber system was originally designed for gasoline applications and was not specifically customized for hydrogen. The PC (M10 thread) is located in a central position in the cylinder head, which is advantageous for the scavenging of the PC and ensures the uniform distance of the jets from the PC to the cylinder walls.

To allow lean operating conditions, a unique low-tumble design is combined with a piston crown, creating a highly turbulent in-cylinder flow close to the firing top dead center (FTDC). More information can be found in [26,34]. An overview of the single-cylinder engine’s specifications is given in

Table 3. Fraunhofer single-cylinder engine specifications.

Bore	[mm]	82.5
Stroke	[mm]	80
Displacement volume	[ccm]	430
Compression ratio	[-]	12.2:1
Injection pressure (H2 PFI)	[bar abs.]	2–8
Max. cylinder pressure	[bar abs.]	180

.

3.3. Operating Points

The 3D-CFD investigations of the combustion process were performed for six operating points (OP) at 3000 RPM, among which three operating points are part of a load variation from 10 to 23 bar IMEP, while the other three operating points represent an air-to-fuel ratio variation for lambda values between 1.5 and 2.5. The intake air temperature T2 was kept constant at 40 °C, and, since there was no exhaust turbocharger within this test bench setup, both the intake pressure p2 and the exhaust pressure p3 were set to the same value. These settings were necessary to ensure that the resulting residual gas content in the cylinder and the pre-chamber was accurately reproduced, as they have a significant influence on the knocking behavior of the real engine. Moreover, a 50% mass fraction burned (MFB50) of 8 °CA after FTDC was aimed across all operating points, and the end of injection (EOI) was kept constant at −280 °CA.

For the load variation, the lambda values were kept constant at around 2.0. This means that the air and fuel demand increased for higher loads and had to be compensated for by a higher boost pressure and a greater amount of injected fuel.

Table 4. Comparison of engine load variation at 3000 RPM between test bench and simulation.

		10 bar IMEP		15 bar IMEP		23 bar IMEP
		TB	Sim	TB	Sim	TB	Sim
Lambda	-	2.00	2.06	2.06	2.07	2.05	2.08
Boost press. (p2)	bar	1.76	1.78	2.58	2.59	3.86	3.84
Air cons.	kg/h	51.8	55.9	80.4	82.8	118.9	121.5
Fuel cons.	kg/h	0.79	0.79	1.13	1.16	1.69	1.70
MFB50%	°CA a. TDC	8.2	7.8	7.9	7.3	8.0	6.8

shows a comparison of the relevant engine values for the test bench and the converged 3D-CFD simulation after multiple consecutive simulated engine cycles.

The lambda variation was implemented at a constant engine load of 15 bar IMEP. Higher boost pressures provide more air mass to increase the lambda values, whereas the fuel mass remains more or less constant (see

Table 5. Comparison of lambda variation at 15 bar IMEP and 3000 RPM between test bench and simulation.

		λ = 1.5		λ = 2.0		λ = 2.5
		TB	Sim	TB	Sim	TB	Sim
Boost press. (p2)	bar	2.21	2.22	2.58	2.59	2.99	3.02
Air cons.	kg/h	64.0	66.3	80.4	82.8	97.1	100.4
Fuel cons.	kg/h	1.20	1.23	1.13	1.16	1.11	1.13
MFB50%	°CA a. TDC	10.6	10.2	7.9	7.3	7.9	6.4
Ind. efficiency	%	40.1	40.5	42.6	42.5	43.4	43.0

). Differences in the engine’s indicated efficiency are the reason for the slightly higher fuel consumption at lower lambda values, because, typically, the wall heat losses are higher here [26]. Moreover, in the case of lambda 1.5, the MFB50 is slightly delayed due to knocking limitations.

The calibrated simulations show good agreement for the different operating points. The geometries of the peripheral components play a decisive role in a single-cylinder engine in particular, as the acoustics inside the intake and exhaust system can lead to a significant difference during engine operation. The different boundary conditions for each load point and lambda affect the behavior of these pressure oscillations and the entire gas exchange process. For this reason, not only the cycle-averaged values are important for the simulation calibration, but also the comparison of the boost pressure, back pressure, and in-cylinder pressure with the indicated pressure curves from the test bench. The accurate reproduction of the geometric properties of the test bench periphery and the calibrated injection model enables the simulation to reproduce the same pressure curves as on the test bench. A more detailed comparison of the calibration process is shown in [26,34]. The analysis of the in-cylinder pressure and the overall combustion characteristics for different reaction mechanisms is given in the following section.

4. 3D-CFD Virtual Engine Results

After comparing different reaction mechanisms in a 0D/1D model, as described in Section 2, the resulting look-up tables for the LFS and the IDT were applied to the 3D-CFD virtual engine test bench. For this analysis, only the three mechanisms of LLNL, Polimi, and Keromnes (see Table 2) were considered. For the following comparisons, the exact same engine cycle for each engine operating point was used to ensure identical conditions in the flow field before the ignition timing. Only the look-up tables were exchanged to clearly highlight the differences in the combustion.

4.1. Engine Load Variation

In order to compare the combustion behavior over a wider engine operating range, while the inner cylinder flow conditions and the mixture formation were as equivalent as possible, a variation in the engine’s IMEP at 3000 RPM with a constant lambda of 2.0 was implemented (see Section 3.2).

Table 6. Simulated conditions for the load variation at 3000 RPM and a lambda of 2.0.

IMEP	[bar]	10	15	23
Ignition point (IP)	[°CA a. FTDC]	−10	−12	−14
Temperature in cyl. @ −20 °CA	[K]	729	735	754
Pressure in cyl. @ −20 °CA	[-]	28.2	41.7	62.5
TKE in cyl. @ −20 °CA	[m2/s2]	30.4	29.9	31.0
Residual gas in cyl. @ IP	[%]	5.7	5.2	4.7
TKE in PC @ −20 °CA	[m2/s2]	51.0	51.1	51.3
Lambda in PC @ IP	[-]	2.04	2.09	2.07
Residual gas in PC @ IP	[%]	9.6	9.2	9.1

shows the relevant conditions inside the cylinder and the pre-chamber (PC) shortly before the combustion started. It shows that, with the developing IMEP, the pressure and temperature increase, whereas the mixture properties are slightly improved with a decrease in residual gas of 1%-pt. The local differences between the cylinder and PC are more distinctive for the three cases. Meanwhile, inside the PC, the residual gas concentrations and the turbulent kinetic energy (TKE) of the fluid are almost doubled, while the air-to-fuel ratio is slightly leaner.

These different conditions of the mixture result in the LFS changing individually for each load point, and different values are imported from the given look-up tables. The local turbulence field, i.e., TKE, plays a crucial role in the overall combustion progress. Depending on the flow conditions, the wrinkling factor K_wrink (see Section 3.1) determines the turbulent flame speed, which is responsible for the global flame progress and the resulting heat release rate (HRR). Figure 6 presents a more detailed depiction of the occurring combustion, as it displays the crank angle resolved pressure and HRR trace inside the cylinder, as well as the LFS and K_wrink as averaged values in the flame front for each load point and reaction mechanism.

It should be mentioned here that the turbulence intensity factor k_turb (see Section 3.1) was adjusted individually for each reaction mechanism so that an equivalent HRR was achieved and the test bench’s cylinder pressure could be accurately reproduced in the 3D-CFD simulations. To realize this, the factor k_turb had to be increased by 7% for Polimi and reduced by 1.5% for Keromnes compared to the case with LLNL. The resulting K_wrink was then determined by Equation (1) in conjunction with the prevailing conditions within the individual flame front. The differences in the respective laminar flame speeds were thus compensated for by the factor k_turb to obtain the same HRR. This behavior could be observed to the same extent for all three load points (see Figure 6), whereby the case with the Polimi mechanism in particular deviated, as its relatively lower LFS level compared to the other mechanisms required increased flame wrinkling in order to create equivalent combustion behavior. Due to the very similar LFS inputs of LLNL and Keromnes, their differences were quite small and they only showed a noticeable variation for the lower load case.

By comparing the load points with each other, it can be seen that, for higher loads, the HRR and pressure peaks increase, while the MFB10–90% is slightly prolonged. All three show relatively uniform flame wrinkling, with K_wrink between 20 and 30 during the main combustion, caused by the similar flow and mixture conditions. The contribution of the LFS, on the other hand, decreases steadily for higher loads. This is mainly induced by the higher pressure levels during the combustion (see Figure 1), even though the temperature increases and the mixture composition is slightly improved. In addition, the local conditions within the PC at the start of combustion (first 5 °CA after IP) lead to a low LFS, which is caused by the increased residual gas concentration. Furthermore, a higher K_wrink appears due to the locally greater turbulence level. In summary, it can be concluded that the decreasing LFS leads, under equivalent flow conditions, to a slightly longer combustion duration.

4.2. Lambda Variation

This section now analyzes the results with different air-to-fuel ratios at the same load point of 15 bar IMEP and 3000 RPM. In

Table 7. Simulated conditions for the lambda variation at 3000 RPM and 15 bar IMEP.

Lambda	[-]	1.5	2.0	2.5
Ignition point (IP)	[°CA a. FTDC]	−3	−12	−19
Temperature in cyl. @ −20 °CA	[K]	744	735	727
Pressure in cyl. @ −20 °CA	[bar]	35.8	41.7	48.5
TKE in cyl. @ −20 °CA	[m2/s2]	30.6	29.9	30.0
Residual gas in cyl. @ IP	[%]	5.5	5.2	5.2
Lambda in PC @ IP	[-]	1.58	2.09	2.59
TKE in PC @ −20 °CA	[m2/s2]	51.1	51.1	51.3
Residual gas in PC @ IP	[%]	9.1	9.2	9.1

, the relevant conditions inside the cylinder and the pre-chamber (PC) are shown for the period shortly before the combustion started. It shows that the mixture properties and flow conditions were relatively constant, but, with an increasing air-to-fuel ratio, the thermodynamic conditions worsened due to the increasing pressure and lower temperatures. Again, as with the load variation, the local differences between the cylinder and PC showed an increase in the TKE, while the mixture composition deteriorated.

The differences in the mixture conditions and particularly the air-to-fuel ratio influenced the combustion progress for the three variations significantly. Figure 7 displays the crank angle resolved pressure and HRR trace inside the cylinder, as well as the LFS and K_wrink as flame front averaged values for each lambda case and reaction mechanism. Due to the varying LFS input tables for each reaction mechanism, the turbulence intensity factor k_turb is similarly adjusted here as described in Section 4.1. The resulting K_wrink is again determined by Equation (1) in correlation with the existing flow and mixture conditions at the flame front. The variation in the respective LFS are compensated for by the factor k_turb, which is adapted to the same extent for all three lambda cases (see Figure 7). These differences in the mechanisms become particularly clear at low lambda values, where the laminar flame component is especially high.

Nonetheless, the variation in the air-to-fuel ratio here also shows that the resulting wrinkling factor K_wrink is at very different levels, even though the turbulent field is comparable. This factor therefore increases from low values of approximately 10 for lambda 1.5 to higher levels of over 60 at lambda 2.5. An explanation for this rapid increase in the turbulent component is that the excess air reduces the local laminar flame speed on the wrinkled flame surface to such an extent that microscopic turbulent effects become more dominant and increasingly enlarge the flame surface, thus further accelerating the turbulent flame’s wrinkling. The formulation of K_wrink used in this work (see Section 3.1) considers this relation and takes into account the resulting increasing behavior of K_wrink with a decreasing LFS. As additional factors, however, it cannot be completely ruled out that all three reaction mechanisms underestimate the LFS due to extensive extrapolation or that the RANS approach used here in combination with a coarse CFD mesh underestimates the mixture’s homogenization.

Overall, Figure 7 shows that, in the leaner cases, the LFS component decreases to a greater extent than the increase in the turbulent proportion, which ultimately leads to a much slower heat release rate. This greatly extends the combustion duration, and the ignition point requires early onset to reach a similar center of combustion.

The detailed evaluation of the 3D-CFD simulation makes it possible to obtain more precise insights into the combustion chamber. In Figure 8, the mixture formation and flame progress are shown from the top view of the cylinder for all three lambda cases. It displays the current flame propagation at the individual center of combustion for the three reaction mechanisms, and they are based on the exact same lambda distribution inside the cylinder. The differences with the varying lambda are clearly visible here. While, with lambda 1.5, the individual flame jets of the pre-chamber are still clearly visible, this texture becomes increasingly blurred for the leaner cases. In addition, all cases show a tendency towards faster propagation in the direction of the exhaust valves, which agrees well with the lambda distributions, as these are slightly inhomogeneous in all three cases and there is a higher fuel concentration on the exhaust side. Distinguishing between the results with different reaction mechanisms, on the other hand, is more challenging. The results only show small differences in the flame front on closer inspection. Due to the alignment of the heat release rate through the flame wrinkling adaption, all three tend to behave similarly. Only the flame shape in the case with lambda 2.5 and the look-up table based on Polimi shows a greater deviation, as instead of being directed towards the intake side, it tends to travel in the orthogonal direction.

4.3. Knock Evaluation

Compared to all other operating points considered in this work, the case at 15 bar IMEP and lambda 1.5 is the only one with a limitation in the spark advance. Indeed, under these conditions, the engine could not be operated at the optimum center of combustion of 8 °CA after FTDC, as knock phenomena occurred and the ignition angle had to be delayed accordingly. An effective way to estimate knocking within a simulation environment is to analyze the actual mass in self-ignition conditions, which can be evaluated by the Livengood–Wu integral (see Section 3.1). Based on our experience with different engines, a value of 5% of the mass in the self-ignition condition, with respect to the total gas mass inside the cylinder, can be set as a virtual limit to identify knocking conditions. This calibrated threshold value is representative of the engine’s knock limit, as also shown in previous publications by the authors (e.g. [26,35,36]). Since, in this approach, the ignition delay time (IDT) is evaluated locally at all simulation steps, taking into account the pressure, temperature, and chemical composition of the fluid, the influence of the reaction mechanism is evident. A detailed comparison of the mechanisms is given in Figure 9, which shows the temporal progression of the mass in the self-ignition state during combustion (left) and a 3D illustration of the mass’s locations at 11 °CA after FTDC (right). The graph reveals that the knock spots are positioned at the edge of the flame front and in the area between the flame jets and therefore are similar for all three cases. They differ, however, in the respective intensity of the critical areas. The case with Keromnes clearly stands out from the other two cases, as it has larger pink-colored areas and also shows an increase in the maximum mass in the self-ignition condition of 1%-point in the quantitative evaluation. The insensitivity to residual gas rates of Keromnes (as shown in Figure 4) could be a possible explanation for this more aggressive behavior.

5. Summary and Conclusions

The current work addresses the holistic approach of modeling hydrogen’s flame speed and self-ignition behavior in a 0D/1D environment with the consideration of various reaction mechanisms and its transition to a complete 3D-CFD virtual testbench.

With the use of adapted zero- and one-dimensional reactor applications, the values of the laminar flame speed and the ignition delay time were calculated with the reaction kinetics solver Cantera under simplified conditions. The following conclusions can be drawn from the 0D/1D analysis.

• The properties of H₂–air mixtures (including residual gas) can be reduced to four main species and calculated in advance for all cases. Five different reaction mechanisms are compared: LLNL, Polimi, Keromnes, Burke, and SanDiego (Section 2.1 and Section 2.2 ).
• Upon comparing the LFS, the mechanisms are very similar under laboratory conditions. For elevated pressures and temperatures, they all show qualitatively similar behavior, but, quantitatively, they exhibit different levels, while the Polimi mechanism in particular has the lowest flame speeds (Section 2.3).
• When comparing the IDT, the mechanisms are again very similar under laboratory conditions. For further extrapolation ranges, a very strong temperature dependence is decisive for all mechanisms, and, with SanDiego, the resulting delays are particularly long (Section 2.4).
• The influence of the residual gas concentration was also investigated. A dilution effect could be observed for all five mechanisms, which reduced the LFS at higher residual gas rates and prolonged the IDT. Only the case with Keromnes showed noticeable insensitivity with increasing residual gas concentrations (Section 2.5).

This practical 0D/1D approach makes it possible to map the LFS and IDT for all relevant operating points and to use them as look-up tables for the much more computationally intensive 3D-CFD engine simulation. Here, the turbulent flow conditions of the single-cylinder engine were simulated with the 3D-CFD tool QuickSim, whereby detailed models for the ignition and flame propagation were implemented. Special attention during this investigation was paid to the influence of different engine loads and varying lambda operations. The following conclusions can be drawn from the 3D-CFD investigation.

• The engine load variation shows that the LFS is slightly reduced with an increasing load due to elevated pressure levels. However, under comparable flow conditions, similar flame wrinkling effects are found. Local mixture differences between the cylinder and the pre-chamber play a visible role and the reaction mechanisms require individual adjustments to the turbulent intensity (Section 4.1).
• The lambda variation shows that the HRR is clearly influenced by the air-to-fuel ratio. The LFS decreases considerably for higher lambda cases, while the flame wrinkling increases significantly despite a comparable turbulence level. As a result, the burning duration is greatly extended, and the 3D-CFD images also show that the shape of the flame clearly adapts to the cylinder’s mixture conditions (Section 4.2).
• In the case of lambda 1.5 and 15 bar IMEP, the engine was operated at the knock limit. The virtual analysis of the knock origins shows similar tendencies between the compared reaction mechanisms, but the knock intensity with Keromnes is estimated to be particularly high (Section 4.3).

In summary, it can be stated that the transfer between the simplified 0D/1D simulation and the 3D-CFD engine simulation leads to acceptable results for all reaction mechanisms under consideration. In each case, the different influences of the pressure, temperature, and mixture composition on the LFS and the IDT could be demonstrated and transferred to different engine operating points and air-to-fuel ratios. However, due to their different extrapolation methods, individual calibration in the virtual engine environment is necessary.