The effects of hydrogen addition on the flame dynamics of a bluff-body stabilized methane–hydrogen turbulent flame are studied with large eddy simulation (LES). The LES is carried out with the thickened flame model and global kinetic mechanisms calibrated for the methane–hydrogen mixtures. Conjugate heat transfer is included in the LES to consider a proper wall temperature while the flame shape changes with hydrogen addition. A data-based calibration of the global mechanisms is done with a methodology based on reproducing the net species production rates computed with a detailed kinetic mechanism. An improvement in this methodology is proposed to increase its accuracy and reliability. The calibrated mechanisms accurately describe the variation of the laminar flame speed and the thermal flame thickness with hydrogen addition and equivalence ratio in a freely propagating premixed flame. The variations of the consumption speed and the thermal flame thickness with the strain rate in a symmetric counterflow premixed flame are also well predicted. The numerical simulations reproduce the transition from V- to M-shape flame induced by hydrogen addition, and the axial distribution of the heat release agrees with the experimental measurements of OH chemiluminescence. The unit impulse response and the flame transfer function are computed from the LES data using system identification (SysID). The flame transfer functions show a remarkable agreement with the experimental data, demonstrating that the LES-SysID approach using properly calibrated global mechanisms can predict the response of turbulent methane–hydrogen flames to velocity fluctuations. A comparison of the unit impulse response for the various hydrogen additions is presented, and the effect of hydrogen in the flow–flame interaction of the burner evaluated is discussed.

The blending of hydrogen with conventional hydrocarbon gaseous fuels has received increased attention for its potential to contribute to CO2 reduction.1 The effect of hydrogen blending on the thermoacoustic response of a lean premixed combustion system depends on the combustion system and its operating conditions.2,3 Most experimental studies show that thermoacoustic instability may be induced by adding hydrogen to methane–air flames.4–6 Conversely, hydrogen addition may also suppress thermoacoustic instability.2,7 The effect of hydrogen on thermoacoustic instabilities is highly related to the impact of hydrogen on flame shape and length.5,6,8 Hydrogen addition reduces the flame length and may trigger a transition from V- to M-shape in non-adiabatic turbulent flames,6,9–11 even if the unstretched laminar flame speed remains constant.9,10

Numerical simulation can be used to determine the frequency response of turbulent flames to study their thermoacoustic behavior.12,13 However, numerical simulation of turbulent combustion of hydrogen blends is challenging due to the properties of hydrogen. Direct numerical simulations of hydrogen-enriched methane–air turbulent flames expose that considering accurate chemistry and diffusion are necessary to correctly describe the flame behavior.14,15 These requirements can be met by performing large eddy simulation using a combustion model based on transported chemistry, where the transport and reaction of the species consider effects, such as heat loss, stretch, and preferential diffusion (non-unity Lewis number and differential diffusion).

The accuracy of any model based on transported chemistry lies, among other aspects, on the chemical kinetic mechanism used to describe the reaction process. Detailed kinetic mechanisms (DM) are impractical in LES due to the large number of species and reaction steps to be considered. Alternatively, reduced versions of these mechanism (the so-called skeletal mechanisms) can be used in LES,16 but still with a relatively high computational cost. Analytically reduced mechanism (ARM)17,18 allows reducing the computational cost by identifying the species that do not need to be transported with the flow. However, different from a skeletal mechanism, the production rate of the transported species in an ARM is not expressed as a combination of elementary reaction in the Arrhenius form but rather as a combination of complex analytical relations.17 This complexity, together with the number of transported species that could still be more than twenty, makes ARM not the first option for applications where lower computational cost is needed due to computational power limitations or a high number of conditions required to be simulated. In such cases, global mechanisms are an alternative that reduces the computational cost.

A global mechanism (GM) is an empirical scheme that only considers the major species involved in the combustion process, now described by a reduced number of global reaction steps. Global mechanisms are calibrated by adjusting the rate parameters in the Arrhenius equation of each reaction step to match some of the main physical properties of the flame, such as laminar flame speed, auto-ignition delay time, or species production rates, obtained from experiments or using a detailed kinetic mechanism for given operation conditions. The calibration is then formulated as an optimization problem, where various methods exist to find a proper set of parameters.19–22 The use of global mechanisms in LES has reproduced important aspects of turbulent combustion.23–26 Although some aspects of the turbulence–chemistry-diffusion interactions are lost if intermediate species are not considered,22,27 the study of global mechanism capabilities to describe different aspects of combustion is still an ongoing task.28–30 

Given that thermoacoustic instabilities represent an important challenge in the design of combustion systems, comprehensive validation of LES models should include the prediction of flame dynamics. Therefore, the present work further tests the capabilities of global mechanisms by considering methane–hydrogen mixtures with highly different chemical-diffusion properties and focusing on the prediction of flame dynamics. Both hydrogen combustion and thermoacoustic are two major research topics on present combustion. Global mechanisms are calibrated for various lean methane–hydrogen–air mixtures based on the methodology proposed by Polifke et al.19 This methodology is improved in the present work to make it more accurate. Canonical 1D laminar flames are used to validate the global mechanisms compared with detailed and skeletal kinetic mechanisms. Then, the global mechanisms are used in the LES of a turbulent premixed burner with various methane–hydrogen mixtures. The results are compared against experimental data. First, predictions of the mean flame shape and length are evaluated, and then, the flame dynamics are analyzed through the unit impulse response and the flame transfer functions computed using system identification.13 

The global mechanism used in this study is a three-step mechanism with an initial step for methane (CH4) breakdown into carbon monoxide (CO) and hydrogen (H2) and two additional steps for CO and H2 oxidation into carbon dioxide (CO2) and water (H2O), respectively, as follows:
CH 4 + 0. 5 O 2 CO + 2 H 2 , ( R 1 ) CO + 0. 5 O 2 CO 2 , ( R 2 ) H 2 + 0. 5 O 2 H 2 O . ( R 3 )

This mechanism is selected due to its simplicity and the good results for wide ranges of hydrogen content. In this mechanism, only the step for CO oxidation is reversible. The rate parameters of this mechanism are calibrated based on the methodology proposed by Polifke et al.,19 which allows calibrating each reaction step independently based on a target reaction progress rate. The methodology is divided in two parts:

