Comparison of Combustion Models for Lifted Hydrogen Flames within RANS Framework

Abstract

Within the framework of a Reynolds averaged numerical simulation (RANS) methodology for modeling turbulence, a comparative numerical study of turbulent lifted H₂/N₂ flames is presented. Three different turbulent combustion models, namely, the eddy dissipation model (EDM), the eddy dissipation concept (EDC), and the composition probability density function (PDF) transport model, are considered in the analysis. A wide range of global and detailed combustion reaction mechanisms are investigated. As turbulence model, the Standard k-ε model is used, which delivered a comparatively good accuracy within an initial validation study, performed for a non-reacting H₂/N₂ jet. The predictions for the lifted H₂/N₂ flame are compared with the published measurements of other authors, and the relative performance of the turbulent combustion models and combustion reaction mechanisms are assessed. The flame lift-off height is taken as the measure of prediction quality. The results show that the latter depends remarkably on the reaction mechanism and the turbulent combustion model applied. It is observed that a substantially better prediction quality for the whole range of experimentally observed lift-off heights is provided by the PDF model, when applied in combination with a detailed reaction mechanism dedicated for hydrogen combustion.

1. Introduction

Current power generation techniques by means of thermal machinery [1] are mainly dependent on combustion. As the endeavors for utilizing new energy sources [2] and recovery technologies [3] continue, combustion is prognosed to remain as an important technology for energy conversion, also for the future. It should be noted that combustion plays an important role also for renewable energies, as it is the conversion technology for biomass [4].

Combustion of gases that contain hydrogen [5,6] occupies a significant role for clean and efficient energy supply. Hydrogen allows a useful way of storing the surplus power generated by photovoltaics and wind [7]. Moreover, in place of combustion [8,9], the gasification of waste, biomass, and coal [10,11] provides a useful procedure for efficient and clean conversion, whereby the product synthesis gas (syngas) contains rather considerable portions of hydrogen. Furthermore, there is a developing tendency in hydrogen production using nuclear energy [12]. Obviously, for the environmental aspect, the subsequent combustion of hydrogen is very welcome as it produces no carbon dioxide.

The use of hydrogen-containing gases as fuel in combustion equipment is a considerably demanding task. Hydrogen is very reactive and has quite different material properties in comparison to the gases that are commonly encountered in combustion applications. Thus, it can remarkably affect the properties of the gaseous mixture, even in small amounts. In premixed combustion, a potential complication related with the use of hydrogen blend fuels is the increased tendency of flashback [13]. The counter piece of flashback is blow-off [14]. The predecessor of blow-off is the lift-off that can occur at large jet speed, causing the flame root to leave the burner rim. After the lift-off, depending on the inlet velocity of the jet, a flame stabilized at a distance downstream the rim, i.e., a lifted flame, can be achieved [15].

Computational modeling of turbulent lifted flames that exhibit intricate interactions of turbulence with chemistry [16] is an ambitious undertaking, as the partial premixing occurring at the flame root leads to complex stabilization mechanisms [17].

For modeling the flow turbulence, the large eddy simulation (LES) technique [18,19] is increasingly being applied. However, its application in the industrial development processes is still very challenging because of its high requirements on computational resources. Therefore, the more traditional Reynolds averaged numerical simulation (RANS) approach [20] is still frequently adopted within this context. Given this perspective, the present investigation is focused on the RANS approach for turbulence modeling. Consequently, in the following review, the RANS-based approaches are addressed.

Experimental and numerical investigations of lifted flames of different fuels have been presented by different authors. For clarity, only the hydrogen flames are considered in the following review. One of the early experimental results for lifted hydrogen flames was provided by Cabra et al. [21], who measured a lifted H₂/N₂ flame in vitiated co-flow. A single flame (with a single co-flow temperature and a single lift-off height) was investigated [21]. Another configuration was experimentally investigated, later on, by Markides and Mastorakos [22]. In Ref. [21], a stable lifted flame was reported, whereas in Ref. [22], different regimes were measured, with emphasis on the unstable regime dominated by auto-ignition. Further experiments on lifted H₂/N₂ flames were performed by Wu et al. [23] and Gordon et al. [24] in a very similar configuration to that of Ref. [21].

The lifted H₂/N₂ flame of Cabra et al. [21] was also investigated numerically by the same authors, applying a composition probability density function (PDF) transport model [25] as well as the eddy dissipation concept (EDC) [26] as turbulent combustion models. The GRI-2.1 [27] was used as the detailed reaction mechanism (stripping out the carbon- and nitrogen-containing species and reactions). The overall performance of both of the combustion models were reported to be rather comparable [21]. The measured lift-off height of H/d = 10 (H: Lift-off height, d: Fuel jet diameter) was under-predicted by PDF and EDC as 7 and 8.5, respectively [21]. Later on, the same flame was computationally investigated by many investigators. Masri et al. [28] and Wu et al. [29] calculated the flame [21], adopting a composition PDF transport method as the turbulent combustion model. As reaction mechanism, the mechanisms of Mueller et al. [30] and GRI-2.1 [27] were used. They concluded that even though the mixing rates are still important, the flame is rather controlled by the chemical kinetics. Using the same basic modeling approach (PDF, GRI-2.1) as in Ref. [21], a larger lift-off height, H/d = 8, was predicted, which is closer to the experimental value. The Mueller et al. [30] mechanism led to an over-prediction, with H/d = 14 [28,29].

Lee and Kim [31] investigated the Cabra et al. [21] flame numerically by a direct-quadrature method of moments (DQMOM)-based PDF transport method as the turbulent combustion model, using the reaction mechanism of Mueller et al. [30]. With H/d = 5.7, their results indicated a stronger under-prediction of the lift-off height compared to the previous calculations of other authors [21,28,29]. They achieved a closer prediction (H/d = 9.8) to experiments by reducing the temperature of the co-flowing vitiated gases by 10 K. They justified this modification of the boundary condition with possible turbulent fluctuations of temperature that were not directly captured by the applied RANS approach. A quite similar modeling approach to the previous work [21,28,29] was applied by Mir Najafizadeh et al. [32] using a composition PDF transport turbulent combustion model based on the reaction mechanism of Li et al. [33], which is an updated version of the mechanism of Mueller et al. [30]. In difference to the previously mentioned studies [21,28,29,31], where the Standard k-ε model [20] was used, a Reynolds stress model (RSM) was applied as the turbulence model by Mir Najafizadeh et al. [32]. They could predict a lift-off height close to the experiments [21] by assuming (similar to the Ref. [31]) a lower temperature (by 5 K) for the co-flow, justifying this by the measurement uncertainty.

