Assessment of the reactingFoam Solver in Simulating the Hydrogen-Methane Bluff-Body-Stabilized Turbulent Diffusion Flame of the University of Sydney

1. Introduction

1.1. Flame Types

A flame is a self-sustaining zone of a combustion and heat release process. Such a flame is produced when a burnable substance (fuel) undergoes exothermic oxidation reactions with an oxidizer [1]. Flames represent a localized phenomenon. Therefore, they occupy a small portion of the combustible substance (fuel) at a given instance instead of taking place uniformly throughout that entire substance [2]. Typically, flames form as a result of reactions between gaseous primary reactants (the fuel and oxidizer) [3,4]. If the fuel is liquid, it undergoes local vaporization and forms a flammable gas before combustion occurs [5]. If the fuel is solid, it undergoes local pyrolysis and breaks down into a flammable gas before combustion occurs [6].

Based on the mixing process between the oxidizer and the fuel, it is possible to classify flames into two main categories [7].

The first flame type according to this classification is the premixed flame. In this type, the fuel and the oxidizer mix and form a homogeneous gas mixture before reaching the combustion zone [8]. Such premixed flames can be visually identified due to being non-luminous and translucent, having pale colors typically varying from yellow to green [9]. The reason that premixed flames tend to be non-luminous is the absence of the incandescent soot particulates [10,11]. This, in turn, is attributed to the consumption of generated soot particles in a premixed flame as they meet the oxidizer available within the flame and react with it [12]. Therefore, premixed flames tend to have a small amount of soot [13]. A famous example of premixed flames is the core of the blue conical flame of the Bunsen burner that is found in laboratory settings [14]. The flame found in most spark-ignited automotive engines is another example of premixed flames [15]. Also, the inner cone of an oxy-acetylene gas welding flame is a good example of premixed flames, where the mixing happens inside the welding torch [16].

The second type of flame, based on the fuel-oxidizer mixing, is the diffusion flame (also called a non-premixed flame). Such a flame is called “non-premixed” because the oxidizer and the fuel are introduced into the combustion zone as separate streams without prior mixing [17]. On the other hand, the name “diffusion” arises from the fact that the chemical reactions take place only where the fuel and the oxidizer are mixed [18] at the molecular level via differential molecular diffusion (DMD) [19]. In non-premixed flames, the time scale of chemical reaction is much smaller than the time scale for molecular diffusion. In other words, reaction is much faster than diffusion. As a result, molecular diffusion becomes the rate-limiting phenomenon [20] because the fuel and oxidizer diffuse into each other relatively at a much slower rate than their chemical reaction. Hence, non-premixed flames are named after the rate-controlling step of their formation, which is diffusion and mixing on a molecular scale [21]. Thus, non-premixed flames usually propagate and burn more slowly than premixed flames [22]. Diffusion (non-premixed) flames are commonly characterized by the release of soot particles [23]. This feature is caused by the insufficient oxidizer within the flame zone, resulting in the famous orange luminous color of such flames. Non-premixed flame burners often have larger flames than premixed flames, and typically exhibit more uniform temperature distribution (lower peak temperature) and heat flow [24]. These burners, known for their high heat production and efficiency, are designed to ensure the complete combustion of the fuel oil, thus reducing emissions [25]. Compared to premixed flames, non-premixed flames are encountered much more frequently in daily life [26]. Examples of non-premixed flames include matchstick flames, candle flames, open fires, the outer cone of the flame of the Bunsen burner [27], and the outer envelope of an oxy-acetylene gas welding flame [28].

Flames can also be categorized based on the regime of their movement, which can be laminar or turbulent [29]. The core of the flame in the case of a Bunsen burner represents a laminar premixed flame [30]. The flame of a candle and the outer cone of the flame in the case of a Bunsen burner are representative of laminar non-premixed flames [31]. Turbulent flames are often found in practical and larger applications [32], such as power generation [33] and thermal industrial processes [34]. Commercial-scale combustion [35] demands high speeds of supply for the reactants [36]. Such speeds exceed the thresholds of a laminar flame remarkably. Laminar flames are of special importance in academic research because they are simpler and easier to characterize and understand [37]. Laminar flames allow for building mathematical models for the flame combustion process, which can be extended afterwards to the more complex turbulent flames [38].

1.2. Turbulent Non-Premixed Flames

Turbulent non-premixed flames (TNF) combine the complexity of turbulence and the challenge of properly mixing the individual reactant streams [39]. Despite this, turbulent non-premixed flames are favored over premixed flames in a large number of industrial systems [40] due to two reasons:

First, non-premixed burners are simpler to design and build than premixed burners [41]. This is due to the higher tolerance of the reactant proportions for achieving a desired level of mixing [42]. There are no premixing requirements for non-premixed burners [43].

Second, non-premixed burners are safer to operate than premixed burners. This is because non-premixed burners do not exhibit flashback or autoignition in undesirable locations [44].

With turbulence introduced in turbulent non-premixed flames, the term “diffusion” that strictly applies to the differential molecular diffusion of chemical species remains the controlling mechanism. Initially, macroscopic turbulent convection is responsible for mixing the fuel and the oxidizer. Eventually, molecular diffusion completes the mixing process at the microscopic scale, such that chemical reaction can happen [45].

1.3. Instability of Turbulent Non-Premixed Flames

High jet speeds in turbulent non-premixed flames seem commercially attractive because they permit a high heat duty [46]. Such high-velocity burners are found in the aluminum, steel, or glass industries. However, such an approach to intensify the flame tends to make the flame unstable, where undesirable phenomena such as blow-off and lift-off appear. Lift-off refers to the phenomenon where the flame is detached from the burner, and an intervening zone of turbulent pre-mixing is established between the flame and the burner [47], as a result of excessive speed of the jet [48]. At larger jet speeds, the flame is extinguished, which is the blow-off phenomenon [49]. Either phenomenon (lift-off or blow-off) is undesirable and is considered a deviation from the stable attached flame [50].

High velocity burners in furnaces have another advantage beyond the intensified heating rates, which is environmental [51]. These burners produce low NOx (nitrogen oxides) emissions, mostly nitric oxide (NO) [52]. This happens because these burners inherently entrain significant amounts of cool furnace gases that are injected into the flame envelope before the complete combustion [53]. This results in reducing the average in-flame temperature, which hinders the high-temperature thermal formation of NOx [54]. This thermal formation of nitrogen oxides requires high temperatures (above 1,200 °C) to allow the molecular nitrogen gas (N₂) and the molecular oxygen gas (O₂) in the air to dissociate [55] into atomic nitrogen gas (N) and atomic oxygen gas (O) [56]. The atomic oxygen reacts with the molecular nitrogen, and the atomic nitrogen reacts with the molecular oxygen and forms nitric oxide (NO) [57].

Alternatively, reducing the environmentally-harmful NOx emissions [58] can be achieved by controlling the combustion conditions [59]. For example, enforcing lean fuel concentration also lowers the temperature and reduces the emissions [60]. However, this approach is considered harsh, and it increases the chance of flame instability [61].

1.4. Stabilizing Turbulent Non-Premixed Flames

Stabilizing high-jet turbulent non-premixed flames (this means anchoring the flames to the burner) in industrial burners can be achieved through three methods [62].

One of these stabilization methods uses a pilot burner (also called a pilot light) [63] as an active stabilization technique [64]. This technique involves adding a small auxiliary flame that burns initially in order to ignite the main burner [65]. The pilot flame can later be deactivated if the main flame does not show signs of instability [66]. However, the pilot flame is kept in continuous operation if this is necessary to avoid extinction of the main flame, thus to maintain a stable burner operation [67]. The auxiliary pilot flame can be implemented in the form of an electric igniter or in the form of a small premixed flame. Pilot-stabilized turbulent non-premixed flames (TNF) are known to be very stable over a wide range of operational conditions [68]. Nevertheless, pilot-stabilized turbulent non-premixed flames imply the penalty of installing and maintaining the additional pilot burner [69]. Also, there is extra energy consumption (either in the form of electricity or fuel) [70] to enable the pilot-flame stabilization method.