First, a detailed mechanism is used to compute a 1D freely propagating laminar flame at a reference operating condition: pressure, temperature, hydrogen fraction, and equivalence ratio. The profiles of temperature T, species mole concentration [ C i ], and species production rates ω ̇ i are obtained from this simulation. Then, the goal is to determine the “target reaction rates” of the GM so that the production rate of the species in the GM matches the ones obtained with the DM along the 1D flame.

For a given detailed kinetic mechanism, the production rate of each species is given by
ω ̇ i = j = 1 M ν i , j r j , with ν i , j = ν i , j ν i , j ,
(1)
where rj is the reaction progress rate of the reaction j and ν i , j and ν i , j are the stoichiometric coefficient of species i in reaction j as reactant and product, respectively. For the whole system of N species and M reaction steps, Eq. (1) can be expressed in matrix form as
p [ N , 1 ] = M [ N , M ] r [ M , 1 ] ,
(2)
where p and r are the vectors of production rates and reaction progress rates, respectively, and M is the matrix of stoichiometric coefficients. The same applies for the GM, yielding
p ̂ [ N ̂ , 1 ] = M ̂ [ N ̂ , M ̂ ] r ̂ [ M ̂ , 1 ] ,
(3)
where . ̂ refers to the global mechanism. Then, the right-hand side of Eq. (3) is equaled to the left-hand side of Eq. (2) ignoring the production rates of the species in the DM not included in the GM, such as
p [ N ̂ , 1 ] = M ̂ [ N ̂ , M ̂ ] r ̂ [ M ̂ , 1 ] .
(4)
Note that p corresponds to the production rates of the species obtained from the 1D flame simulation using the DM. Then, Eq. (4) can be solved to obtain the reaction progress rates r ̂ of the GM that produce the same species production rates that the DM. This, however, is not possible in general, because the matrix M ̂ is not square ( N ̂> M ̂). Polifke et al.19 addressed this issue by selecting a number M ̂ of species in the GM, denoted “principal species,” to build a square matrix M ̂. The drawback of this approach is that the overall solution depends on the non-unique choice of these principal species. To overcome this limitation, a different approach is proposed in this work. It consists of solving Eq. (4) using the method of least squares so that all species production rates are considered for the reaction progress rates of the GM. The vector r ̂ * of the target progress rates is, thus, defined as
r ̂ * = ( M ̂ T M ̂ ) 1 M ̂ T p .
(5)

This enhancement improves the accuracy of the calibration methodology and makes it more physics-driven.

The target rates for an example of GM calibration are shown in Fig. 1(a) with solid lines. The GM with the reaction rates r ̂ * closely reproduces the species production rates obtained with the detailed mechanism, as shown in Fig. 1(b). Note that only a single 1D laminar flame calculation is used to obtain the target reaction rates for each calibration, and no more are needed from now on.

FIG. 1.

Example of GM calibration, X H 2 = 0 and ϕ = 0.7. (a) Reaction progress rates of GM. (b) Species production rates ω ̇ i from the 1D laminar flame, the GM target reaction rates, and the GM reaction rates after calibration.

FIG. 1.

Example of GM calibration, X H 2 = 0 and ϕ = 0.7. (a) Reaction progress rates of GM. (b) Species production rates ω ̇ i from the 1D laminar flame, the GM target reaction rates, and the GM reaction rates after calibration.

Close modal
The second part of the calibration is to adjust the rate parameters of each reaction step to match the target reaction progress rates r ̂ j *. The reaction progress rate of each reaction step j is given by
r ̂ j = k j i = 1 N ̂ [ C i ] ν i , j non - reversible k j K eqj i = 1 N ̂ [ C i ] ν i , j reversible , with k j = A j T b j exp ( E a j R T ) ,
(6)
where A is the pre-exponential factor, b is the temperature exponent, Ea is the activation energy, R is the gas constant, and Keq is the equilibrium constant as a function of T.

The calibration of each reaction step is done using a genetic algorithm31 to minimize the cost-function19 relating r ̂ j and r ̂ j *, where the values of T, [ C i ], and r ̂ j * along the 1D laminar flame are the reference data and Aj, bj, and E a j are the parameters to be calibrated. For the nonreversible reaction steps, ν i , j is also calibrated but restricted to be positive to increase the robustness of the GM.

The reaction progress rates obtained after the calibration are shown in Fig. 1(a) with dashed lines and agree reasonably well with the target reaction progress rates. The same applies to the species progress rates obtained with these reaction rates, as seen in Fig. 1(b). The higher difference is obtained for the H2 oxidation step in the high-temperature region corresponding to the oxidation layer.

The reference operating condition for the calibration of the mechanism must agree with the operating condition in the LES, which, in the present work, corresponds to an atmospheric premixed CH4–H2–air flame with an equivalence ratio ϕ = 0.7 and hydrogen fraction by volume of fuel X H 2 = 0.0, 0.25, 0.57, and 0.67. The GM is calibrated for two conditions, X H 2 = 0.0 and X H 2 = 0.57, leading to two sets of parameters, referred to as 3S-00H and 3S-57H, respectively. The detailed kinetic mechanism UC San Diego,32 referred to as UC-SD, is used to compute the reference laminar flame for the calibration. The calibrated parameters for the three-step GM are presented in Table I.

TABLE I.

Global kinetic mechanism calibrated for CH4–H2–air. ϕ = 0.7. 3S-00H X H 2 = 0.0 and 3S-57H X H 2 = 0.57. Units: cal, K, mol, cm3, s.

A b Ea ν CH 4 ν H 2 ν O 2
3S-00H             
(R1)  4.65 × 107  0.0  32 082  0.64  …  0.17 
(R2)  1.51 × 109  0.33  24 334  …  …  … 
(R3)  3.55 × 1013  1.0  20 847  …  1.99  0.19 
3S-57H             
(R1)  2.52 × 109  0.0  21 482  0.49  …  0.82 
(R2)  1.68 × 109  0.0  15 694  …  …  … 
(R3)  9.38 × 1011  1.0  20 938  …  1.91  0.10 
A b Ea ν CH 4 ν H 2 ν O 2
3S-00H             
(R1)  4.65 × 107  0.0  32 082  0.64  …  0.17 
(R2)  1.51 × 109  0.33  24 334  …  …  … 
(R3)  3.55 × 1013  1.0  20 847  …  1.99  0.19 
3S-57H             
(R1)  2.52 × 109  0.0  21 482  0.49  …  0.82 
(R2)  1.68 × 109  0.0  15 694  …  …  … 
(R3)  9.38 × 1011  1.0  20 938  …  1.91  0.10 