Mouangue et al. [34] applied Lagrangian intermittent modeling (MIL, Modèle Intermittent Lagrangien) [35] as turbulent combustion model to simulate the flame of Cabra et al. [21], in combination with the GRI-2.1 [27] as the reaction mechanism. Naud et al. [36] analyzed the Cabra et al. [21] flame numerically, using an unsteady flamelet progress variable (UFPV) approach [37] as turbulent combustion model based on the detailed reaction mechanism of Saxena and Williams [38]. As RANS turbulence model, RSM was used. For the experimental conditions, the lift-off height was substantially over-predicted (H/d ≈ 16). A close prediction to the experimental value (H/d = 10, [21]) could be obtained by increasing the co-flow temperature by 8 K and 17 K (depending on the alternative definitions of the progress variable within the UFPV model). The conditional moment closure (CMC) [39] was applied as the turbulent combustion model by De et al. [40] to simulate the flame of Cabra et al. [21], using the Mueller et al. [30] reaction mechanism along with the Standard k-ε as the RANS turbulence model. The lift-off height was under-predicted by H/d = 6.8. In a quite recent study of Larbi et al. [41], Eulerian (DQMOM) and Lagrangian composition PDF transport models were used as turbulent combustion models, while only Lagrangian PDF models were adopted in the majority of the above-stated studies. The effect of model constants was studied. As turbulence model, a modified k-ε model [41] and RSM was applied, where the former was performing better than the latter. As chemical reaction mechanism, GRI-2.1 [27] was used. A remarkable under-prediction of the lift-off height was reported without an explicit quantification.

The lifted H₂/N₂ flame measured by Wu et al. [23], which is very similar to that of Cabra et al. [21], was computationally analyzed by Cao et al. [42]. As turbulent combustion model, the joint velocity-turbulence frequency-composition PDF method was applied. A detailed study of the mixing models was performed. As RANS turbulence model, the Standard k-ε model was used. As chemical reaction mechanism, the mechanisms of Mueller et al. [30] and Li et al. [33] were applied, concluding that the latter predicts shorter lift-off heights than the former, especially for low co-flow temperatures. The sensitivity of the results onto details of boundary conditions (profile shapes, turbulence) was found to be low, except the strong sensitivity to the co-flow temperature. A quite good overall agreement of the predicted lift-off heights with the measured values was observed.

In the above overview, it can be seen that various and quite mature models are applied to calculate turbulent lifted hydrogen flames. Nonetheless, to the assessment of the present authors, there is still demand for further analysis. The goal of the current investigation is the presentation of a “coherent” validation study, for a series of different turbulent combustion modeling approaches within a wide range. Using the same methods for all other aspects of mathematical and numerical modeling, a more representative and direct comparison of the combustion modeling approaches is aimed for, unmasked by the differences in other aspects of the overall mathematical and numerical modeling, such as the numerical grids, discretization schemes. This coherent and comprehensive validation study on the combustion modes, based on the same approaches in all other aspects of the mathematical and numerical modeling, on the same and widely used software platform is assumed to be additionally valuable for the researchers in the field.

2. Problem Definition

2.1. Non-Reacting H₂/N₂ Jet

For determining the adequate turbulence model, a preliminary analysis was performed on a non-reacting test case, which has allowed the validation of the turbulence models isolated from the turbulent combustion models. Sautet and Stepowski [43] performed measurements on non-reacting H₂/N₂ jets in co-flow of air, in a configuration similar to that of the flame considered in the main body of the present study. The flow at ambient conditions is isothermal, but exhibits density variations due to the changes in the composition. A sketch of the configuration is depicted in Figure 1.

The considered experimental conditions are listed in

Table 1. Experimental conditions for non-reacting H₂/N₂ jet [43] (V: Velocity; T: Temperature; d, D: Diameters of central and annular jets; Re: Reynolds number; X_i: Mole fractions).

	Central Jet	Co-flow
V (m/s)	45	4.5
T (K)	300	300
d, D (mm)	10	112
Re	8700	27,000
XH2	0.50	0
XO2	0	0.21

.

2.2. Lifted H₂/N₂ Jet Flame

As the main test case, the atmospheric, lifted flame of a H₂/N₂ jet in vitiated co-flow was considered, which was investigated experimentally by Wu et al. [23]. A sketch of the configuration is presented in Figure 2.

The fuel jet is positioned at the center of a disk with perforations, which supports a large number of lean premixed flames that provide a hot co-flowing stream. The central fuel jet protrudes about 70 mm over the perforated disk. Thus, it can be assumed that the central fuel jet issues into a homogeneous composition of the vitiated co-flow. The latter consists of the combustion products of an H₂/Air mixture, which delivers the inlet boundary conditions for the calculation of the lifted flame. The test conditions for two typical cases for different co-flow conditions (Flame A and Flame B) with different H₂/Air ratios are represented in

Table 2. Experimental conditions for lifted H₂/N₂ jet flame [23] (V: Velocity; T: Temperature; d, D: Diameters of central and annular jets; Re: Reynolds number; Q_i: Volume flow rates).

	Central Jet	Co-flow
V (m/s)	110	4
d, D (mm)	4.57	190
		Flame A	Flame B
Re		18,965	19,200
T (K)	317.5	1044	1013
QH2 (slm)	25	190.1	185.2
QN2 (slm)	75	-	-
QAIR (slm)	-	1720	1720

.

