Laminar premixed hydrogen flame transfer functions: Comparison between flow solvers and transport models

Abstract

The computation of flame transfer functions (FTFs), which characterize the flame response to acoustic fluctuations, is of central importance to the assessment of thermoacoustic stability. Hydrogen is an important carbon-free fuel to ensure sustainable power production in gas turbines. Since hydrogen possesses different reactivity and thermo-diffusive characteristics in comparison to natural gas, the nominal flame structure is different from natural gas, and therefore it also exhibits different dynamics. In the context of computing hydrogen flames, the non-unity Lewis numbers call for transport models that consider the enhanced molecular diffusion of hydrogen. This article answers an open question of how present-day computational fluid dynamics (CFD) solvers and transport models perform in the calculation of flame dynamics. Computations of FTFs of laminar premixed hydrogen flames are performed using OpenFOAM, ANSYS Fluent, S3D DNS code, and the AVBP solver. These solvers encompass the spectrum of numerical schemes and transport models used in computational combustion. Transport models taking into account the enhanced molecular diffusion of hydrogen are used in each solver. The results show that despite using similar transport models and identical chemical mechanisms, quantitative differences in the laminar flame speed, mean flame shapes, and the FTFs are seen between the solvers. However, the qualitative features of the FTF remain the same and the quantitative differences in terms of the root-mean-square error metric are within $15 %$ in gain and $6 %$ in phase. However, with transport models, larger quantitative differences in the phase of the FTF are observed.

Keywords

Hydrogen flames flame transfer functions thermo-diffusive effects computational fluid dynamics laminar premixed flame dynamics transport models CFD solvers

Introduction

Combustion instability, an important practical problem in gas turbine and rocket engines, is caused due to the coupled interaction between the flame, flow, and the resonant acoustic modes of a combustion chamber.¹ Flames respond to acoustic waves by modulating their heat release rate, position, and shape² in time. Unsteady flames, in turn, act as a source of sound creating acoustic waves.³ This coupled feedback loop involves an interplay between the intricate combustion, hydrodynamic, and acoustic phenomena and is generally non-trivial to predict. Insight into the thermoacoustic behavior of a given flame–combustion chamber system can be obtained by computing the open-loop response of the flame to acoustic perturbations. For this purpose, the flame is subjected to acoustic waves either experimentally using a loudspeaker/siren^4,5 or numerically by imposing acoustic waves at the inlet/exit of the computational domain.^6,7 The resulting global heat release response of the flame is measured/calculated to derive a flame transfer function (FTF) which can be mathematically written as

F (ω) = \frac{\hat{{\dot{Q}}^{'}} / {\dot{Q}}_{0}}{\hat{u^{'}} / u_{0}}

(1)

where

\dot{Q}

represents the integrated heat release rate,

u

is the flow velocity at a reference location,

()_{0}

represents a time-averaged quantity,