Another stabilization method for high-jet TNF utilizes swirl or recirculation, which is a passive technique (this means it does not consume external energy) [71]. In this method, a recirculation zone is carefully introduced such that it provides favorable conditions for the chemical reactions to take place. However, this swirling technique can be accompanied by unfavorable effects because a precessing vortex core can lead to combustion instability. Furthermore, the introduced swirl affects the unsteady location where the flame is anchored. Therefore, further combustion instabilities may occur, and this was observed experimentally for a gas turbine combustor [72].

The third TNF stabilization method is the method we explore in the current study. It is also a passive one. It. This method is the utilization of a bluff-body burner with a central fuel jet. This burner type is also called a flame holder [73]. This burner is characterized by a flat or a block-shaped object, which is the bluff body [74]. This bluff body helps in establishing a zone of low-speed hot recirculation zone [75]. When this bluff body is placed in the fast jet streams of the reactants, the flow is slowed down and heated, and this enables the chemical reactions to proceed toward ignition. Bluff-body non-premixed burners are particularly important in several industrial applications. Bluff body burners are also used in premixed flames [76,77]. Therefore, it is very useful to have simulation and design tools that are reliable for analyzing and predicting the flow fields and compositional profiles of the produced flames [78]. Bluff-body turbulent non-premixed flames burners serve as a good benchmarking case [79] for computational combustion modeling because they form a reasonable trade-off between the complication level associated with real, more sophisticated combustors and the relative simplicity and clear boundary conditions for modeling purposes.

1.5. HM1 Flame

The bluff-body stabilized flame adopted here as a benchmarking case for modeling is called the HM1 flame. It is an axisymmetric turbulent non-premixed flame that was studied experimentally at the University of Sydney [80]. It has been referenced online by the TNF Workshop Series [81,82,83,84,85] of the International Workshop on Measurement and Computation of Turbulent Non-premixed Flames [86].

This flame is assigned the name “HM1” because of the chemical composition of the gaseous fuel stream, which is composed of molecular hydrogen (H₂) and methane (CH₄) with equal volume fractions of 50% each (1:1 molar ratio). The burner is surrounded by a coflow of air. The HM flame family has three members: (1) HM1 flame, which has a fuel jet speed of 118 m/s, (2) HM2 flame, which has a faster fuel jet of 178 m/s, and (3) HM4 flame, which has a faster fuel jet of 214 m/s [87].

Given the 50%/50% volume composition (mole fractions) of CH₄ and H₂ in the fuel, the corresponding mass fractions are 8/9 or 0.8889 for methane (CH₄) and 1/9 or 0.1111 for molecular hydrogen (H₂) [88]. These mass fractions are based on molecular weights of 16 kg/kmol for methane [89] and 2 kg/kmol for molecular hydrogen [90]. The mixture’s molecular weight can be computed based on these data as the arithmetic mean of the two molecular weights, which is 9 kg/kmol.

The HM1 flame’s experiment is also assigned the code name “B4F3” by the experimental investigators at the University of Sydney, who are members of the Clean Combustion Research Group at the School of Aerospace, Mechanical and Mechatronic Engineering (AMME) [91].

The HM1 flame is an attractive benchmarking problem for computational fluid dynamics (CFD) solvers. This is due to the practicality of its high-velocity bluff-body burner as well as the detailed information about the burner, the flow settings, and the measurement process. This information allows for establishing a reasonably equivalent computational model with less deviation from the actual geometric or input conditions. In addition, the detailed experimental information for the HM1 flame is useful for making a non-ambiguous interpretation of the experimental data. These features of the HM1 flame consequently lead to meaningful and consistent assessment of the similarity between the CFD results and the experimental results.

1.6. Some Computational Studies About the HM1 Flame

The HM1 bluff-body stabilized turbulent non-premixed benchmarking flame considered in the current study has been utilized in earlier computational works.

For example, Yan et al. [92] solved the axisymmetric governing equations for this flame using a general-purpose computational fluid dynamics (CFD) computer code that employed a multi-block structured grid that was constructed based on general non-orthogonal coordinates and the finite volume methods (FVM) [93]. The non-uniform grid consisted of three blocks with a total of approximately 18,000 cells, which were distributed more densely near the bluff-body burner. They used the flamelet model for representing the turbulence-chemistry interaction (TCI) [94]. Their primary goal was to examine the influence of the turbulence model on the overall numerical predictions. They explored three different turbulence models, which are: (1) the basic k-epsilon ($k-ϵ$) model, (2) a modified k-epsilon model ($k-ϵ$) with a varied anisotropy parameter, and (3) the explicit algebraic stress model (EASM) [95,96]. It is worth mentioning that the coefficient $C_{ϵ2}$ in their basic k-epsilon model was 1.83, which is slightly below the standard value of 1.92 [97]. The researchers found that the numerical predictions depend on the selected turbulence model. Their numerical predictions were generally satisfactory, despite notable deviations in the radial profile of the root mean square (RMS) of the axial velocity fluctuations, the root mean square (RMS) of the mixture fraction variance, and the average mass fraction of carbon dioxide (CO₂) for each of the three turbulence models attempted.

Another research team, Liu et al. [98], treated the problem as axisymmetric also. However, they employed a single-block structured grid with a total of 6,912 cells (72 cells along the axial direction and 96 cells along the radial direction). They modeled the turbulence-chemistry interaction (TCI) through the joint probability density function (PDF) technique [99]. This research team investigated the sensitivity of their numerical predictions to the mixing model for the probability density function calculations. This mixing model represents the influences of molecular diffusion. The research team tried two different mixing models, which are (1) the interaction by exchange with the mean (IEM) model [100,101], and (2) the Euclidean minimum spanning tree (EMST) [102]. Some numerical results of the radial profiles of the mean mixture fraction and its root mean square (RMS) values showed strong deviations near the burner’s centerline. It was observed that the model gave much lower values than the experimental values.

In another study, Odedra and Malalasekera [103] again treated the HM1 flame as an axisymmetric problem in their numerical modeling. They performed Favre-averaged Navier Stokes (FANS) simulations using steady and unsteady flamelet models. The time-dependence was captured through a post-processing technique using the Eulerian particle flamelet model (EPFM) [104]. They implemented their simulations using the commercial computational fluid dynamics software Ansys Fluent [105]. In their simulations, they selected the Reynolds stress model (RSM) [106] for representing the turbulence effects. Excluding a small upstream extension, the main computational domain was discretized into 44,200 (260 axially × 170 radially) quadrilateral cells. The computed radial profiles of the axial velocity were found to be in good agreement with the experimental data. However, the radial velocity profiles showed big deviations from the experimental data. The computed mass fraction of molecular hydrogen (H₂) as a reactant component was over-predicted near the centerline. Nevertheless, the computed results quickly became aligned with the experimental data as the distance from the centerlines increased.

Wang et al. [107] used a modified version of the HM1 flame as a representative test case to examine the finite-rate chemical reaction in the case of a specialized MILD (moderate or intense low oxygen dilution) combustion [108], with a Jet-into-Hot-Coflow (JHC) burner. They investigated the differences between the environmentally-attractive [109] reduced-temperature low-pollution (low-NOx) MILD combustion and the conventional hot-temperature combustion. These differences include the dissimilarities in the flame characteristics. They used joint velocity-frequency-composition transported probability density function (TPDF) simulations, combined with a sensitivity analysis. It should be noted that their modeling does not correspond to the original HM1 flame we refer to here, but to a variant of it, since the burner they modeled has a MILD JHC [110] design that is different from the original experimental configuration, and also the coflow was oxygen-enriched. The coflow is a hot exhaust gas with oxygen and nitrogen. The chemical composition of that oxidant coflow, in the form of mass fractions, was 0.065 H₂O, 0.055 CO₂, 0.030 O₂, and 0.850 N₂ [111]. The research team attempted to maintain some features of the HM1 flame, such as the 1:1 volumetric ratio of the molecular hydrogen (H₂) and methane (CH₄).