In the present study, seven cases with seven different co-flow temperatures were considered (1010, 1013, 1020, 1030, 1040, 1044, and 1050 K). The different co-flow temperatures result from different H₂/Air ratios of the co-flow. The corresponding compositions of the vitiated co-flow were calculated using the chemical kinetics code Cantera [44]. In doing so, the following procedure was applied: With the experimentally given flow rates of H₂ and Air (Table 2), the complete combustion of the mixture was calculated assuming a well-stirred reactor. First, an adiabatic reactor was considered. In such calculations, the resulting exhaust gas temperatures turned out to be larger than the measured ones (Table 2). This is attributed to the heat loss in the experimental combustion device. Thus, in the well-stirred reactor calculations, a sufficient amount of heat loss was considered to obtain the experimental exhaust gas temperatures. The resulting exhaust gas composition was then taken as the composition to be used as the inlet boundary condition for the subsequent (main) calculation of the lifted jet flame. The considered co-flow temperatures and the corresponding gas compositions (calculated as mentioned above) are displayed in

Table 3. Considered co-flow temperatures and compositions for lifted H₂/N₂ jet flame [23] (T: Temperature; X_i: Mole fractions).

T (K)	XO2	XH2O
1010	0.1489	0.1005
1013	0.1487	0.1008
1020	0.1483	0.1016
1030	0.1477	0.1028
1040	0.1470	0.1039
1044	0.1468	0.1043
1050	0.1464	0.1050

.

In the study of Cabra et al. [21], a single flame for a single co-flow temperature (1045 K) was published that was rather high to lead to a rather moderate lift-off height (H/d = 10). A different and interesting feature of the presently considered experiments of Wu et al. [23] is the investigation of a number of flames for a range of co-flow temperatures (1010–1050 K). The lift-off height increases with decreasing co-flow temperatures and the focus of the present study is the predictability of this relationship. For the present class of flames, autoignition and flame propagation, both, were reported to play a role [21,22,23,24]. For the flame of Ref. [21] with a moderate lift-off height, it was postulated that the transient startup was governed by autoignition, where the subsequent flame stabilization was affected by flame propagation in partial premixing. However, for varying co-flow temperatures, a shift in the relative importance of the mechanisms can be expected. Obviously, both processes, i.e., the autoignition and flame propagation (laminar flame speed), are strongly temperature-dependent, correlating with the observed [23] influence of the co-flow temperature (where the differential/preferential diffusion can also play role for the laminar flame speed of hydrogen flames). On the other hand, the temperature and composition of the burnable mixture that influence the chemical kinetics depend on the mixing that is strongly affected by the turbulent flow field.

3. Modeling

In this section, the modeling details are described in two sub-sections. In the first sub-section, an outline of the complete mathematical and numerical modeling is given. Combustion modeling is described separately in the subsequent sub-section.

3.1. Outline of the Mathematical and Numerical Modeling

An ideal gas mixture with Newtonian behavior was assumed. Possible buoyancy effects were omitted. This was justified by the prevailing high Froude numbers [45], indicated by the inlet conditions. Heat transfer due to radiation [46] was also omitted, assuming an adiabatic system. The molecular properties were modeled accurately. The specific heat capacity of each species was calculated via a pair of fourth order polynomials of temperature [47] (for low and high ranges of temperature). The molecular transport properties such as the viscosity, thermal conductivity, and the binary diffusion coefficients were obtained from the kinetic theory [14,48]. The multi-component diffusion in the mixture was treated by assuming Fickian diffusion, assuming “effective” binary diffusion coefficients [14]. This was considered to be a sufficiently accurate approach for mixtures having a species in excess [49], which was nitrogen in the present case. For some of the turbulent combustion models, an equality of the molecular Schmidt numbers was assumed, as will be indicated below.

As discussed above, turbulence was modeled within a RANS framework. Since the considered geometries and boundary conditions were axi-symmetric, the sought steady-state RANS solution also needed to be axi-symmetric. Thus, a two-dimensional axi-symmetric formulation was adapted.

As RANS turbulence models, two-equation turbulent viscosity models were considered [20]. Although higher order closures such as RSM [20] are potentially more accurate, they were not considered, since their application in the industrially relevant problems remained restricted due to the increased computational overhead and convergence problems encountered frequently. Among the turbulent viscosity-based turbulence models, those relying on the specific dissipation rate (omega, ω) are being increasingly used [50]. Nonetheless, since the present flow was of free-shear type, the turbulent viscosity models that are based on the dissipation rate (epsilon, ε) were considered only. The Standard k-ε [51], the RNG k-ε [52] (RNG: Renormalization group theory), and the Realizable k-ε models [53] were used. For modeling the turbulent fluxes of scalar quantities, the gradient–diffusion assumption was employed. In doing so, constant turbulent Prandtl-Schmidt numbers were assumed, assigning the value of 0.85 for the energy equation, and the value 0.7 for the remaining scalar transport equations. At the inlet boundaries, the boundary conditions of the turbulence variables were always determined by assuming a turbulence intensity of 4% (along with a locally isotropic turbulence) and a macro length scale based on the jet diameters [20].

The SIMPLEC pressure correction method [54] was used to treat the coupling of the velocity and pressure fields. For the discretization of the convective terms, the QUICK algorithm [55] was employed with a third order formal accuracy. The multi-purpose, finite volume method-based CFD code ANSYS Fluent 18.0 [56] was employed in the computational analysis. No special ad-hoc modifications tailored to the present type of problem were applied in the turbulence and turbulence combustion models. Thus, the standard settings [56] should be assumed for any model, unless the contrary is explicitly stated.

3.2. Combustion Modeling

3.2.1. Reaction Mechanisms

Two single-step global reaction mechanisms were considered (that include only the main species, H₂, O₂, H₂O). These are the mechanism of Kudriakov et al. [57] (KU) and the mechanism of Marinov et al. [58] (MA). The former and latter assume an irreversible and a reversible reaction, respectively.

As the detailed reaction mechanisms, four mechanisms were considered. Each of them comprise about 20 elementary reactions (without counting the reverse ones) between the eight main and intermediate species, namely, H₂, O₂, H₂O, H, O, OH, HO₂, and H₂O₂. The considered detailed mechanisms are the Gri-Mech 2.11 [27] (G2), Gri-Mech 3.0 [59] (G3), the mechanisms of Conaire et al. [60] (CO), Keromnes et al. [61] (KE), and Li et al. [33] (LI). Gri-Mech 2.11 and Gri-Mech 3.0 are mechanisms that are optimized to model natural gas combustion, whereas the mechanism of Keromnes et al. is designed for syngas combustion. In using these mechanisms, the carbon- and nitrogen-containing species and reactions are stripped off. The mechanisms of Conaire et al. and Li et al. were developed specifically for hydrogen combustion.