To verify that the GM can be used in the LES of the CH4–H2 flames, an evaluation of the global mechanism is done in the frame of 1D premixed laminar flames using the software Cantera.33 Full multicomponent mass diffusion and Soret effect are considered. The reference pressure and temperature of the reactants are equal to 101.3 kPa and 300 K, respectively. The global mechanism is compared with the detailed mechanism UC-SD used as a reference for the calibration and the well-known mechanism DRM19.34 The latter as a reference for the accuracy of a skeletal mechanism.

The freely propagating premixed flame (free-flame) is used to compute the unstretched laminar flame speed S L 0 and the unstretched thermal flame thickness defined by δ T = ( T b T u ) / max ( d T / d x ), where subscripts u and b denote the unburnt and burnt sides of the flame, respectively. Figure 2 shows the results for ϕ = 0.7 and CH4–H2 mixtures ranging from pure methane to pure hydrogen. Both global mechanisms can predict the values of S L 0 for a wide range of hydrogen concentrations although they are calibrated for only one specific condition. The GM calibrated for X H 2 = 0.0 shows better accuracy for low hydrogen additions, while the one calibrated for X H 2 = 0.57 has better accuracy for high hydrogen additions. In these ranges, the error in S L 0 with the GM is less than 10% respect to the DM.

FIG. 2.

Variation of unstretched laminar flame speed, S L 0, and thermal flame thickness, δ T 0, with hydrogen addition, X H 2. ϕ = 0.7.

FIG. 2.

Variation of unstretched laminar flame speed, S L 0, and thermal flame thickness, δ T 0, with hydrogen addition, X H 2. ϕ = 0.7.

Close modal

Regarding the flame thickness, the GM results are in good agreement with the detailed and skeletal mechanisms. However, the flame thickness is underpredicted for mixtures close to pure hydrogen where the GM effectively reduces to a one-step reaction.

Although only perfectly premixed flames are considered in this study, it is relevant to evaluate the accuracy of the global mechanism for various equivalence ratios ϕ since the latter may vary along the flame front due to preferential diffusion of hydrogen. Figure 3 shows the evolution of S L 0 as a function of ϕ. The GM 3S-00H is used for the cases with X H 2 = 0.0 and 0.25, while the GM 3S-57H is used for X H 2 = 0.57 and 0.67, both will be referred to as GM hereafter. In all cases, the GM archives a good accuracy for a wide range of equivalence ratios compared with the detailed and skeletal mechanisms. The DRM19 mechanism consistently over-predicts S L 0, especially for higher equivalence ratio and hydrogen addition. The accuracy of the GM diminishes when the equivalence ratio is far from the calibration condition. The GM could be tailored for rich and very lean mixtures by calibrating the mechanism for multiple equivalences ratios following the presented methodology and letting one of the Arrhenius coefficients vary with the equivalence ratio [e.g., A j = A j ( ϕ )], as can be found in the literature for CH4–air combustion.24,35

FIG. 3.

Variation of unstretched laminar flame speed, S L 0, with equivalence ratio, ϕ.

FIG. 3.

Variation of unstretched laminar flame speed, S L 0, with equivalence ratio, ϕ.

Close modal

In the present work, the global mechanism is calibrated for a reference temperature of 300 K. However, the GM is evaluated for various unburnt gas temperatures Tu ranging from 250 to 450 K, as it may be relevant to accurately represent flames with heat loss or gain. Figure 4 shows the evolution of S L 0 as a function of Tu. For all cases, the effect of the unburnt gas temperature on the laminar flame speed is well captured by the GM. The major discrepancy with the detailed mechanism is presented at high Tu by the case with X H 2 = 0.67, a condition distinct from the one for the calibration. Nevertheless, the error is still in the magnitude range of the skeletal mechanism, showing the strength of the present calibration methodology.

FIG. 4.

Variation of unstretched laminar flame speed, S L 0, with unburnt gas temperature, Tu.

FIG. 4.

Variation of unstretched laminar flame speed, S L 0, with unburnt gas temperature, Tu.

Close modal

A global mechanism calibrated with the present methodology can also resolve the internal structure of the laminar flame with reasonable accuracy, as shown in Fig. 5 for X H 2 = 0.0 and 0.57. The position across the free-flame x is normalized by the flame thickness δ T o REF computed with the detailed mechanism, and the origin is located at max ( d T / d x ). The agreement between the GM and the DM is remarkable in the inner layer and around the peak value. The higher difference corresponds to the production rate of H2 in agreement with the difference in the reaction rate of the H2 oxidation step shown in Fig. 1. This difference may be related to the fact that in the GM chosen, all the hydrogen in the methane molecule has to pass through H2 to oxidize to H2O.

FIG. 5.

Profiles of heat release rate, Q ̇, (×108 W/m3) and production rate, ω ̇, of CH4, CO, and H2 (kmol/m3/s) across the freely propagating laminar flame, for cases with X H 2 = 0.0 (left) and 0.57 (right).

FIG. 5.

Profiles of heat release rate, Q ̇, (×108 W/m3) and production rate, ω ̇, of CH4, CO, and H2 (kmol/m3/s) across the freely propagating laminar flame, for cases with X H 2 = 0.0 (left) and 0.57 (right).

Close modal

Nevertheless, the species mole fractions and temperature are also in good agreement with the DM, as shown in Fig. 6. The peak concentration of the intermediate species CO and H2 in the case of pure methane is well predicted by the GM. The GM reaches the equilibrium condition sooner, which is inherent to the absence of more intermediate species.

FIG. 6.

Profiles of temperature and mole fractions of CH4, O2, H2O, CO2, CO, and H2 across the freely propagating laminar flame, for cases with X H 2 = 0.0 (left) and 0.57 (right).

FIG. 6.

Profiles of temperature and mole fractions of CH4, O2, H2O, CO2, CO, and H2 across the freely propagating laminar flame, for cases with X H 2 = 0.0 (left) and 0.57 (right).

