Comparison of Combustion Models for Lifted Hydrogen Flames within RANS Framework

Within the framework of a Reynolds averaged numerical simulation (RANS) methodology for modeling turbulence, a comparative numerical study of turbulent lifted H2/N2 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 H2/N2 jet. The predictions for the lifted H2/N2 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.


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 considered experimental conditions are listed in Table 1.

Lifted H2/N2 Jet Flame
As the main test case, the atmospheric, lifted flame of a H2/N2 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.  [19,23]. The considered experimental conditions are listed in Table 1.

Lifted H 2 /N 2 Jet Flame
As the main test case, the atmospheric, lifted flame of a H 2 /N 2 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. Stepowski [43].
The considered experimental conditions are listed in Table 1.

Lifted H2/N2 Jet Flame
As the main test case, the atmospheric, lifted flame of a H2/N2 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.  [19,23]. Figure 2. Sketch of experimental set-up of Wu et al. [19,23]. 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 2 /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 2 /Air ratios are represented in Table 2. 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 2 /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 2 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. 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.

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.

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 Energies 2020, 13, 152 7 of 24 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.

Reaction Mechanisms
Two single-step global reaction mechanisms were considered (that include only the main species, H 2 , O 2 , H 2 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.  [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: where M i and ν 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]: where C i stay for the species concentrations. The exponents n i are equivalent to ν 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.
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.

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: where ρ, k, ε 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 α, β 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: 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: In Equation (5), the variables ξ * and τ * 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 τ * 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: The functional relationships φ (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 ( φ) are obtained by integrating the relationships over the Z and C spaces, with the help of presumed probability density functions (Beta functions are used): For the construction of the beta functions, modeled transport equations are solved for mean values ( Z, C) and variances ( Z 2 , C 2 ) of Z, C [56,66,67]. The source term of the modeled 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: 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 Ψ and J i,k denote the composition space and molecular diffusion flux vectors, respectively. The terms A|B 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.

On Solution Domains, Boundary Conditions, Grids
The two-dimensional axisymmetric solution domains of the non-reacting H 2 /N 2 jet (Section 2.1) and the lifted H 2 /N 2 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. 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.

On Solution Domains, Boundary Conditions, Grids
The two-dimensional axisymmetric solution domains of the non-reacting H2/N2 jet (Section 2.1) and the lifted H2/N2 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 (Tables 1-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.

Non-Reacting H2/N2 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. 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 (Tables 1-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.

Non-Reacting H 2 /N 2 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 Energies 2020, 13, 152 11 of 24 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 H2/N2 flame.

Lifted H2/N2 Jet Flame
A grid independence study was conducted also for the lifted H2/N2 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. 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 2 /N 2 flame. 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 H2/N2 flame.

Lifted H2/N2 Jet Flame
A grid independence study was conducted also for the lifted H2/N2 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.

Lifted H 2 /N 2 Jet Flame
A grid independence study was conducted also for the lifted H 2 /N 2 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, 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 TCO = 1050 K (left) and TCO = 1010 K (right). In the figure, a sub-region in the nearfield of the central jet is displayed, with a size of 30dx10d. 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 TCO = 1050 K and TCO = 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. Not all turbulent combustion models were used in combination with all reaction mechanisms. An overview of the investigated combinations is provided in Table 4, where a cross (X) means that the combination was investigated, and a minus (−) means that the combination was not considered. Table 4. Investigated combinations of turbulent combustion models and reaction mechanisms. 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.     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.
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 TCO 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 TCO, 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. 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 TCO (TCO > 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 TCO (even showing a slight reverse trend for some values of TCO). 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 TCO. 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 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. 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. 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 H2-related reactions in between the versions seem to imply a strong difference for H2 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. 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 2 -related reactions in between the versions seem to imply a strong difference for H 2 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. 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 TCO (Figure 10). CO-EDC exhibits two regions with different qualitative behaviors ( Figure 11). For high TCO (TCO > 1040 K), i.e., low values of H, CO-EDC shows a slowly increasing H with decreasing TCO, and a good quantitative agreement with the measurements. However, for lower TCO values (TCO < 1400), a very rapid, abrupt increase of H is observed with decreasing TCO (with a remarkable over-prediction of the measured values), which quickly leads to blow-off beyond TCO < 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 TCO > 1030 K, the predicted values by KE-PDF agree rather well with the experimental ones. For TCO < 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 TCO, followed by a rapid increase and blow-off by a further decrease of TCO. Compared to CO-EDC (Figure 11), the characteristic TCO values are lower for KE-EDC ( Figure  12). The border between the two regimes is about TCO ≈ 1030 K and the blow-off occurs for TCO < 1020 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 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 TCO 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 SL (undisturbed, planar flame) were calculated using the abovementioned 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 TCO = 1010 K and TCO = 1050 K, spanning the whole considered range. Note that for a given coflow 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). 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. 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 (Figures 9-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.   (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 (Figures 9-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. 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.

Conclusions
Lifted H2/N2 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. 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.

Conclusions
Lifted H 2 /N 2 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 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.