In the studies of Raman and Pitsch [112], Kempf et al. [113], James et al. [114], and Chen et al. [115], the HM1 benchmarking flame was modeled using the computationally demanding [116] large eddy simulation (LES) method with a three-dimensional geometry. Given the large computational resources and time demands involved in this modeling method, these studies are considered less practical and more academic than strategies aimed at general engineering and industry problems [117].

In our computational fluid dynamics (CFD) simulations, we report here for the HM1 flame, we use a simplified axisymmetric geometry, which is adequate for this axisymmetric flame. We solve the FANS (Favre-averaged Navier Stokes) governing equations using the reactingFoam solver, which is part of the open-source software OpenFOAM (Open Field Operation And Manipulation) [118]. Our assessment of the open-source reactingFoam solver provides the interested researchers with a source of information to assess the solver and make a guided decision about the solver’s capabilities and limitations when modeling turbulent non-premixed flames (TNF) or related problems involving reactive gases and combustion. This assessment is achieved both qualitatively by visually comparing axial and radial profiles of velocity, chemical, or thermal field variables along different positions within the flame region, and quantitatively through computed aggregate deviation scalar quantities, such as the mean absolute deviation (MAD).

2. Experimental Geometry of the HM1 Flame

Figure 1 shows a schematic of the experimental cylindrical HM1 burner. The burner is situated inside a wind tunnel whose cross-section is a square. Despite the enclosing wind tunnel walls, the separation between the burner’s centerline and the wind tunnel walls is large enough to make the HM1 flame effectively open (not confined) [119]. The burner has a flat annular bluff-body face, with a fuel jet nozzle located in the center.

The fuel nozzle (orifice) has a diameter of 3.6 mm, while the outer burner diameter is 50 mm. The diagonal distance between the two opposite corners of the wind tunnel is 431.3 mm, while the perpendicular face-to-face distance is 305 mm. The cross-sectional area of the wind tunnel is 93,025 mm², while the cross-sectional area of the burner’s circular cross-section (including the fuel nozzle) is 1.963.5 mm², and the fuel nozzle’s area is 10.179 mm². Thus, the area ratio of the wind tunnel to the burner to fuel nozzle is 47.3772:1.0000:0.0052 or 9139.13:192.90:1.00. These ratios indicate the big differences between the three areas.

The relatively large annular area of the bluff body (193 times the area of the fuel jet) results in forming a recirculation zone in which high-temperature gases are located, which are sufficiently hot to stabilize the flame by keeping it attached to the burner.

Table 1 lists several properties of the HM1 flame (experiment code B4F3, dated 1995) [120]:

In a turbulent non-premixed (diffusion) flame (TNF), if the fuel jet speed increases and reaches a threshold limit (the blow-off speed), then the flame can be extinguished. This flame extinction happens downstream of the recirculation zone. The extinguished flame may reignite at a further downstream point. From the data given in Table 1, it can be seen that the HM1 flame is approximately at 50% of the blow-off threshold. Therefore, this blow-off instability is not encountered.

Figure 2 is a photo of an HM flame [121]. It is noticeable in the photo that this turbulent non-premixed flame (TNF) has a blue color rather than the more expected luminous orange color for non-premixed flames due to soot incandesce. The explanation of the uncommon blue color is also the explanation for blending methane (CH₄) with molecular hydrogen (H₂) in the HM family of flames. The presence of hydrogen in this fuel blend reduces the soot formation compared to burning pure methane or natural gas [122,123]. The suppression of soot formation in the synthetic CH₄/H₂ fuel is aligned with the formation of a recirculation zone that is nearly clean from soot [124]. The addition of hydrogen results in a carbon dilution effect (changing the concentration of carbon per unit mass of the fuel), which in turn causes a change in the distribution of the temperature field within the flame [125]. The addition of hydrogen also inhibits soot precursor species, such as acetylene or ethyne (C₂H₂) and polycyclic aromatic hydrocarbons (PAHs) [126]. The flames of pure methane or natural gas (without hydrogen blending) generally exhibit the formation of soot in the recirculation zone. As this formed soot is carried downstream (upward) [127], it undesirably interferes with the Raman signal. This interference harmfully affects the accuracy of the measured species mass fractions. Adding hydrogen to methane satisfactorily frees the recirculation zone from the undesirable soot.

In the experimental setting, the measurements of the axial velocity component and the radial velocity component were made concurrently using the non-intrusive optical technique of laser Doppler velocimetry (LDV) [128], with a two-color argon-ion laser beam having a power of 5 W [129]. This laser type has two major emission lines, one at 488.0 nm (blue) and another at 514.5 nm (green) [130,131].

Raman spectroscopy [132], Rayleigh spectroscopy [133], and laser-induced fluorescence (LIF) [134] were used to measure instantaneous values of the temperature and the species mass fractions at different spatial locations.

3. Computational reactingFoam Model

3.1. Overview of OpenFOAM and reactingFoam

Our computational modeling for the HM1 flame was conducted using the computational fluid dynamics (CFD) solver reactingFoam [135]. This solver supports multi-species finite-rate chemistry, thermodynamics, flow turbulence, and turbulent combustion models. The solver is widely used in the combustion community [136]. The solver reactingFoam is one of several CFD solvers distributed under the OpenFOAM open-source software, which is a collection of computer codes written using the object-oriented programming (OOP) concepts [137] and the principles of the C++ programming language. The computer codes are augmented with various text files that serve as conveniently-editable datasets and lookup input files [138,139] for describing the computational geometry, defining the boundary conditions, setting the initial conditions, determining the physical modeling parameters, selecting the desired numerical solver settings, computing needed fluid properties, and listing the chemical reaction steps.

3.2. Main Governing Equations

The scalar mass conservation equation, the vector momentum conservation equation, the scalar species mass conservation equation (expressed in terms of the mass fraction), and the scalar energy conservation equation (expressed in terms of the static specific enthalpy) in the reactingFoam solver, respectively, are [140,141]

$${\displaystyle \frac{\partial \rho }{\partial t}}\mathrm{ }+\nabla \cdot \left(\rho \mathrm{ }\overrightarrow{U}\right)\mathrm{ }\mathrm{ }$$

(1)

$${\displaystyle \frac{\partial \left(\rho \overrightarrow{U} \right)}{\partial t}}+\nabla \cdot \left(\rho \overrightarrow{\overrightarrow{UU}}\right)-\nabla \cdot \left({\mu }_{eff} \overrightarrow{\overrightarrow{\nabla U}}\right)=-\overrightarrow{\nabla p}+\nabla \cdot \left({\mu }_{eff}\left[{\overrightarrow{\overrightarrow{\nabla U}}}^T-\overrightarrow{\overrightarrow{I}}{\displaystyle \frac{2}{2}}tr\left(\overrightarrow{\overrightarrow{\nabla U}}\right)\right]-\overrightarrow{\overrightarrow{I}}{\displaystyle \frac{2}{3}}\rho k\right)+\rho \stackrel{⃑}{g}$$

(2)

$${\displaystyle \frac{\partial \left(\rho Y_i\right)}{\partial t}}+\nabla \cdot \left(\rho \overrightarrow{U}{Y}_i\right)-\nabla \cdot \left(\rho D_{i,eff} \overrightarrow{\nabla Y_i}\right)=\kappa {\mathbb{R}}_i$$

(3)