Close modal
The transition from V- to M-shape flame induced by hydrogen addition has been related to the different response of the flame to stretch.5,10 Therefore, the symmetric counterflow premixed flame (twin-flame) is used to evaluate the capability of the GM to describe the response of the laminar flame to positive strain. Figure 7 shows the variation of the consumption speed with the increase in the strain rate, which is expressed in a dimensionless form using the Karlovitz number Ka. The results for X H 2 = 0.57 are not presented for the sake of conciseness, but the trend is the same. The consumption speed for the multicomponent fuel is defined in the present work as
S c i = 1 N f η i ω ̇ i d x ρ u i = 1 N f η i ( Y i , b Y i , u ) ,
(7)
where ρ is the density, Yi is the mass fraction of the species i, η i = Y i , u / Y f , u, and f denotes the fuel mixture. The GM can describe the increase in the flame speed with a strain rate produced by the effect of preferential diffusion. This result is especially relevant for hydrogen-enriched turbulent flames where the stretch induced by the turbulence produces local changes of the flame speed.14,15 The maximum Sc predicted by the GM in the twin-flame differs from the DM by less than 10%.
FIG. 7.

Variation of consumption speed, Sc, with a strain rate. Extinction strain rate shown by the marker (×).

FIG. 7.

Variation of consumption speed, Sc, with a strain rate. Extinction strain rate shown by the marker (×).

Close modal

As shown in Fig. 8, the decrease in the flame thickness due to positive strain is also well described by the GM up to the maximum strain rate supported by the laminar twin-flame, known as the extinction strain rate. The values of the extinction strain rate computed with the GM are consistent with those obtained with the detailed and skeletal mechanisms. The skeletal mechanism consistently under-predicts the extinction strain rate in 17% compared to the DM, while the GM is in better agreement with the DM.

FIG. 8.

Variation of thermal flame thickness, δT, with a strain rate. Extinction strain rate shown by the marker (×).

FIG. 8.

Variation of thermal flame thickness, δT, with a strain rate. Extinction strain rate shown by the marker (×).

Close modal

Another effect of preferential diffusion is the appearance of thermo-diffusive instabilities exhibited as cellular instabilities in the flame front. Such thermo-diffusive instabilities increase the flame surface area, increasing the flame propagation speed even in turbulent flames.36 An evaluation of the capability of the GM to capture thermo-diffusive instabilities is presented in the  Appendix. However, it is not clear if this type of instability plays an important role in flame dynamics.

The burner considered in the LES is a bluff body stabilized turbulent flame confined in a cylindrical combustion chamber open to the atmosphere, investigated extensively at NTNU.6,8,37 The numerical domain is shown in Fig. 9. The solid regions corresponding to the combustion chamber, the burner body, and the bluff body are included in the LES through conjugate heat transfer formulation (CHT). The combustion chamber is made of quartz, and all the other parts are made of SAE-316L stainless steel. The thermal power is fixed at 7 kW, and four hydrogen concentrations are evaluated, X H 2 = 0.0, 0.25, 0.57, and 0.67. Fuel and air are perfectly premixed, with a fixed equivalence ratio of ϕ = 0.7. The Reynolds number and Karlovitz number for all flames at the dump plane are around 38 300 and 5, respectively.

FIG. 9.

Schematic of the computational domain used in the LES. The white solid regions are not considered for the CHT. Dimensions are given in mm.

FIG. 9.

Schematic of the computational domain used in the LES. The white solid regions are not considered for the CHT. Dimensions are given in mm.

Close modal

The computational domain is discretized with a structured mesh of 7.82 M hexahedral cells, of which 0.34 M correspond to the solid regions. The mean cell size is 0.16 mm in the flame area and 0.25 mm in the space between the grub screws and the dump plane. The mesh is refined near the walls yielding a y + 1. The software Simcenter STAR-CCM+38 is used to solve the filtered incompressible multi-species Navier–Stokes equations. The sub-grid stress tensor is modeled by the WALE39 model without using a wall function, and the PISO algorithm40 is used for the pressure–velocity coupling.

Turbulence–flame interaction is considered with the thickened flame model (TFM). The flame thickness is discretized by at least five grid points, and the power law is used to model the sub-grid flame interactions and wrinkling.41 The reference flame thickness is computed as a function of T and S L 0. As for the 1D flames, GM 3S-00H is used for the cases with X H 2 = 0.0 and 0.25, while GM 3S-57H is used for X H 2 = 0.57 and 0.67. The dynamic viscosity is computed as a function of the temperature using the power-law, with the coefficients for air. The mixture thermal conductivity and the mass diffusion coefficient of each species are computed from the viscosity using constant Prandtl and species Schmidt numbers, respectively. These numbers are obtained from the mean values across the 1D laminar free-flame for each fuel mixture. Turbulent Prandtl and Schmidt numbers are set to 0.7. A second-order bounded-central difference scheme42 is used for the convective term of the momentum equation, while the second-order upwind scheme is used for enthalpy and mass species equations. The backward time difference scheme is used with a 2 μs time step to achieve a CFL below 0.9.

The energy conservation equation is solved in the solid regions of the burner to consider a proper wall temperature while the flame shape is changing with hydrogen addition. Wall heat loss by radiation and convection are considered with constant emissivity ϵ, and convective heat transfer coefficient h reported in Fig. 9 for the combustion chamber and the burner body. The period of transient heating of the solid regions is computed using the multi-timescale workflow43 to reduce computational time. The LES statistics are calculated during six flow-through times (0.12 s) after the solid region has reached the steady-state thermal condition.

The flame shape is one of the main aspects to predict with the LES, since it plays an important role in the thermoacoustic instabilities.5,6 Figure 10 shows a comparison of the flame shape obtained from the experiments and the LES. The heat release rate Q ̇ is used to characterize the flame shape in the LES, while the OH* chemiluminescence is used in the experiments.6,8 The LES with the global mechanism is able to reproduce the change of the flame shape with hydrogen addition. For pure methane, the flame is attached to the inner shear layer with large portions of the flame brush in contact with the combustion chamber wall. The flame becomes more compact with hydrogen addition, and a transition from V to M shape occurs for X H 2 between 0.25 and 0.57.

FIG. 10.

Contours of heat release rate from LES and experimental OH* chemiluminescence. (a) Instantaneous Q ̇ and Abel deconvolution of mean OH*. (b) Line-of-sight of mean Q ̇ and line-of-sight of mean OH*.

FIG. 10.

