Pressurized turbulent premixed CH4/H2/air flame validation using OpenFOAM

Guerrero, J.

doi:10.1063/5.0098715

In this study, we use the Algebraic Flame Surface Wrinkling (AFSW) model to conduct numerical simulations of the PSI (Paul Scherrer Institute) high pressure, turbulent premixed Bunsen flame experiments. The AFSW model was implemented in the open-source computational fluid dynamics solver OpenFOAM, and the numerical simulations were performed using a two-dimensional axial-symmetric model with the standard k–ɛ turbulence model with wall functions. The numerical simulations were performed for two different fuel/air mixtures, namely, 100% CH₄ volumetric ratio and 60% CH₄ + 40% H₂ volumetric ratio. The thermophysical and transport properties of the mixture were calculated as a function of temperature using the library Cantera (an open-source suite of tools for problems involving chemical kinetics, thermodynamics, and transport processes), together with the GRI-Mech 3.0 chemical mechanism. It was found that the outcome of the AFSW model implemented in OpenFOAM was in good agreement with the experimental results, quantitatively and qualitatively. Further assessment of the results showed that, as much as the chemistry, the turbulence model and turbulent boundary/initial conditions significantly impact the flame shape and height.

I. INTRODUCTION

The study of premixed turbulent combustion is an area of active research as mastering this technology can directly translate into increased efficiency and reduced NOx and other pollutant emissions. However, modeling premixed combustion phenomena is a non-trivial task because of the interaction between turbulence and chemical reactions. To overcome this challenge, many combustion models have been developed and tested, to name a few, the eddy breakup model,¹ turbulent flame speed closure,² the coherent flamelet model,³ the level set approach for corrugated flamelet regimes (G-equation models),^4,5 and the algebraic flame surface wrinkling model.⁶

Hereafter, we perform numerical simulations of the well-known Paul Scherrer Institute (PSI) high pressure, turbulent premixed Bunsen CH₄/H₂/air flame experiments.^7,8 The combustion process is modeled using the Algebraic Flame Surface Wrinkling (AFSW) model. The AFSW model is an algebraic model originally derived by Muppala et al.⁶ through curve fitting of the Kobayashi experiments on turbulent flame speed measurements for methane and propane flames.⁹ The model was further improved by Dinkelacker et al.,¹⁰ where the authors included an effective Lewis number term to extend its applicability to blended hydrogen mixtures.

We aim to reproduce the experimental results presented in Ref. 10 by using the open-source numerical library OpenFOAM,^11–13 where the AFSW model has been implemented by modifying the solver XiFOAM distributed with OpenFOAM. The implementation of the AFSW model in OpenFOAM’s XiFOAM solver requires the modification and renaming of two files in the source code of the numerical library (XiFoam.C → myXiFoam.C and bEqn.H → mybEqn.H) where the algebraic model is implemented. It also requires the creation of two extra files (AFSW.H and myCreateFields.H) where the new variables are declared. In the software repository,¹⁴ the interested reader can find the modified solver and the input files needed to reproduce the results presented in this manuscript. The same repository serves as a revision control system.

II. THEORETICAL BACKGROUND

A. Premixed combustion model

Under the assumption of simple one-step chemistry with the unity Lewis number and adiabatic conditions, the species transport equations can be reduced to a single combustion progress variable equation, as follows:^2,15,16

\frac{\partial (\bar{ρ} \tilde{c})}{\partial t} + \frac{\partial (\bar{ρ} {\tilde{u}}_{i} \tilde{c})}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} [(\bar{ρ} α + \frac{μ_{t}}{S c_{t}}) \frac{\partial \tilde{c}}{\partial x_{i}}] + {\bar{\dot{ω}}}_{c} .

(1)

In Eq. (1), the transported quantity c refers to the normalized mass fraction of the products,² which we call the progress variable. The combustion progress variable describes the thermochemical state of the mixture at any point in space and time (products and reactants). A detailed derivation of Eq. (1) can be found in Refs. 5 and 15–18.

The solution of Eq. (1), together with additional closure models (turbulence, reaction rate source term, and turbulent flame speed) and the thermophysical and transported properties of the unburnt/burnt mixture and flame, gives the propagation of the premixed flame.

It is worth mentioning that OpenFOAM, instead of solving for the progress variable c, solves for the regress variable b,¹⁹ where

b = 1 - c .

(2)

Henceforth, Eq. (1) can be written in terms of the regress variable b as follows:

\frac{\partial (\bar{ρ} \tilde{b})}{\partial t} + \frac{\partial (\bar{ρ} {\tilde{u}}_{i} \tilde{b})}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} [(\bar{ρ} α + \frac{μ_{t}}{S c_{t}}) \frac{\partial \tilde{b}}{\partial x_{i}}] - {\bar{\dot{ω}}}_{b} .

(3)

In Eqs. (1) and (3), the overbar $̄$ and the tilde $̃$ refer to the Reynolds and Favre averaging, respectively. The molecular thermal diffusivity α appearing in Eqs. (1) and (3) is usually much smaller than the turbulent diffusivity μ_t/Sc_t and can simply be neglected in RANS/URANS simulations.¹⁵

In the equation governing the regress variable b [Eq. (3)], b is bounded between zero and one, where b = 0 represents the products and b = 1 represents the reactants. Equation (3) is the equation used in OpenFOAM and, in particular, in the solver XiFoam (the solver used in this study). OpenFOAM’s XiFoam solver is a transient solver for compressible premixed and partially premixed combustion with turbulence modeling.

In Eq. (3), evaluating the mean reaction rate source term $\bar{\dot{ω}}$ is the main problem in modeling premixed turbulent combustion. This term can be modeled using algebraic methods or methods based on additional transport equations. In this study, we use an algebraic model, namely, the AFSW model.^6,10,20,21 By using this model, the reaction rate source term for the regress variable can be expressed as follows:

{\bar{\dot{ω}}}_{b} = ρ_{u} S_{t} | \nabla \tilde{b} |,

(4)

where $| \nabla \tilde{b} |$ represents the magnitude of the gradient and can be evaluated as follows:

| \nabla \tilde{b} | \equiv {[\sum_{i = 1}^{3} {(\frac{\partial \tilde{b}}{\partial x_{i}})}^{2}]}^{1 / 2} .

(5)

The turbulent flame speed S_t appearing in Eq. (4) is then modeled using the AFSW model^6,10,20,21 as follows:

S_{t} = S_{L 0} + S_{L 0} \frac{0.46}{L e_{e f f}} {R e_{t}}^{0.25} {(\frac{u^{'}}{S_{L 0}})}^{0.3} {(\frac{P}{1 a t m})}^{0.2},

(6)

where S_L0 is the laminar unstrained flame speed (flame without strain or stretch) and is dependent on the chemistry. It can be calculated using one-dimensional chemical kinetics solvers with detailed chemistry. In order to have a closed-form of Eq. (6), the following closure relationships are used:

\frac{1}{L e_{e f f}} = \frac{D}{α} = \frac{X_{{C H}_{4}} D_{{C H}_{4}}}{α} + \frac{X_{H_{2}} D_{H_{2}}}{α} = \frac{X_{{C H}_{4}}}{L e_{{C H}_{4}}} + \frac{X_{H_{2}}}{L e_{H_{2}}},

(7)

R e_{t} = \frac{u^{'} l_{t}}{ν},

(8)

l_{t} = 0.0 9^{3 / 4} \frac{k^{3 / 2}}{ε},

(9)

u^{'} = \sqrt{\frac{2}{3} k} .

(10)

The effective Lewis number Le_eff appearing in Eq. (7) models the preferential diffusion effect of the species H₂ in the unburnt mixture on the turbulent flame speed [Eq. (6)].¹⁰ In Eq. (7), the mixture molecular thermal diffusivity α is computed using the following equations:

α = \frac{λ}{ρ C_{p}},

(11)

where α is computed at the temperature at which the maximum heat is released (T_mhrr in Table I).

TABLE I.

Premixed mixture properties. The sub-index u refers to unburnt conditions. The sub-index ad refers to adiabatic conditions. The sub-index mhrr refers to the temperature at which the maximum heat is released. The sub-indices CH₄ and H₂ refer to the species. The variable ϕ is the equivalence ratio and is computed using Eq. (23).

$X_{{C H}_{4}}$ (−)	$X_{H_{2}}$ (−)	ϕ (−)	T_u (K)	P (atm)	ρ_u (kg/m³)	μ_u (m²/s)	$D_{{C H}_{4}}$ (m²/s)	$D_{H_{2}}$ (m²/s)	α (m²/s)	S_L0 (m/s)	T_mhrr (K)	T_ad (K)
1.0	0.0	0.5	673	5	2.52	3.25 × 10⁻⁵	7.93 × 10⁻⁵	2.60 × 10⁻⁴	7.56 × 10⁻⁵	0.232	1606	1777
0.6	0.4	0.5	673	5	2.46	3.25 × 10⁻⁵	7.79 × 10⁻⁵	2.55 × 10⁻⁴	8.05 × 10⁻⁵	0.334	1589	1803

$X_{{C H}_{4}}$ (−)	$X_{H_{2}}$ (−)	ϕ (−)	T_u (K)	P (atm)	ρ_u (kg/m³)	μ_u (m²/s)	$D_{{C H}_{4}}$ (m²/s)	$D_{H_{2}}$ (m²/s)	α (m²/s)	S_L0 (m/s)	T_mhrr (K)	T_ad (K)
1.0	0.0	0.5	673	5	2.52	3.25 × 10⁻⁵	7.93 × 10⁻⁵	2.60 × 10⁻⁴	7.56 × 10⁻⁵	0.232	1606	1777
0.6	0.4	0.5	673	5	2.46	3.25 × 10⁻⁵	7.79 × 10⁻⁵	2.55 × 10⁻⁴	8.05 × 10⁻⁵	0.334	1589	1803

In this study, we assumed the turbulent Lewis number Le_t to be equal to one. This implies that the turbulent Schmidt number Sc_t and the turbulent Prandtl number Pr_t are the same. The main consequence of this assumption is that the eddy mass diffusivity and eddy thermal diffusivity are the same. In turbulent combustion modeling, this hypothesis is widely used and often necessary when conducting RANS/URANS simulations.^15,17,18 For completeness, recall that the Le_t, the Sc_t, and the Pr_t numbers are defined as follows:

L e_{t} = \frac{S c_{t}}{\Pr_{t}}, S c_{t} = \frac{ν_{t}}{D_{t}}, P r_{t} = \frac{ν_{t}}{α_{t}} .

(12)

By combining these expressions, we can derive a relation for the turbulent diffusivity D_t [Eq. (13)], where the turbulent Schmidt number Sc_t is assumed to be constant,

D_{t} = \frac{ν_{t}}{S c_{t}} .

(13)

It is worth noting that the evaluation of D_t is not simple. The closure given in Eq. (14) is often used together with the k–ɛ turbulence model, where μ_t is computed from the transported turbulent quantities k and ɛ. However, while it is commonly accepted that the value of the coefficient C_μ is equal to 0.09, different authors have reported different values of Sc_t. In this study, we assumed Sc_t = Pr_t = 0.7,^17,22

D_{t} = \frac{C_{μ}}{S c_{t}} \frac{\bar{k^{2}}}{\bar{ε}} .

(14)

Finally, by substituting Eq. (13) into Eq. (3) and by neglecting the molecular thermal diffusivity α, we obtain the following solvable equation for the regress variable b:

\frac{\partial (\bar{ρ} \tilde{b})}{\partial t} + \frac{\partial (\bar{ρ} {\tilde{u}}_{i} \tilde{b})}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} (\bar{ρ} D_{t} \frac{\partial \tilde{b}}{\partial x_{i}}) - {\bar{\dot{ω}}}_{b} .

(15)

Equation (15), together with the compressible RANS equations, the k–ɛ closure turbulence model, and proper boundary conditions and initial conditions, is solved using OpenFOAM’s XiFoam solver.

B. Energy equation treatment

In OpenFOAM’s XiFoam solver, the energy equation is treated in terms of the absolute enthalpy h_a formulation (the sensible enthalpy plus the chemical enthalpy) or in terms of the internal energy e_i formulation.²³ In this study, we used the absolute enthalpy formulation in OpenFOAM, that is, Eq. (16). In this equation, K is the specific kinetic energy and P is the thermodynamic pressure (absolute pressure). Since h_a already contains the chemical enthalpy, there is no chemical heating source in the following equation:

\begin{align} \frac{\partial (\bar{ρ} {\tilde{h}}_{a})}{\partial t} + \frac{\partial (\bar{ρ} {\tilde{u}}_{i} {\tilde{h}}_{a})}{\partial x_{i}} + \frac{\partial (\bar{ρ} \tilde{K})}{\partial t} + \frac{\partial (\bar{ρ} {\tilde{u}}_{i} \tilde{K})}{\partial x_{i}} \\ = \dots \frac{\partial \bar{P}}{\partial t} + \frac{\partial}{\partial x_{i}} [(\frac{\tilde{μ}}{\Pr} + \frac{{\tilde{μ}}_{t}}{\Pr_{t}}) \frac{\partial {\tilde{h}}_{a}}{\partial x_{i}}] . \end{align}

(16)

C. Thermophysical and transport properties

The thermophysical and transport properties of the unburnt/burnt mixtures and flame were calculated as a function of the temperature using the library Cantera 2.4.0.²⁴ Cantera is an open-source suite of tools for problems involving chemical kinetics, thermodynamics, and transport processes. For the 1D chemistry calculations, we used the library GRI-Mech 3.0.²⁵ GRI-Mech is a database of chemical reactions and associated rate constant expressions capable of the best representation of natural gas flames and ignition. Previous studies have shown that the calculated S_L0 for CH₄/H₂/air mixtures using this mechanism is in good agreement with experimental data.^26,27

In Table I, we list the calculated properties. The variables $α, D_{{C H}_{4}}$ ⁠, and $D_{H_{2}}$ are calculated at the maximum heat release rate temperature (T_mhrr). As can be seen in this table, the addition of H₂ to the mixture increases the laminar flame speed S_L0. This is attributed to higher H₂ reactivity with respect to CH₄.²⁸

The solver XiFoam in OpenFOAM requires the definition of the JANAF table coefficients for C_p and the definition of the coefficients of the Sutherland law for the computation of the dynamic viscosity. The thermal conductivity λ is then computed from the definition of the molecular Prandtl number Pr,

P r = \frac{C_{p} μ}{λ} .

(17)

The JANAF table coefficients were calculated using Eqs. (18)–(20).²⁹ In these equations, the superscript ^o refers to standard state pressure conditions or P = 1 atm. In Eqs. (18)–(20), R is the universal gas constant and H and S refer to the standard state enthalpy and entropy, respectively; both defined at a reference temperature of T_ref = 298 K,

\frac{C_{p}^{o}}{R} = a_{0} + a_{1} T + a_{2} T^{2} + a_{3} T^{3} + a_{4} T^{4},

(18)

\frac{H^{o}}{R T_{r e f}} = a_{0} + \frac{a_{1}}{2} T_{r e f} + \frac{a_{2}}{3} T_{r e f}^{2} + \frac{a_{3}}{4} T_{r e f}^{3} + \frac{a_{4}}{5} T_{r e f}^{4} + \frac{a_{5}}{T_{r e f}},

(19)

\frac{S^{o}}{R} = a_{0} \ln T_{r e f} + a_{1} T_{r e f} + \frac{a_{2}}{2} T_{r e f}^{2} + \frac{a_{3}}{3} T_{r e f}^{3} + \frac{a_{4}}{4} T_{r e f}^{4} + a_{6} .

(20)

The coefficients a₀, a₁, a₂, a₃, and a₄ in Eq. (18) were obtained by curve fitting at different temperature values (T). Once a₀, …, a₄ have been determined, they are substituted in Eqs. (19) and (20) at T_ref. At this point, the coefficient a₅ appearing in Eq. (19) (enthalpy jump) and the coefficient a₆ appearing in Eq. (20) (entropy jump) can be computed. These coefficients were calculated for low- (200 to 1000 K) and high- (1000 to 6000 K) temperature ranges for both unburnt and burnt mixtures.

The dynamic viscosity of the mixture was computed using the Sutherland law,³⁰ which is expressed as follows:

μ = A_{s} \frac{T^{0.5}}{1 + T_{s} / T},

(21)

where μ is the temperature-dependent mixture dynamic viscosity, A_s and T_s are the model coefficients, and T is the varying mixture temperature. Curve fitting was performed under varying temperatures at the operating pressure. In the cases presented in this study, the operating pressure corresponded to the experimental conditions and was equal to 5 atm.^7,8,10

In XiFoam, the coefficients C_p and μ must be defined separately for the unburnt and the burnt mixtures. Additionally, the molar weights of the unburnt and the burnt mixture must be defined. A Python-based Cantera script was written to achieve this, and the required properties were calculated with the GRI-Mech 3.0 mechanism. For completeness, the script is provided in the software repository.¹⁴

The final two parameters, stoichiometric air-to-fuel mass ratio AFR_st and equivalence ratio ϕ, can be calculated using Eqs. (22) and (23), respectively. In Eq. (23), the sub-index fuel refers to the fuel mixture. In this study, we simulated two different fuel mixture cases: one case corresponding to a species volumetric ratio equal to 100% CH₄ (case 1) and the other case corresponding to a species volumetric equivalent to 60% CH₄ + 40% H₂ (case 2),

A F R_{s t} = {(\frac{m_{a i r}}{m_{fu e l}})}_{s t},

(22)

ϕ = \frac{\frac{m_{fu e l}}{m_{o x i d i z e r}}}{{(\frac{m_{fu e l}}{m_{o x i d i z e r}})}_{s t}} = \frac{\frac{n_{fu e l}}{n_{o x i d i z e r}}}{{(\frac{n_{fu e l}}{n_{o x i d i z e r}})}_{s t}} .

(23)

III. DESCRIPTION OF THE EXPERIMENTAL CONFIGURATION AND NUMERICAL BACKGROUND

A. Experimental setup and data

The experimental data referred to in this study were taken from Refs. 7, 8, and 10. In Fig. 1, the Bunsen-type high-pressure generic burner used in the experimental setup is depicted. Turbulence was created using a grid placed at 50 mm upstream of the inlet in the experimental test rig. The turbulent fluctuations u′ and length scales l_t were measured at the inlet using the Particle Image Velocimetry (PIV) method. In Table II, the operating conditions of the experiments are listed. Further information about the experimental setup can be found in Refs. 7, 8, and 10.

FIG. 1.

View large Download slide

Schematic view of the PSI generic, high pressure Bunsen burner. Image taken from Ref. 8. P. Siewert, Ph.D. thesis, ETH Zurich, 2006; licensed under a Creative Commons Attribution (CC BY) license.

TABLE II.

Experiment conditions.^7,8,10

T $(K)$	P $(a t m)$	ϕ $(-)$	u_inlet $(m / s)$	u′ $(m / s)$	L_t $(m m)$
673	5	0.5	40	6.16	2.2

B. Computational domain and mesh

In Figs. 2 and 3, we depict the computational domain and the mesh, respectively. The domain simulated corresponds to the region with dimension quotes in Fig. 1. The two-dimensional mesh is made up of 76 220 hexahedral elements. The minimum cell side length of the hexahedral cells was equal to 0.25 mm. This mesh criterion was taken from Ref. 10, where a mesh independence study was carried out. To conduct the simulations, we used two-dimensional axial symmetry.

FIG. 2.

View large Download slide

Computational domain and dimensions. The figure is not to scale.

FIG. 3.

View large Download slide

Visualization of the two-dimensional mesh.

C. Boundary and initial conditions

In this study, we used the standard k–ɛ model with wall functions;^31,32 therefore, we need to provide proper boundary conditions and initial conditions for the turbulent variables. In Table III, we list the boundary conditions used in OpenFOAM for all field variables (including the wall functions). All the field variables were initialized using uniform values. In all simulations, the flow was assumed to be adiabatic. That is, there is no heat or mass transfer between the computational domain and the external environment. The initial estimates of the turbulent quantities k and ɛ were obtained using Eqs. (24) and (25), and the values of u′ and l_t are listed in Table II. A fully developed velocity profile was defined at the inlet boundary (plotted in Fig. 4), which was obtained from a separate two-dimensional axisymmetric pipe flow precursor simulation,

k = \frac{3}{2} u^{' 2},

(24)

ε = 0.0 9^{3 / 4} \frac{k^{3 / 2}}{l_{t}} .

(25)

TABLE III.

Boundary conditions used in OpenFOAM. Note that the type of the boundary conditions is specific to OpenFOAM. FV = fixedValue, C = Calculated, ZG = zeroGradient, NS = noSlip, UDP = User-Defined velocity profile (plotted in Fig. 4), KWF = kqRWallFunction, EWF = epsilonWallFunction, NKWF = nutkWallFunction.

Variable	Inlet	Outlet	Wall
k	FV	ZG	KWF
ɛ	FV	ZG	EWF
ν_t	C	ZG	NKWF
U	UDP	ZG	NS
p	ZG	FV	ZG
T	FV	ZG	ZG
b	FV	ZG	ZG

FIG. 4.

View large Download slide

Inlet velocity profile computed using a precursor simulation. The inlet extreme is located at 12.5 mm from the pipe axis line.

D. Numerical setup and discretization schemes

The governing equations were solved using the collocated finite volume method in OpenFOAM.^12,13 As the solution takes place in collocated meshes, the Rhie–Chow interpolation method is used to prevent the pressure checkerboard instability.³³ OpenFOAM’s XiFoam solver is a transient solver for compressible premixed and partially premixed combustion with turbulence modeling. To conduct the numerical simulations with XiFoam, we used the least-squares cell-based method for gradient discretization. In order to prevent spurious oscillations, a multi-dimensional gradient limiter was used.³⁴ The diffusive terms were discretized using a second-order centered differences scheme with non-orthogonality corrections. This corrected numerical scheme for the diffusive terms takes into account mesh non-orthogonality and mesh stretching. The convective terms appearing in the governing equations were discretized using a limited linear scheme. This scheme tends toward upwind in regions of rapidly changing gradients. In this discretization scheme, the amount of upwinding can be controlled by using a blending coefficient bounded between zero and one, where one means aggressive limiting (upwinding) and zero means pure linear scheme (very accurate but oscillatory).

For the time discretization, we used OpenFOAM’s Crank–Nicolson formulation. This formulation uses a blending coefficient bounded between zero and one. A value of zero is equivalent to the first-order Euler scheme (stable but diffusive), and a value of one is equivalent to the pure Crank–Nicolson scheme (second-order accurate but slightly oscillatory). By setting the blending coefficient to a value between 0.7 and 0.9, we can obtain a good compromise between stability and accuracy. The time step was chosen in such a way so that the Courant number does not exceed a value of three, which was found to be a good compromise among stability, accuracy, boundedness, and convergence rate. To obtain extra stability, we used explicit and implicit under-relaxation, where all the Under-Relaxation Factors (URF) were set to 0.9.

The pressure–velocity coupling was achieved by means of the iterative PISO algorithm³⁵ (PIMPLE method in OpenFOAM), where we used at least five corrector steps and one outer-corrector step. During the iterative marching, the number of correction steps was increased in order to stabilize the solution and improve the accuracy of the results. This added extra stability and accuracy comes at a higher computational cost because we need to solve the pressure, momentum, energy, and additional transport equations (turbulence equations, regress variable equations, etc.) several times. As a side note, all the simulations were started using a stable method (first-order accuracy), and, as the solution progressed, we switched to a more accurate but potentially unbounded method (second-order accuracy).

Finally, the convergence of the solution was judged based on the shape of the flame (the brush). As soon as a steady flame shape (or brush) was achieved, the simulation was stopped. This requires the computation and continuous monitoring of the average field values during the simulation.

IV. RESULTS AND DISCUSSION

In this section, we present the numerical results obtained using the implementation of the AFSW model in OpenFOAM. In addition, we also compare OpenFOAM’s results with the results obtained using a similar implementation of the AFSW model in Ansys Fluent.³⁶ In order to make a fair comparison between both solvers, we used the same mesh, with equivalent boundary conditions and initial conditions, and similar numerical discretization schemes and iterative marching techniques. For a more complete discussion regarding the implementation of the AFSW model in OpenFOAM and Fluent, and the main differences between both computational fluid dynamics (CFD) solvers, the interested reader is invited to consult Ref. 37. All the necessary files needed to reproduce the numerical results presented in this work can be found in Ref. 14.

Figure 5 shows the comparison of the flame shapes. The experimental flame shapes were obtained using OH* chemiluminescence imaging, as described in Ref. 8. The numerical flame shapes were obtained using the scalar variable $\tilde{c}$ (Favre-averaged progress variable), which was converted to the scalar variable $\bar{c}$ (Reynolds-averaged progress variable) by means of Eq. (26).¹⁵ This approach was adopted since Reynolds averaging provides the mean value over a time sequence, similar to the averaging in the experiments, where simultaneous images are averaged to obtain the mean flame shape,

τ = \frac{ρ_{u}}{ρ_{b}} - 1, \bar{c} = \frac{(1 + τ) \tilde{c}}{1 + τ \tilde{c}} .

(26)

FIG. 5.

View large Download slide

Comparison of flame shapes (brush). Top row: experimental normalized OH* chemiluminescence.⁸ P. Siewert, Ph.D. thesis, ETH Zurich, 2006; licensed under a Creative Commons Attribution (CC BY) license. Middle row: numerical flame shapes visualization using $\bar{c}$ (Ansys Fluent results). Bottom row: numerical flame shapes using $\bar{c}$ (OpenFOAM results). The left column of all rows corresponds to the case 100% CH₄. The right column of all rows corresponds to the case of 60% CH₄ + 40% H₂.

From Fig. 5, we can observe that the AFSW model implemented in OpenFOAM provides a reasonably good qualitative agreement with the experimental results and the numerical results obtained with Ansys Fluent. The AFSW model captures the flame shortening effect due to the Le_eff and S_L0 terms appearing in Eq. (6). As the H₂ content increased, S_L0 increased as well (refer to Table I). Le_eff decreased due to high mass diffusivity $(D_{H_{2}})$ ⁠; consequently, S_t increased.

In Fig. 6, we plot the $\bar{c}$ distribution along the axisymmetric axis. In the figure, a value of $\bar{c} = 0.5$ (the horizontal green line) corresponds to the mid-isosurface between the unburned state (reactants) and the burned state (products).

FIG. 6.

View large Download slide

Plot of the $\bar{c}$ distribution over the axisymmetric axis. In the figure, the horizontal green line represents $\bar{c} = 0.5$ ⁠. A value of $\bar{c} = 0.5$ corresponds to the mid-isosurface between the unburned state and the burned state. The intersection of this line with the axial $\bar{c}$ distributions (the numerical results) was used to compute the flame height.

In Table IV, we list the values of the flame height at $\bar{c} = 0.5$ $(h_{\bar{c} = 0.5})$ for the different conditions studied. The flame height was calculated from Fig. 6 by projecting the intersection points of the line $\bar{c} = 0.5$ (the horizontal green line in Fig. 6) and the numerical $\bar{c}$ distributions over the abscissa. Note that the values shown in Table IV correspond to the values sampled along the axisymmetric axis (refer to Fig. 2).

TABLE IV.

Numerical and experimental flame height values $h_{\bar{c} = 0.5}$ ⁠. The numerical values are estimated at $\bar{c} = 0.5$ ⁠.

Case	Volumetric ratio $(%)$	OpenFOAM $(m m)$	Fluent $(m m)$	Experiment $(m m)$
1	100% CH₄ + 0% H₂	103.74	86.35	140.36
2	60% CH₄ + 40% H₂	58.98	51.56	50.16

Case	Volumetric ratio $(%)$	OpenFOAM $(m m)$	Fluent $(m m)$	Experiment $(m m)$
1	100% CH₄ + 0% H₂	103.74	86.35	140.36
2	60% CH₄ + 40% H₂	58.98	51.56	50.16

Comparing OpenFOAM’s XiFoam solver results against the outcome of Ansys Fluent, we can observe that OpenFOAM slightly over-predicts the flame height in reference to Ansys Fluent. For the case of 100% CH₄, both solvers predicted shorter flames than the experiments. For the case of 60% CH₄ + 40% H₂, OpenFOAM slightly over-predicted the flame height, whereas Ansys Fluent better predicted the flame height (refer to Table IV), always in reference to the experiments.

In Table V, we list the values of the flame speed at $\bar{c} = 0.5$ $(S_{t, \bar{c} = 0.5})$ ⁠. The flame speed values shown in Table V were calculated using Kobayashi’s method.⁶ Comparing Tables IV and V, it can be evidenced that $h_{\bar{c} = 0.5}$ and $S_{t, \bar{c} = 0.5}$ are inversely proportional. That is, larger flame heights results in smaller flame speeds, and vice versa.

TABLE V.

Numerical and experimental flame speed values $S_{t, \bar{c} = 0.5}$ ⁠. The flame speed was calculated using Kobayashi’s method.⁶

Case	Volumetric ratio $(%)$	OpenFOAM $(m / s)$	Fluent $(m / s)$	Experiment $(m / s)$
1	100% CH₄ + 0% H₂	4.78	5.73	3.54
2	60% CH₄ + 40% H₂	8.29	9.42	9.67

Case	Volumetric ratio $(%)$	OpenFOAM $(m / s)$	Fluent $(m / s)$	Experiment $(m / s)$
1	100% CH₄ + 0% H₂	4.78	5.73	3.54
2	60% CH₄ + 40% H₂	8.29	9.42	9.67

In Ref. 38, it is stated that the first term on the RHS of Eqs. (1) and (3) controls the flame brush thickness and, the second term (which is a function of $S_{t, \bar{c} = 0.5}$ ⁠) controls the flame length. Based on this statement, the good agreement of the flame brush thicknesses with the experiments confirmed the validity of the turbulent Schmidt number used in the CFD simulations (Sc_t = Pr_t = 0.7).

As seen in Eqs. (1) and (3), the RHS depends on the turbulence model via the scalar variables μ_t and Sc_t. Therefore, discrepancies between Ansys Fluent and OpenFOAM can be attributed to the specific implementation of the turbulence model and the additional corrections implemented. On the other hand, the turbulence flame speed S_t, responsible for the turbulence chemistry interaction, strongly depends on the turbulence model through the u′ and l_t terms. Hence, the turbulence model is of paramount importance when modeling combustion.

Finally, it is important to note that it is difficult to do an honest comparison of the results shown in Fig. 5, mainly because the experimental image is not 100% equivalent to the numerical images. We can do a similar statement regarding the results presented in Fig. 6, where we used the line $\bar{c} = 0.5$ to compute the flame height $h_{\bar{c} = 0.5}$ ⁠. Slightly smaller or larger values of the reference $\bar{c}$ value can result in different flame heights.

V. CONCLUSIONS AND PERSPECTIVES

In this study, we implemented the Algebraic Flame Surface Wrinkling (AFSW) model in the CFD solver OpenFOAM. We also validated the performance of the model using the Paul Scherrer Institute (PSI) high-pressure, turbulent premixed Bunsen flame experiments.⁸ It was found that the outcome of the numerical simulations obtained with OpenFOAM was in good agreement with the experimental results described in Ref. 8, quantitatively and qualitatively speaking. Further assessment of the results confirms that, as much as the chemistry, the turbulence model and turbulent conditions (boundary and initial conditions) significantly impact the flame height and speed.

To further improve the results presented in this manuscript and to reduce the uncertainty related to the numerical simulations and physical models, we envisage extending the numerical simulations to three-dimensional domains (abandoning in this way the axial symmetry hypothesis), a better estimation of the turbulent conditions (boundary and initial conditions), the inclusion of the turbulence grid upstream of the domain (as illustrated in Fig. 1), and the use of non-linear eddy viscosity models. To facilitate the monitoring of the solution, we also expect the use of image and pattern recognition methods to compare in real time a target image (experimental results) with a source image (numerical results), similar to the workflow presented in Ref. 39.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

DATA AVAILABILITY

The data that support the findings of this study are openly available in figshare at http://doi.org/10.6084/m9.figshare.16774336, Ref. 14.

NOMENCLATURE

Notation

a: JANAF table coefficient, non-dimensional
AFR: air-to-fuel ratio, non-dimensional
A_s: Sutherland model coefficient, measured in $kg/m s \sqrt{K}$
b: regress variable, non-dimensional
c: progress variable, non-dimensional
CH₄: methane species
C_p: specific heat, measured in J/kg K
C_p^o: standard state specific heat, measured in J/kg K
D: binary mass diffusivity, measured in m²/s
h: enthalpy (specific), measured in J/kg
$h_{\bar{c} = 0.5}$: flame height at $\bar{c} = 0.5$ ⁠, measured in mm
H₂: hydrogen species
H^o: standard state enthalpy (specific), measured in J/kg
k: turbulent kinetic energy, measured in m²/s²
K: specific kinetic energy, measured in m²/s²
Le: Lewis number, non-dimensional
l_t: turbulence length scale, measured in m
m: mass, measured in kg
n: mole number, measured in mol
N₂: nitrogen species
NOx: nitrogen oxide
O₂: oxygen species
P: absolute pressure, measured in Pa
Pr: Prandtl number, non-dimensional
R: universal gas constant, measured in J/kg K
Re: Reynolds number, non-dimensional
Sc: Schmidt number, non-dimensional
S_L0: unstrained adiabatic laminar flame speed, measured in m/s
S^o: standard state entropy (specific), measured in J/kg K
S_t: turbulent flame speed, measured in m/s
t: time, measured in s
T: temperature, measured in K
T_ref: standard state temperature, measured in K
T_s: Sutherland model coefficient, measured in K
u: velocity, measured in m/s
u′: turbulent fluctuating velocity, measured in m/s
x: spatial coordinates, measured in m
X: molar ratio, non-dimensional
Y_f: fuel mixture mass fraction, non-dimensional
¯: Reynolds average
˜: Favre average
$Δ H_{c}^{o}$: lower heating value of combustion, measured in J/kg
α: thermal diffusivity (molecular), measured in m²/s
ε: turbulence dissipation rate, measured in m²/s³
λ: thermal conductivity, measured in W/m K
μ: dynamic viscosity, measured in kg/m s
ν: kinematic viscosity, m²/s
ρ: density, measured in kg/m³
τ: heat release factor, non-dimensional
ϕ: equivalence ratio, non-dimensional
$\bar{\dot{ω}}$: reaction rate source term, measured in kg/m³ s
${\bar{\dot{ω}}}_{b}$: $\bar{\dot{ω}}$ in the regress variable b equation, measured in kg/m³ s
${\bar{\dot{ω}}}_{c}$: $\bar{\dot{ω}}$ in the progress variable c equation, measured in kg/m³ s
∆: gradient operator

Sub-indices

_a: absolute
_ad: adiabatic conditions
_air: mix of O₂ and N₂ species
_b: burnt conditions
_chem: chemical
${C H}_{4}$: methane species
_eff: effective
_f: fuel
_fuel: 100% CH₄ or 60% CH₄ + 40% H₂
_i: index of the spatial dimension, where i = 1, 2, 3
$H_{2}$: hydrogen species
_t: turbulence related quantity
_u: unburnt conditions

Acronyms

AFSW: algebraic flame surface wrinkling
CFD: computational fluid dynamics
JANAF: joint army, navy, and air force
k–ε: standard k–ε turbulence model
PDF: probability density function
PSI: Paul Scherrer Institute
RANS: Reynolds-averaged Navier–Stokes
RHS: right-hand side
SIMPLE: semi-implicit method for pressure linked equations
UDF: user-defined function
URANS: unsteady Reynolds-averaged Navier–Stokes

REFERENCES

1.

D. B.

Spalding

, “

Development of the eddy-break-up model of turbulent combustion

,”

Symp. (Int.) Combust., [Proc.]

16

(

1

),

1657

–

1663

(

1977

).

https://doi.org/10.1016/s0082-0784(77)80444-x

Google Scholar

Crossref

2.

V.

Zimont

,

W.

Polifke

,

M.

Bettelini

, and

W.

Weisenstein

, “An efficient computational model for premixed turbulent combustion at high Reynolds numbers based on a turbulent flame speed closure,”

J. Eng. Gas Turbines Power

120

,

526

–

532

(

1998

).

https://doi.org/10.1115/1.2818178

Google Scholar

Crossref

3.

T.

Poinsot

,

D.

Veynante

, and

S.

Candel

, “

Diagrams of premixed turbulent combustion based on direct simulation

,”

Symp. (Int.) Combust., [Proc.]

23

(

1

),

613

–

619

(

1991

).

https://doi.org/10.1016/s0082-0784(06)80308-5

Google Scholar

Crossref

4.

N.

Peters

, “A spectral closure for premixed turbulent combustion in the flamelet regime,”

J. Fluid Mech.

242

,

611

–

629

(

1992

).

https://doi.org/10.1017/s0022112092002519

Google Scholar

Crossref

5.

N.

Peters

,

Turbulent Combustion

(

Cambridge University Press

,

2000

).

Google Scholar

Crossref

6.

S. P. R.

Muppala

,

N. K.

Aluri

,

F.

Dinkelacker

, and

A.

Leipertz

, “Development of an algebraic reaction rate closure for the numerical calculation of turbulent premixed methane, ethylene, and propane/air flames for pressures up to 1.0 MPa,”

Combust. Flame

140

,

257

–

266

(

2005

).

https://doi.org/10.1016/j.combustflame.2004.11.005

Google Scholar

Crossref

7.

P.

Griebel

,

E.

Boschek

, and

P.

Jansohn

, “Lean blowout limits and nox emissions of turbulent, lean premixed, hydrogen-enriched methane/air flames at high pressure,”

J. Eng. Gas Turbines Power

129

,

404

–

410

(

2007

).

https://doi.org/10.1115/1.2436568

Google Scholar

Crossref

8.

P.

Siewert

, “Flame front characteristics of turbulent lean premixed methane/air flames at high-pressure,” Ph.D. thesis,

Swiss Federal Institute of Technology Zurich (ETH Zurich)

,

2006

.

9.

H.

Kobayashi

,

T.

Tamura

,

K.

Maruta

,

T.

Niioka

, and

F. A.

Williams

, “Burning velocity of turbulent premixed flames in a high-pressure environment,”

Proc. Combust. Inst.

26

,

389

–

396

(

1996

).

https://doi.org/10.1016/s0082-0784(96)80240-2

Google Scholar

Crossref

10.

F.

Dinkelacker

,

B.

Manickam

, and

S. P. R.

Muppala

, “Modelling and simulation of lean premixed turbulent methane/hydrogen/air flames with an effective Lewis number approach,”

Combust. Flame

158

,

1742

–

1749

(

2011

).

https://doi.org/10.1016/j.combustflame.2010.12.003

Google Scholar

Crossref

11.

See http://www.openfoam.org for more information about the openfoam foundation; accessed May 2022.

12.

H. G.

Weller

,

G.

Tabor

,

H.

Jasak

, and

C.

Fureby

, “A tensorial approach to computational continuum mechanics using object-oriented techniques,”

Comput. Phys.

12

,

620

–

631

(

1998

).

https://doi.org/10.1063/1.168744

Google Scholar

Crossref

13.

C.

Greenshields

and

H.

Weller

,

Notes on Computational Fluid Dynamics: General Principles

(

CFD Direct Ltd.

,

Reading, UK

,

2022

).

Google Scholar

14.

J.

Guerrero

, “

Case repository - Pressurized turbulent premixed CH₄/H₂/air flame validation using OpenFOAM

,”

Figshare

(

2022

).

https://doi.org/10.6084/m9.figshare.16774336.v15

Google Scholar

15.

T.

Poinsot

and

D.

Veynante

,

Theoretical and Numerical Combustion

, 3rd ed. (

Institut de Mécanique des Fluides de Toulouse

,

2012

), http://elearning.cerfacs.fr/combustion/onlinePoinsotBook/buythirdedition/index.php.

Google Scholar

16.

P. A.

Libby

and

F. A.

Williams

,

Turbulent Reacting Flows

(

Springer-Verlag

,

1980

).

Google Scholar

Crossref

17.

A.

Lipatnikov

,

Fundamentals of Premixed Turbulent Combustion

, 1st ed. (

CRC Press

,

2012

).

Google Scholar

Crossref

18.

K. K.

Kuo

and

R.

Acharya

,

Fundamentals of Turbulent and Multiphase Combustion

, 1st ed. (

Wiley

,

2012

).

Google Scholar

19.

H. G.

Weller

,

G.

Tabor

,

A. D.

Gosman

, and

C.

Fureby

, “Application of a flame-wrinkling les combustion model to a turbulent mixing layer,”

Proc. Combust. Inst.

27

,

899

–

907

(

1998

).

https://doi.org/10.1016/s0082-0784(98)80487-6

Google Scholar

Crossref

20.

N. K.

Aluri

,

S.

Muppala

, and

F.

Dinkelacker

, “Large-eddy simulation of lean premixed turbulent flames of three different combustion configurations using a novel reaction closure,”

Flow, Turbul. Combust.

80

,

207

–

224

(

2008

).

https://doi.org/10.1007/s10494-007-9114-2

Google Scholar

Crossref

21.

B.

Manickam

,

J.

Franke

,

S.

Muppala

, and

F.

Dinkelacker

, “Large-eddy simulation of triangular-stabilized lean premixed turbulent flames: Quality and error assessment,”

Flow, Turbul. Combust.

88

,

1573

–

1987

(

2012

).

https://doi.org/10.1007/s10494-011-9385-5

Google Scholar

Crossref

22.

P. K.

Yeung

, “Lagrangian investigations of turbulence,”

Annu. Rev. Fluid Mech.

34

,

115

–

142

(

2002

).

https://doi.org/10.1146/annurev.fluid.34.082101.170725

Google Scholar

Crossref

23.

D.

Iurashev

, Numerical and analytical study of combustion instabilities in industrial gas turbines, Ph.D. thesis,

University of Genoa

,

Italy

,

2012

.

24.

D. G.

Goodwin

,

R. L.

Speth

,

H. K.

Moffat

, and

B. W.

Weber

, “Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes,” version 2.5.1, https://www.cantera.org, accessed May 2022.

25.

G. P.

Smith

,

D. M.

Golden

,

M.

Frenklach

,

N. W.

Moriarty

,

B.

Eiteneer

,

M.

Goldenberg

,

C. T.

Bowman

,

R. K.

Hanson

,

S.

Song

,

W. C.

Gardiner

,

V. V.

Lissianski

, and

Z.

Qin

, Gri-mech 3.0, http://combustion.berkeley.edu/gri-mech/version30/text30.html; accessed May 2022.

26.

N.

Donohoe

,

A.

Heufer

,

W. K.

Metcalfe

,

H. J.

Curran

,

M. L.

Davis

,

O.

Mathieu

,

D.

Plichta

,

A.

Morones

,

E. L.

Petersen

, and

F.

Güthe

, “Ignition delay times, laminar flame speeds, and mechanism validation for natural gas/hydrogen blends at elevated pressures,”

Combust. Flame

161

,

1432

–

1443

(

2014

).

https://doi.org/10.1016/j.combustflame.2013.12.005

Google Scholar

Crossref

27.

C.

Ji

,

D.

Wang

,

J.

Yang

, and

S.

Wang

, “A comprehensive study of light hydrocarbon mechanisms performance in predicting methane/hydrogen/air laminar burning velocities,”

Int. J. Hydrogen Energy

42

,

17260

–

17274

(

2017

).

https://doi.org/10.1016/j.ijhydene.2017.05.203

Google Scholar

Crossref

28.

H.

Kutkan

,

A.

Amato

,

G.

Campa

,

G.

Ghirardo

,

L.

Tay Wo Chong

, and

E.

Asoy

, “Modelling of turbulent premixed CH₄/H₂/air flames including the influence of stretch and heat losses,”

J. Eng. Gas Turbines Power

144

,

011020

(

2021

).

https://doi.org/10.1115/1.4051989

Google Scholar

Crossref

29.

S.

Gordon

and

B.

McBride

, “Computer program for calculation of complex chemical equilibrium compositions and applications,” NASA Reference Publication 1311,

1994

, https://shepherd.caltech.edu/EDL/PublicResources/sdt/refs/NASA-RP-1311-1.pdf.

30.

J. J.

Bertin

and

R. M.

Cummings

,

Aerodynamics for Engineers

, 6th ed. (

Cambridge University Press

,

2021

).

Google Scholar

Crossref

31.

B. E.

Launder

and

D. B.

Spalding

, “The numerical computation of turbulent flows,”

Comput. Methods Appl. Mech. Eng.

3

,

269

–

289

(

1974

).

https://doi.org/10.1016/0045-7825(74)90029-2

Google Scholar

Crossref

32.

K.

Hanjalic

, “Two-dimensional asymmetric turbulent flow in ducts,” Ph.D. thesis,

Faculty of Engineering, University of London

,

1970

.

33.

C. M.

Rhie

and

W. L.

Chow

, “Numerical study of the turbulent flow past an airfoil with trailing edge separation,”

AIAA J.

21

,

1525

–

1532

(

1983

).

https://doi.org/10.2514/3.8284

Google Scholar

Crossref

34.

S.

Kim

,

B.

Makarov

, and

D.

Caraeni

, “A multidimensional linear reconstruction scheme for arbitrary unstructured grid,” in

AIAA 16th Computational Fluid Dynamics Conference

,

2003

.

35.

R. I.

Issa

, “Solution of the implicitly discretized fluid flow equations by operator-splitting,”

J. Comput. Phys.

62

,

40

–

65

(

1985

).

https://doi.org/10.1016/0021-9991(86)90099-9

Google Scholar

Crossref

36.

Ansys Fluent Academic Research, Release 19, Help System

(

Ansys Fluent Theory Guide, Ansys, Inc

).

37.

H.

Kutkan

and

J.

Guerrero

, “Turbulent premixed flame modeling using the algebraic flame surface wrinkling model: A comparative study between OpenFOAM and ansys fluent,”

Fluids

6

,

462

(

2021

).

https://doi.org/10.3390/fluids6120462

Google Scholar

Crossref

38.

A. N.

Lipatnikov

and

J.

Chomiak

, “Turbulent flame speed and thickness: Phenomenology, evaluation, and application in multi-dimensional simulations,”

Prog. Energy Combust. Sci.

28

,

1

–

74

(

2002

).

https://doi.org/10.1016/s0360-1285(01)00007-7

Google Scholar

Crossref

39.

J.

Guerrero

,

L.

Mantelli

, and

S. B.

Naqvi

, “Cloud-based cad parametrization for design space exploration and design optimization in numerical simulations,”

Fluids

5

,

36

(

2020

).

https://doi.org/10.3390/fluids5010036

Google Scholar

Crossref

2022

Author(s)

Pressurized turbulent premixed CH₄/H₂/air flame validation using OpenFOAM

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Premixed combustion model

B. Energy equation treatment

C. Thermophysical and transport properties

III. DESCRIPTION OF THE EXPERIMENTAL CONFIGURATION AND NUMERICAL BACKGROUND

A. Experimental setup and data

B. Computational domain and mesh

C. Boundary and initial conditions

D. Numerical setup and discretization schemes

IV. RESULTS AND DISCUSSION

V. CONCLUSIONS AND PERSPECTIVES

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

NOMENCLATURE

Notation

Sub-indices

Acronyms

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Pressurized turbulent premixed CH4/H2/air flame validation using OpenFOAM

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Premixed combustion model

B. Energy equation treatment

C. Thermophysical and transport properties

III. DESCRIPTION OF THE EXPERIMENTAL CONFIGURATION AND NUMERICAL BACKGROUND

A. Experimental setup and data

B. Computational domain and mesh

C. Boundary and initial conditions

D. Numerical setup and discretization schemes

IV. RESULTS AND DISCUSSION

V. CONCLUSIONS AND PERSPECTIVES

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

NOMENCLATURE

Notation

Sub-indices

Acronyms

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Pressurized turbulent premixed CH₄/H₂/air flame validation using OpenFOAM