$${\displaystyle \frac{\partial \left(\rho h\right)}{\partial t}}+\nabla \cdot \left(\rho \overrightarrow{U}h\right)-\nabla \cdot \left({\alpha }_{eff} \overrightarrow{\nabla h}\right)={\displaystyle \frac{Dp}{Dt}}+\sum _{i=1}^{N_s}\left(\nabla \cdot 〈h_i\left[\rho D_{i,eff}-{\alpha }_{eff}\right]〉 \overrightarrow{\nabla Y_i}\right)$$

(4)

These equations apply to a compressible turbulent gaseous mixture. They are the Favre-averaged Navier Stokes formulation (FANS), which is theoretically derived based on density-weighted averaging of field variables [142]. The FANS formulation here is extended to a multi-species reacting gas. They are also called compressible Reynolds-averaged Navier-Stokes (RANS) equations [143].

In the above partial differential equations, a single arrow above a variable’s symbol means that it is a vector field, while a double arrow means that it is a tensor field.

The density is denoted by the symbol $\rho$, the velocity vector is denoted by $\overrightarrow{U}$, the effective dynamic viscosity (sum of molecular viscosity and eddy viscosity) [144] is denoted by ${\mu }_{eff}$, the static pressure is denoted by $p$, the turbulent kinetic energy per unit mass is denoted by $k$, the effective kinematic mass diffusivity of the ith gaseous species in the mixture is denoted by $D_{i,eff}$, the mass fraction of the ith gaseous species in the mixture is denoted by $Y_i$, the reactive volume fraction [145] (which is a positive dimensionless scaling multiplier bounded by an upper limit of 1.0) is denoted by $\kappa$, the unscaled chemical reaction source [146] of the ith gaseous species is denoted by ${\mathbb{R}}_i$, the mixture’s static specific enthalpy is denoted by $h$, the static specific enthalpy of the ith gaseous species is denoted by $h_i$, the effective dynamic thermal diffusivity is denoted by ${\alpha }_{eff}$, and the total number of gaseous species in the mixture is denoted by $N_s$.

In the momentum equation, Equation (2), the operator ($\mathrm{t}\mathrm{r})$ is the trace of a tensor. The identity tensor in that equation is denoted by $\overrightarrow{\overrightarrow{I}}$ . The superscript ($T$) is the transpose operator.

In the static-enthalpic energy equation, Equation (4), the term ${\displaystyle \frac{D}{Dt}}p$ is the substantial time derivative of the static pressure [147]. The last term of that equation was added as a differential diffusion term. It represents the difference between the species’ effective mass diffusivity and the mixture’s effective thermal diffusivity. This term implies the assumption of a unity molecular Lewis number and unity turbulent Lewis number for each gaseous species in the mixture [148]. The Lewis number is the ratio of the thermal diffusivity to the mass diffusivity [149]. It is admitted that this differential diffusion effect is generally neglected when modeling turbulent reacting flows [150].

All flow variables are stored at the cell centers (the center of the control volumes), rather than the cell’s face centers. However, spatial interpolation is heavily implemented during the simulation in order to compute the needed cell’s face values of the variables as temporary quantities. With this spatial scheme, the computational mesh (grid) is described here as having a collocated, rather than staggered, layout [151]. Collocated-mesh (CM) layouts are often favored over staggered-mesh (SM) layouts in the case of turbulence simulations [152]. Compared to staggered grids, collocated grids have the advantages of being easier to extend for three-dimensional nonorthogonal grids and allowing for better treatment of the domain boundaries in the case of body-fitted grids [153].

3.3. Turbulence-Chemistry Interaction

The turbulence-chemistry interaction (TCI) in our simulations is based on the Chalmers’ partially stirred reactor (CPaSR) model [154]. We adopted a minor variant of the standard k-epsilon model for modeling the unresolved turbulent scale. The minor improvement change adopted here is increasing the value of the model parameter ($C_{ϵ1}$) from the standard level of 1.44 to a slightly higher level of 1.60. This 11.1% increase is introduced because it better suits problems with self-similar round jets as the one encountered here [155]. In addition, this adjustment was found to yield favorable numerical stability [156].

With the simplifying assumption of unity laminar Schmidt numbers and unity turbulent Schmidt numbers for each species, $\rho D_{i,eff}$ can be replaced by ${\mu }_{eff}$in Equation (3) and Equation (4). The modified versions of these two equations, respectively, are

$${\displaystyle \frac{\partial \left(\rho Y_i\right)}{\partial t}}+\nabla \cdot \left(\rho \overrightarrow{U}{Y}_i\right)-\nabla \cdot \left({\mu }_{eff} \overrightarrow{\nabla Y_i}\right)=\kappa {\mathbb{R}}_i$$

(5)

$${\displaystyle \frac{\partial \left(\rho h\right)}{\partial t}}+\nabla \cdot \left(\rho \overrightarrow{U}h\right)-\nabla \cdot \left({\alpha }_{eff} \overrightarrow{\nabla h}\right)={\displaystyle \frac{Dp}{Dt}}+\sum _{i=1}^{N_s}\left(\nabla \cdot 〈h_i\left[{\mu }_{eff}-{\alpha }_{eff}\right]〉 \overrightarrow{\nabla Y_i}\right)$$

(6)

The chemical kinetics are described by the Jones-Lindstedt (JL) reaction mechanism [157], which involves six species (CH₄, H₂, O₂, CO₂, H₂O, CO) [158] and four reaction steps [159]. Nitrogen in the air is treated as an inert gas [160,161].

3.4. Pressure-Velocity Coupling

Our simulations adopt the PISO (Pressure Implicit Splitting of Operators) coupling scheme for the pressure and velocity [162]. In the PISO pressure-velocity coupling scheme, the primitive form of the continuity equation, Equation (1), is not directly solved. Instead, mathematical manipulation is performed to give two derived equations. The first derived equation is a Poisson-type elliptic partial differential equation (PDE) for the pressure. The source term in this derived PDE depends on the velocity vector. The second derived equation is an explicit expression for the velocity vector, which relates the velocity vector to the gradient vector of the pressure.

To better explain the PISO pressure-velocity coupling scheme, we consider the case of an incompressible fluid [163]. Writing a semi-discretized version of the vector momentum equation, Equation (2), gives [164]

$$a_P\overrightarrow{U_P}=\overrightarrow{H}\left(\overrightarrow{U}\right)-\overrightarrow{\nabla p}$$

(7)

where $a_p$ is a fixed-value coefficient for each cell, and it does not change over time. This coefficient depends on the geometry and discretization scheme used for the momentum equation [165]. This coefficient is the multiplier quantity for the cell’s own velocity vector $\overrightarrow{U_P}$. The source vector function $\overrightarrow{H}\left(\overrightarrow{U}\right)$ is vector that accounts for the discretized terms of the momentum equation while excluding the pressure gradient. This means that this vector source function accounts only for the transport part of the momentum equation.

The derived implicit elliptic equation for the pressure and the companion derived explicit equation for the velocity vector correction can be obtained directly from Equation (7). Respectively, these two derived equations of the PISO scheme are

$$\nabla \cdot \left({\displaystyle \frac{1}{a_p}} \overrightarrow{\nabla p}\right)=\nabla \cdot \left({\displaystyle \frac{1}{a_p}}\overrightarrow{H}\left(\overrightarrow{U^{*}}\right)\right)$$

(8)

$$\overrightarrow{U_P^{corr}}={\displaystyle \frac{1}{a_p}}\overrightarrow{H}\left(\overrightarrow{U^{*}}\right)-{\displaystyle \frac{1}{a_p}}\overrightarrow{\nabla p}$$

(9)