Contours of heat release rate from LES and experimental OH* chemiluminescence. (a) Instantaneous Q ̇ and Abel deconvolution of mean OH*. (b) Line-of-sight of mean Q ̇ and line-of-sight of mean OH*.

Close modal

A comparison of the flame length using the stream-wise profiles of heat release integrated along the transverse plane, q x ¯ ( x ) = Q ̇ ( x , y , z ) d y d z, is presented in Fig. 11. The OH* data are normalized such that the volume integral is equal to 7 kW. The general distribution of q x ¯ and the location of its maximum value are well reproduced by the LES for all methane–hydrogen mixtures. For the case of pure methane, the flame length is under-predicted by the LES. Good agreement is obtained near the dump plane so that the difference may be associated with the flame–wall interaction. A LES including radiation heat transfer with a gray gas model does not significantly improve the flame length discrepancy for pure methane, neither does a higher convective heat transfer coefficient outside of the combustion chamber. (These results are not presented for conciseness.) The difference is then attributed to a limitation of the global mechanism to describe direct wall–flame interaction. The lack of intermediate species may result in a flame less sensitive to direct contact with the wall. This is not a problem for cases where considerable portions of the flame are not attached to the chamber walls.

FIG. 11.

Stream-wise heat release distribution, q x ¯, for various hydrogen additions. The origin is located at the dump plane.

FIG. 11.

Stream-wise heat release distribution, q x ¯, for various hydrogen additions. The origin is located at the dump plane.

Close modal

The dynamic response of a flame to upstream flow perturbations can be represented by the unit impulse response in the time domain and the flame transfer function (FTF) in the frequency domain, which are important to predict the acoustic response of a combustion system. Computing the flame dynamics with LES is challenging for the combustion model, so it constitutes an excellent way to examine the performance of the LES based on global mechanisms. The flame dynamics are obtained from the LES data using a method based on system identification (LES-SysID).13 For this purpose, the flow is forced with a Daubechies wavelet-based broadband signal44 superimposed on the mean flow at the inlet. The length of the signal is 0.2 s, and its amplitude is 4% of the mean value. This amplitude is the same one used in the experiments to ensure a linear response of the flame. A reference velocity ur is measured at the dump plane as the area-weighted average axial velocity. At the same time, the heat release rate Q ̇ is integrated over the volume of the computational domain. The unit impulse response is determined by the optimal linear least squares estimation between the auto-correlation matrix and the cross correlation vector of the time series data ur and Q ̇. Then, the FTF is obtained by a z-transformation of the unit impulse response.

It is important to note that the flame dynamics in thermo-acoustics instabilities follow different mechanisms of the acoustics–flame interaction. In the case of perfectly premixed low Mach deflagration flames, heat release fluctuations induced by flow disturbances, such as reactant mass flow fluctuations and hydrodynamic fluctuations, are expected to be the most important since the effect of pressure, temperature, and strain rate perturbations directly accompanying the acoustic waves are pretty weak.45,46 The latter justifies the use of incompressible LES to characterize the flame dynamics of these types of flames. In the context of LES-SysID, the assumption of an incompressible flow not only decreases the computational cost of the LES but also breaks the coupling from the flame back to acoustics, increasing the signal-to-noise ratio of the time series data ur and Q ̇ since self-excited acoustics cannot occur.

The respective unit impulse response of the turbulent flames for the four hydrogen additions is presented in Fig. 12. All cases present a similar initial response to a sudden increase in axial velocity, which produces an increase in the mixture supplied to the flame. The extra mixture is transported through the flame as it burns and perturbs the flame surface area. This results in an overall increase in heat release, with a temporal distribution similar to the spatial distribution of heat release in the axial direction,47 seen in Fig. 11. The increase in the heat release is followed by a decrease while the flame returns to its initial state. In the cases of high hydrogen addition, the heat release increases again when the vortices produced by the grub screws reach the flame. These vortex waves increase the flame surface area and then the heat release rate. Since the vortices do not supply more mixture to the flame, the increase in heat release is followed by a similar decrease due to mass and energy conservation.

FIG. 12.

Unit impulse responses for various hydrogen additions. τ1: time-delay from dump plane to mean flame length (). τ2: time-delay from screws to mean flame length ().

FIG. 12.

Unit impulse responses for various hydrogen additions. τ1: time-delay from dump plane to mean flame length (). τ2: time-delay from screws to mean flame length ().

Close modal
The described flow-flame interaction is characterized by two different time-delays.8 The first time-delay from the dump plane to the mean flame length is
τ 1 H / u p ,
(8)
where up is the bulk velocity at the dump plane and H is the mean flame length defined by H = x q x ¯ d x / q x ¯ d x. The second time-delay from the screws to the mean flame length is
τ 2 L / u b + H / u p ,
(9)
where ub is the bulk velocity after the screws and L is the distance from the screws to the dump plane. Both time-delays are retrieved by the LES-SysID.

In the cases of low hydrogen addition, the effect of the vortices produced by the grub screws on the flame response is not significant. This happens when the frequency of the vortex waves is close or higher than the cutoff frequency of the flame response.37 The cutoff frequency f c 1 / τ 1 defines the band of frequencies in which the flame is more sensitive to upstream flow fluctuations and acts as an amplifier.

For the cases with low hydrogen addition, X H 2 = 0.0 and 0.25, the unit impulse response presents an initial decrease in the heat release before its main increase. This behavior may be related to a downstream movement of the leading edge of the flame closer to the dump plane produced by the higher axial velocity, momentarily decreasing the flame surface area. This initial decrease in the heat release does not happen for cases with high hydrogen addition, X H 2 = 0.57 and 0.67, due to the higher flame speed that produces a more robust flame close to the dump plane, as shown in Figs. 10 and 11.

The flame transfer function is expressed as the ratio of the total heat release rate fluctuations to the fluctuations of the reference velocity in the frequency domain as
FTF ( f ) = Q ̇ ( f ) / Q ̇ ¯ u r ( f ) / u r ¯ .
(10)