()^{'}

represents a fluctuating quantity, and

(\hat{})

denotes the Fourier transform of the associated quantities.

A FTF has a characteristic frequency $(ω)$ -dependent behavior and sheds crucial insight into the combustion instability phenomenon since it gives information on the time lag(s) characteristic to the flame and the frequency range over which the flame exhibits a significant response to acoustic perturbations.⁸ Thus, it can be used both for prediction and control of combustion instabilities. Indeed, an FTF can be used, in combination with an acoustic network model of the combustion chamber, to predict the linear thermoacoustic frequencies and growth rates.^9,10 More recently, it has been shown that self-excited thermoacoustic instabilities can be efficiently suppressed by targeted modifications to FTFs.^11,12

Hydrogen is regarded as a crucial carbon-free fuel that will play a key role in enabling power generation in gas turbine engines with low greenhouse gas emissions. However, it possesses different physical and chemical properties to that of natural gas, and therefore, must be re-assessed from the standpoint of combustion instability. In other words, the combustion (in)stability characteristics of existing burners must be re-evaluated when using hydrogen since it results in flames with nominally different shapes and heat release distributions. Two unique properties of hydrogen are particularly important in this regard. First, hydrogen exhibits an increased reactivity (increased laminar flame speed) in comparison to natural gas. For a fixed flow velocity and flame temperature, hydrogen flames are much shorter than natural gas flames, and this reflects in FTFs as a reduced phase slope and increased cutoff frequency.^13–15 Additionally, hydrogen also exhibits enhanced diffusivities of mass in relation to heat, and therefore has a strong stretch/curvature sensitivity of the local consumption speed.^16,17 For a fixed unstretched laminar flame speed and flow velocity, hydrogen flames can exhibit shorter or longer flame lengths in comparison to natural gas depending on the flame front curvature and the Markstein length.¹⁸ Therefore, it is natural to expect, even under conditions where the flow velocity and the unstretched laminar flame speed are maintained constant, differences in the FTF of hydrogen and natural gas flames. This warrants the need for tools to computationally determine the FTFs of hydrogen flames to get insight into the combustion instability characteristics of these carbon-free fuels.

Since the enhanced mass diffusivity of hydrogen is a key characteristic of the fuel, computational approaches to extract the FTF of hydrogen flames must capture the non-unity Lewis numbers of hydrogen. More specifically, transport models that set all species Lewis numbers to unity, while working well for natural gas flames, are not expected to capture the physics of hydrogen flames. To consistently capture the enhanced mass diffusivity of hydrogen in a flame computation, three approaches are possible. First, a “multi-component” transport model where the species diffusional velocities are computed by solving a system of linear equations involving the binary diffusion coefficients $(D_{i j})$ of each species pair, mass fractions of the species $(X_{i})$ , and pressure. This entails inverting a system of equations of size $N_{s}^{2}$ , where $N_{s}$ is the number of species, along each direction, at each spatial point and at each time instant, therefore leading to significant computational costs.¹⁹ The second approach, commonly christened as “mixture-averaged” transport model, employs the Hirschfelder and Curtiss approximation where an effective diffusion coefficient $(D_{i})$ of any species into the mixture is determined from the binary diffusion coefficients and the species mass fractions, and does not involve the solution of system of equations.²⁰ Such an approach accurately captures the enhanced mass diffusion of hydrogen, where the non-unity Lewis numbers arise from the mixture-averaged mass diffusion coefficient rather than being prescribed as constants. Finally, in the third approach, it is also possible to prescribe constant non-unity values to the individual species Lewis numbers, characteristic to the state of the unburnt mixture, for example. In this paper, the latter two approaches will be used to consistently capture thermo-diffusive effects of hydrogen flames.

The importance of correctly capturing the enhanced mass diffusivity of hydrogen in computations of the static and dynamic flame characteristics prompts us to the following key question. How do different flow solvers and transport models compare in computing key combustion characteristics (static and dynamic) of hydrogen flames? Furthermore, we also intend to get quantitative insight into the relative differences in FTFs obtained when using different solvers and transport models. Answering the preceding research questions will be the main aim of this paper. To accomplish this, we perform computations of one-dimensional (1D) laminar premixed flames and two-dimensional (2D) slit-stabilized laminar premixed flame dynamics of hydrogen flames using different computational fluid dynamics (CFD) solvers and different transport models. While we acknowledge that the present study only focuses on laminar premixed flames, whereas real-world gas turbine combustors mostly feature turbulent flames, we emphasize that the current work represents a first step to understand how different flow solvers and transport models behave in computing hydrogen flames, and can be easily extended to turbulent premixed flames in the future. The computations in this work are performed using the open-source flow solvers OpenFOAM,²¹ ANSYS Fluent v2024 R1, S3D DNS code of Sandia National Laboratories,²² and the AVBP solver from CERFACS.²³ The next section gives details of the geometrical configuration and the computational methods employed.

Geometrical configuration and computational methodology

Both 2D and 1D computations of laminar premixed flames are performed in this work. While the 2D computations are used to characterize the flame shapes, global consumption speeds, and the flame dynamics, the 1D computations (see Figure 1) are intended to compute the unstretched laminar flame speeds. The 1D computations feature a premixed propagating flame and are simulated computationally using a three-dimensional domain with an axial extent of $20$ mm and transverse dimensions of $2$ mm. The domain is discretized with $15,000$ points in the axial direction and with one grid point in the transverse $(y, z)$ directions. This grid results in $400$ grid points across the flame thickness and is more than sufficient to capture the rapid variations in flow quantities. Symmetry plane boundary conditions are imposed on the top and bottom and the front and back boundaries. At the inlet, a preheated premixed reacting mixture is supplied at a flow velocity $(u_{0})$ identical to the unstretched laminar flame speed computed from the chemical kinetics code Cantera.²⁴ The computation is initialized with an initial flame solution positioned at the midpoint of the domain obtained from Cantera. This solution is then progressed until a physical time of 200 $δ_{l}^{0} / S_{l}^{0}$ , where $δ_{l}^{0}$ is the flame thickness and $S_{l}^{0}$ is the unstretched laminar flame speed. The actual value of the unstretched laminar flame speed predicted by a given CFD solver can then be computed from

S_{l}^{0} = u_{0} - u_{f}

(2)

where

u_{f}

is the steady velocity of the flame in the laboratory reference frame computed as

d X_{f} / d t

, where

X_{f}

is the flame position given as the spatial location of the maximum heat release rate in the domain.

Figure 1.

One-dimensional geometrical configuration used for the unstretched laminar flame speed computations.

For the 2D computations, a slit-flame combustor (see Figure 2), originally introduced by Meindl et al.,²⁵ is employed. This configuration consists of a straight channel with a length of $20$ mm and a width $(b)$ of $2$ mm. Flame stabilization is achieved by utilizing a composite wall. The first $7$ mm of this wall is maintained to be isothermal at the reactants’ temperature, while the rest $13$ mm is maintained to be adiabatic, thus, anchoring the flame in the transition point between the isothermal and adiabatic boundaries. Preheated reactants at a temperature of $750$ K and a mean velocity $u_{0}$ are introduced at the inlet of the domain. The transverse variation of the inlet velocity is given by the classical plane channel velocity profile

u (y) = \frac{3}{2} u_{0} (1 - \frac{y^{2}}{b^{2}})

(3)

For studying the flame dynamics, acoustic velocity fluctuations with an amplitude of

4 %

of the mean value are imposed. To reduce computational effort, only half of the channel is computed by imposing a symmetry boundary condition. A structured grid of cell size 20

μ m

is used for discretizing the domain, which ensures a minimum of

24

grid points in the flame thickness across all cases. Next, we present the CFD solvers employed in this work.

Figure 2.

Slit-flame configuration featuring a premixed propagation-stabilized flame. This configuration is used for characterizing the stretch-affected laminar flame speed and the flame dynamics.

OpenFOAM computational methodology

The reactingFOAM compressible solver from the OpenFOAM v10^21,26 toolbox is used for computing the laminar premixed flame and its dynamics. Viscosity and thermal conductivity of each species are specified as polynomial functions of temperature. From these individual species viscosities, thermal conductivities, and the species mass fractions, the mixture-averaged viscosity and thermal conductivity are then calculated using standard mixture-averaging rules.²⁰ The mixture-specific heats and enthalpies are also computed from the corresponding species values, which are given in the form of NASA polynomials²⁷ as a function of temperature. The reacting mixture in OpenFOAM is specified as a “coefficientWilkeMultiComponentMixture” which implicitly implements the mixture-averaging rules.

A vital aspect in the computation of hydrogen flames is the representation of the enhanced mass diffusivities of hydrogen. This is done using a mixture-averaged (MA) transport model approach, which is incorporated using the “FickianFourier” utility in OpenFOAM. In this framework, the species diffusional velocities are calculated using the Fickian law,²⁸ which involves the spatial gradients of species mass (or mole) fractions and the mixture-averaged diffusion coefficients $(D_{i}^{mix})$ , where $i$ represents the species index. The $D_{i}^{mix}$ is, in turn, calculated from the binary diffusion coefficients $D_{i j}$ of each species pair.²² The binary diffusion coefficients of each species pair are specified as tabulated functions of pressure and temperature precomputed using Cantera. The mathematical expressions used for the determination of transport properties are listed in Table 1. Furthermore, in the calculation of the species diffusion velocities, thermal, and barodiffusion effects are neglected. In addition to the mixture-averaged approach, computations are also performed using an unity-Lewis (Le $_{1}$ ) transport model wherein all species Lewis numbers are set to unity, that is, $D_{i}^{mix}$ is calculated as $k / ρ C_{p}$ .

Table 1.

Method used for computation of transport properties for the MA transport model in the various solvers.

AVBP
Species properties	$μ_{i} = μ_{i} (M W_{i}, T, σ, k_{B}, Ω_{i i}^{(2, 2)})$ ,
	$k_{i} = k_{i} (μ_{i}, C_{p, rot}, C_{p, trans}, C_{p, vib}, M W_{i}, R, D_{i i})$ ,
	$D_{i j} = D_{i j} (T, σ, k_{B}, M W_{i}, M W_{j}, Ω_{i j}^{(1, 1)})$
Mixture properties	$μ = \sum_{i} \frac{X_{i} μ_{i}}{X_{i} + \sum_{j \neq i} X_{j} Φ_{i j}}$ ,
	$Φ_{i j} = [1 + (μ_{i} / μ_{j})^{1 / 2} (M W_{j} / M W_{i})^{1 / 4}]^{2} [8 (1 + M W_{i} / M W_{j})]^{- 1 / 2}$
	$k = \frac{1}{2} (\sum_{i} X_{i} k_{i} + \frac{1}{\sum_{i} X_{i} / k_{i}})$
	$D_{i}^{mix} = (1 - Y_{i}) / (\sum_{j \neq i} X_{j} / D_{i j})$
S3D
Species properties	Same as AVBP
Mixture properties	Same as AVBP with additional Soret effect included for calculation of species diffusion
	$D_{i}^{mix} = (1 - X_{i}) / (\sum_{j \neq i} X_{j} / D_{i j})$
Fluent
Species properties	Same as AVBP
Mixture properties	Same as AVBP except $k$ and $D_{i}^{mix}$
	$k = \sum_{i} (X_{i} k_{i} / \sum_{j} X_{j} Φ_{i j})$ , $D_{i}^{mix} = (1 - X_{i}) / (\sum_{j \neq i} X_{j} / D_{i j})$
OpenFOAM
Species properties	$μ_{i} = \sum_{n = 1}^{8} a_{i, n} (\ln T)^{n - 1}$ ,
	$k_{i} = \sum_{n = 1}^{8} b_{i, n} (\ln T)^{n - 1}$ ,
	$D_{i j} :$ Function of $T, p$ in tabular form
Mixture properties	Same as Fluent

$M W_{i}$ denotes the molecular weight of species $i$ , $σ$ denotes the Lennard–Jones length parameter, $ϵ / k_{B}$ represents the Lennard–Jones energy parameter where $k_{B}$ is the Boltzmann constant, $R$ is the universal gas constant, $Ω_{i j}^{(l, s)}$ denotes the collision integrals, $X_{i}$ is the mole fraction of species $i$ , and all the summations over indices $i$ and $j$ are performed over the number of species.^{20,22,49,51,52} MA: mixture-averaged.

With regard to the boundary conditions imposed for the computations, at the inlet, the temperature and species mass fractions are specified as fixed values using the “fixedValue” type condition, and a “zeroGradient” condition is imposed for pressure. The parabolic velocity profile at the inlet is specified using the “codedFixedValue” condition. A “waveTransmissive” boundary condition is employed for the pressure at the outlet boundary, which imposes a partially non-reflective outflow in accordance with the Navier–Stokes characteristic boundary condition (NSCBC) methodology with a specified relaxation coefficient.^29,30 This relaxation coefficient is specified as a length parameter (lInf) which is kept constant at 1 m for all computations reported in this paper. A “pressureInletOutletVelocity” condition is imposed for velocity at the outlet. No-slip isothermal/adiabatic condition is imposed at the walls, while the “symmetry” condition is employed for the bottom boundary. Except for the inlet, all the other boundaries employ a “zeroGradient” condition for species mass fractions. A PIMPLE algorithm with 15 outer iterations and 4 inner iterations is used for the computations. The numerical schemes employed for the temporal and spatial discretization are listed in Table 2. The maximum convective Courant–Friedrichs–Lewy (CFL) number is maintained as 0.4 for all the computations.

Table 2.

Numerical schemes used for time and space discretization in OpenFOAM.

Discretization variable/operator	Scheme
Time	backward Euler
Gradient	Gauss linear
Divergence	Gauss limitedCubicV
Laplacian	Gauss linear orthogonal
Interpolation (default)	cubic
Interpolation (species)	limitedCubic01

AVBP computational methodology

AVBP is a massively parallel, fully compressible, reactive turbulent flow solver designed for both Large eddy simulations (LES) and direct numerical simulation (DNS) on GPU and CPU architectures. It solves the multi-species reactive Navier–Stokes equations in conservative form on both structured and unstructured grids, using a cell-vertex finite-volume method.²³ The viscous–convective part of the flux tensor is discretized with a revisited version of the multi-parameter family of a two-step Taylor–Galerkin (TTG) scheme,³¹ denoted TTGC,³² which is the third-order accurate both in space and time. Since AVBP is extensively employed for LES of industrial configurations, the scheme is designed to handle hybrid and irregular meshes without exhibiting over-dissipative behavior at high frequencies, which would otherwise damp important turbulent structures at practical mesh resolutions. At the same time, it achieves a near fourth-order accuracy in some norms for relatively regular meshes,³³ as is the case in the present work. On the other hand, the inviscid-diffusive part of the flux tensor is discretized with a centered second-order scheme. The acoustic CFL number, set to 0.7 in this work, is chosen by the user and kept below the critical CFL threshold of the convection scheme. This corresponds to a time-dependent time step of approximately 10 ns across all simulated cases.

The thermodynamic properties of the mixture are derived from the NIST–JANAF tables²⁷ or from the CHEMKIN database.³⁴ As far as the transport model is concerned, to adequately capture the enhanced mass diffusivity of hydrogen, a mixture-averaged approach is employed. This method provides a more accurate determination of transport coefficients by accounting for the composition-dependent contributions of each species in the mixture. The resulting transport properties correspond to those of a hypothetical mixture with averaged characteristics. According to the Chapman–Enskog theory,^35,36 the evaluation of the mixture viscosity, thermal conductivity, and diffusivities requires first the computation of the viscosity and thermal conductivity of each species, as well as the binary diffusion coefficients for each species pair. These quantities depend on several parameters (temperature, Boltzmann constant, molecular mass, collision diameter, and collision integrals) and are, therefore, fitted against temperature. The mixture viscosity is then obtained using Wilke’s law,³⁷ while individual species viscosities are evaluated through kinetic theory³⁸; species’ thermal conductivities are determined using the approach proposed by Eucken,³⁶ and the thermal conductivity of the mixture is computed following the method of Mathur et al.³⁹ Finally, species diffusivities are evaluated using the mixing law of Hirschfelder and Curtiss.⁴⁰ Currently, Soret effect, that is, the influence of thermal diffusivity on species diffusion velocities, is not accounted for in AVBP. This procedure is repeated at every iteration and for every computational node. The equations mentioned above can be found in Table 1.

Since AVBP is a fully compressible code, the implementation of boundary conditions (BC) is of paramount importance. Historically, the formalism introduced by Poinsot and Lele⁴¹ has been used in AVBP in which all variations of physical quantities at the boundaries are decomposed into incoming and outgoing characteristic waves. The key concept underlying this approach is the following: the outgoing characteristic waves leaving the domain are accurately computed by the numerical scheme and must therefore be left unchanged, whereas the incoming characteristic waves entering the domain cannot be computed by the numerical scheme and are replaced by user-defined values determined by the physics of the boundary condition. For an inlet-type boundary condition, four waves must be specified, determined by velocity, temperature, and species mass fractions, whereas for an outlet, only one wave, determined by pressure, needs to be prescribed. In this study, partially non-reflective boundary conditions with specified relaxation coefficients ( $K_{i}$ corresponding to each variable $i$ ) are imposed. The magnitude of the relaxation coefficients denotes the degree of non-reflectiveness of the boundary, with $K_{i} = 0$ ensuring a fully non-reflective BC. In this study, $K_{u} = 1000$ s⁻¹, $K_{T} = 100$ s⁻¹, and $K_{Y} = 100$ s⁻¹ were employed at the inlet boundary, while $K_{p} = 50$ s⁻¹ was used for the outflow boundary. The chosen value of $K_{p}$ is sufficiently low to be nearly numerically anechoic, yet high enough to prevent pressure drift.

S3D computational methodology

S3D is a compressible, reactive, DNS flow solver, developed at the Sandia National Laboratories,²² specifically designed to yield highly accurate solutions of reacting flows in simple geometrical configurations. The Navier–Stokes equations for a compressible fluid in the conservative form are solved on structured, Cartesian meshes. For spatial discretization, an eighth-order finite difference scheme is used in combination with a tenth-order, explicit, spatial filter to remove spurious numerically induced high-wavenumber errors.⁴² A fourth-order, six-stage, explicit Runge–Kutta scheme is used for time integration with the time step set to $5$ ns.⁴³ The thermodynamic properties of the mixture are determined from those of the individual species, which are calculated from polynomial functions of temperature using the CHEMKIN package.³⁴ The mixture-averaged transport model described in detail in Kee et al.²⁰ and in Chen et al.²² is used for computing the transport properties and the species diffusional fluxes. For the calculation of the species diffusional velocities, thermal (Soret) diffusion effects are taken into account, while barodiffusion is neglected. The methodology and the equations employed for this transport model are very similar to those used in AVBP and are listed in Table 1.

In addition to the mixture-averaged transport model, S3D also features a “LEWIS” transport model where individual species Lewis numbers can be specified. In this method, the thermal conductivity is calculated using the relation proposed by Smooke et al.⁴⁴

λ = A C_{p} {(\frac{T}{T_{0}})}^{0.7}

(4)

where

T_{0} = 298

K and

A = 2.58 \times 10^{- 4}

g (cm s)⁻¹, and

C_{p}

is the constant pressure specific heat. The mixture viscosity is calculated using a specified Prandtl number and the following relation

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

(5)

where a Prandtl number of

0.60

is used for all fuel computations reported in this paper. The species diffusion coefficients are then calculated using

D_{i} = \frac{λ}{ρ C_{p} L e_{i}}

(6)

where

L e_{i}

represents the Lewis numbers of the individual species, computed for each operating condition using the unburnt state of the premixed mixture.

Two crucial deficiencies associated with the “LEWIS” transport model must be mentioned at this point. First, the power-law expression of equation (4) has been verified only for premixed natural gas flames under stoichiometric and lean conditions.⁴⁴ For hydrogen flames, a revision to the form and coefficients of equation (4) may be required. However, this is not attempted in this work, since the aim is to compare the existing state of the “LEWIS” model in S3D with other detailed approaches. Furthermore, a constant value of Prandtl number of $0.60$ is assumed for all fuel computations in this work. Cantera calculations of Prandtl numbers of hydrogen flames considered in this work revealed that the values of this constant varied between 0.56 and 0.61 in the unburnt side. Additionally, in the burnt side of the flame, the Prandtl number increases to a value of around 0.70–0.71. Therefore, we chose a constant value of $0.60$ . The incorporation of a better model for the thermal conductivity and a more accurate specification of temperature-dependent Prandtl numbers will be implemented in the flow solver in future works.

With regard to the boundary conditions, open (inlets, outlets) and closed (walls, symmetry plane) boundaries are treated differently in S3D. For open boundaries, the NSCBC methodology^41,45 as previously explained is used to impose boundary conditions. A “non-reflecting” outflow with a relaxation coefficient of $K_{p} = 930 s^{- 1}$ is used, while a hard-inflow condition, where all flow quantities are specified in a time-dependent manner, as described in Sutherland and Kennedy⁴⁵ is employed. In contrast to open boundaries, for closed boundaries, most of the values of the solution vector $[ρ u, ρ v, ρ w, ρ, ρ e_{t}, ρ Y_{i}]^{T}$ , which consists of the velocities $(u, v, w)$ , density $(ρ)$ , total internal energy $(e_{t})$ , and species mass fractions $(Y_{i})$ , at the boundary points are directly specified. This is done by overwriting the values of the conservative variables obtained from the time-marching procedure when needed. For instance, in the case of an isothermal wall, the value of $ρ e_{t}$ is overwritten by calculating it using the wall temperature, while for an adiabatic wall $e_{t}$ is obtained by naturally solving the energy equation by setting the normal heat flux to zero.

In the following, we detail the implementation of the symmetry plane and the composite wall in S3D, since both of these were newly added to the solver framework. The symmetry plane was modeled as an adiabatic slip wall wherein the wall-normal velocity is zero. Therefore

\begin{aligned} u_{n} = 0 \\ RHS ({Momentum}_{n}) = 0 \end{aligned}

(7)

where

u_{n}

is the normal velocity at the boundary point and

RHS ({Momentum}_{n})

represents the right-hand side of the normal momentum equation. The tangential velocities are computed by solving the tangential momentum equations. In accordance with Poinsot and Veynante,¹⁹ the tangential shear stresses in the wall-normal plane are set to zero, that is, for a

y

-normal symmetry plane

\begin{aligned} τ_{y x} = & τ_{x y} = 0 \\ τ_{y z} = & τ_{z y} = 0 \end{aligned}

(8)

Furthermore, due to the non-porous nature, the species diffusional fluxes normal to the wall are zero, resulting in

\frac{\partial Y_{i}}{\partial n} = 0

(9)

which can be inverted to determine the species mass fraction at the walls. The total internal energy, and therefore the temperature at the symmetry plane, is computed by solving the energy equation with the boundary-normal heat flux set to zero

q_{n} = 0

(10)

The final variable needed to be specified is the density, which is computed using the equation of state

\begin{aligned} ρ_{w} = & p_{w} / R T_{w} \\ RHS (Continuity) = & 0 \end{aligned}

(11)

where

p_{w}

is the pressure at the symmetry plane, and

T_{w}

is the temperature at the symmetry plane. Following Gruber,⁴⁶ the pressure at the symmetry plane is deduced by consideration of the wall-normal momentum equation. Again, considering a

y

-normal symmetry plane, the momentum equation in the wall-normal direction can be written as

\begin{aligned} ρ \frac{\partial v}{\partial t} + ρ u \frac{\partial v}{\partial x} + ρ v \frac{\partial v}{\partial y} + ρ w \frac{\partial v}{\partial z} \\ = \frac{\partial}{\partial x} (τ_{x y}) + \frac{\partial}{\partial x} (τ_{y y}) + \frac{\partial}{\partial x} (τ_{z y}) - \frac{\partial p}{\partial y} \end{aligned}

(12)

For a symmetry plane, the fact that the normal velocity is zero at all points results in the left-hand side of the above equation becoming zero. Therefore, at the symmetry plane, the following condition is satisfied

\frac{\partial p}{\partial y} = \frac{\partial}{\partial x} (τ_{x y}) + \frac{\partial}{\partial x} (τ_{y y}) + \frac{\partial}{\partial x} (τ_{z y})

(13)

which can be solved to determine the pressure.

The specification of the solution vector for a composite wall follows similarly to the implementation for a symmetry plane. In the following, we only highlight the differences between the implementation of a composite wall in relation to the symmetry plane. For velocities, in addition to equation (7), the tangential velocity components and the RHS of the momentum equations in the tangential direction are set to zero. The total internal energy for the adiabatic portion of the wall is computed by solving the energy equation with the normal heat flux set to zero, while for the isothermal portion of the wall it is determined using the wall temperature $(T_{w})$ and species mass fractions at the wall $(Y_{i, w})$

e_{t, wall} = e_{t} (T_{w}, Y_{i, w})

(14)

Simultaneously, for the isothermal portion of the wall, the RHS of the energy equation is set to zero. The specification of density for a composite wall follows the identical procedure for a symmetry plane, where equation (13) is satisfied due to the no-slip condition.

Fluent computational methodology

The incompressible reactive flow solver in ANSYS Fluent v2024 R1 is used for the laminar premixed flame computations reported in this paper. In this framework, variations in density are only due to temperature changes and the pressure is assumed to be spatio-temporally constant. We use the incompressible flow solver since prior works have shown that for propagating flames stabilized at low Mach numbers, the resolution of acoustic waves is not strictly essential for the correct determination of the flame dynamics. At low Mach numbers, the magnitude of pressure and temperature fluctuations imposed by an acoustic wave are weak, and heat release fluctuations caused by an acoustic wave is mainly due to velocity fluctuations and the hydrodynamic processes induced by these oscillations, which are well captured by an incompressible framework.^47,48 Furthermore, an incompressible approach precludes the possibility of acoustic resonances in the system and does not demand carefully tailored boundary conditions and input signal amplitude selection to avoid non-linear effects.⁴⁷ The validity of the incompressible approach is also asserted by comparison of the FTFs with a fully compressible approach in this paper (see Figure 10).

The thermodynamic properties of the reactive mixture are computed using the “ideal-gas-mixing-law” in Fluent. This framework requires the specification of the species viscosity, thermal conductivity, specific heats, and enthalpies. The species-specific heats and enthalpies are specified as polynomial functions of temperature, while the viscosity and thermal conductivity are computed using the “kinetic-theory” model, which uses the Lennard–Jones characteristic length and energy parameters.⁴⁹ The mass diffusion coefficients of each species are also determined using the “kinetic-theory” model, which implements the mixture-averaged transport model approach. The binary diffusion coefficients of each species pair is calculated using the Lennard–Jones characteristic length and energy parameter. The equations involved in the implementation of the mixture-averaged transport model are listed in Table 1. In addition to the mixture-averaged approach, computations are also performed using the unity-Lewis approach, where the species diffusion coefficients are computed by assuming a Lewis number of unity.

With regard to the boundary conditions, a “velocity-inlet” condition is imposed at the inlet which specifies the flow velocity, temperature, and species mass fractions. For the acoustic forcing at the inlet, a time-varying inlet velocity profile is specified using a user-defined function. Isothermal/adiabatic wall boundary conditions are imposed at the top boundary, while a “Symmetry” condition is imposed for the bottom boundary. At the outlet, the “pressure-outlet” condition is imposed which specifies a gauge pressure. The “Finite-Rate” chemistry model is used with a stiff chemistry solver utilizing the “Direct Integration” methodology. The SIMPLE scheme⁵⁰ is employed for the pressure–velocity coupling. Under-relaxation factors of $0.3$ for pressure and $0.7$ for momentum are employed. Table 3 lists the numerical schemes used for spatial and temporal discretization. A constant time step of $5 \times 10^{- 7}$ s, corresponding to a convective CFL of 0.32 is used for all computations.

Table 3.

Numerical schemes used for time and space discretization in Fluent.

Discretization variable/operator	Scheme
Time	Second-order implicit
Gradient	Least squares cell based
Pressure	Second-order
Species	Second-order upwind
Energy	Second-order upwind

Results

This section presents the key results of this article, and is organized into two key subsections. The first subsection presents the static and dynamic flame characteristics obtained from the different flow solvers with the mixture-averaged transport model, while the second subsection discusses the impact of the transport models on the hydrogen flame characteristics.

Comparison of flow solvers

First, it is instructive to compare the unstretched laminar flame speed obtained from the various solvers employing the mixture-averaged (MA) transport model. Figure 3 plots the values of the unstretched laminar flame speeds obtained from the 1D computations of hydrogen flames employing different solvers (the tables are given in the Appendix). The operating conditions are chosen in the low equivalence ratio regime where the thermo-diffusive effects are prominent.¹⁷ The reduced version of the San-Diego mechanism designed for hydrogen and ammonia-hydrogen combustion under preheated conditions proposed by Jiang et al.⁵³ is used for the chemical kinetics description in all solvers employed in this study. First, all solvers predict the unstretched laminar flame speed with a good quantitative accuracy in comparison to Cantera. Second, the influence of Soret diffusion has a noticeable impact on the unstretched laminar flame speed of the hydrogen flames. Indeed, the $S_{l}^{0}$ obtained from Cantera shows a $6 %$ reduction when incorporating the Soret effect. This change also explains the relative differences in the $S_{l}^{0}$ values between the solvers. S3D, with the inclusion of Soret diffusion, exhibits a slightly lower flame speed than the OpenFOAM and AVBP solvers which neglect this effect. Encouragingly, all these solvers (S3D, AVBP, OpenFOAM) agree very well with the corresponding Cantera counterparts to within $4 %$ . However, the Fluent results, though computed with neglecting Soret diffusion, exhibit a slightly higher difference of $11 %$ in relation to the corresponding Cantera values. In summary, it is encouraging to observe that all flow solvers predict the laminar flame speed to within $11 %$ of the corresponding Cantera values, with Fluent exhibiting the maximal differences. Furthermore, Soret diffusion has a non-negligible impact on the unstretched laminar flame speed. Next, we examine the mean flame structure obtained with the various flow solvers.

Figure 3.

Unstretched laminar flame speeds computed for the one-dimensional (1D) premixed hydrogen flame at equivalence ratios of $0.20, 0.25, 0.30$ , and $0.35$ from various flow solvers. The mixture-averaged (MA) transport model is used for all the computations.

Figures 4 and 5 depict the flame structure (contours of the heat release rate) for the H $_{2}$ flame at an equivalence ratio of $0.20$ computed with OpenFOAM, Fluent, and S3D. The laminar flame speed predicted by Fluent is quantitatively lower than that predicted with OpenFOAM by $18 %$ at this equivalence ratio, which results in an increased flame length for a fixed flow velocity, as observed in Figure 4. On the other hand, the $S_{l}^{0}$ predicted by S3D and OpenFOAM are in better agreement ( $7.5 %$ ) at this condition, resulting in identical flame lengths as in Figure 5. However, from a flame dynamics perspective, the distance between the flame anchoring point and the centroid of the heat release distribution is the key length scale that governs the time lag of the FTF.⁴ The spatial location of the centroid of the heat release distribution, plotted with the red “+” in Figures 4 and 5 computed with the different solvers are in close agreement. This observation suggests that despite the different $S_{l}^{0}$ (and consequently, the flame length) predicted by the solvers, the phase of the FTF would be quantitatively the same. An additional aspect noticeable in Figures 4 and 5 is the variations in the heat release rate along the flame sheet. The positively curved regions near the flame root exhibit higher local flame speeds (and therefore, higher heat release rates) in comparison to the negatively curved regions at the flame tip. The low local consumption speeds at the tip of the hydrogen flame is a well-known consequence of thermo-diffusive effects,⁵⁴ which is captured by all solvers.

Figure 4.

Heat release rate contours of hydrogen flame at an equivalence ratio of $0.20$ computed using OpenFOAM (top) and Fluent (bottom). The red “+” denotes the location of the centroid of the heat release distribution. MA transport model is used in both solvers. MA: mixture-averaged.

Figure 5.

Heat release rate contours of hydrogen flame at an equivalence ratio of $0.20$ computed using OpenFOAM (top) and S3D (bottom). The red “+” denotes the location of the centroid of the heat release distribution. MA transport model is used in both solvers. MA: mixture-averaged.

A key parameter that quantifies the flame speed enhancement due to thermo-diffusive effects is the ratio of the global consumption speed ( $S_{c}^{g}$ ) to the unstretched laminar flame speed, also popularly known as the stretch factor ( $I_{0}$ ) in combustion literature.⁵⁵ The global consumption speed can be obtained by applying mass balance and is given by

S_{c}^{g} = \frac{{\dot{m}}_{R}}{ρ_{R} A_{c = 0.8}}

(15)

where

{\dot{m}}_{R}

is the reactant mass flow rate,

ρ_{R}

is the density of reactants, and

A_{c = 0.8}

is the area of the flame sheet represented by the surface along which the progress variable is

0.8

.⁵⁶ The progress variable is defined as

c = 1 - \frac{Y_{H_{2}}}{Y_{H_{2}}^{u}}

(16)

where

Y_{H_{2}}

is the mass fraction of hydrogen, and the superscript “

u

” represents the mass fraction of hydrogen in the unburnt mixture. For the present 2D slit-flame setup, equation (15) reduces to

S_{c}^{g} = \frac{b u_{0}}{L_{c = 0.8}} ⟹ \frac{S_{c}^{g}}{S_{l}^{0}} = \frac{b}{L_{c = 0.8}} \cdot \frac{u_{0}}{S_{l}^{0}}

(17)

where

L_{c = 0.8}

is the length of the

c = 0.8

curve.

Figure 6 plots the ratio $S_{c}^{g} / S_{l}^{0}$ for different equivalence ratios obtained with the flow solvers (the tables are given in the Appendix). At low equivalence ratios, thermo-diffusive effects associated with the hydrogen flame result in global consumption speeds higher than the unstretched laminar flame speed. As the equivalence ratio increases, thermo-diffusive effects diminish, and the ratio $S_{c}^{g} / S_{l}^{0}$ converges to unity. It is also interesting to note that the intricate differences in the implementation of the mixture-averaged transport model and the numerical methods used to solve the governing equations in the flow solvers produce variations in the values of the stretch factor. Across all conditions, S3D and Fluent predict higher flame speed enhancement compared to the other flow solvers. Thus, the differences in the implementation of the transport models, combined with the different numerical methods, cause non-negligible variations in the stretch factor of laminar premixed hydrogen flames when computed using different flow solvers.

Figure 6.

Comparison of the ratio of the global consumption speed to the unstretched laminar flame speed corresponding to the two-dimensional (2D) slit-stabilized flame obtained with different flow solvers. MA transport model is used in all solvers. MA: mixture-averaged.

Next, the response of the laminar premixed hydrogen flames computed using the different solvers is presented and compared. The flame dynamics are quantified via the heat release rate transfer function (equation (1)), with the reference location for the velocity fluctuations being the inlet of the domain. Two methods are used to compute the FTF. First, by imposing velocity fluctuations at the inlet which are purely harmonic in nature, that is

u^{'} (x = 0, y, t) = u_{a} \cdot P (y) \cdot \sin ω t

(18)

where

u_{a}

denotes the amplitude of velocity fluctuations,

P (y)

is the kernel which describes the variation of velocity fluctuations in the transverse direction, and

ω

is the angular frequency of the imposed disturbances. The kernel function should ensure that the velocity fluctuations go to zero at the walls, and is represented by a hyperbolic tangent profile in S3D, Fluent, and OpenFOAM, and is given by

P (y) = 0.5 (1 + \tanh (\frac{y - y_{0}}{w}))

(19)

where

y_{0}

denotes the transverse location of the variation (here, close to the wall), and

w

is the characteristic thickness of the hyperbolic tangent profile, taken as

0.08

mm in this work. In AVBP, the kernel function is identically equal to unity across the inflow patch. From a numerical standpoint, this implies that the temperature fluctuations induced by the incoming acoustic wave,

T_{a}^{'}

, must generate an entropy disturbance, represented by temperature fluctuations

T_{s}^{'}

, confined to the wall in order to satisfy the iso-thermal wall condition

T_{w}^{'} = 0

. This entropy perturbation is small in amplitude, does not generate sound, and is dissipated downstream of the flame anchoring point, where an adiabatic wall is imposed. Overall, it is considered to have a negligible influence on the flame dynamics.

Alternatively, the imposed time variation of the velocity fluctuations can also take the form of a discrete random binary signal (DRBS), which essentially represents a signal containing a superposition of many frequencies. In this work, a low-pass filtered DRBS signal which takes binary values corresponding to the maximum and minimum velocity fluctuations is used following Huber.⁵⁷ The complete impulse response and the FTF can then be obtained with one computation using system identification procedures as explained in detail in Polifke et al.⁵⁸ and Polifke.⁸ The optimal length of the impulse response is decided by gradually increasing it from $100$ until a value where the FTF shows less than $5 %$ change in the gain and phase. Furthermore, the standard deviation is computed using a bootstrapping method where the FTF computed by randomly drawing several time series samples from the original data set is statistically analyzed.⁵⁹ Figure 7 shows a representative impulse response obtained for the hydrogen flame at $ϕ = 0.20$ using the AVBP solver. After the application of the input signal, the flame responds after a characteristic time delay (corresponding to the time at which the impulse response is maximum) and is distributed over a characteristic time interval—the spread of the impulse response (the range of time over which the impulse response is non-zero).

Figure 7.

Unit impulse response of the hydrogen flame at an equivalence ratio of $ϕ = 0.2$ , computed using AVBP. The shaded area indicates the standard deviation, $σ$ , associated with each coefficient of the impulse response.

In an open-loop framework, it is advisable to impose non-reflecting boundary conditions at the outlet to prevent possible reflected waves from altering the flame dynamics, and hence the FTF. Reflected waves at the outflow would introduce pressure and velocity perturbations. Since acoustically compact premixed flames are more sensitive to velocity disturbances than to pressure disturbances,¹⁷ and because the reflection delay time ( $\approx 0.035$ ms) is much shorter than the characteristic flame time scale ( $\approx 0.6$ ms, Figure 7), any reflections, if present, will manifest themselves as a slight modulation of the reference velocity, and thus, are unlikely to have a significant impact on the FTF.

The quantitative differences in the FTFs obtained from the various solvers and transport models are quantified in terms of the root-mean-square-difference (RMSD) defined as⁶⁰

RMSD = \frac{{[(1 / N) \sum_{i = 1}^{N} {(y_{i} - y_{i}^{ref})}^{2}]}^{1 / 2}}{max (y_{i}^{ref}) - min (y_{i}^{ref})}

(20)

where

y_{i}

represents the computed FTF metric (real/imaginary part or gain/phase) and

y_{i}^{ref}

is the computed FTF metric using the reference solver. The denominator of the RMSD expression represents the range of values that the considered FTF metric can take, and the summation in the numerator is performed over all discrete frequency points. For the gain, the denominator is usually close to unity, while for the phase it is

π

. Since premixed FTFs exhibit low-pass filter behavior, the phase of the FTF has no significance at high frequencies. To prevent the non-physical values of phase at high frequencies from contaminating the RMSD metric, the data at frequencies wherein the gain is lower than

0.05

is discarded.

Figure 8 plots the FTFs computed for the hydrogen flames under the conditions considered in this work using the different solvers. The inlet mean velocity is maintained to be thrice the unstretched laminar flame speed obtained from Cantera, and is therefore fixed for a given equivalence ratio across solvers. The FTFs are computed by imposing a time-harmonic signal at multiple discrete frequencies in S3D, and using the DRBS combined with system identification for OpenFOAM, Fluent, and AVBP. To verify the correctness of the implemented system identification techniques to reconstruct the FTF when using the DRBS signal, companion harmonically forced computations were also performed in OpenFOAM (not shown here). For all the equivalence ratios, an excellent match between the two approaches was observed, thus verifying the implemented system identification procedure.

Figure 8.

Flame transfer functions of hydrogen flames at equivalence ratios of 0.20 (a), 0.25 (b), 0.30 (c), and 0.35 (d) computed with different solvers. Mean inlet velocities of $4$ , $7$ , $10$ , and $13$ m s⁻¹ are imposed for the equivalence ratios of $0.20, 0.25, 0.30$ , and $0.35$ , respectively. The symbols denote the results from harmonically forced computations, while the lines denote the FTF reconstructed using a DRBS signal and system identification procedures. FTF: flame transfer function; DRBS: discrete random binary signal.

Two key observations can be deduced by comparing the FTFs computed with the different solvers in Figure 8. First, the phase of the FTF, which contains the information about the time lag(s) characteristic of the flame, specifically, to the fundamental mechanism that causes the flame response, agrees extremely well across the solvers. Discrepancies, if any, between the solvers are only observed at frequencies where the gain is low. Thus, the essential physics of the physical mechanism(s) resulting in the flame response is captured by all solvers. Second, with regard to the gain, all the solvers still agree well with each other, but noticeable quantitative differences are observed, especially between OpenFOAM and the other solvers. S3D, AVBP, and Fluent agree very closely in gain with each other. Quantitative estimates on the differences between the solvers can be obtained by considering the RMSD metric (equation (20)) tabulated in Table 4. The RMSD of the phase is within $10 %$ across all solvers, while for the gain the RMSD is within $15 %$ . The OpenFOAM FTFs show maximal differences to other solvers, especially with regard to the gain. To summarize, despite the differences in the numerical methods, transport model implementation and boundary treatment between the CFD solvers, all of them show a good agreement in the FTF, which, in terms of the RMSD metric, agrees to within $15 %$ in gain and $10 %$ in phase.

Table 4.

Maximum and minimum values of the RMSD (equation (20)) of the (real, imaginary) part and (gain, phase) of the transfer function compared for different solvers and transport models.

Quantity	Fuel	Solver/transport model combination	Max ( $%$ )	Min ( $%$ )
Real(FTF)	H $_{2}$	OpenFOAM–Fluent	$9.69$	$4.10$
		Fluent–AVBP	$6.66$	$5.25$
		OpenFOAM–AVBP	$12.36$	$6.93$
		MA–Le $_{1}$ (OpenFOAM)	$22.17$	$10.43$
Imag(FTF)	H $_{2}$	OpenFOAM–Fluent	$9.69$	$3.79$
		Fluent–AVBP	$7.00$	$6.54$
		OpenFOAM–AVBP	$11.55$	$6.54$
		MA–Le $_{1}$ (OpenFOAM)	$21.90$	$10.38$
Gain(FTF)	H $_{2}$	OpenFOAM–Fluent	$12.70$	$4.62$
		Fluent–AVBP	$5.18$	$4.35$
		OpenFOAM–AVBP	$13.45$	$8.37$
		MA–Le $_{1}$ (OpenFOAM)	$9.57$	$7.56$
Phase(FTF)	H $_{2}$	OpenFOAM–Fluent	$3.23$	$0.94$
		Fluent–AVBP	$9.10$	$4.45$
		OpenFOAM–AVBP	$5.71$	$4.05$
		MA–Le $_{1}$ (OpenFOAM)	$24.70$	$17.54$

The comparison between transport models is performed for a specific solver listed in parentheses, while the comparison between solvers is done for a particular transport model (MA). For a particular solver/transport model combination, the reference is always the first listed value, i.e. for the OpenFOAM-Fluent pair, the reference is OpenFOAM. RMSD: root-mean-square-difference; MA: mixture-averaged; FTF: flame transfer function.

Next, we examine the influence of variations in the mean inlet velocity profile and the role of pressure fluctuations (compressibility) on the FTFs. To get insight into the former, a plug flow inlet velocity profile and a parabolically varying inlet velocity profile for the identical area-averaged mean velocity are considered. These are shown in Figure 9, where it is to be noted that both profiles have identical area-averaged velocity magnitude. While the parabolic profile is the exact solution of the Navier–Stokes equations in this configuration, the plug flow profile is imposed by default in OpenFOAM and Fluent. To get insight into the role of compressibility, we compute the FTFs of the hydrogen flame in Fluent using the “incompressible-ideal-gas” (weakly compressible) and “ideal-gas” (fully compressible) models. In the weakly compressible approach, density fluctuations are only due to changes in temperature, and the pressure in the domain is assumed to be constant. Figure 10 compares the FTFs obtained using these approaches at two equivalence ratios. In accordance with previous works,^47,48 the role of pressure fluctuations on the flame dynamics of the H $_{2}$ flame is indeed minor. On the other hand, noticeable differences, especially in the phase of the FTF, are obtained with the different inlet velocity profiles.

Figure 9.

Mean velocity profile at a location just upstream of the flame front (6 mm). The velocity profiles are plotted for a plug flow and a parabolic channel flow inlet velocity profile. The area-averaged inlet velocity is identical ( $= 10$ m s⁻¹) for both cases. FTF: flame transfer function.

Figure 10.

Comparison of the flame transfer functions computed with a fully compressible approach and a weakly compressible approach in Fluent at an equivalence ratio of $0.25$ (a) and $0.35$ (b). Also shown is the comparison between the FTFs obtained with a plug flow and a parabolic inlet velocity profile for the same area-averaged inlet velocity. FTF: flame transfer function.

The fact that the FTF phase is sensitive to the transverse variation of the velocity profile can be explained by the time-lag characteristic to the flame, which, in this case, is given by the time taken for a flame sheet disturbance to convect from the root to the centroid of the flame with a speed corresponding to the local flow velocity along the flame-tangential direction.⁶¹ This time delay can be mathematically expressed as

τ = \int_{0}^{S_{f}} \frac{d S}{u_{t}}

(21)

where the integration is performed along the flame contour (

S

S_{f}

is the location of the flame centroid, and

u_{t}

is the flame-tangential component of the local mean flow velocity. The key difference comes from the tangential velocity along the flame sheet. Although both the profiles have identical area-averaged inlet velocity, the velocity “seen” by the flame sheet disturbance during convection is different, since the transverse distribution of the velocity varies, which results in different convection times of the flame sheet disturbance from the root to the centroid.

Finally, we conclude this section by summarizing the computational times taken for the FTF calculations using the various solvers and present an analysis for the optimal grid to use in accurately computing the FTFs of the H $_{2}$ laminar slit-stabilized flame. Table 5 lists the wall-clock time taken to compute the FTF using the different solvers. All the computations were carried out for a flow time of $0.1$ s, which was sufficient to accurately reconstruct the FTF using the DRBS signal and system identification technique. The explicit codes are observed to be more expensive than the semi-implicit codes with the exception of AVBP, which gives comparable estimates in computational time to that of OpenFOAM. However, it is important to keep in mind that the computational times reported in Table 5 are very specific to the solver settings used in this study, which have not been optimized for each solver. Furthermore, the computational times are also reported “as-is” and not penalized taking into account the accuracy of the solution. For instance, the S3D solver features high-order numerical methods that are much more accurate than the others, which automatically implies longer compute times. Therefore, the figures given in Table 5 should be regarded as specific estimates, particular to the current case and solver settings, and are not intended to benchmark or assess the computational performance of the solvers.

Table 5.

Wall-clock time taken to compute the FTF using the various CFD solvers.

Solver	Wall-clock time (h)
Fluent	$93.75$
AVBP	$130$
OpenFOAM	$136.5$
S3D	$1152.4$

The computations are performed for a time of $0.1$ s using 20 cores of AMD EPYC 7763 processor. FTF: flame transfer function; CFD: computational fluid dynamics.

Furthermore, to determine the “cheapest” setup that yields a good prediction of the FTF, computations of the hydrogen flame dynamics at $ϕ = 0.20$ with varying number of grid points ( $N_{f}$ ) across the flame thickness were performed in OpenFOAM. The results, shown in Figure 11, demonstrate that with increasing the number of grid points the FTFs progressively converge. Treating the $N_{f} = 30$ case, which is the grid used for all computations reported in this paper, as the reference, the root-mean-square error (RMSE) in the gain and phase vary from $3 %$ for the $N_{f} = 20$ grid to $15 %$ for the $N_{f} = 8$ grid. The RMSE remained within $10 %$ for all grids that had at least $10$ points within the flame thickness. Thus, $N_{f} = 10$ is the optimal number of points within the flame thickness to ensure an accurate determination of the FTF in the current configuration. Identical results were observed when repeating this exercise in AVBP as well.

Figure 11.

FTFs of the hydrogen flame at $ϕ = 0.20$ computed with varying numbers of grid points ( $N_{f}$ ) within the flame thickness. All computations were performed with the mixture-averaged transport model in OpenFOAM. FTF: flame transfer function.

Comparison of transport models

In this section, a comparison of the flame dynamics with different transport models is presented. Three transport models are considered: the mixture-averaged transport model (designated as “MA” henceforth), the transport model with specification of individual species Lewis numbers corresponding to the unburnt state of the reactant mixture (designated as “Le $_{i}$ ” henceforth), and the transport model where all species Lewis numbers are set to unity (designated as “Le $_{1}$ ”).

First, it is instructive to compare the unstretched laminar flame speeds for a given solver (S3D) employing different transport models (see Figure 12, values are listed in the Appendix). It is evident that the unity-Lewis approach, due to neglect of the fast molecular diffusion of hydrogen, gives an incorrect estimate of $S_{l}^{0}$ with differences as high as $30 %$ in comparison to the mixture-averaged transport model for hydrogen. With the Le $_{i}$ transport model, where the individual species numbers are specified as non-unity constants, the prediction of $S_{l}^{0}$ greatly improves with maximal differences of only $9 %$ compared to a mixture-averaged approach. These observations can now be used to interpret the structure of the mean flames obtained using the various transport models.

Figure 12.

Comparison of the unstretched laminar flame speed obtained in S3D using the MA, $L e_{i}$ , and Le $_{1}$ transport models. MA: mixture-averaged.

Figures 13 and 14 depict the mean flame structure in terms of the normalized heat release rate ( $\dot{Q} / {\dot{Q}}_{max}$ ) obtained when using the transport models accounting for and neglecting the thermo-diffusive effects, for an equivalence ratio of $0.20$ and $0.35$ , respectively. First, consider Figure 13. It is evident that the Le $_{i}$ transport model captures the variation in the heat release distribution along the flame sheet attributed to thermo-diffusive effects and shows excellent agreement in both the flame length and the flame structure in comparison to the MA transport model (in Figure 5). On the other hand, the hydrogen flame computed using the Le $_{1}$ transport model shows a uniformly distributed heat release rate along the flame sheet (except close to the root). As a consequence of this, the locations of the centroid of the individual heat release distributions are also different. Next, consider the flame structure for the H $_{2}$ flame at $ϕ = 0.35$ computed with the MA and Le $_{1}$ approaches (Figure 14). Under this condition, while the distribution of the heat release rate along the flame sheet is uniform even for the MA model (due to reduced propensity of thermo-diffusive effects), the flame lengths are drastically different. This is due to the Le $_{1}$ flame exhibiting a reduced flame speed compared to the MA flame (see Figure 12).

Figure 13.

Normalized heat release rate contours of hydrogen flame at an equivalence ratio of $0.20$ computed using S3D with $L e_{i}$ transport model (top) and Le $_{1}$ transport model (bottom). The red “+” denotes the location of the centroid of the heat release distribution.

Figure 14.

Normalized heat release rate contours of hydrogen flame at an equivalence ratio of $0.35$ computed using OpenFOAM with MA transport model (top) and Le $_{1}$ transport model (bottom). The red “+” denotes the location of the centroid of the heat release distribution. MA: mixture-averaged.

In summary, when using a unity-Lewis approach to compute the flame structure, drastic differences may be observed, first due to neglect of thermo-diffusive effects resulting in an incorrect heat release distribution, and second due to the different $S_{l}^{0}$ predicted by the Le $= 1$ approach resulting in different flame lengths. Encouragingly, the Le $_{i}$ approach captures key features associated with the MA flame such as the heat release distribution and flame length. Indeed, Figure 15 plots the ratio of the global consumption speed to the unstretched laminar flame speed obtained for the H $_{2}$ flame in S3D using different transport models. The Le $_{1}$ approach, as expected, predicts no flame speed enhancement across all conditions, but the Le $_{i}$ transport model does predict the ratio $S_{c}^{g} / S_{l}^{0}$ to be greater than unity at low equivalence ratios. While the value predicted by the Le $_{i}$ model does not agree quantitatively with the MA model, it is nevertheless able to capture the fact that the H $_{2}$ flame exhibits an increased flame speed at low equivalence ratios. Next, we compare the FTFs obtained with the different transport models.

Figure 15.

Comparison of the ratio of the global consumption speed to the unstretched laminar flame speed computed from the S3D solver with different transport models.

Figure 16 compares the FTFs obtained for hydrogen flames using the mixture-averaged and the unity-Lewis number transport model in OpenFOAM. Noticeable differences not only in the gain, but also in the phase of the transfer functions are observed. The reason for this can be easily deduced from Figure 12. For a given solver, the unstretched laminar flame speed computed with a unity-Lewis transport model is different from the corresponding mixture-averaged flame speed. For a fixed flow velocity, this results in a different flame length, which shifts the centroid of the heat release distribution, resulting in different time delays (see e.g. Figures 13 and 14). Furthermore, with regard to the gain, the unity-Lewis transport model does not capture the thermo-diffusive effects. Therefore, it does not capture the influence of flame speed fluctuations on the heat release dynamics. In addition, due to the different unstretched laminar flame speed and the global consumption speed between the unity-Lewis and mixture-averaged approaches, the flame shapes are different resulting in different area responses. It is, thus, no surprise that the gain of the FTF changes when using mixture-averaged and unity-Lewis transport models. In terms of the quantitative RMSD metric, the phase of the FTFs shows differences as high as $25 %$ between the mixture-averaged and the unity-Lewis transport models (see Table 4), while the differences in the gain between the two approaches are within $10 %$ .

Figure 16.

Flame transfer functions of hydrogen flames at equivalence ratios of $0.20, 0.25, 0.30$ , and $0.35$ (top to bottom) computed with MA and Le $_{1}$ transport models in OpenFOAM. The mean inlet velocity imposed for each equivalence ratio is the same as that used for Figure 8.

Finally, we assess the effectiveness of the Le $_{i}$ transport model, where the individual species Lewis numbers are specified as non-unity constants, in the computation of the FTFs of hydrogen flames. Figure 17 plots the FTFs obtained from the MA, Le $_{1}$ , and Le $_{i}$ transport models using the S3D solver. While noticeable deviations in the gain and phase are observed between the MA and Le $_{1}$ transport model, the Le $_{i}$ transport model yields results that are very identical to the more detailed mixture-averaged approach. This result gives impetus for the use of the Le $_{i}$ transport model in yielding quantitatively accurate estimates of the acoustic response of hydrogen flames. Furthermore, the Le $_{i}$ transport model is also computationally efficient and easy to implement in flow solvers. In computational terms, it was observed that the Le $_{i}$ transport model exhibits a $50 %$ speedup in computational time for a given FTF calculation in comparison to the MA transport model.

Figure 17.

Flame transfer functions of H $_{2}$ flames at equivalence ratios of $0.20$ , $0.25$ , and $0.30$ (top to bottom) computed with MA, Le $_{1}$ , and $L e_{i}$ transport models in S3D.

In summary, the results of this section demonstrate that the proper selection of transport model is crucial in the computation of FTFs of hydrogen flames. First, using the unity-Lewis transport model for hydrogen can result in inconsistent results due to the inability of this model in capturing the thermo-diffusive effects. Indeed, for the current operating conditions, the quantitative deviations in the FTF of the unity-Lewis approach in relation to the more detailed mixture-averaged transport model were about $25 %$ in phase. The Le $_{i}$ transport model is an ideal alternative, both in terms of accuracy and cost, since the FTFs predicted by this approach match well with the MA model and offers a $50 %$ computational speedup.

Conclusions

The advent of carbon-free fuels, with different reactivity and transport properties in relation to conventional fuels, for sustainable power production calls for robust computational tools to accurately predict the combustion characteristics of these fuels. Therefore, it is important to compare existing CFD solvers from the standpoint of static and dynamic flame characteristics of these carbon-free fuels. In this paper, computations of laminar premixed flame shapes and dynamics (response to acoustic perturbations) with four CFD solvers (ANSYS Fluent, OpenFOAM, AVBP, and S3D) were performed and compared. It was first shown that, despite using the identical chemical kinetic mechanism and similar transport models, quantitative differences up to $18 %$ were observed in the prediction of the most fundamental quantity associated with laminar premixed flames—the unstretched laminar flame speed $(S_{l}^{0})$ . These differences consequently translate to differences in the time-averaged flame structure, which was found to be rather minimal. In terms of the flame dynamics, all solvers gave similar quantitative and qualitative FTF characteristics. The quantitative differences were characterized in terms of the RMSD metric which was within $10 %$ in phase and within $15 %$ in gain between the different solvers. Furthermore, the role of pressure fluctuations on the flame response was found to be rather minimal in this configuration, but the transverse variation of the inlet velocity had a noticeable effect on the dynamics.

The second key contribution of this work is the observation that the incorporation of non-unity-Lewis numbers of hydrogen is highly essential in capturing the correct flame dynamics. Using a unity-Lewis approach gives inaccurate estimates of the unstretched laminar flame speeds, and does not capture key processes characteristic of hydrogen flames such as flame speed enhancement due to stretch, non-uniform heat release distributions along the flame sheet and tip-opening. Furthermore, assuming all species Lewis numbers to unity results in errors as high as $25 %$ in the calculation of the hydrogen flame FTFs. Incorporation of the thermo-diffusive effects can be done via a mixture-averaged transport model where the mass diffusion coefficient of each species is determined from the binary diffusion coefficients. Alternatively, the diffusion coefficients can also be computed by specifying a non-unity constant Lewis number of each species. In the current paper, it is shown that the latter approach, despite its simplicity, captures important thermo-diffusive effects associated with hydrogen flames and yields quantitatively similar estimates of FTFs in relation to the more detailed mixture-averaged approach. Therefore, this approach can be used in flow solvers to obtain flame dynamics of hydrogen flames with time-efficient and easy computational implementation.

Footnotes

ORCID iDs

Harish S Gopalakrishnan

Mirko R Bothien

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Harish S. Gopalakrishnan was funded by the Swiss National Science Foundation (SNSF) under the project ADONIS (Ammonia hydrogen combustion in micro gas turbines) with the grant number 206244.

Declaration of conflicting interests

The authors declare no potential conflicts of interest with respect to research, authorship, and publication of this article.

Appendix

This section presents the key tables of the unstretched laminar flame speed computed from the various solvers (Table 6) and transport models (Table 7) and the global consumption speed computed using the various solvers and transport models (Table 8).

References

Lieuwen

Yang

. Combustion instabilities in gas turbine engines: operational experience, fundamental mechanisms, and modeling. Reston, VA: American Institute of Aeronautics and Astronautics, 2005.

Schuller

Durox

Candel

. A unified model for the prediction of laminar flame transfer functions: comparisons between conical and V-flame dynamics. Combust Flame 2003; 134: 21–34.

Strahle

. Combustion noise. Prog Energy Combust Sci 1978; 4: 157–176.

Ducruix

Durox

Candel

. Theoretical and experimental determinations of the transfer function of a laminar premixed flame. Proc Combust Inst 2000; 28: 765–773.

Freitag

Konle

Lauer

, et al.

Pressure influence on the flame transfer function of a premixed swirling flame

Turbo Expo: power for land, sea, and air. American Society of Mechanical Engineers (ASME), Vol. 42363, 477–486.

Daviller

Oztarlik

Poinsot

. A generalized non-reflecting inlet boundary condition for steady and forced compressible flows with injection of vortical and acoustic waves. Comput Fluids 2019; 190: 503–513.

Polifke

Wall

Moin

. Partially reflecting and non-reflecting boundary conditions for simulation of compressible viscous flow. J Comput Phys 2006; 213: 437–449.

Polifke

. Black-box system identification for reduced order model construction. Ann Nucl Energy 2014; 67: 109–128.

Bothien

Moeck

Lacarelle

, et al. Time domain modelling and stability analysis of complex thermoacoustic systems. Proc Inst Mech Eng Part A: J Power Energy 2007; 221: 657–668.

10.

Dowling

. The calculation of thermoacoustic oscillations. J Sound Vib 1995; 180: 557–581.

11.

Æsøy

Dawson

Moeck

. Systematic modulation of the flame transfer function and its effect on thermoacoustic stability. Combust Flame 2024; 265: 113494.

12.

Æsøy

Jankee

Yadala

, et al. Suppression of self-excited thermoacoustic instabilities by convective-acoustic interference. Proc Combust Inst 2023; 39: 4611–4620.

13.

Æsøy

Aguilar

Wiseman

, et al. Scaling and prediction of transfer functions in lean premixed H

_{2}

/CH

_{4}

-flames. Combust Flame 2020; 215: 269–282.

14.

Blondé

Schuermans

Pandey

, et al. Effect of hydrogen enrichment on transfer matrices of fully and technically premixed swirled flames. J Eng Gas Turbine Power 2023; 145: 121009 1–11.

15.

Ghani

Polifke

. Control of intrinsic thermoacoustic instabilities using hydrogen fuel. Proc Combust Inst 2021; 38: 6077–6084.

16.

Kido

Nakahara

Nakashima

, et al. Influence of local flame displacement velocity on turbulent burning velocity. Proc Combust Inst 2002; 29: 1855–1861.

17.

Lieuwen

. Unsteady combustor physics. New York, NY: Cambridge University Press, 2021.

18.

Wiseman

Gruber

Dawson

. Flame transfer functions for turbulent, premixed, ammonia-hydrogen-nitrogen-air flames. J Eng Gas Turbine Power 2023; 145: 031015.

19.

Poinsot

Veynante

. Theoretical and numerical combustion. Philadelphia, PA: RT Edwards, Inc., 2005.

20.

Kee

Dixon-Lewis

Warnatz

, et al. Transport: a software package for the evaluation of gas-phase, multicomponent transport properties. Chemkin Collect 1999; 1: 12–28.

21.

Greenshields

. OpenFOAM v10 user guide. London, UK: The OpenFOAM Foundation, 2022.

22.

Chen

Choudhary

De Supinski

, et al. Terascale direct numerical simulations of turbulent combustion using S3D. Comput Sci Discov 2009; 2: 015001.

23.

Schönfeld

Rudgyard

. Steady and unsteady flow simulations using the hybrid flow solver AVBP. AIAA J 1999; 37: 1378–1385.

24.

Goodwin

Moffat

Schoegl

, et al. Cantera: an object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, Version 3.0.0, 2023. https://www.cantera.org. DOI: 10.5281/zenodo.8137090.

25.

Meindl

Silva

Polifke

. On the spurious entropy generation encountered in hybrid linear thermoacoustic models. Combust Flame 2021; 223: 525–540.

26.

Greenshields

Weller

. Notes on computational fluid dynamics: general principles. Reading, UK: CFD Direct Ltd, 2022.

27.

Burcat

Ruscic

. Third millennium ideal gas and condensed phase thermochemical database for combustion (with update from active thermochemical tables). Technical report, Argonne National Lab. (ANL), Argonne, IL (USA), 2005.

28.

Turns

. Introduction to combustion. Vol. 287. New York, NY, USA: McGraw-Hill Companies, 1996.

29.

Lucchese

Nilsson

. Proceedings of CFD with OpenSource software. Chapter: Implementation of non-reflecting boundary conditions in OpenFOAM. Chalmers University of Technology, 2022.

30.

Selle

Nicoud

Poinsot

. Actual impedance of nonreflecting boundary conditions: implications for computation of resonators. AIAA J 2004; 42: 958–964.

31.

Quartapelle

Selmin

. High-order Taylor–Galerkin methods for nonlinear multidimensional problems. Finite Elements in Fluids 1993; 76: 46.

32.

Colin

Rudgyard

. Development of high-order Taylor–Galerkin schemes for LES. J Comput Phys 2000; 162: 338–371.

33.

Martin

Cheng

Rudgyard

, et al. A new class of Taylor–Petrov–Galerkin schemes using polynomial reconstruction. J Comput Phys 2025; 527: 113797.

34.

Kee

Rupley

Meeks

, et al. Chemkin-III: a Fortran chemical kinetics package for the analysis of gas-phase chemical and plasma kinetics. Technical report, Sandia National Lab. (SNL-CA), Livermore, CA (USA), 1996.

35.

Chapman

. The kinetic theory of simple and composite monatomic gases: viscosity, thermal conduction, and diffusion. Proc R Soc London Ser A, Contain Papers a Mathemat Phys Character 1916; 93: 1–20.

36.

Eucken

. Über das wärmeleitvermögen, die spezifische wärme und die innere reibung der gase. Physikalische Zeitschrift 1913; 14: 324–332.

37.

Wilke

. A viscosity equation for gas mixtures. J Chem Phys 1950; 18: 517–519.

38.

Bird

Stewart

Lightfoot

. Transport phenomena. New York, NY: Wiley, 2006.

39.

Mathur

Tondon

Saxena

. Thermal conductivity of binary, ternary and quaternary mixtures of rare gases. Mol Phys 1967; 12: 569–579.

40.

Hirschfelder

Curtiss

Bird

. The molecular theory of gases and liquids. New York, NY: Wiley, 1964.

41.

Poinsot

Lele

. Boundary conditions for direct simulations of compressible viscous flows. J Comput Phys 1992; 101: 104–129.

42.

Kennedy

Carpenter

. Several new numerical methods for compressible shear-layer simulations. Appl Numer Math 1994; 14: 397–433.

43.

Kennedy

Carpenter

Lewis

. Low-storage, explicit Runge–Kutta schemes for the compressible Navier–Stokes equations. Appl Numer Math 2000; 35: 177–219.

44.

Smooke

. Reduced kinetic mechanisms and asymptotic approximations for methane-air flames: a topical volume. Berlin Heidelberg: Springer-Verlag, 1991.

45.

Sutherland

Kennedy

. Improved boundary conditions for viscous, reacting, compressible flows. J Comput Phys 2003; 191: 502–524.

46.

Gruber

. Direct numerical simulation of turbulent combustion near solid surfaces [Doctoral Thesis]. Trondheim, Norway: Norwegian University of Science and Technology, 2006.

47.

Eder

Silva

Haeringer

, et al. Incompressible versus compressible large eddy simulation for the identification of premixed flame dynamics. Int J Spray Combust Dynam 2023; 15: 16–32.

48.

Han

Morgans

. Simulation of the flame describing function of a turbulent premixed flame using an open-source LES solver. Combust Flame 2015; 162: 1778–1792.

49.

ANSYS, Inc. ANSYS Fluent theory guide, Release 2024, 2024. ANSYS, Inc. https://www.ansys.com/products/fluids/ansys-fluent.

50.

Ferziger

Perić

. Computational methods for fluid dynamics. 2nd ed. Berlin: Springer-Verlag, 1999.

51.

Ahtye

. Thermal conductivity in vibrationally excited gases. J Chem Phys 1972; 57: 5542–5555.

52.

Chapman

Cowling

. The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases. Cambridge, UK: Cambridge University Press, 1990.

53.

Jiang

Gruber

Seshadri

, et al. An updated short chemical-kinetic nitrogen mechanism for carbon-free combustion applications. Int J Energy Res 2020; 44: 795–810.

54.

Law

. Combustion physics. 1st ed. Cambridge: Cambridge University Press, 2006.

55.

Driscoll

. Turbulent premixed combustion: Flamelet structure and its effect on turbulent burning velocities. Prog Energy Combust Sci 2008; 34: 91–134.

56.

Berger

Attili

Pitsch

. Synergistic interactions of thermodiffusive instabilities and turbulence in lean hydrogen flames. Combust Flame 2022; 244: 112254.

57.

Huber

. Impact of fuel supply impedance and fuel staging on gas turbine combustion stability. PhD Thesis, Technical University of Munich, 2009.

58.

Polifke

Poncet

Paschereit

, et al. Reconstruction of acoustic transfer matrices by instationary computational fluid dynamics. J Sound Vib 2001; 245: 483–510.

59.

Tay-Wo-Chong

Bomberg

Ulhaq

, et al. Comparative validation study on identification of premixed flame transfer function. J Eng Gas Turbine Power 2012; 134: 021502 1–8.

60.

Guo

Silva

Polifke

. Robust identification of flame frequency response via multi-fidelity Gaussian process approach. J Sound Vib 2021; 502: 116083.

61.

Boyer

Quinard

. On the dynamics of anchored flames. Combust Flame 1990; 82: 51–65.