It can be noticed that Equation (9) is very similar to Equation (7). Both sides of Equation (7) were divided by $a_p$ to reach Equation (9), which strengthens its nature as an explicit equation. Also, the right-hand side velocity vector in Equation (9) has a temporary (predicted) velocity vector $\overrightarrow{U^{*}}$ that needs a corrector step in order to become a corrected velocity vector $\overrightarrow{U_P^{corr}}$ in the left-hand side.

From Equation (9), it is possible to derive Equation (8) after applying the divergence operation to both sides of the equation, and then enforcing the condition that the corrected velocity field has zero divergence, which is mathematically expressed as

$$\nabla \cdot \overrightarrow{U^{corr}}=0$$

(10)

The vector momentum equation, Equation (2), is solved initially using the latest existing pressure field. This yields a temporary velocity vector field that does not necessarily satisfy the continuity equation.

This predictor stage is followed by a three-step PISO corrector stage, which is described as follows:

• First, the elliptic pressure equation, Equation (8), is solved to give a pressure field.
• Second, the mass fluxes at the cell faces are updated using the obtained pressure field.
• Third, the temporary velocity vector field obtained during the predictor stage is corrected using the explicit expression in Equation (9). This takes into account the newly resolved pressure field.

The above sequence represents a single PISO correction step [166]. It can be repeated for improved performance.

For the overall solution procedure, the governing equations are corrected five times per time step in a sequential manner [167,168]. Within each of these “outer” corrections, two PISO “inner” corrections are performed.

Therefore, the following governing equations are solved five times in each time step:

• Momentum equation
• Species mass-fraction equations
• Energy equation

On the other hand, the following governing equations are solved ten times in each time step:

• Elliptic pressure equation
• Explicit velocity vector correction equation

The following solution order takes place in the computational algorithm during a single “outer” loop:

4.

The momentum equation, Equation (2), is solved.

5.

The species equations, Equations (5), are solved.

6.

The energy equation, Equation (6), is solved.

7.

The PISO correction “inner” loop is performed.

For each implicit equation (all field equations except the velocity correction equation), an algebraic system of equations is obtained due to the finite volume discretization at the cell level [169]. These algebraic systems are solved using iterative solvers. The preconditioned bi-conjugate gradient (PBiCG) solver is used [170,171], with a diagonal-based incomplete lower-upper LU (DILU) preconditioner [172,173]. A tolerance of 10^–5 is set for the convergence criterion for this iterative algebraic-system solver [174,175].

3.5. Computational Domain and Mesh

The reactingFoam solver and other OpenFOAM solvers always adopt three-dimensional meshes even when modeling two-dimensional planar or two-dimensional axisymmetric problems [176]. An appropriate set of boundary conditions is then applied to the three-dimensional domain to account for the reduced dimensionality [177]. For planar problems, the equivalent three-dimensional domain formed by OpenFOAM is an extruded version of the planar area [178], with a one-cell layer thickness [179]. For axisymmetric problems, the equivalent three-dimensional domain formed by OpenFOAM is a wedge-like body [180].

In our current axisymmetric HM1 flame configuration, the three-dimensional computational domain we constructed through OpenFOAM is the wedge illustrated in Figure 3. It is full-wedge-angle is 5°. The domain extends axially for 400 mm downstream the burner face, with an intentionally added auxiliary part upstream the inlet of the fuel jet (7.2 mm extension) and upstream the inlet air coflow (25 mm extension) [181]. These additional zones allow the inlet flow to develop a more natural turbulent profile before they are introduced into the main flame zone [182].

To allow the incoming fuel stream and the air stream to build a realistic non-uniform profile at the inlet of the computational domain, the computational domain is extended in the upstream direction by four jet radii for the fuel jet and by one bluff-body radius for the air coflow. The figure shows three blocks that represent a structured non-uniform (stretched) multi-block mesh [183]. Cells near the centerline and near the burner are smaller to allow for better capturing of stepper changes.

As shown in the figure, the three blocks are identified by 12 grid points (numbered from 0 to 11). OpenFOAM has a built-in mesh generator utility that understands this numbering and properly constructs a wedge geometry and its mesh [184].

There is one limitation in the described axisymmetric geometric model, which is that it is unable to precisely replicate the square wind-tunnel boundary around the air coflow. To circumvent this limitation, the coflow passage is approximated here as an annulus (which is compatible with the axisymmetric geometric modeling) with a radius of 172 mm. This radius is computed based on an area equivalence condition, as shown in the following equation:

$$305\times 305\approx \pi {\left(172\right)}^2\approx \mathrm{93,000}$$

(11)

It is useful to add here that the impact of this geometric simplification is considered insignificant because the flow variable variations largely decay well before the outer boundary is reached. The HM1 flame is essentially an open flame, and the shape of the side geometric boundary practically acts as a far-stream condition.

We solved the problem with two mesh resolutions. We used a coarse mesh with 12,836 cells and a fine mesh with 37,202 cells. The ratio of the number of cells is 2.898 (nearly triple the cell count). The use of two levels of resolution helps in testing the convergence of the solution and the adequacy of the grid. In the coarse resolution, the axial cell size starts from 1.89 mm at the burner face and increases to 9.44 mm at the downstream exit boundary. In the fine resolution, the axial cell size starts from 0.4 mm (21% of the coarse-mesh value) at the burner face and increases to 8.20 mm (87% of the coarse-mesh value) at the downstream exit boundary. In the coarse resolution, the radial cell size at the centerline is 0.223 mm. In the fine resolution, the radial cell size at the centerline is 0.126 mm (57% of the coarse-mesh value).

At the downstream exit boundaries (the outflow face of the computational domain), a pressure boundary condition is used. At the flow inlet boundaries, a zero-gradient pressure condition is applied. The wall-function treatment is employed at the walls [185].

The chemical composition of the fuel stream and the air stream is specified in terms of the mass fractions (0.8889 for CH₄ and 0.1111 for H₂), consistent with experimental data of the HM1 flame.

The temperature and species mass fractions are fixed at the inlet boundaries. The fuel jet temperature and the coflow temperature are set to 298 K and 300 K, respectively.

A zero-gradient condition is specified at the outflow boundaries. At the front and back faces of the wedge, an axisymmetric boundary condition is applied. In OpenFOAM, this condition type is called “wedge”.

3.6. Mesh and Boundaries Script

The input data file for constructing the three-block computational domain and setting the boundary conditions is listed below.

FoamFile

{

version 2.0;

format ascii;

class volVectorField;

object U;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -1 0 0 0 0];

internalField uniform (1e-4 0 0);

boundaryField

{

primaryJet

{

type swirlMassFlowRateInletVelocity;

massFlowRate 0.0001375;

rpm 0;

value uniform (1 0 0);

}

secondaryJet

{

type fixedValue;

value nonuniform List<vector> 8

(

(18.27047713 0 27.748)

(18.26797534 0 29.893)

(18.26618485 0 32.173)

(18.26529552 0 34.596)

(18.26483522 0 37.117)

(18.2642634 0 39.570)

(18.26454143 0 41.848)

(18.26228408 0 43.963)

);

}

centerline

{

type fixedValue;

value uniform (0 0 0);

}

outlet

{

type zeroGradient;

}

horizAndVerticWalls

{

type fixedValue;

value uniform (0 0 0);

}

frontAndBack_pos

{

type wedge;

}

frontAndBack_neg

{

type wedge;

}

}

// *********************************************************************** //

3.7.reactingFoam Code

The main code file of the used reactingFoam solver is listed below.

The outer solver loop starts at the line “for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)”, while the inner PISO loop within that outer loop starts with the line “for (int corr=1; corr<=nCorr; corr++)”.

\*---------------------------------------------------------------------------*/

#include “fvCFD.H”

#include “hCombustionThermo.H”

#include “compressible/RASModel/RASModel.H”

#include “chemistryModel.H”

#include “chemistrySolver.H”

#include “multivariateScheme.H”

#include “ignition.H” // added

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])