Figures 13 and 14 show the flame frequency response as FTF(f) = G ( f ) exp ( i θ ( f ) ), where G is the gain and θ is the phase. The confidence interval of the SysID is shown with two standard deviations σ, which represents the aleatory uncertainty caused by the resolved turbulent structures in the LES. The experimental FTF is measured at each frequency value as described by Æsøy et al.6,8 A remarkable agreement between the experiments and the LES-SysID is obtained, for both gain and phase, in particular for high hydrogen content. The latter is particularly important for thermoacoustic instability because it gives the time-delay, τ ( f ) = θ / ( 2 π f ), between the upstream flow perturbations and the heat release fluctuations, where the average τ ( f ) is close to τ1. The LES-SysID method is able to reproduce well the modulation in the gain of the FTF, produced by the interaction between the flame and the vortex shedding from the grub screws.37 The addition of hydrogen produces an increase in the cutoff frequency and a decrease in the time-delay. Both effects are related to the flame length and are captured very well by the LES using the global mechanism. In the case of pure methane, the time-delay is under-predicted in agreement with the shorter flame in the LES compared with the experiments.

FIG. 13.

Gain of the FTF for various hydrogen additions. Confidence interval from SysID (shaded area).

FIG. 13.

Gain of the FTF for various hydrogen additions. Confidence interval from SysID (shaded area).

Close modal
FIG. 14.

Phase of the FTF for various hydrogen additions. Confidence interval from SysID (shaded area).

FIG. 14.

Phase of the FTF for various hydrogen additions. Confidence interval from SysID (shaded area).

Close modal

The result of the FTF shows that the LES-SysID approach can predict the dynamics of turbulent CH4–H2 flames considering the different mechanisms of flow-flame interaction.

This work presents the large eddy simulation of a bluff-body stabilized turbulent flame with various methane–hydrogen mixtures using calibrated global mechanisms. A three-step mechanism with a H2 oxidation step was calibrated for various methane–hydrogen mixtures. A calibration methodology based on the computation of a single 1D laminar flame per calibration has been selected for its robustness and its time-efficiency. The calibration methodology has been enhanced in this work in order to make it more accurate and less user dependent. The calibrated global mechanisms accurately describe global flame parameters, such as the variation of the laminar flame speed and the thermal flame thickness with hydrogen addition and equivalence ratio, but also the internal flame structure and the response to positive strain rate.

Good prediction of the flame shape is obtained with the LES using thickened flame mode with GM, describing the transition from V- to M-shape flame induced by hydrogen addition. The axial distribution of the heat release and the flame length from the LES compares well with the experimental data. For the case of pure methane, where large portions of the flame are in contact with the walls, the flame length is under-predicted. Future work may focus on the effect of the GM and the reference conditions used for its calibration on the flame–wall interaction and the impact of the near-wall turbulence modeling.

The flame transfer functions computed from the LES data using system identification are in good agreement with the experiments. The effects of hydrogen addition on the flame frequency response are captured by the LES with the GM, including the modulation of the gain and phase produced by the flow–flame interaction. These results show the potential of global kinetic mechanisms and the calibration methodology in large-eddy simulation of hydrogen blends for thermoacoustics applications. The same methodology can be applied for various fuel mixtures and operating conditions.

This work was part of the Marie Skłodowska–Curie Innovative Training Network Pollution Know-How and Abatement (POLKA). We gratefully acknowledge the financial support from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie Grant Agreement No. 813367. The authors would also like to thank Eirik Æsøy and co-workers at NTNU for the valuable information and experimental data provided.

The authors have no conflicts to disclose.

Alex Mauricio García Vergara: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Writing – original draft (lead); Writing – review and editing (lead). Sophie Le Bras: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (supporting); Writing – review and editing (equal). Jens Prager: Conceptualization (equal); Investigation (supporting); Methodology (equal); Software (equal); Writing – review and editing (supporting). Matthias Häringer: Investigation (equal); Methodology (equal); Writing – original draft (equal). Wolfgang Polifke: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review and editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The cellular instabilities in the flame front induced by thermo-diffusive instabilities are more easily observed in the absence of turbulence. Therefore, to evaluate the capability of the GM to capture thermo-diffusive instabilities, the simulation of the two-dimensional (2D) version of the symmetric counterflow premixed flame is performed. This configuration consists of two opposing streams of reactants, creating a stagnation point flow with a laminar flame located at each side of the stagnation point. The numerical domain with a symmetric boundary condition at the stagnation point is shown in Fig. 15. The domain has a width of 50 mm and a distance between the flow inlet and the stagnation point of 30 mm. It is discretized with a structured mesh of 0.11 M quadrilateral cells with a cell size of 0.1 mm in the flame area. The laminar flame is located where the flame speed balances the local axial velocity; thus, an inlet velocity of 2.5 times the unstretched laminar flame speed is specified to get a similar flame location for the various hydrogen fuel additions. The numerical setup is the same as the one used for the LES of the turbulent flame but with no turbulent model. The results obtained from the global mechanisms are compared with those obtained from the skeletal mechanism.

FIG. 15.

Schematic of the computational domain used in the simulation of the 2D symmetric counterflow premixed flame.

FIG. 15.

Schematic of the computational domain used in the simulation of the 2D symmetric counterflow premixed flame.

Close modal

The results of the 2D simulation are presented in Fig. 16 using the heat release rate to characterize the flame shape. For high hydrogen content, the flame front exhibits cellular structures. These cellular structures are induced mainly by a combination of thermo-diffusive and hydrodynamic instabilities. While the hydrogen addition increases, the flame becomes more thermo-diffusive unstable, increasing the number of cellular structures per unit of area. The heat release rate decreases in concavely curved flame regions due to preferential diffusion. In contrast, the opposite occurs in convexly curved regions but to a lesser extent due to the lower curvature. The GM seems to overpredict this effect. It may be due to the lack of intermediate species, particularly the radical H, which diffuses tangential to the flame front,48 slightly attenuating the impact of the hydrogen differential diffusion. Nevertheless, the GM is able to capture the thermo-diffusive instabilities with good accuracy compared with the skeletal mechanism.

FIG. 16.

Contours of instantaneous heat release rate.

FIG. 16.

Contours of instantaneous heat release rate.