A reversible reaction with N species can be represented as:

$${\displaystyle \sum _{i=1}^N{\nu }_i^{\prime }{M}_i \underset{k_b}{\overset{k_f}{\rightleftarrows }}} {\displaystyle \sum _{i=1}^N{\nu }_i^{″}{M}_i},$$

(1)

where $M_i$ and ${\nu }_i$ stay for the species and the related stoichiometric coefficients, while k_f and k_b represent rate coefficients of the forward and backward reactions, respectively.

According to the law of mass action, using the rate coefficients provided by the respective chemical reaction mechanism, the molar conversion rate of a species i according to chemical kinetics (R_K,i) can be calculated as [14]:

$$R_{K,i}=\left({\nu }_i^{″}-{\nu }_i^{\prime }\right)\left(k_f{\displaystyle \prod _{j=1}^N{C}_j^{n_j^{\prime }}-k_b{\displaystyle \prod _{j=1}^N{C}_j^{n_j^{″}}}}\right),$$

(2a)

where C_i stay for the species concentrations. The exponents n_i are equivalent to ${\nu }_i$ for elementary reactions, as they may be different for empirical global reaction mechanisms. The rate coefficients are obtained by the Arrhenius equation [14], which is expressed for k_f below.

$$k_f=A{T}^m\exp \left(-\frac{E_a}{R_{u}T}\right),$$

(2b)

where A, E_a, R_u, m denote the pre-exponential factor, activation energy, universal gas constant, temperature exponent, respectively, while T stays for the thermodynamic temperature.

Source terms of species transport equations are given by above expressions (Equation (2)). For a turbulent flow, within a RANS formulation, where time-averaged equations are solved, it is obvious that severe closure problems emerge due to the strong non-linearity of the source terms. Therefore, the so-called “turbulent combustion” models are needed that take the interaction between turbulence and chemistry, in an adequate manner into account. The presently considered turbulence combustion models are outlined below.

3.2.2. Turbulent Combustion Models

The presently considered turbulent combustion models are outlined in this subsection. Note that such models are, in general, quite complex, and details in their implementation may be influential on the results. The considered models refer to the versions as they are implemented in the used software ANSYS Fluent 18.0 [56].

Eddy Dissipation Model with Kinetics (EDM+K)

The time-averaged transport equations for each species are solved. The eddy dissipation model (EDM) assumes a mixing controlled combustion and calculates time-averaged conversion rates of the species based on the dissipation rate of turbulence eddies [62]. Here, the time-averaged molar conversion rate (R_EDM,i) is calculated from:

$$R_{EDM,i}={\nu }_i\alpha \rho \frac{\epsilon }{k}\min \left(\min \left(\frac{Y_R}{{\nu }_R^{\prime }{M}_R}\right),\beta \frac{{\displaystyle \sum _P{Y}_P}}{{\displaystyle \sum _j^N{\nu }_j^{″}{M}_j}}\right),$$

(3)

where $\rho$, k, $\epsilon$ stay for the averaged mixture density, turbulence kinetic energy, and its dissipation rate, respectively. The terms Y_i and M_i denote the averaged mass fraction of species i and its molar mass, respectively. The subscript R denotes each of the reacting species, whereas the subscript P refer to the sum of the product species. Coefficients $\alpha$, $\beta$ are empirical model constants. Standard values are used [62].

The chemical kinetics effects (K) are considered approximately, obtaining the conversion rates (R_K,i) from the Arrhenius expressions of the underlying chemical reaction scheme by neglecting the fluctuations of the thermochemical variables, comparing the both rates and assigning the smaller one as the effective conversion rate [56]. Thus, the resulting time-averaged conversion rate is obtained from:

$$R_{EDM+K,i}=\min \left(R_{EDM,i},R_{K,i}\right).$$

(4)

Due to the rather ad-hoc consideration of the chemical kinetics, the model makes sense only in combination with global reaction schemes. Therefore, this method is applied to Kudriakov et al. [57] (KU) and Marinov et al. [58] (MA) schemes, only. It is also to be noted that the MA scheme comprises a backward reaction, which is deactivated when applied in combination with this (EDM+K) turbulent combustion model.

Eddy Dissipation Concept (EDC)

The eddy dissipation concept (EDC) has evolved from the above-mentioned eddy dissipation idea [63] and exhibits slightly different versions that emerged in time. It can model the combined mixing and kinetics effects in a sounder manner compared to EDM+K. As the main modeling concept, it is assumed that reactions occur within small structures that are called fine-scales. The time-averaged source term (S_EDC,i) of the species mass fraction transport equation is obtained from:

$$S_{EDC,i}=\rho \frac{1}{{\tau }^{*}}\frac{{\xi }^{{*}^2}}{1-{\xi }^{{*}^3}}\left(Y_i^{*}-Y_i\right).$$

(5)

In Equation (5), the variables ${\xi }^{*}$ and ${\tau }^{*}$ represent the length fraction and time scale of the fine structures that relate to the volume and residence time of the assumed fine-scale reactors. These quantities are estimated from the velocity and length scales of turbulence. The variable $Y_i^{*}$ denotes the fine-scale species mass fractions that are obtained after a reaction time of ${\tau }^{*}$ by solving chemical reactor equations. The version of the model that underlies the used software [56] corresponds to that of Gran and Magnussen [63]. Still, there are some deviations from the original model [63] in the implementation [56], such as the use of batch reactors to model the fine structures instead of perfectly stirred ones. In the present calculations, the direct integration method was applied in combination with the stiff chemistry solver [56] instead of making use of in situ adaptive tabulation (ISAT) [56]. The time-averaged transport equation was solved for each species, closing the source term as outlined above. Due to the rather sound modeling of turbulence-chemistry interaction, the model is generally observed to lead to good results in combination also with detailed reaction schemes.

Flamelet Generated Manifold (FGM)