{

# include “setRootCase.H”

# include “createTime.H”

# include “createMesh.H”

# include “readChemistryProperties.H”

# include “readCombustionProperties.H”

# include “readEnvironmentalProperties.H”

# include “createFields.H”

# include “createAdditionalFields.H”

# include “initContinuityErrs.H”

# include “readTimeControls.H”

# include “compressibleCourantNo.H”

# include “setInitialDeltaT.H”

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

while (runTime.run())

{

# include “readTimeControls.H”

# include “readPISOControls.H”

# include “compressibleCourantNo.H”

# include “setDeltaT.H”

runTime++;

Info<< “Time = “ << runTime.timeName() << nl << endl;

# include “chemistry.H”

# include “rhoEqn.H”

# include “UEqn.H”

# include “ignite.H”

for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)

{

# include “YEqn.H”

# define Db turbulence->alphaEff()

# include “hEqn.H”

for (int corr=1; corr<=nCorr; corr++)

{

# include “pEqn.H”

}

}

turbulence->correct();

rho = thermo->rho();

runTime.write();

}

return(0);

}

4. Results

4.1. Selected Data

We present our HM1 simulation results in the form of:

• radial profiles of the axial velocity component
• radial profiles of the radial velocity components
• radial profiles and axial profiles of the reactant methane (mass fraction)
• radial profiles and axial profiles of the product water vapor (mass fraction)
• axial profile of the temperature

Despite the large number of results available, we selected the above-listed set. It allows a reasonable level of diversity, and it is also effective in showing the spatial variations of the flow field within the flame. It also enables convenient judgment of the level of agreement between the simulation results and the measurements, as well as between the coarse-mesh solution and the fine-mesh solution.

4.2. Reynolds Averaging Versus Favre Averaging

Before comparing the simulation results from the reactingFoam CFD modeling of the HM1 TNF, under two mesh resolutions, with the experimental data, it is important to clarify here that the reported experimental data for the scalar variables (the temperature and the species mass fractions) are available in two forms. These are (1) ensemble average (Reynolds mean) and Favre average (density-weighted mean) [186].

Because our computational fluid dynamics (CFD) model handles a compressible fluid and inherently solves for variables in a Favre-averaged type, we select that type of experimental data for the comparison.

On the other hand, only the Reynolds-averaged type is reported in the experimental data for the velocity fields, but we compare it with the CFD Favre average. While this mismatch introduces errors, our examination of the difference in the experimental data between the Favre and Reynolds averages for some scalar quantities shows that the differences are generally small. This mitigates the impact of this inconsistency. For example, at an axial distance of 13 mm downstream the burner and a radial distance of 1.1 mm from the centerline, the Reynolds-average temperature is 369 K (with RMS 45 K) and the Favre-average value is 365 K (with RMS 43 K). The difference is very small. At the same location, the Reynolds-average mass fraction of methane is 0.81216 (with RMS 0.05715) and the Favre-average value is 0.81633 (with RMS 0.05476). Again, the difference is very small. As another example, at an axial distance of 65 mm downstream the burner and a radial distance of 1.1 mm from the centerline, the Reynolds-average temperature is 831 K (with RMS 93 K) and the Favre-average value is 826 K (with RMS 91 K). At that location, the Reynolds-average mass fraction of methane is 0.45185 (with RMS 0.05997) and the Favre-average value is 0.45444 (with RMS 0.05921). The differences are very small.

Also, it is worth mentioning that while our reactingFoam modeling is unsteady, we wait for a sufficient time until a steady state (temporal changes decay to trivial levels) is reached before sampling the desired fields and exporting the computational data as external data files for visualization and post-processing.

4.3. Radial Profiles of the Axial Velocity

The axial velocity component at increasing axial locations beyond the jet inlet is presented in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8. The corresponding axial distances are 5 mm (10% of the bluff-body burner diameter of 50 mm), 13 mm (26% of the burner diameter), 30 mm (60% of the burner diameter), 45 mm (90% of the burner diameter), and 65 mm (130% of the burner diameter), respectively.

For better visualization purposes and to emphasize the flame region near the burner, the span of the radial coordinate is limited to 30 mm (120% of the burner radius) in these figures. The experimental values are compared with the reactingFoam results for both the coarse mesh and the fine mesh.

Due to diffusion, the variation in the axial velocity is sharper near the burner, while it decays downstream. At the first axial station (5 mm after the bluff-body face), the coflow entrance velocity of 40 m/s can be identified at radial coordinates exceeding the burner radius of 25 mm.

We point out that although the fuel jet velocity in Table 1 is 118 m/s, a higher value of about 160 m/s is reached at the centerline in the case of the simulations. The reason for this apparent discrepancy is that the lower tabulated fuel jet velocity is a bulk (average) speed, while the actual centerline speed should exceed it, and the near-nozzle speeds should be below that average value.

Although the axial station of 5 mm only downstream of the burners is a challenging station for the reactingFoam CFD solver because the velocity variation is steep, the computational results are in good overall agreement with the experimental results. The numerical model was able to capture well the steep velocity drop at the edge of the jet core.

At the other four axial stations, there is a clear zone of negative axial velocity (reversed flow) between the jet core and the surrounding coflow. This negative axial velocity reflects the presence of a recirculation zone that is established behind the bluff-body face. The peak reversal velocity resolved in the simulation is –14.35191717 m/s (at a radial distance of 10.62124 mm) in the case of the coarse mesh (12,836 cells), and –14.27113179 m/s (at a radial distance of 10.92147 mm) in the case of the coarse mesh (37,202 cells). The two peak values are close to each other.

The computational fluid dynamics model performs very well in capturing the axial decay of the velocity, which is indicated by the agreement with the experimental results near the centerline, where the velocity decreases as the axial distance from the burner increases, particularly at the axial distance of 65 mm.

The predictions of the two mesh resolutions are very similar. This suggests that a mesh-independent solution is reached. This also shows that the coarse mesh (12,836 cells) is actually enough to model the HM1 flame.

4.4. Radial Profiles of the Radial Velocity

Next, we present the profiles of the radial velocity. These profiles are shown for five successive axial stations after the jet inlet in Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13. Respectively, they correspond to the same axial distances considered in the previous section for the sampled axial velocity component. Again, these five axial distances after the burner are 5 mm (10% of the bluff-body burner diameter of 50 mm), 13 mm (26% of the burner diameter), 30 mm (60% of the burner diameter), 45 mm (90% of the burner diameter), and 65 mm (130% of the burner diameter).

Similar to the axial velocity component discussed in the previous subsection, the variation in the radial velocity component is steeper near the bluff-body face.

In the immediate vicinity of the burner’s centerline, the radial velocity is positive (pointing outward, away from the centerline). However, radial velocity becomes negative as the radial distance increases (inward motion), and then its sign is reversed again (from negative to positive), and it is reversed once again (from positive to negative) when the coflow zone is reached. This is a characteristic profile attributed to the recirculation phenomenon caused by the bluff body and also attributed to the entrainment phenomenon [187] caused by the core jet.

At the smallest axial distance of 5 mm, the peak negative radial velocity is –2.383640222 m/s, which is located at a radial distance of 3.02089 mm in the case of the coarse grid. This becomes –2.628699787 m/s at a radial distance of 2.81983 mm in the case of the fine grid. On the other hand, the peak positive radial velocity is 3.384016841 m/s at a radial distance of 17.36856 mm in the case of the coarse grid. This becomes 3.179678864 m/s at a radial distance of 17.62448 mm in the case of the fine grid.

At the next axial distance of 13 mm, the peak negative radial velocity is –2.246834789 m/s at a radial distance of 4.07027 mm in the case of the coarse grid. This becomes –2.45293912 m/s at a radial distance of 3.83079 mm in the case of the fine grid. On the other hand, the peak positive radial velocity is 1.785740987 m/s at a radial distance of 17.36856 mm in the case of the coarse grid. This becomes 1.803396426 m/s at a radial distance of 17.29624 mm in the case of the fine grid. The present small changes in the results due to the mesh resolution change suggest that adopting the fine grid is justified despite the additional computational work and time.