Close modal
1.
K.-K.
Lam
,
P.
Geipel
, and
J.
Larfeldt
, “
Hydrogen enriched combustion testing of siemens Industrial SGT-400 at atmospheric conditions
,”
J. Eng. Gas Turbines Power
137
,
021502
(
2015
).
2.
G.
Oztarlik
,
L.
Selle
,
T.
Poinsot
, and
T.
Schuller
, “
Suppression of instabilities of swirled premixed flames with minimal secondary hydrogen injection
,”
Combust. Flame
214
,
266
276
(
2020
).
3.
I.
Chterev
and
I.
Boxx
, “
Effect of hydrogen enrichment on the dynamics of a lean technically premixed elevated pressure flame
,”
Combust. Flame
225
,
149
159
(
2021
).
4.
J.
Zhang
and
A.
Ratner
, “
Experimental study on the excitation of thermoacoustic instability of hydrogen-methane/air premixed flames under atmospheric and elevated pressure conditions
,”
Int. J. Hydrogen Energy
44
,
21324
21335
(
2019
).
5.
S.
Shanbhogue
,
Y.
Sanusi
,
S.
Taamallah
,
M.
Habib
,
E.
Mokheimer
, and
A.
Ghoniem
, “
Flame macrostructures, combustion instability and extinction strain scaling in swirl-stabilized premixed CH4/H2 combustion
,”
Combust. Flame
163
,
494
507
(
2016
).
6.
J. G.
Aguilar
,
E.
Æsøy
, and
J. R.
Dawson
, “
Predicting the influence of hydrogen in combustion instabilities
,” in
Symposium on Thermoacoustics in Combustion: Industry Meets Academia (SoTiC 2021)
(
Munich
,
Germany
,
2021
), p.
9
.
7.
A.
Ghani
and
W.
Polifke
, “
Control of intrinsic thermoacoustic instabilities using hydrogen fuel
,”
Proc. Comb. Inst.
38
,
6077
6084
(
2021
).
8.
E.
Æsøy
,
J. G.
Aguilar
,
S.
Wiseman
,
M. R.
Bothien
,
N. A.
Worth
, and
J. R.
Dawson
, “
Scaling and prediction of transfer functions in lean premixed H2/CH4-Flames
,”
Combust. Flame
215
,
269
282
(
2020
).
9.
T. F.
Guiberti
,
D.
Durox
,
L.
Zimmer
, and
T.
Schuller
, “
Analysis of topology transitions of swirl flames interacting with the combustor side wall
,”
Combust. Flame
162
,
4342
4357
(
2015
).
10.
S.
Taamallah
,
S. J.
Shanbhogue
, and
A. F.
Ghoniem
, “
Turbulent flame stabilization modes in premixed swirl combustion: Physical mechanism and Karlovitz number-based criterion
,”
Combust. Flame
166
,
19
33
(
2016
).
11.
R.
Mao
,
J.
Wang
,
W.
Zhang
,
Z.
An
,
W.
Lin
,
M.
Zhang
, and
Z.
Huang
, “
Effect of high hydrogen enrichment on the outer-shear-layer flame of confined lean premixed CH4/H2/air swirl flames
,”
Int. J. Hydrogen Energy
46
,
17969
17981
(
2021
).
12.
A.
Gentemann
,
C.
Hirsch
,
K.
Kunze
,
F.
Kiesewetter
,
T.
Sattelmayer
, and
W.
Polifke
, “
Validation of flame transfer function reconstruction for perfectly premixed swirl flames
,” in
ASME Turbo Expo 2004: Power for Land, Sea, and Air
(
ASME
,
Vienna, Austria
,
2004
).
13.
L.
Tay-Wo-Chong
,
T.
Komarek
,
R.
Kaess
,
S.
Föller
, and
W.
Polifke
, “
Identification of flame transfer functions from LES of a premixed swirl burner
,” in
Proceedings of ASME Turbo Expo 2010
(
ASME
,
Glasgow, UK
,
2010
), pp.
623
635
.
14.
E. R.
Hawkes
and
J. H.
Chen
, “
Direct numerical simulation of hydrogen-enriched lean premixed methane–air flames
,”
Combust. Flame
138
,
242
258
(
2004
).
15.
M. S.
Day
,
X.
Gao
, and
J. B.
Bell
, “
Properties of lean turbulent methane-air flames with significant hydrogen addition
,”
Proc. Combust. Inst.
33
,
1601
1608
(
2011
).
16.
N.
Zettervall
,
K.
Nordin-Bates
,
E.
Nilsson
, and
C.
Fureby
, “
Large eddy simulation of a premixed bluff body stabilized flame using global and skeletal reaction mechanisms
,”
Combust. Flame
179
,
1
22
(
2017
).
17.
A.
Felden
,
P.
Pepiot
,
L.
Esclapez
,
E.
Riber
, and
B.
Cuenot
, “
Including analytically reduced chemistry (ARC) in CFD applications
,”
Acta Astronaut.
158
,
444
459
(
2019
).
18.
D.
Laera
,
P.
Agostinelli
,
L.
Selle
,
Q.
Cazères
,
G.
Oztarlik
,
T.
Schuller
,
L.
Gicquel
, and
T.
Poinsot
, “
Stabilization mechanisms of CH4 premixed swirled flame enriched with a non-premixed hydrogen injection
,”
Proc. Combust. Inst.
38
,
6355
6363
(
2021
).
19.
W.
Polifke
,
W.
Geng
, and
K.
Döbbeling
, “
Optimization of rate coefficients for simplified reaction mechanisms with genetic algorithms
,”
Combust. Flame
113
,
119
134
(
1998
).
20.
A.
Abou-Taouk
,
B.
Farcy
,
P.
Domingo
,
L.
Vervisch
,
S.
Sadasivuni
, and
L.-E.
Eriksson
, “
Optimized reduced chemistry and molecular transport for large eddy simulation of partially premixed combustion in a gas turbine
,”
Combust. Sci. Technol.
188
,
21
39
(
2016
).
21.
J.
Si
,
G.
Wang
,
X.
Liu
,
M.
Wu
, and
J.
Mi
, “
A new global mechanism for MILD combustion using artificial-neural-network-based optimization
,”
Energy Fuels
35
,
14941
14953
(
2021
).
22.
J.
Armengol
,
O. L.
Maitre
, and
R.
Vicquelin
, “
Bayesian calibration of a methane-air global scheme and uncertainty propagation to flame-vortex interactions
,”
Combust. Flame
234
,
111642
(
2021
).
23.
G.
Staffelbach
,
L. Y. M.
Gicquel
,
G.
Boudier
, and
T.
Poinsot
, “
Large Eddy simulation of self excited azimuthal modes in annular combustors
,”
Proc. Combust. Inst.
32
,
2909
2916
(
2009
).
24.
B.
Franzelli
,
E.
Riber
,
L. Y.
Gicquel
, and
T.
Poinsot
, “
Large Eddy simulation of combustion instabilities in a lean partially premixed swirled flame
,”
Combust. Flame
159
,
621
637
(
2012
).
25.
S.
Hermeth
,
G.
Staffelbach
,
L. Y. M.
Gicquel
, and
T.
Poinsot
, “
LES evaluation of the effects of equivalence ratio fluctuations on the dynamic flame response in a real gas turbine combustion chamber
,”
Proc. Combust. Inst.
34
,
3165
3173
(
2013
).
26.
C.
Kraus
,
L.
Selle
, and
T.
Poinsot
, “
Coupling heat transfer and large eddy simulation for combustion instability prediction in a swirl burner
,”
Combust. Flame
191
,
239
251
(
2018
).
27.
A.
Ghani
and
T.
Poinsot
, “
Flame quenching at walls: A source of sound generation
,”
Flow, Turbul. Combust.
99
,
173
184
(
2017
).
28.
B.
Rochette
,
F.
Collin-Bastiani
,
L.
Gicquel
,
O.
Vermorel
,
D.
Veynante
, and
T.
Poinsot
, “
Influence of chemical schemes, numerical method and dynamic turbulent combustion modeling on LES of premixed turbulent flames
,”
Combust. Flame
191
,
417
430
(
2018
).
29.
Y.
Dagan
,
N. W.
Chakroun
,
S. J.
Shanbhogue
, and
A. F.
Ghoniem
, “
Role of intermediate temperature kinetics and radical transport in the prediction of leading edge structure of turbulent lean premixed flames
,”
Combust. Flame
207
,
368
378
(
2019
).
30.
D.
Brouzet
,
M.
Talei
,
M.
Brear
, and
B.
Cuenot
, “
The impact of chemical modelling on turbulent premixed flame acoustics
,”
J. Fluid Mech.
915
,
A3
(
2021
).
31.
R. M.
Solgi
, see https://pypi.org/project/geneticalgorithm/ for “
Geneticalgorithm 1.0.2
” (
2020
).
32.
University of California at San Diego
, see http://ucsd.edu for “
Chemical-Kinetic Mechanisms for Combustion Applications
” (
2016
).
33.
D. G.
Goodwin
,
H. K.
Moffat
, and
R. L.
Speth
,
Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes
, version 2.4.0, Caltech, Pasadena, CA, 2009.
34.
A.
Kazakov
and
M.
Frenklach
,
Reduced Reaction Sets Based on GRI-Mech 1.2
(
University of California at Berkeley
,
Berkeley, CA
,
1994
).
35.
A.
Abou-Taouk
,
S.
Sadasivuni
,
D.
Lörstad
, and
L.-E.
Eriksson
, “
Evaluation of global mechanisms for LES analysis of SGT-100 DLE combustion system
,” in
Volume 1B: Combustion, Fuels and Emissions
(
American Society of Mechanical Engineers
,
San Antonio, Texas
,
2013
), p.
V01BT04A036
.
36.
L.
Berger
,
A.
Attili
, and
H.
Pitsch
, “
Synergistic interactions of thermodiffusive instabilities and turbulence in lean hydrogen flames
,”
Combust. Flame
244
,
112254
(
2022
).
37.
E.
Æsøy
,
J. G.
Aguilar
,
M. R.
Bothien
,
N. A.
Worth
, and
J. R.
Dawson
, “
Acoustic-convective interference in transfer functions of methane/hydrogen and pure hydrogen flames
,”
j. Eng. Gas Turbines Power
143
,
121017
(
2021
).
38.
SID Software
,
Simcenter STAR-CCM+
(
Siemens
,
2021
).
39.
F.
Nicoud
and
F.
Ducros
, “
Subgrid-scale stress modelling based on the square of the velocity gradient tensor
,”
Flow Turbul. Combust.
62
,
183
200
(
1999
).
40.
N. W.
Bressloff
, “
A parallel pressure implicit splitting of operators algorithm applied to flows at all speeds
,”
Int. J. Numer. Methods Fluids
36
,
497
518
(
2001
).
41.
F.
Charlette
,
C.
Meneveau
, and
D.
Veynante
, “
A power-law flame wrinkling model for LES of premixed turbulent combustion Part I: Non-dynamic formulation and initial tests
,”
Combust. Flame
131
,
159
180
(
2002
).
42.
M. S.
Darwish
and
F. H.
Moukalled
, “
Normalized variable and space formulation methodology for high-resolution schemes
,”
Numer. Heat Transfer
26
,
79
96
(
1994
).
43.
M.
Karalus
,
D.
Brandt
,
A.
Brown
, and
V.
Lister
, “
A multi-timescale approach for the prediction of temperatures in a gas turbine combustor liner
,” in
Turbo Expo: Power for Land, Sea, and Air
(
Virtual
,
Online
,
2020
), Vol.
84164
, p.
11
.
44.
S.
Föller
and
W.
Polifke
, “
Advances in identification techniques for aero-acoustic scattering coefficients from large eddy simulation
,” in
18th International Congress on Sound and Vibration
(
Rio de Janeiro
,
Brazil
,
2011
), Vol.
4
, pp.
3122
3129
.
45.
T.
Lieuwen
, “
Modeling premixed combustion—Acoustic wave interactions: A review
,”
J. Propul. Power
19
,
765
781
(
2003
).
46.
S.
Ducruix
,
T.
Schuller
,
D.
Durox
, and
B.
Candel
, “
Combustion dynamics and instabilities: Elementary coupling and driving mechanisms
,”
J. Propul. Power
19
,
722
734
(
2003
).
47.
W.
Polifke
, “
Modeling and analysis of premixed flame dynamics by means of distributed time delays
,”
Prog. Energy Combust. Sci.
79
,
100845
(
2020
).
48.
X.
Wen
,
T.
Zirwes
,
A.
Scholtissek
,
H.
Böttler
,
F.
Zhang
,
H.
Bockhorn
, and
C.
Hasse
, “
Flame structure analysis and composition space modeling of thermodiffusively unstable premixed hydrogen flames—Part I: Atmospheric pressure
,”
Combust. Flame
238
,
111815
(
2022
).
Published open access through an agreement with Technische Universität München Fakultät für Maschinenwesen