2.1. 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,
20,
21],
In Equation (
1), the transported quantity
c refers to the normalized mass fraction of the products [
2], 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 Equation (
1) can be found in references [
7,
20,
21,
22,
23].
The solution of Equation (
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.
To conduct the numerical simulations, we used the CFD solvers OpenFOAM and Ansys Fluent. It is worth mentioning that Ansys Fluent [
17] solves Equation (
1). The CFD solver OpenFOAM [
15,
16], instead of solving for the progress variable
c, solves for the regress variable
b [
24], where,
Henceforth, Equation (
1) can be written in terms of the regress variable
b as follows,
In Equations (
1) and (
3), the overbar
and the tilde
refer to the Reynolds and Favre averaging, respectively. In these equations, the molecular thermal diffusivity
is usually much smaller than the turbulent diffusivity
and can simply be neglected in RANS/URANS [
20]. Recall that the molecular thermal diffusivity
is computed using the following relation,
In Equation (
1) (the equation used in Ansys Fluent), the progress variable
c is bounded between zero and one, where
represents the reactants (unburned state) and
represents the products (burned state). Any intermediate values between
denote the gas mixture with temperatures and composition between those of the reactants and products.
In the equation governing the regress variable
b (Equation (
3)),
b is also bounded between zero and one. However, in this case,
represents the products and
represents the reactants. Equation (
3) is the equation used in OpenFOAM and, in particular, in the solver XiFoam. OpenFOAM’s XiFoam solver is a transient solver for compressible premixed and partially premixed combustion with turbulence modeling. XiFoam is the specific OpenFOAM solver that was used in this study.
In Equations (
1) and (
3), evaluating the mean reaction rate source term
is the central 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 used an algebraic model, in particular the AFSW model [
8,
11,
13,
14]. By using this model, the reaction rate source term for the progress variable and the regress variable can be expressed as follows,
where the magnitude of the gradient term in Equation (
5) can be evaluated as follows,
The turbulent flame speed
appearing in Equation (
5) is then modeled using the AFSW model [
8,
11,
13,
14], as follows,
In this equation,
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 Equation (
7), the following closure relationships are used,
The effective Lewis number
appearing in Equation (
8) models the preferential diffusion effect of the species
in the unburnt mixture on the turbulent flame speed (Equation (
7)) [
11]. In Equation (
8), the mixture molecular thermal diffusivity
is computed using Equation (
4), where
is computed at the temperature at which the maximum heat is released (
in
Table 1).
In this study, we assumed the turbulent Lewis number
to be equal to one. This implies that the turbulent Schmidt number
and the turbulent Prandtl number
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 [
20,
22,
23]. For completeness, recall that the
, the
, and the
numbers are defined as follows,
by combining these expressions, we can derive a relation for the turbulent diffusivity
(Equation (
13)), where the turbulent Schmidt number
is assumed to be constant.
It is worth noting that the evaluation of
is not simple. The closure given in Equation (
14) is often used together with the
turbulence model, where
is computed from the transported turbulent quantities
k and
. However, while it is commonly accepted that the value of the coefficient
is equal to 0.09, different authors have reported different values of
. In this study, we assumed
[
22,
25].
Finally, by substituting Equation (
13) into Equations (
1) and (
3) and by neglecting the molecular thermal diffusivity
, we obtain the following solvable equations for the progress variable
c and the regress variable
b,
Equation (
15) or Equation (
16), together with the compressible RANS equations, the
closure turbulence model, and proper boundary conditions and initial conditions, are solved using OpenFOAM’s XiFoam solver and Ansys Fluent.
2.2. Energy Equation Treatment
In OpenFOAM’s XiFoam solver, the energy equation is treated either in terms of the absolute enthalpy
formulation (the sensible enthalpy plus the chemical enthalpy) or in terms of the internal energy
formulation [
26]. In this study, we used the absolute enthalpy formulation in OpenFOAM, namely Equation (
17). In this equation,
K is the specific kinetic energy and
P is the thermodynamic pressure (absolute pressure). Since
already contains the chemical enthalpy, there is no chemical heating source in Equation (
17).
In contrast, Ansys Fluent uses the sensible enthalpy
formulation when dealing with premixed and partially premixed combustion [
17]. In the formulation used hereafter (Equation (
18)), the thermodynamic pressure, the specific kinetic energy, and the viscous dissipation (or viscous heating) terms are not included (these terms are often negligible in incompressible flows). Instead, it includes the chemical enthalpy
as a source term on the RHS of Equation (
18), which needs to account for the conservation of energy in the sensible enthalpy formulation [
17,
26].
In Ansys Fluent, the partially premixed setup was used, but with the PDF equations disabled, which is identical to perfectly premixed combustion. Therefore, all compressibility effects were taken into account using the PDF look-up tables computed by Ansys Fluent. Finally, in Equation (
18), the chemical enthalpy
is defined as,
where
is the reaction rate source term and is defined in Equation (
5),
is the lower heating value of combustion, and
is the fuel mixture mass fraction.
2.3. 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 [
27]. 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 [
28]. GRI-Mech is a database of chemical reactions and associated rate constant expressions capable of the best representation of natural gas flames and ignition. It is a compilation of 325 chemical reactions and related rate coefficient expressions and thermochemical parameters for the 53 species involved in them. Previous studies have shown that the calculated
for
/
/air mixtures using this mechanism is in good agreement with experimental data [
29,
30].
In
Table 1, we show the calculated properties. In this table, the variables
, and
are calculated at the maximum heat release rate temperature (
). As can be seen in this table, the addition of
to the mixture increases the laminar flame speed
. This is attributed to higher
reactivity with respect to
[
31].
In Ansys Fluent, the temperature-dependent mixture properties, namely, specific heat
, thermal conductivity
, and dynamic viscosity
, are calculated using polynomial models. To compute the coefficients of the polynomials, the physical properties
,
, and
are calculated across the flame at certain temperature values (equally spaced between unburnt and burnt temperature) using the FreeFlame configuration in Cantera. Then, by using polynomial fitting, the coefficients of the polynomial model are obtained. In addition, the mixture law model is used to model the composition-dependent mixture specific heat [
17].
On the other hand, the solver XiFoam in OpenFOAM requires the definition of the JANAF table coefficients for
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
,
The JANAF table coefficients were calculated using Equations (
21)–(
23) [
32]. In these equations, the superscript
o refers to standard state pressure conditions or
. In Equations (
21)–(
23),
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
K.
The coefficients
, and
in Equation (
21) were obtained by curve fitting at different temperature values (
T). Once
have been determined, they are substituted in Equations (
22) and (
23) at
. At this point, the coefficient
appearing in Equation (
22) (enthalpy jump) and the coefficient
appearing in Equation (
23) (entropy jump) can be computed. These coefficients were calculated for low- (
) and high- (
) temperature ranges for both unburnt and burnt mixtures.
The dynamic viscosity of the mixture was computed using the Sutherland law [
33], which is expressed as follows,
where
is the temperature-dependent mixture dynamic viscosity,
and
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
[
9,
10,
11].
In XiFoam, the coefficients
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 [
34].
The final two parameters, stoichiometric air-to-fuel mass ratio
and equivalence ratio
, can be calculated using Equations (
25) and (
26), respectively. In Equation (
26), the sub-index
refers to the fuel mixture. In this study, we simulated two cases corresponding to two different fuel mixtures. The species volumetric ratio of each fuel mixture is given in
Table 2.