At the smaller axial distance of 5 mm from the burner bluff-body face, the overall radial profile is predicted correctly by the computational model, with a good qualitative agreement with the experimental results. However, the coarse mesh is accompanied by a solution that exhibits small oscillations aligned with the radial edge of the bluff body. These oscillations occur at the interface between the coflow and the recirculation zone, and they rapidly decay and disappear at a radial distance of 29 mm (thus, they are present over a radial range of only 4 mm). These spurious oscillations [188] disappear in the fine mesh. Other than these numerical wiggles, the solution with the fine mesh is close to the solution with the coarse mesh. It can be said that the coarse mesh is satisfactory for modeling the HM1 benchmarking flame.

At the larger axial distance of 13 mm from the burner bluff-body face, it is noticed that the experimental profile of the radial velocity component exhibits discernible irregularities, although the profile was smoother and more aligned with the recirculation and entrainment phenomena in the upstream station at 5 mm. The computational results at the axial distance of 13 mm are smoother and more consistent with the recirculation and entrainment, as well as consistent with the upstream profile at 5 mm. Upon inspecting the experimental dataset, we found that the RMS of the radial velocity component is relatively large compared to the magnitudes of the mean radial velocities. For example, at the axial location of 13 mm, the measured mean radial velocity jumps sharply from–0.03 m/s at a radial distance of 12 mm to 1.29 m/s at a radial distance of 14 mm. However, the RMS values for these two locations are identically 2.32 m/s, which is nearly twice the mean value. This strange variability in the radial velocity, as well as the deviation from the expected logical development of the profile of the radial velocity, suggests that the computational results are more reliable than the experimental results. With this, we consider the computational results valuable here, as they give useful insights about the development of the velocity field within the HM1 flame, with a high resolution and reproducibility that are not attainable in experimental settings. At the three later axial distances (30 mm, 45 mm, and 65 mm), the deviation between the computational results and the experimental results is more obvious, although there is still a weak qualitative similarity. The high scatter in the experimental results favors the computational results due to perceived high uncertainty in the measurements. Unlike the axial velocity component, the radial profiles of the radial velocity component are not simply decaying away from the burner, but it is remarkably reshaped. For example, at the axial station of 65 mm downstream the burner, the radial velocity is positive near the centerline, while it changes sign once and becomes negative approximately at a radial distance of 8.6 mm (17.2% of the burner diameter) rather than changing the sign three times at the earlier axial stations at 5 mm, 13 mm, and 30 mm.

The spurious wiggles for the coarse mesh at the axial distance of 13 mm are much weaker than those observed at the axial distance of 5 mm. As observed in the results at the axial distance of 5 mm, these fictitious oscillations in the radial velocity are eliminated when the mesh resolution is upgraded. These numerical ripples with the coarse mesh are very faint at the axial station of 30 mm, while they disappear at the farther axial stations of 45 mm and 65 mm.

4.5. Radial Profiles of the Methane Mass Fraction

Moving from velocity variables to chemical variables, we present in the current subsection profiles of the reactant methane (CH₄) at two selected intermediate axial stations (13 mm and 30 mm from the burner).

At 13 mm after the burner face (thus, downstream the jet inlet by a distance that is 26% of the diameter of the bluff-body burner), the radial profile for the mass fractions of methane is shown in Figure 14.

It can be noticed that at this early axial station, the center of the jet core (the fuel jet) still maintains the specified inlet mass fraction of 0.8889 (corresponding to a mole fraction of 50%). This initial high concentration drops steeply and monotonically as a result of the combined effect of chemical reaction and spatial diffusion. The experimental data points are in satisfactory agreement with the computational results. The mesh resolution is nearly insignificant, where the coarse resolution and the fine resolution give very similar profiles. Both experimental and computational profiles show that methane is confined to the burner’s radial extent of 25 mm at this axial station.

At the farther downstream location of 30 mm (60% of the burner diameter), the developed radial profile of the methane mass fraction is shown in Figure 15. This profile is less steep, reflecting a diffusive process. The computational results are in reasonable agreement with the experimental data points, although the computational model predicts a steeper profile (weaker spreading) than the experimental profile. When compared to the earlier station at 13 mm, it is observed that methane is confined to a smaller radial extent of nearly 20 mm at this axial station.

4.6. Radial Profiles of Water Vapor Mass Fraction

As a representative product species, water vapor (H₂O) variations in the HM1 flame are presented in the current subsection for the same two intermediate axial stations considered in the previous subsection (which was dedicated to the reactant species of methane, CH₄). We deliberately keep the same scale of the mass fraction here as the one used to fit the range of mass fraction for methane in the previous subsection (this unified mass fraction range is from 0 to 0.9) for effective visualization and contrast of the concentration for the two gaseous species.

At the earlier axial station located 13 mm after the burner face, the experimental and computational mass fractions of water vapor (H₂O) are compared in Figure 16. Water vapor is absent in the fuel jet core area, which is still occupied by the fuel stream at this axial station. A distinct hat-shaped profile is visible in the figure. This feature strongly manifests the role of the recirculation zone that is established behind the bluff-body face. In that zone, hot combustion product gases are present, and they are able to sustain the HM1 turbulent nonpremixed flame without undesirable instabilities.

There is an excellent agreement between the computational profiles and the experimental data points. Furthermore, the two computational profiles corresponding to two mesh resolutions give almost identical solutions that are not easily distinguishable. This testifies to a successful state of the mesh independence [189].

The mass fraction profiles of H₂O at the farther downstream location of 30 mm are compared in Figure 17. The diffusive and mixing action is discernible in the smoothened hat-shaped portion of the profile, which develops more gradual edges accompanied by a convex inner distribution rather than being nearly flat. The computational predictions reasonably capture the changes in the water vapor profile at this axial stage when compared with the experimental data points. As in the case of the earlier axial stage of 13 mm, the computational predictions are practically mesh-independent, which testifies to the adequacy of the selected mesh resolutions.

4.7. Near-Centerline Axial Profiles

The investigated radial profiles earlier are augmented in the current subsection with additional profiles that are axial (showing variations versus the axial distance from the burner).

These axial profiles are located near the burner centerline. No experimental data were available exactly along the centerline. We select here a close location that is 1.1 mm from the centerline. This location is within the projection of the fuel jet nozzle (whose radial range is from 0 to 1.8 mm), and it corresponds to a relative coordinate of 0.022 in terms of the bluff-body burner diameter of 50 mm (0.044 in terms of the bluff-body burner radius of 25 mm).

We select three key variables, which are the mass fractions of methane (CH₄) as a representative reactant species, the mass fraction of water vapor as a representative product species (H₂O), and the temperature. All experimental axial profiles shown here are Favre-averaged data of the HM1 flame.

For these axial profiles, the experimental data are limited to six axial positions (downstream distances from the burner), which are 13 mm, 30 mm, 45 mm, 65 mm, 90 mm, and 120 mm. This range extends to 2.4 times the burner diameter, which is not sufficient to show the full changes within the HM1 flame near the burner centerline. For the computational results, we extend the visualized profiles to 400 mm (eight times the burner diameter), which is the entire axial extent of the computational domain after the burner. We found that this range is large enough to show the full evolution within the HM1 flame, where the changes nearly vanish at the end of that large range.