According to the flamelet generated manifold (FGM) method [56,64], it is assumed that the realized trajectories on the thermochemical manifold of a turbulent flame can approximately be represented by the scalar evolution similar to that of a laminar flame. Under similar assumptions to the laminar flamelet method [65], the elemental mass fractions are associated with a single scalar quantity called mixture fraction (Z) that is governed by a source-free differential transport equation. The state of the reaction is described by an additional scalar variable, the so-called reaction progress variable (C) [64]. Thus, the thermochemical dependent variables can be expressed as unique functions of Z and C:

$$\varphi =\varphi \left(Z,C\right),$$

(6)

The functional relationships $\varphi$ (Equation (6)) are obtained from one-dimensional laminar flame calculations that are performed in advance, and prepared in tabulated form for the subsequent CFD analysis. In the present application, the one-dimensional laminar flame calculations were performed based on a non-premixed (counterflow) configuration. About 100 different values of strain are considered, covering the whole spectrum from negligibly low strain up to extinction. Assuming that Z and C are statistically independent, the Favre-averaged values of the thermo-chemical variables ($\tilde{\varphi }$) are obtained by integrating the relationships over the Z and C spaces, with the help of presumed probability density functions (Beta functions are used):

$$\tilde{\varphi }={\displaystyle \iint \varphi \left(Z,C\right)} \beta \left(Z;\tilde{Z},\widetilde{Z{″}^2}\right) \beta \left(C;\tilde{C},\widetilde{C{″}^2}\right)dZdC.$$

(7)

For the construction of the beta functions, modeled transport equations are solved for mean values ($\tilde{Z}$, $\tilde{C}$) and variances ($\widetilde{Z{″}^2}$, $\widetilde{C{″}^2}$) of Z, C [56,66,67]. The source term of the modeled $\tilde{C}$ equation is extracted from the one-dimensional laminar flame data like any other thermochemical variable. In constructing the beta functions, the mixture fraction and reaction progress variable space, each, are resolved by 41 points, while 21 points are used for the variances. The attractiveness of the method lies therein that it can incorporate detailed chemistry in a CFD calculation of turbulent combustion, with comparably quite low computational costs, since it requires the solution of only four differential transport equations (mean values and variances of Z and C), irrespective of the complexity of the underlying chemical reaction scheme. On the other hand, an assumption underlying the model is the equality of the molecular Schmidt numbers of the species transport equations (which is not necessary for EDM and EDC). Thus, the effects of differential and preferential diffusion [13,14] that are known to be important for flames of hydrogen containing fuels cannot be considered.

Composition PDF Transport (PDF)

To determine the averaged values of the thermochemical quantities, a single-point, joint probability density function (PDF) can be obtained from its transport equation. The joint probability density function represents the time fraction that the mixture spends at each species, temperature, and pressure state, leading to an (N + 2) dimensional joint probability density function for N species, temperature, and pressure spaces. The transport equation for the Favre composition probability density function (P) can be derived from the governing equations as:

$$\frac{\partial }{\partial t}\left(\rho P\right)+\frac{\partial }{\partial x_i}\left(\rho u_{i}P\right)+\frac{\partial }{\partial {\Psi }_k}\left(\rho S_{k}P\right)=-\frac{\partial }{\partial x_i}\left[\rho \langle u_i^{"}|\Psi \rangle P\right]+\frac{\partial }{\partial {\Psi }_k}\left[\rho \langle \frac{1}{\rho }\frac{\partial J_{i,k}}{\partial x_i}|\Psi \rangle P\right].$$

(8)

In Equation (8), $u_i$ and ${u_i}^{″}$ denote the Favre mean and fluctuation values of the velocity vector. $S_k$ is the reaction rate for species k, while $\Psi$ and $J_{i,k}$ denote the composition space and molecular diffusion flux vectors, respectively. The terms $\langle A|B\rangle$ express the conditional probability of event A, given that event B occurs. The terms on the left-hand side are closed. The third term (which is closed) is the change due to chemical reactions, which makes up the strength of the method [25,56]. The right-hand side terms need to be modeled. The first term on the right-hand side, i.e., the turbulent flux term, is approximated by a gradient-diffusion approximation. The second right-hand side term is due to molecular mixing/diffusion. For its modeling, the Euclidian minimum spanning tree (EMST) model is applied [68], which is known to be more accurate than the alternatives models [56] such as the modified curl (MC) and interaction by exchange with the mean (IEM). Reaction source terms are obtained by direct integration, which was affordable in the present two-dimensional application. A Lagrangian Monte-Carlo method was adopted to solve the PDF transport equation. The number of particles per cell was chosen to be 50, which was shown by Cao et al. [42] to lead to a deviation of around 1% for the main species compared to the calculation with 100 particles per cell.

4. Results

4.1. On Solution Domains, Boundary Conditions, Grids

The two-dimensional axisymmetric solution domains of the non-reacting H₂/N₂ jet (Section 2.1) and the lifted H₂/N₂ flame (Section 2.2) have the same topology. The domains have a rectangular shape in the plane of axial (x) and radial (r) coordinates, which is depicted in Figure 3, with indications of the domain size and boundary types.

At the outlet boundary, a constant static pressure is prescribed, along with zero-gradient conditions for the remaining variables. For the lifted flame (Section 2.2), the left boarder of the domain boundary is completely covered by the central jet and co-flow inlet boundaries. For the isothermal case (Section 2.1), the co-flow inlet extends up to a radial position of approximately 6d (indicated by a dashed line in the figure), where the remaining outer part (ambient inlet, Figure 3) is defined to be a pressure boundary that allows an inflow, i.e., suction of ambient air by the ejector effect.

At the central jet and co-flow inlet boundaries, uniform profiles are prescribed in accordance with the measured values (Table 1, Table 2 and Table 3). A minor influence of the profile shapes was reported by Cao et al. [42] in a similar study. Computational grids were generated as structured, rectangular grids with axial and radial clustering of the cells near the jet inlet.

4.2. Non-Reacting H₂/N₂ Jet

For obtaining a sufficiently accurate grid resolution, a grid independence study was conducted. In doing so, the Standard k-ε model was employed as the turbulence model. Figure 4 shows the predicted length (L) of the potential core (non-dimensionalized by d) as a function of grid resolution, where N represents the total number of grid nodes. It can be observed that a satisfactory grid independence was achieved for N ≥ 5000. In the subsequent computations for validating turbulence models, the grid with the highest resolution, i.e., the one with approximately 16,000 nodes, was used.

The profiles of the calculated half value radius (δ) along the jet axis were compared with the measurements [43] in Figure 5. It can be observed that the results obtained by the Standard k-ε model show a quite well agreement with the measured values, which is better than those of the Realizable and RNG versions, for the current isothermal, variable density jet system. Thus, the Standard k-ε model was selected as the turbulence model for the subsequent study of the lifted H₂/N₂ flame.

4.3. Lifted H₂/N₂ Jet Flame

A grid independence study was conducted also for the lifted H₂/N₂ jet flame. In Figure 6, the profile of the calculated centerline temperature (T) at ten jet diameters downstream the jet inlet (x/d = 10) for six different grids with varying number of total nodes (N) (obtained by the KU-EDM+K) approach) is presented. It can be observed that a satisfactory grid independence is obtained for the finer grids. In the further calculations, the finest grid, having approximately 16,000 nodes, was used. This complies well with the previous studies, where a grid independence was reported even for coarser grids with about 10,000 nodes.

Not all turbulent combustion models were used in combination with all reaction mechanisms. An overview of the investigated combinations is provided in

Table 4. Investigated combinations of turbulent combustion models and reaction mechanisms.

	KU	MA	G2	G3	CO	KE	LI
EDM+K	X	X	−	−	−	−	−
EDC	X	X	X	X	X	X	X
FGM	−	−	X	X	X	X	X
PDF	−	X	X	X	X	X	X

Note: KU: Kudriakov et al. [57]; MA: Marinov et al. [58]; G2: Gri-Mech 2.11 [27]; G3: Gri-Mech 3.0 [59]; CO: Conaire et al. [60]; KE: Keromnes et al. [61]; LI: Li et al. [33].

, where a cross (X) means that the combination was investigated, and a minus (−) means that the combination was not considered.

As mentioned before, the models were assessed based on their prediction accuracy of the lift-off height as a function of co-flow temperature. In many of the previous computational investigations (e.g., [42]) of similar cases, the flame root was marked in terms of the OH mass fraction. Since this cannot be applied while using a global/reduced reaction mechanism without comprising OH species, a temperature-based criterion was used in the present study, which is more general, in this sense. The flame root location, i.e., the lift-off height (H) was defined to be given by the position, where the local Favre-averaged temperature exceeds the co-flow temperature by 5 K (at any radius). This is exemplarily illustrated in Figure 7, where the predicted temperature fields by the G2-PDF approach are displayed for two co-flow temperatures, namely, for T_CO = 1050 K (left) and T_CO = 1010 K (right). In the figure, a sub-region in the nearfield of the central jet is displayed, with a size of 30d × 10d. According to the present approach, the flame root location is marked by the attainment of the local temperature values of 1055 K and 1015 K, for T_CO = 1050 K and T_CO = 1010 K, respectively. The corresponding isotherms T = 1055 K (left) and T = 1015 K (right) are also indicated in the figure. Thus, the lift-off height (H) is defined to be the distance between the inlet boundary and the nearest point of the corresponding isotherm to the inlet, as illustrated in Figure 7.

No flame lift-off could be predicted by the global reaction mechanism of Kudriakov et al. In combination with both turbulent combustion models (KU-EDM+K, KU-EDC), the predicted flame remained attached to the burner rim for all values of the co-flow temperature within the range 1010–1050 K.

The results obtained by the global reaction mechanism of Marinov et al. in combination with different turbulent combustion models are shown in Figure 8, where the lift-off height (H) is presented as a function of the co-flow temperature. The measured values by Wu et al. are also shown in the figure.

One can see that the strong dependence of the experimental lift-off height on the co-flow temperature is under-predicted (Figure 8). In combination with all three combustion models, a quite weak increase of H with decreasing T_CO is predicted. The slope of the MA-EDC curve is closer to that of the experimental curve, compared to MA-EDM+K and MA-PDF curves, especially for the higher values of T_CO, where a rather good quantitative agreement with the experiments and MA-EDC curve is observed.

The lift-off heights predicted by the detailed reaction mechanism GRI-Mech 3.0 in combination with different turbulent combustion models are compared with the measurements of Wu et al. in Figure 9.

Using the eddy dissipation concept as the turbulent combustion model (G3-EDC), a stable lifted flame could not be obtained for the whole range of co-flow temperatures. A blow-off was always observed. Thus, a G3-EDC curve is missing in the figure. Applying the composition transport PDF method (G3-PDF), a lifted flame with a very large lift-off height of about H/d ≈ 55 is predicted for the highest co-flow temperature of 1050 K, which increases to H/d ≈ 75 for T_CO = 1040 K. For lower values of T_CO, the flame is blown-off (or, moves out of the solution domain). The predicted values (H/d > 55) for T_CO > 1040 K are obviously very much higher than the measured values.

Applying the flamelet generated manifold method as the turbulent combustion model (G3-FGM), the predicted value agrees very well with the measured value (H/d ≈ 5) for high T_CO (T_CO > 1040 K) (Figure 9). However, the increase of the lift-off height with decreasing co-flow temperature is not predicted, as the lift-off height remains rather unaffected by T_CO (even showing a slight reverse trend for some values of T_CO). A similar behavior is observed also for the remaining detailed reaction mechanisms using the FGM as the turbulent combustion model, where the model cannot predict the dependence of H onto T_CO. This may be seen to imply that the assumptions underlying FGM lose their validity for conditions prevailing for decreasing co-flow temperatures. The assumption of species Schmidt numbers to be equal (and Lewis number being unity), which is an essential part of FGM, but not required for the other turbulent combustion models, may also be a cause of the inaccuracy. Furthermore, the present approach of generating the one-dimensional laminar flame data based on a non-premixed configuration may be less adequate for large lift-off heights, with increasing degree of premixing. Moreover, the lack of ability of the presently applied procedure (FGM) for considering the autoignition effects can be seen as a further reason for the observed inferior performance.

The results obtained by the detailed reaction mechanism GRI-Mech 2.11 in combination with different turbulent combustion models are presented in Figure 10, where the experimental results of Wu et al. are also displayed.

G2-EDC and G2-PDF show a qualitative agreement with the measurements, where they under-predict the measured values quantitatively. The under-prediction increases with decreasing co-flow temperatures, i.e., increasing lift-off heights. The G2-PDF results agree quite well with the measurements for high T_CO (T_CO > 1040), also showing a quantitatively better agreement with the experiments compared to the G2-EDC ones.

The observed difference in the predictions of the Gri-Mech 3.0 (G3) (Figure 9) and the Gri-Mech 2.11 (G2) (Figure 10) is quite substantial. This, and an even better performance of the earlier version of the mechanism (Gri Mech 2.11), were not necessarily expected. The modifications in the H₂-related reactions in between the versions seem to imply a strong difference for H₂ combustion. The observed relative performance of the schemes (Figure 9 vs. Figure 10) can be seen to be in line with the published comparisons [69], in which the premixed hydrogen-air flame speeds by Gri-Mech 2.11 were found to be larger than those of Gri-Mech 3.0. On the other hand, it should be noted that the observed performance of Gri-Mech 3.0 (which was optimized for natural gas) for hydrogen flames cannot fully be attributed to the mechanism itself, but should be seen in combination with the applied RANS framework, which is intrinsically less accurate in resolving flame structure/heat release patterns compared to higher order methods such as LES.

The predictions based on the detailed reaction mechanism of Conaire et al., in combination with various turbulent combustion models, are compared with the experimental results of Wu et al. in Figure 11.

The results obtained by the composition PDF transport model (CO-PDF) as the turbulent combustion model show a qualitative agreement with the measured values. Nonetheless, quantitatively, the experimental values are highly over-predicted (Figure 11).

The results obtained by the eddy dissipation concept (CO-EDC) (Figure 11) show a qualitatively different behavior compared to that observed for Gri-Mech 2.11 (Figure 10). The G2-EDC curve shows a continuous increase of H with decreasing T_CO (Figure 10). CO-EDC exhibits two regions with different qualitative behaviors (Figure 11). For high T_CO (T_CO > 1040 K), i.e., low values of H, CO-EDC shows a slowly increasing H with decreasing T_CO, and a good quantitative agreement with the measurements. However, for lower T_CO values (T_CO < 1400), a very rapid, abrupt increase of H is observed with decreasing T_CO (with a remarkable over-prediction of the measured values), which quickly leads to blow-off beyond T_CO < 1030 K (Figure 11).

The lift-off heights that are predicted by the reaction mechanism of Keromnes et al. using different turbulent combustion models are compared with the experimental results of Wu et al. in Figure 12.

Qualitatively, the predictions obtained by the Keromnes et al. mechanism (Figure 12) show a similar behavior to that observed for the Conaire et al. mechanism (Figure 11). Quantitatively, the predicted values by EDC (KE-EDC) and PDF (KE-PDF) turbulent combustion models agree better with the experiments (Figure 12) than the CO-EDC and CO-PDF do (Figure 11). The KE-PDF curve predicts lower H values compared to the CO-PDF curve and is much closer to the experimental curve (Figure 12). For T_CO > 1030 K, the predicted values by KE-PDF agree rather well with the experimental ones. For T_CO < 1030 K, a remarkable over-prediction by KE-PDF is still observed. The results by EDC model (KE-EDC) (Figure 12) are qualitatively similar to those of CO-EDC (Figure 11), showing first a shallow increase with decreasing T_CO, followed by a rapid increase and blow-off by a further decrease of T_CO. Compared to CO-EDC (Figure 11), the characteristic T_CO values are lower for KE-EDC (Figure 12). The border between the two regimes is about T_CO ≈ 1030 K and the blow-off occurs for T_CO < 1020 K (the values were 1040 K and 1030 K for CO-EDC).

The results obtained by the detailed reaction mechanism of Li et al. in combination with different turbulent combustion models are displayed in Figure 13, where the lift-off height (H) is presented as a function of the co-flow temperature. The measured values by Wu et al. are also displayed in the figure.

The predictions obtained by the mechanism of Li et al. (Figure 13) using different turbulent combustion models are qualitatively similar to those obtained by the Keromnes et al. mechanism (Figure 12). LI-FGM and LI-EDC curves (Figure 13) are also quantitatively very close to KE-FGM and KE-EDC (Figure 12), respectively. The results obtained by the composition PDF transport model (LI-PDF) are (Figure 13), quantitatively, closer to the experimental curve compared to the results by KE-PDF (Figure 12). Among all of modeling approaches considered, the Li et al. mechanism [33] combined with the composition PDF transport model [25] provides the best overall agreement with the experiments for the lift-off height (Figure 13). However, the prediction accuracy by LI-PDF still leaves much to be desired, as measured lift-off height is remarkably over-predicted, especially for rather low T_CO values.

As an attempt for a better understanding of the mechanisms and performance of the modeling approaches, the important parameters of chemical kinetics were also analyzed. For different mixtures of the central jet and co-flow streams at different equivalence ratios (ϕ), the ignition delay time τ_IGN, and the laminar flame speed S_L (undisturbed, planar flame) were calculated using the above-mentioned detailed reaction mechanisms, using the chemical kinetics code Cantera [44]. Note that the calculations hold for premixed conditions for each equivalence ratio, and, thus, the results (Figure 14) can only be used as indications for the present partially premixed situation.

Figure 14 illustrates the calculated ignition times (Figure 14a) and laminar flame speeds (Figure 14b) as a function of the equivalence ratio, for two different temperatures of the co-flow stream, i.e., for T_CO = 1010 K and T_CO = 1050 K, spanning the whole considered range. Note that for a given co-flow temperature, the temperature of the mixture changes remarkably with changing equivalence ratio, due to the large difference between the jet and co-flow temperatures.

One can see that the ignition delay times are quite long for ~0.2 < ϕ and ϕ < ~0.005, due to the too low mixture temperatures for the former, and too low fuel mass fractions for the latter (Figure 14a). Within the region ~0.005 < ϕ < ~0.2, finite ignition delay times are predicted in all cases, mainly due to the high mixture temperature although the mixture is very fuel-lean. A rough estimation of the convective residence time can be done based on the co-flow inlet velocity and domain length, which gives a value of about 100 ms. This is within the range of the ignition delay times that occur in Figure 14a, and indicates the likeliness of autoignition for the present problem.

The laminar flame speeds are also quite high for the fuel-lean values of the equivalence ratio, due to the high mixture temperatures (Figure 14b). For some cases, values for very fuel-lean mixtures are not plotted due to the occurrence of autoignition events, hindering the determination of flame speed in classical sense. Note that the calculated laminar flame speed values (Figure 14b) cope well with the flow velocity, which may roughly be characterized by the co-flow inlet velocity of 4 m/s.

The calculated ignition delay times and laminar flame speeds (Figure 14) agree qualitatively well with the predicted performance of the detailed reaction mechanisms shown in the previous figures (Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13). The high values of the lift-off height correlate with the high values of the ignition delay time and low values of the laminar flame speed, as one would expect.

An inspection of the ignition delay times for all considered detailed reaction mechanisms show that the shortest ignition delay times occur for equivalence ratios (ϕ) within the range 0.02–0.05. Variation of the radial coordinate as a function of the axial coordinate along these isolines of the equivalence ratio (ϕ = 0.02, ϕ = 0.05) are plotted in Figure 15. These isolines are obtained from a non-combusting calculation (pure mixing) for T_CO = 1010 K. Additionally, the coordinates of the flame roots (the flame anchoring positions used in determining the lift-off heights, H, as illustrated in Figure 7) are plotted in the figure, as predicted by different detailed reaction mechanisms using the PDF (Figure 15a) and EDC (Figure 15b) turbulent combustion models.

It is interesting to observe that the predicted flame root positions lie very close to the region of shortest autoignition times (Figure 14a) for all cases (Figure 15). This implies that autoignition is very likely to play the major role in anchoring the flame. On the other hand, the very rapid increase of the autoignition time with the variation of the equivalence ratio indicates that the flame propagation phenomenon is likely to be dominant in the other regions of the flame front, determining its shape.

It is also interesting to note that the predictions obtained by PDF lie quite close to the ϕ = 0.02 curve (Figure 15a), whereas the EDC results turn out to be closer to the ϕ = 0.05 curve (Figure 15b).

For a fuller assessment of the overall performance of the different modeling approaches, the Central Processing Unit (CPU) times were compared. The mechanisms of Marinov et al. [58] and Li et al. [33] were taken as the representative global and detailed reaction mechanisms, respectively. The CPU times required for the same number of iterations are compared in Figure 16. The CPU times are scaled in such a way that the largest CPU time (which is for the detailed reaction mechanism in combination with the PDF) takes the value of 1000.

Obviously, the PDF turbulent combustion model leads to very much larger CPU times in comparison to the alternative modeling approaches, even in combination with a global reaction mechanism, which is the main obstacle for the use of this potentially most accurate model in complex, three-dimensional problems. The lowest CPU time is demanded by the FGM turbulent combustion model, although it integrates a detailed reaction mechanism. Unfortunately, the model turns out to be quite inaccurate for the present problem.

5. Conclusions

Lifted H₂/N₂ flame in vitiated co-flow was computationally investigated within a RANS approach for turbulence modeling. Global and detailed reaction mechanisms were applied in combination with different models for turbulent combustion. The focus of the investigation was the prediction of the lift-off height that varies with varying co-flow temperature. Predictions were compared with published measurements of other authors.

The results indicate that an accurate numerical prediction of the lifted hydrogen flame is, generally, a very challenging task, especially for low co-flow temperatures, i.e., for large lift-off heights. For rather low experimental lift-off heights (H/d < 10), some modeling approaches deliver fairly accurate results. For higher values of the experimental lift-off height, all modeling approaches show a rather large deviation from the experimental value.

Among the considered global reaction mechanisms, the mechanism of Marinov et al. [58] is observed to be able to predict the lift-off phenomenon when applied in combination with the eddy dissipation model (EDM+K) [56,62], eddy dissipation concept (EDC) [56,63], and composition PDF transport (PDF) [25,56] as the turbulent combustion model. However, the sensitivity of the lift-off height to the co-flow temperature is under-predicted. For low H/d (H/d < 10), EDM+K, and EDC, the results show a fair agreement with the experiments.

It is observed that Gri-Mech 3.0 [59] and Gri-Mech 2.11 [27] perform very differently. Gri-Mech 3.0 predicts complete blow-off or extremely large lift-off (at high co-flow temperatures in combination with PDF), whereas Gri-Mech 2.11 generally under-predicts the lift-off height, while delivering a quite good agreement with the measurements for H/d < 10, combined with the PDF turbulent combustion model.

The flamelet generated manifold (FGM) [56,63] turbulent combustion modeling approach, when used with a detailed reaction mechanism, predicts a rather small lift-off that may show a fair agreement with the experimental value only for high co-flow temperatures for some reaction mechanisms. However, the sensitivity of the lift-off height to the co-flow temperature is not predicted, leading to a strong under-prediction for all higher values of the co-flow temperature.

The results obtained by the detailed reaction mechanisms of Conaire et al. [60], Keromnes et al. [61], and Li et al. [33] show a similar qualitative behavior. The results obtained by the EDC turbulent combustion model exhibit two different characteristics for high and low co-flow temperature ranges. The predicted lift-off height varies rather weakly in the high co-flow temperature range. In the low co-flow temperature range, a quite steep increase of the lift-off height is observed, leading to a rapid blow-off. The predictions delivered by PDF as the turbulent combustion model agree qualitatively well with the measurements. As the results obtained by the Conaire et al. mechanism strongly over-predict the experimental values, the predictions of the Li et al. mechanism show the closest agreement with the experiments. Although the Li et al. mechanism combined with the PDF turbulent combustion model delivers the best agreement with the measurements, its accuracy still leaves much to be desired for low co-flow temperatures, where the increasing lift-off heights are remarkably over-predicted.