The axial profiles of the methane (CH₄) mass fraction as obtained from our two computational solutions and as recorded from the third-party benchmarking experimental measurements are compared in Figure 18. There is good agreement between the experimental data and the computational results, which demonstrate the ability of the solver in capturing the decline in the methane concentration downstream of the inlet nozzle. At this location, methane extends axially to about 350 mm (seven burner diameters). Initially, the methane leaving the burner nozzle at a mass fraction of 0.8889 is consumed (in terms of the mass proportions) at a roughly uniform rate. However, this rate decelerates further downstream, which is a logical consequence of the diminishing concentration of reactants there. Both mesh resolutions have similar results, which is a favorable observation.

The axial profiles of the water vapor (H₂O) mass fraction as obtained from our two computational solutions and as recorded from the third-party benchmarking experimental measurements are compared in Figure 19. We adopted the same mass fraction range of (0–0.9) in this figure as was done for the methane (CH₄) mass fraction in the previous figure. This unification of the scale helps in establishing an effective visual comparison between the two species. Starting at a zero fraction, water vapor is produced gradually due to the combustion process, until its mass fraction approximately saturates near 0.15 at an axial distance of about 200 mm. Like the case of methane, the change in the water vapor mass fraction is roughly linear at the beginning (although it is a positive change for water vapor but a negative change for methane), but then it decays following a nonlinear decelerating trend. As in the case of methane, the agreement between the experimental data points and the computational profiles is reasonable.

The axial profiles of the temperature as obtained from our two computational solutions and as recorded from the third-party benchmarking experimental measurements are compared in Figure 20. The temperature increases from 25 °C (298 K), which is compatible with the inlet boundary condition for the fuel jet, to about 1,625 °C (1,908 K) at the outlet boundary. There is an overall consistency between the profiles obtained under the two resolutions. However, a mild underprediction of the computational models is observed between 30 mm and 120 mm. There are no sufficient measurements to decide if this deviation is suppressed at subsequent axial stations, although visual inspection suggests this because the underprediction at the axial distance of 65 mm is about 167 °C, which decreases to 155 °C at 90 mm, and decreases again to 141 °C at 120 mm. This deviation is ameliorated when knowing that the RMS for the Favre-average temperature measurements at these three spatial points are 91 °C, 89 °C, and 125 °C, respectively. These spread levels in the experimental data are comparable with the underprediction deviations themselves.

4.8. Quantified Deviations

The one-dimensional profiles displayed earlier are useful in identifying not only the agreement level between the experimental data and either of the computational solutions (corresponding to either of the two mesh resolutions), but also in explaining the behavior of different variables within the domain of the HM1 flame, such as the recirculation induced successfully by the bluff-body design.

In the current subsection, we aim to improve the assessment of the deviation between the experimental data and the computational data by listing scalar statistical values that aggregate the deviations between individual data points for a given profile. Such a single-value representation of the deviations permits quicker and easier evaluation of the reactingFoam solver in handling the HM1 flame, as well as more succinct evaluation of the gain achieved when the resolution is improved by nearly tripling the number of cells in the fine mesh (37,202 cells) compared to the coarse mesh (12,836 cells).

We use two deviation metrics that are commonly used for quantifying spread within the same dataset or deviations between two data points, which are the root mean square (RMS) [190] and the mean absolute deviation (MAD) [191]. The RMS deviation metric may also be called root mean square deviation (RMSD) [192] or root mean square error (RMSE) [193].

For two sets of data points (one being the experimental results and the other being computational results with either the coarse mesh or the fine mesh), the MAD metric is computed as the arithmetic mean of the absolute deviation values, or

$$\mathrm{R}\mathrm{M}\mathrm{S}={\displaystyle \frac{1}{\mathrm{n}}}\sum _{i=1}^n\left|d_i\right|$$

(12)

where $d_i$ represents the individual deviation values, $\left|d_i\right|$ represents the absolute value of these individual deviations, and $n$ represents the count of deviation values.

The RMS metric is computed as the square root of the mean squared deviation values, or

$$\mathrm{R}\mathrm{M}\mathrm{S}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{d}_i^2}$$

(13)

For fair and meaningful evaluation of the deviations, the computational data are interpolated using the cubic spline interpolation technique [194] such that their interpolated values corresponding exactly to the location (radial or axial) of the reported experimental data are used in computing the deviation (rather than using the closest available computational value, for example).

Table 2 summarizes the MAD and RMS metrics for the 17 profiles covered in the current study. Because there are two deviation metrics (MAD and RMS) per profile, and there are two solutions (coarse and fine) per profile, the total number of deviation metrics in the table is 68 (17 × 2 × 2).

While the listed deviation metrics in the previous table are informative, it might be better to normalize them such that they convey relative deviations rather than absolute deviations. With this scaling, the influence of the size of the variable is suppressed, and the normalized deviations become of comparable magnitudes. To achieve this, we use the maximum value of the experimental data point of the reach profile as the reference value for normalization. These reference values for all 17 profiles are always positive numbers here.

The 17 reference scaling values for the 17 profiles, and the corresponding 78 transformed (normalized) deviations are listed in Table 3. The %MAD in the table is the ratio of the MAD metric to the maximum experimental value for the given profile. Similarly, the %RMS in the table is the ratio of the RMS metric to the maximum experimental value for the given profile. With the normalization process, it becomes easier to judge the overall deviation for each profile.

It seems valid to exclude the radial velocity while assessing the overall performance of a computational solver, given the large uncertainty and the detected discrepancy in their measurements. If the radial velocity is excluded, then more than half of the relative deviations (either %MAD or %RMS) are within 10%.

Furthermore, if the normalized deviation metrics are averaged over the 12 profiles (the five radial velocity profiles are excluded from the 17 profiles), then the %MAD metric is approximately 8%, and the %RMS deviation is approximately 10%. These are small levels of deviation.

For the %MAD, using the fine grid (37,202 cells) results in a minor improvement with a reduction in the %MAD by only 0.05% percentage points compared to the coarse grid (12,836 cells). For the %RMS, using the fine grid results also in a minor improvement with a reduction in the %RMS by only 0.10% percent points compared to the coarse grid. Thus, the small gain made by increasing the mesh resolution is not attractive enough to justify the incurred increase in the computational demands. In other words, the coarse grid is considered sufficient.

In summary, the reactingFoam solver is considered successful in handling turbulent non-premixed flames.

5. Conclusion

Using the OpenFOAM computational fluid dynamics (CFD) tools and the reactingFoam solver, we constructed a computational model for turbulent non-premixed flames (TNF). We validated this solver through qualitative and quantitative comparisons of line profiles corresponding to the HM1 bluff-body high-speed benchmarking flame. This HM1 flame is an air-fuel open diffusion flame with the fuel being a gaseous mixture of methane and molecular hydrogen with equal volumetric proportions. The flame is axisymmetric, and it is stabilized using a bluff-body burner.

In the computational model, the unsteady Favre-averaged Navier-Stokes equations were integrated in time and space until a steady-state solution was reached. The solution gave the predicted distribution of the various scalar and vector fields. These fields were sampled and post-processed to obtain radial and axial profiles that were compared with the experimental results for assessing the accuracy of the modeling framework.

The simulations were performed twice, with two different resolutions of the multi-block structured meshes. This allowed for estimating the mesh-sensitivity of the computational results. A total of 17 profiles were investigated, covering five variables (axial velocity, radial velocity, methane mass fraction, water vapor mass fraction, and temperature). Our evaluative comparisons suggest that the open-access open-source reactingFoam solver is an adequate tool for modeling combustion problems and diffusion flames. For the HM1 flame, a grid with about 12,000 cells is sufficient to perform reasonable axisymmetric simulations. Although a finer resolution is recommended for better results for the radial velocity in particular.

Based on the findings of the current study, the reactingFoam solver is suitable for use in designing burners and analyzing combustion problems, in research or industrial applications. Also, the findings of the current study show that HM1 flame is an excellent test case for tuning or evaluating CFD models with reactive flows, although special attention should be paid to the potential inaccuracy in the radial velocity.