Thermoacoustic finite element method modeling of partly autoignition and propagation-stabilized flames

Abstract

This paper describes a methodology to thermoacoustically account for different flame types in a single Finite-Element-computation. To do so, the flame segmentation mechanism presented in previous studies is used to characterize the different flame zones. Differentiation is done between propagation-stabilized shear-layer flames and autoignition flames that respond to acoustic perturbations very differently. While autoignition flames mainly respond to acoustic pressure and temperature fluctuations, propagation-stabilized flames respond to acoustic velocity perturbations. Thus, flame transfer functions specific to each flame type are analytically implemented within the frequency domain Finite-Element-computation. Using the novel framework, it is shown that the global flame transfer function obtained from Computational Fluid Dynamics (CFD) simulations of a backward facing step reheat combustor can be reproduced accurately in the low-frequency regime for an operating point where the flame is forced in a planar manner. The investigated operating point operates on hydrogen fully premixed at lean and autoignitive conditions. The autoignition framework is validated by comparison to one-dimensional direct numerical simulation. The time-averaged CFD heat release rate result is validated with large eddy simulation data.

Keywords

Thermoacoustic autoignition reheat-combustor finite element method

Novelty and significance statement

The novelty of this work lies in the incorporation of two different types of flame in a single thermoacoustic finite element method (FEM) calculation. To the best of our knowledge, this has not yet been shown in the literature before. The significance of this work lies within industrial applications. In industrial reheat combustion chambers, the autoignition flame is typically not purely autoignition stabilized, and flame anchoring regions are present which contribute to the stability of the flame.

Introduction

To fulfill the Paris agreement,¹ a large-scale integration of renewable energy systems in the energy landscape is essential.^2–4 To stabilize and balance the grid, gas turbines running on carbon-free fuels can be an ideal asset.^3,5–9 However, these non-conventional fuels possess very different combustion properties when compared to natural gas. Thus, novel gas turbine combustor development is faced with new challenges to enable the fuel flexible firing capability of an engine. The thermoacoustic analysis to aid the mitigation of combustion instabilities remains a key component in the development phase.¹⁰

There are a variety of methods to perform thermoacoustic analysis of combustion chambers. Direct numerical simulation (DNS)^11,12 and large eddy simulation (LES)¹³ are considered high-fidelity analysis tools, which are capable of analyzing certain thermoacoustic effects in the highest detail. However, they come with a high computational cost and are unfeasible for entire industrial combustor geometries. Alternatively, Helmholtz solvers in combination with flame transfer/describing functions extracted from computations or experiments have been shown to accurately capture the thermoacoustics of combustion chambers.^14–20 Nicoud et al.¹⁴ showed that using a finite element method (FEM) Helmholtz solver the stability of combustion chambers with propagation-stabilized flames could be captured accurately. To do so, the flame transfer function (FTF) was described by an $n - τ$ model.²¹ Laera et al.¹⁵ also used a Helmholtz solver together with a flame describing function to correctly capture the stability of a laboratory scale annular combustion chamber comprised of multiple matrix burners. Silva et al.¹⁹ combined a Helmholtz solver with a flame describing function for a swirl-stabilized combustion chamber to predict the acoustic pressure amplitudes of limit-cycles. Generally, the application of Helmholtz solvers to perform thermoacoustic analysis significantly reduces computational time compared to DNS and LES. Alternatively, network-models have shown to be able to reproduce the combustion chamber acoustics and flame interaction with the highest computational efficiency.^22–26

In this paper, we present an FEM-based approach to thermoacoustically model flames that are partly stabilized by autoignition. The overall flame, which consists of a propagation- and an autoignition part, is analytically described in FEM. The segmentation mechanism from prior studies^17,18 is leveraged to define the different flame types. The segmentation mechanism^17,27 was initially designed to thermoacoustically model non-compact reheat flames under the influence of transverse combustor eigenmodes. In the case of transverse eigenmodes, the acoustic wavelength is of similar length as the flame, thus, resulting in an acoustically non-compact flame. To deal with this, the segmentation mechanism divides the flame into multiple subflames in transverse direction. Each of the very thin subflames is then considered acoustically compact without variation of the fluctuation quantities across its height. Thus, single values for the pressure, temperature and velocity fluctuations act on each subflame. For each of the subflames, a transfer function can be assigned. In the study presented here, the flame is acoustically compact and as such no segmentation is needed. However, the segmentation mechanism is leveraged to assign a specific FTF to each heat release rate (HRR) region depending on its flame stabilization type. To obtain the flame characterization parameters, Reynolds Averaged Navier–Stokes (RANS) simulations are done using Fluent with an in-house reheat flame model, which is an adaptation and extension of the model derived by Kulkarni et al.^28,29 The thermoacoustic analysis is done for a planarly forced backward facing step (BFS) burning hydrogen at elevated pressure.

The presented research is relevant, because in industrial reheat combustors the autoignition flame is never only stabilized by autoignition. Typically, the reheat flame consists of propagation-stabilized HRR regions as well as regions predominantly stabilized by autoignition. Propagation-stabilized flames occur within the shear layers of recirculation zones. Their stabilization is governed by the balance of the flame consumption speed and the flow velocity. With respect to thermoacoustics, the acoustic waves affect the propagation-stabilized flame regions by velocity fluctuations, which modulate the equivalence ratio. In contrast, autoignition stabilized flames are governed by the balance of the autoignition time scale and the flow residence time scale. The fluctuating pressure and especially the isentropic temperature fluctuations induced by acoustic waves affect the local reaction rates and thus the ignition chemistry. This can strongly alter the ignition time and therefore the ignition length, which results in HRR fluctuations. Thereby, the history effect of the fluctuations during the ignition phase matters for reheat flames and is accounted for in current state-of-the-art flame response modeling tools, as presented in Refs. Heinzmann et al.,¹⁷ Heinzmann et al.,²⁷ Zellhuber et al.,³⁰ Gant et al.,³¹ Gopalakrishnan et al.,³² Gopalakrishnan et al.³³

Validation data

The segmentation mechanism developed by Heinzmann et al.^17,18 was already validated in prior studies and initially designed to numerically deal with acoustically non-compact flames. This segmentation mechanism can also be efficiently used for the planar wave case to incorporate different flame types. The model is compared to one-dimensional (1D) DNS³² for a flame solely stabilized by autoignition, and to RANS computations for the partly autoignition- and propagation-stabilized flame in a BFS combustion chamber.

1D DNS simulation

The computational fluid dynamics (CFD) study used for validation in Gopalakrishnan et al.³⁴ is a compressible 1D DNS computation burning hydrogen fully premixed at lean conditions. The operating pressure is atmospheric with an inlet temperature of 1100 K and an inlet flow speed of 200 m/s. The duct of length 0.3 m with non-reflecting boundary conditions was forced with acoustic and entropic waves at the inlet. For more details on the numerical setup the reader is referred to Gopalakrishnan et al.³⁴ The duct is modeled in two-dimensional (2D) in FEM with a length of 0.3 m and a height of 0.05 m. The boundary conditions were set to non-reflecting and the ${\hat{f}}_{1}$ and ${\hat{g}}_{1}$ acoustic forcing state from the CFD initialized by forcing with ${\hat{f}}_{1, b}$ , ${\hat{h}}_{1, b}$ and ${\hat{g}}_{2, b}$ waves (Figure 1), where the subscript $b$ stands for boundary. Because the Helmholtz equation does not account for any mean flow effects, entropic waves incident at the inlet were forced directly within the mathematical flame formulation via a harmonic source term. No additional entropy transport equation is solved in the present Helmholtz FEM formulation. Entropy disturbances (if present in the validation case) enter the model only through the prescribed input ${\hat{h}}_{1, b}$ in the analytical autoignition flame response law.

Figure 1.

1D DNS configuration schematic of the modeled duct in the FE-solver. The boundaries are forced with ${\hat{f}}_{1, b}$ , ${\hat{h}}_{1, b}$ and ${\hat{g}}_{2, b}$ waves. DNS: direct numerical simulation; FE: finite element; 1D: one-dimensional.

Autoignition flame RANS computation

The 2D RANS CFD computations are performed for a BFS Figure 2 with a total length of 0.3 m. The BFS is located at 0.15 m, with a symmetric area expansion from 0.01 m to 0.02 m. The investigated operating point burns hydrogen and air fully premixed at lean and autoignitive conditions (mass fractions: $H_{2}$ (0.0053), $O_{2}$ (0.1842), $H_{2} O$ (0.0521), $N_{2}$ (0.7584)). The inlet velocity measures 150 m/s, the mean operating pressure is 12.9 bar, and the inlet temperature is 1143 K. The simulation is done using ANSYS Fluent 2024 R1³⁵ with an in-house combustion model which is based on the work of Kulkarni et al.^28,29 In the combustion model, the radical buildup during the induction phase is represented by a normalized progress variable (PV). The PV and its source term is obtained by tracking representative intermediate and product species (here $H O_{2}$ and $H_{2} O$ ) in zero-dimensional reactors. For normalization, the sum of the tracked species mass fractions is then divided by its sum at equilibrium. It can be shown, that for fully premixed conditions, the PV source term only depends on the temperature, pressure and the current ignition progress.³⁶

{\dot{ω}}_{PV} = {\dot{ω}}_{PV} (T, p, PV)

(1)

Figure 2.

Schematic representation of the backward facing step (BFS) geometry.

Assuming, that the temperature does not change significantly during the radical build up, it can be divided into two parts: The initial temperature $T_{0, init}$ and the fluctuating temperature $T^{'}$ due to the acoustic forcing. Similarly, the pressure is divided into the operating pressure $p_{op}$ and the fluctuating part $p^{'}$ .

{\dot{ω}}_{PV} = {\dot{ω}}_{PV} (T_{0, init}, T^{'}, p_{op}, p^{'}, PV)

(2)

By assuming isentropic correlation of pressure and temperature for acoustic waves, $T^{'}$ can be expressed using $p^{'}$ , $p_{op}$ and $T_{0, init}$ :

T^{'} = T_{0, init} {(\frac{p^{'}}{p_{op}} + 1)}^{\frac{γ - 1}{γ}} - T_{0, init}

(3)

As $p_{op}$ and $T_{0, init}$ are constant and known $a$ $p r i o r i$ , the PV source term is defined by:

{\dot{ω}}_{PV} = {\dot{ω}}_{PV} (p^{'}, PV)

(4)

For the turbulence chemistry interaction a presumed beta PDF³⁷ approach is used. The mean turbulent PV source term is then obtained by folding the source term over a probability distribution $P$ (PV) which is defined by the PV’s mean and variance:

{\bar{\dot{ω}}}_{PV} = \int_{0}^{1} {\dot{ω}}_{PV} (p^{'}, PV) P (PV) d PV

(5)

The reaction rate is then computed by multiplying the Eddy Dissipation Model’s³⁸ reaction rate with the value of the PV in order to delay it until the residence time of a fluid particle is equal to the ignition delay time.

The advantage of this approach is that the correct ignition delay time can be achieved without transporting all the intermediate species and computing their reaction rates. Instead, the PV source terms are computed $a$ $p r i o r i$ using a detailed mechanism (SanDiego³⁹) and stored in a lookup table.

For the CFD simulations the realizable $k - ϵ$ turbulence model and second order temporal and spatial discretization is used with a timestep of 6 $\times 10^{- 6}$ s to obtain an acoustic CFL of 0.96. The structured uniform mesh consists of hexahedral elements with a size of 0.5 mm. A mesh independence study has been made comparing element sizes of 0.2 and 0.125 mm. The RANS mean field is validated by comparison to data from a compressible LES computation of the same geometry and at the same operating point. The OpenFoam code⁴⁰ is used to solve the Navier–Stokes equations. The partially stirred reactor turbulent combustion model is employed. The LES simulation is set up identically to the shorter BFS shown in Gopalakrishnan et al.³³

Methodology

Flame response and segmentation

The flame response to planar acoustic waves is characterized differently for the propagation-stabilized shear layer flames and the autoignition flame. For the shear layer flame, an $n - τ - σ$ model is implemented, similarly to Franke et al.,⁴¹ using Equation 6. The FTF thereof ( ${FTF}_{p s}$ Equation 7) relates the HRR perturbations to the velocity fluctuations at the dump plane, and is derived from the initial $n - τ$ model from Crocco.²¹ The magnitude $n$ has been chosen equal to 1 to satisfy the low frequency FTF limit for fully premixed propagation-stabilized flames (gain equal to 1).⁴² The second exponential term $\exp (- σ_{p s}^{2} τ_{p s}^{2} / 2)$ is a measure for the delay spread across the flame and acts as a low-pass filter at higher frequencies. They are calculated by fitting Gaussian kernels to the time-averaged HRR field and evaluating the location and spread of ignition lengths.

\frac{{\hat{\dot{q}}}_{a, p s}}{{\dot{q}}_{0, a, p s}} = \exp (i ω τ_{p s}) e x p (- σ_{p s}^{2} τ_{p s}^{2} / 2) \frac{{\hat{u}}_{x} (x_{ref, y})}{{\bar{u}}_{0}}

(6)

{FTF}_{p s} = \frac{\frac{{\hat{\dot{q}}}_{a, p s}}{{\dot{q}}_{0, a, p s}}}{\frac{{\hat{u}}_{x} (x_{ref, y})}{{\bar{u}}_{0}}}

(7)

For the autoignition flame, the framework from Gopalakrishnan et al.³² is used to quantify the flame response. The total autoignition HRR response of the flame is assumed to be a linear superposition of the individual flame responses (FTFs) to the acoustic and entropic waves.³¹ Equation 8 represents the HRR FTFs $F_{i} (ω)$ and Equation 9 characterizes the flame movement using the $G_{i} (ω)$ FTFs. ${\hat{\dot{q}}}_{a, a i} / {\dot{q}}_{0, a, a i}$ is the normalized HRR perturbation, ${\hat{x}}_{i g} / x_{i g, 0}$ the normalized flame movement, ${\dot{q}}_{0, a, a i}$ (W/m $^{2}$ ) is the mean HRR per unit area and $x_{i g, 0}$ the mean ignition length.³² ${\hat{f}}_{1}$ is the downstream traveling acoustic wave, ${\hat{g}}_{1}$ the upstream traveling acoustic wave and ${\hat{h}}_{1}$ the with the flow transported entropic wave.³² Typically, during thermoacoustic stability analyses using Helmholtz solvers the interaction between acoustic resonance modes and the flame is investigated. In this case, the resulting fluctuations on the flame are solely of acoustic nature and $F_{3}$ as well as $G_{3}$ are set to zero.

\frac{{\hat{\dot{q}}}_{a, a i}}{{\dot{q}}_{0, a, a i}} = F_{1} (ω) \frac{{\hat{f}}_{1}}{p_{0}} + F_{2} (ω) \frac{{\hat{g}}_{1}}{p_{0}} + F_{3} (ω) \frac{{\hat{h}}_{1}}{p_{0}}

(8)

\frac{{\hat{x}}_{i g}}{x_{i g, 0}} = G_{1} (ω) \frac{{\hat{f}}_{1}}{p_{0}} + G_{2} (ω) \frac{{\hat{g}}_{1}}{p_{0}} + G_{3} (ω) \frac{{\hat{h}}_{1}}{p_{0}}

(9)

The local acoustic field is evaluated from the FEM solution. For the propagation-stabilized flame regions, the local velocity fluctuation is used together with the prescribed ${FTF}_{p s}$ flame model to compute the corresponding HRR fluctuation. For the autoignition-stabilized HRR regions, the Riemann invariants are used. They can be computed from the acoustic quantities using the following relation:³⁴

{\hat{f}}_{1} = ({\hat{p}}_{1} + ρ_{0} c_{0} {\hat{u}}_{1}) / 2

(10)

{\hat{g}}_{1} = ({\hat{p}}_{1} - ρ_{0} c_{0} {\hat{u}}_{1}) / 2

(11)

The HRR fluctuation of the autoignition region is then computed by the multiplication with the transfer functions

F_{1}

F_{2}

, and

F_{3}

, as given by Equation 8.

FEM implementation

To characterize the flame in the FEM simulation and solve the inhomogeneous Helmholtz equation (Helmholtz equation with source terms), the numerical segmentation mechanism derived by Heinzmann et al.^17,18 is used. A main benefit of the segmentation is, that simulations can be done for non-compact flames, that is, flames that are under the influence of transverse eigenmodes. Figure 3 schematically visualizes the segmentation mechanism. Differentiation is made between propagation-stabilized subflames (blue circles) and the autoignition subflames (magenta circles). Due to the fine segmentation, it can be assumed that for each subflame a single value for the pressure and velocity perturbation is present. Thus, the flame also reacts with a single value for the HRR and the flame movement. For each of the subflames, specific FTFs can be implemented to capture the different flame responses of propagation-stabilized and autoignition flames. Hence, the segmentation mechanism allows different flame types in the same FEM computation.

Figure 3.

Numerical segmentation method for a perturbed partly autoignition and propagation-stabilized flame in two-dimensional (2D). The blue circles indicate the propagation-stabilized flame parts, and the magenta circles represent the autoignition flame.

For the propagation-stabilized flame parts (Equation 12) as well as the autoignition flame regions (Equation 13), Gaussian kernels are used to analytically position the flame in the computational domain. $σ_{\dot{q}, p s}$ and $σ_{\dot{q}, a i}$ are measures for the HRR width, $x_{f, p s} (y)$ and $x_{i g, 0} (y)$ the mean flame positions.

{\dot{q}}_{p s} (x, y) = \frac{{\dot{q}}_{0, a, p s} (y)}{σ_{\dot{q}, p s} (y) \sqrt{2 π}} \exp (- \frac{1}{2} {(\frac{x - x_{f, p s} (y)}{σ_{\dot{q}, p s} (y)})}^{2})

(12)

{\dot{q}}_{a i} (x, y) = \frac{{\dot{q}}_{0, a, a i} (y)}{σ_{\dot{q}, a i} (y) \sqrt{2 π}} \exp (- \frac{1}{2} {(\frac{x - x_{i g, 0} (y)}{σ_{\dot{q}, a i} (y)})}^{2})

(13)

The parameters for the HRR width, the mean flame position and the integrated HRR can be calculated from either CFD mean field solutions or experimental chemiluminescence imaging if available. Using these quantities and the assumption of the Gaussian kernel in x-direction (multiple studies used Gaussian distributions^{17,30,31,34,43}), the average flame can be rebuilt fully in FEM. Whilst the procedure here is done in 2D on the $x y$ -plane, the previous study²⁷ derived the identical procedure in three-dimensional (3D). Hence, the new method shown here can also be used in 3D.

At this point, the average flame is fully characterized in the FEM simulation and what remains to be characterized is the representation of the fluctuating HRR perturbations. This is achieved by obtaining the instantaneous HRR after multiplying the mean HRR by the FTF Equation 8 (for the propagation-stabilized flame Equation 6, respectively) and accounting for the flame motion in the Gaussian distribution (shift of the center of the Gaussian kernel) using Equation 9 to obtain Equations 14. For the propagation-stabilized flame region, no flame motion is introduced due to the anchored nature of the flame.

\begin{aligned} {\hat{\dot{q}}}_{a, p s} (y) = & {\dot{q}}_{0, a, p s} (y) {FTF}_{p s} \frac{{\hat{u}}_{x} (x_{ref}, y)}{{\bar{u}}_{0}} \\ {\hat{\dot{q}}}_{a, a i} (y) = & {\dot{q}}_{0, a, a i} (y) (F_{1} (ω) \frac{{\hat{f}}_{1}}{p_{0}} + F_{2} (ω) \frac{{\hat{g}}_{1}}{p_{0}} + F_{3} (ω) \frac{{\hat{h}}_{1}}{p_{0}}) \end{aligned}

(14)

Using Gaussian kernel distributions for both flames the instantaneous volumetric HRR is obtained (Equation 15).

\begin{aligned} {\dot{q}}_{(} x, y) = \frac{{\dot{q}}_{0, a, p s} (y)}{σ_{\dot{q}, p s} (y) \sqrt{2 π}} \exp (- \frac{1}{2} {(\frac{x - x_{f, p s} (y)}{σ_{\dot{q}} (y)})}^{2}) \dots \\ (\frac{{\hat{\dot{q}}}_{a, p s} (y)}{{\dot{q}}_{0, a, p s} (y)} + 1) + \frac{{\dot{q}}_{0, a, a i} (y)}{σ_{\dot{q}, a i} (y) \sqrt{2 π}} \dots \\ \exp (- \frac{1}{2} {(\frac{x - x_{i g, 0} (y) - {\hat{x}}_{i g} (y)}{σ_{\dot{q}, a i} (y)})}^{2}) (\frac{{\hat{\dot{q}}}_{a, a i} (y)}{{\dot{q}}_{0, a, a i} (y)} + 1) \end{aligned}

(15)

The fluctuating HRR $\hat{\dot{q}} (x, y)$ , which is the source term for the Helmholtz equation, is consequently computed by subtracting the mean HRR from the instantaneous counterpart (Equation 16).

\hat{\dot{q}} (x, y) = \dot{q} (x, y) - ({\dot{q}}_{p s} (x, y) + {\dot{q}}_{a i} (x, y))

(16)

The equation for the inhomogeneous Helmholtz equation, which is solved in FEM (in the frequency domain), thus becomes Equation 17 by adding the source term to the right-hand side of the homogeneous Helmholtz equation. From the computations the eigenfrequencies and linear growth rates are extracted. In this study COMSOL 6.2 is used.

\nabla \cdot (\frac{1}{ρ_{0}} \nabla \hat{p}) + \frac{ω^{2}}{γ p_{0}} \hat{p} = - i ω \frac{γ - 1}{γ p_{0}} \hat{\dot{q}} (x, y)

(17)

Results

Results obtained by implementing the proposed methodology are shown in the following section.

1D DNS forced

Figure 4 shows the intermediate result obtained in comparison with the DNS simulation to validate the FEM implementation of a flame solely stabilized by autoignition. An excellent match is achieved in comparison with the DNS result in both phase and magnitude. The small deviation in phase could result from neglecting mean flow effects by employing the Helmholtz equation.

Figure 4.

Comparison of the HRR fluctuations for the 1D DNS configuration forced at 100 Hz. The solid line represents the DNS data and the dashed line the FEM result. FEM: finite element method; 1D: one-dimensional; DNS: direct numerical simulation.

Mean fields of the BFS burning hydrogen

Figure 5 shows the time-averaged mean fields from the RANS simulation. Figure 5(a) shows the $y$ -velocity. The recirculation zone is visible by the blue and red areas. Figure 5(b) displays the $x$ -velocity. The recirculation zone can again be nicely observed by the negative velocity areas in blue. The temperature distribution is given in Figure 5(c) and matches well with the results from CANTERA simulations. Lastly, the pressure is shown in Figure 5(d), where a pressure drop downstream of the area jump is observed and a pressure increase downstream of the flame.

Figure 5.

Mean field solutions of the Reynolds averaged Navier–Stokes (RANS) simulation: (a) $y$ -velocity, (b) $x$ -velocity, (c) temperature and (d) pressure.

The time-averaged flame shape is computed using the in-house autoignition reheat flame model described in the CFD section. To validate the mean field qualitatively, comparisons to center-plane data from the described LES simulation are made. The comparison between the results from the LES and the RANS simulations shows that the model can replicate the flame to a certain extent. The mean ignition-length and overall flame shape is captured well (length scales are small). The RANS simulation however predicts the propagation-stabilized shear layer flames to be further upstream compared to the validation data. Experience shows, that the occurrence of an earlier timed ignition is linked to the RANS solver used in combination with the reheat flame model.

Further comparisons can be made by analyzing the volumetric HRR of Figure 6 for two different sections in Figure 7. The three dash-dotted blue lines show the normalized volumetric HRR along a section ( $y = 0.001$ m) across the shear layer flame. The diamonds represent the LES data, the circles the RANS data and the stars the FEM data, respectively. As was already indicated before, the HRR peak of the shear layer flame occurs further upstream in the RANS simulation at 0.158 m compared to 0.169 m. The FEM value matches the one of the RANS closely with 0.157 m. The peak value is over-predicted by the RANS compared to the LES simulation. The solid magenta lines provide the same comparison for a section at a $y$ -coordinate of 0.001 m, close to the symmetry line. In this region the flame is stabilized by autoignition. The peak of the RANS simulation is at an $x$ -coordinate of 0.188 m. It is predicted further downstream compared to the LES at around $\sim$ 0.18 m. The mean ignition length of the FEM matches the RANS closely with 0.186 m. The HRR distribution is spread more symmetrically around the peak compared to the RANS, which is a prerequisite of the introduced mathematical model (Equation 15). More generally, the flame appears more compact in the LES compared to the RANS computation. Due to the match of the overall flame shape, the RANS mean field can still be used for the subsequent study. Simulations with harmonic acoustic forcing at the inlet with a pressure amplitude of 0.5% are performed using ANSYS Fluent. The values of the $n - τ - σ$ model are computed to be $n = 1$ (low frequency limit for fully premixed propagation stabilized flames⁴²), $σ = 0.099$ ms and $τ = 0.12$ ms using the method described in the methodology section.

Figure 6.

Comparison of the mean HRR between LES (top) RANS (middle) and FEM (bottom). The backward facing step is located at $x = 0.15$ m. FEM: finite element method; LES: large eddy simulation; RANS: Reynolds averaged Navier–Stokes; HRR: heat release rate.

Figure 7.

Comparison of the volumetric HRR of different sections through the shear layer flame (blue dash-dotted lines.) and the autoignition stabilized part (solid lines). The comparison is made between the LES (diamonds), RANS (circles) and the FEM (stars). RANS: Reynolds averaged Navier–Stokes; FEM: finite element method; LES: large eddy simulation.

The forced FEM simulation is set up identically to the RANS computation by forcing harmonically at the inlet with a pressure amplitude of 0.5%. The outlet is set to fully anechoic to mitigate reflections and upstream traveling acoustic $g$ -waves.

Figure 8 shows the FTFs $({\dot{q}}^{'} / {\dot{q}}_{0}) / (u^{'} / u_{0})$ of both the RANS simulations and the FEM simulations. The 1D transfer functions assigned to the different flame regions are the described analytical and numerical models. The quantity denoted here as the global FEM FTF is reconstructed a posteriori from the forced FEM solution and is shown for validation against the forced RANS result. The FEM computation (circles) is able to capture the phase very accurately. The gain of the FTF shows slight differences between the CFD simulation (diamonds) and the FEM computation. The overall agreement of the FTF is good, indicating that the newly introduced segmentation model can combine reheat flame regions with different dynamic responses in a single FEM computation. Further, the comparison to the dashed blue line, which represents the $n - τ - σ$ model, indicates that the investigated flame response is dominated by the propagation-stabilized flames in the shear layer. Consequently, the proposed methodology can be used in a further step for stability calculations of longitudinal eigenmodes.

Figure 8.

FTF Comparisons between the CFD results (diamonds) and the FEM computation (circles). The blue dashed line represents the used $n - τ - σ$ model. FTF: flame transfer function; CFD: computational fluid dynamics; FEM: finite element method.

For completeness, Figure 9 shows the used FTFs (with respect to the velocity fluctuations $({\dot{q}}^{'} / {\dot{q}}_{0}) / (u^{'} / u_{0})$ ) for a larger frequency range. At higher frequencies, the assumption of no mean flow velocity breaks down in the FEM. The low-frequency forced BFS comparison is available only up to 500 Hz. The extended curves shown in Figure 9 therefore represent model-based trends of the prescribed component flame-response models and should not be interpreted as direct BFS CFD validation over the full frequency range. However, the FTFs can still be analyzed to understand the flame response at higher frequencies. Looking at the $n - τ - σ$ FTF (blue dashed line) of the shear layer flame zones, the low-pass filter behavior can clearly be observed. A negligible HRR response contribution from the shear layer flames from 4000 Hz onwards is indicated. The FTF from the solely autoignition stabilized central HRR zone is computed using the Lagrangian framework³² and is shown by the magenta dash-dotted line. The FTF indicates that at higher frequencies the autoignition flame zone contributes to a large extent to the overall flame response. However, no low-pass filter is contained in the model and thus the gain is estimated to be slightly lower in reality. The Lagrangian framework³² was validated with 1D DNS simulations up to 4000 Hz, where deviations of +20% were observed in the gain between DNS and the model. The black solid line represents a superposition of both the shear layer flame FTF and the autoignition flame FTF with an overall autoignition weighting of 0.32 and a shear layer flame weighting of 0.68, respectively. This gives an indication of how the flame behaves globally at higher frequencies, given that the assumption of a superposition of the HRR response of a propagation-stabilized flame region and an autoignition flame region is valid.

Figure 9.

FTF Comparisons between the $n - τ - σ$ model (blue dashed) and the autoignition FTF (magenta dash-dotted line). The solid black line is the combination of both FTFs. The circles represent the FTF from the FEM, and the diamonds the FTF from the RANS simulations. FTF: flame transfer function; RANS: Reynolds averaged Navier–Stokes; FEM: finite element method.

Conclusion

With reheat combustors becoming of increasing interest in industrial gas turbines, autoignition flames have shifted into the scope of current research. Especially in industrial reheat combustion chambers, an autoignition flame will likely have certain flame regions that are propagation-stabilized due to either flow recirculation or vortex breakdown. So far, research on thermoacoustics, especially the flame dynamics, has been performed separately for propagation-stabilized flames and autoignition-stabilized flames. In this paper, we propose an analytical model to predict the response of a composite flame partly stabilized by propagation and partly by autoignition. We leverage a previously developed segmentation method to incorporate both flame types in the same FEM computation. Using the segmentation mechanism, different FTFs are prescribed to the propagation- and autoignition-stabilized flame regions. For the propagation-stabilized HRR regions, a classical $n - τ - σ$ model is used. For the autoignition-stabilized HRR regions, the FTF obtained from a Lagrangian particle tracking based approach is used. Results showed a very good fit in comparison to forced CFD data. Consequently, the framework can be used to compute longitudinal eigenmode stability predictions of partly autoignition and propagation-stabilized flames and can aid the development of novel emissions-free industrial gas turbine combustors.

Footnotes

ORCID iDs

Simon M Heinzmann

Nino A Leder

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: The research work presented in this manuscript was carried out within the FLEX4H2 project. The FLEX4H2 project is supported by the Clean Hydrogen Partnership and its members European Union, Hydrogen Europe and Hydrogen Europe Research (GA 101101427), and the Swiss Federal Department of Economic Affairs, Education and Research, State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or any other granting authority. Neither of them is liable for any use that may be made of the information contained therein.

Declaration of Conflicting Interests

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

Appendix Notation

References

United Nations. Paris agreement, 2015.

IEA. World energy outlook 2022, 2022.

Bundesamt für Energie. Energieperspektiven 2050+, 2021.

Expert Group Security of Supply ETH Zurich. Energy security in a net zero emissions future for switzerland, 2023.

Ciani

Bothien

Bunkute

, et al. Superior fuel and operational flexibility of sequential combustion in ansaldo energia gas turbines. J Glob Power Propuls Soc 2019; 3: 630–638.

Bothien

Ciani

Wood

, et al. Toward decarbonized power generation with gas turbines by using sequential combustion for burning hydrogen. J Eng Gas Turbines Power 2019; 141: 121013.

Farhat

Salvini

. Novel gas turbine challenges to support the clean energy transition. Energies 2022; 15: 5474.

Goldmeer

Buck

Suleiman

, et al. Techno-economic analysis of hydrogen and ammonia as low carbon fuels for power generation. In: Proceedings of ASME Turbo Expo 2023, 2023.

Kloess

Zach

. Bulk electricity storage technologies for load-leveling operation - an economic assessment for the austrian and german power market. Int J Elect Power Energy Syst 2014; 59: 111–122.

10.

Bothien

Pennell

Zajadatz

, et al. On key features of the aev burner engine implementation for operational flexibility. Proceedings of ASME Turbo Expo 2013.

11.

Gruber

Bothien

Ciani

, et al. Direct numerical simulation of hydrogen combustion at auto-ignitive conditions: Ignition, stability and turbulent reaction-front velocity. Combust Flame 2021; 229: 111385.

12.

Casel

Ghani

. Analysis of the flame dynamics in methane/hydrogen fuel blends at elevated pressures. Proc Combust Inst 2023; 39: 4631–4640.

13.

Gruber

Meyer

OHH

Heggset

, et al. Numerical investigation of reheat hydrogen flames in the sequential-combustion stage of a heavy-duty gas turbine. In: Proceedings of ASME Turbo Expo 2024 turbomachinery technical conference and exposition GT2024.

14.

Nicoud

Benoit

Sensiau

, et al. Acoustic modes in combustors with complex impedances and multidimensional active flames. AIAA J 2007; 45: 426–441.

15.

Laera

Schuller

Prieur

, et al. Flame describing function analysis of spinning and standing modes in an annular combustor and comparison with experiments. Combust Flame 2017; 184: 136–152.

16.

Campa

Camporeale

. Eigenmode analysis of the thermoacoustic combustion instabilities using a hybrid technique based on the finite element method and the transfer matrix method. Adv Appl Acoust 2012; 1: 1–14.

17.

Heinzmann

Gopalakrishnan

Gant

, et al. Development and validation of a framework to predict the linear stability of transverse thermoacoustic modes of a reheat combustor. Combust Flame 2025; 275: 114010.

18.

Heinzmann

Gopalakrishnan

Wood

, et al. Linear stability analysis of transversal thermoacoustic modes in reheat combustors. Proceedings of GT2025 ASME Turbo Expo 2025,GT2025-152879, 2025.

19.

Silva

Nicoud

Schuller

, et al. Combining a helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor. Combust Flame 2013; 160: 1743–1754.

20.

Silva

Magri

Runte

, et al. Uncertainty quantification of growth rates of thermoacoustic instability by an adjoint helmholtz solver. J Eng Gas Turbine Power 2017; 139. DOI: 10.1115/1.4034203

21.

Crocco

. Aspects of combustion stability in liquid propellant rocket motors part i: Fundamentals. Low frequency instability with monopropellants. J Am Rocket Soc 1951; 21: 163–178.

22.

Schuermans

. Modeling and control of thermoacoustic instabilities. PhD Thesis, EPFL, 2003.

23.

Schuermans

Bellucci

Guethe

, et al. A detailed analysis of thermoacoustic interaction mechanisms in a turbulent premixed flame. In: Proceedings of ASME Turbo Expo 2004 Power for Land, Sea, and Air.

24.

Noiray

Bothien

Schuermans

. Analytical and numerical analysis of staging concepts in annular gas turbines, 2010, 2010.

25.

Laurent

Bauerheim

Poinsot

, et al. A novel modal expansion method for low-order modeling of thermoacoustic instabilities in complex geometries. Combust Flame 2019; 206: 206.

26.

Bothien

Noiray

Schuermans

. A novel damping device for broadband attenuation of low-frequency combustion pulsations in gas turbines. J Eng Gas Turbines Power 2014; 136: 041504.

27.

Heinzmann

Gopalakrishnan

Wood

, et al. Linear stability analysis of transversal thermoacoustic modes in reheat combustors. J Eng Gas Turbines Power 2025; 147: 121019.

28.

Kulkarni

. Large Eddy Simulation of Autoignition in Turbulent Flows. PhD Thesis, Technische Universität München, 2013.

29.

Kulkarni

Bunkute

Biagioli

, et al. Large eddy simulation of alstom’s reheat combustor using tabulated chemistry and stochastic fields-combustion model. DOI: 10.1115/GT2014-26053.

30.

Zellhuber

Bellucci

Schuermans

, et al. Modelling the impact of acoustic pressure waves on auto-ignition flame dynamics. In: Proceedings of the European combustion meeting 2011.

31.

Gant

Gruber

Bothien

. Development and validation study of a 1D analytical model for the response of reheat flames to entropy waves. Combust Flame 2020; 222: 305–316.

32.

Gopalakrishnan

Gruber

Moeck

. Response of auto-ignition-stabilized flames to one-dimensional disturbances: Intrinsic response. J Eng Gas Turbines Power 2021; 143: 121011.

33.

Gopalakrishnan

Heggset

Gruber

, et al. Prediction of autoignition-stabilized flame dynamics in a backward-facing step reheat combustor configuration. In: Symposium on thermoacoustics in combustion: Industry meets Academia (SoTiC 2023), 2023.

34.

Gopalakrishnan

Gruber

Moeck

. Computation and prediction of intrinsic thermoacoustic oscillations associated with autoignition fronts. Combust Flame 2023; 254: 112844.

35.

ANSYS. Ansys fluent user’s guide, 2024. http://www.ansys.com.

36.

Brandt

. Beschreibung der Selbstzündung in turbulenter Strömung unter Einbeziehung ternärer Mischvorgänge. PhD Thesis, Technische Universität München, 2005. https://mediatum.ub.tum.de/601964.

37.

deBruynKops

Riley

. Large-eddy simulation of a reacting scalar mixing layer with arrhenius chemistry. Comput Math Appl 2003; 46: 547–569.

38.

Magnussen

Hjertager

. On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion. Symp (Int) Combus 1977; 16: 719–729.

39.

Mechanical and Aerospace Engineering (Combustion Research). Chemical-kinetic mechanisms for combustion applications. http://combustion.ucsd.edu, Accessed 2025. San Diego Mechanism web page, University of California at San Diego.

40.

Weller

Tabor

Jasak

, et al. A tensorial approach to computational continuum mechanics using object orientated techniques. Comput Phys 1998; 12: 620–631.

41.

Franke

Bothien

Sattelmayer

. Flame transfer function measurements of hydrogen-enriched reheat flames. Combust Flame 2025; 273: 113949.

42.

Polifke

Lawn

. On the low-frequency limit of flame transfer functions. Combust Flame 2007; 151: 437–451.

43.

Rosenkranz

Sattelmayer

. Experimental investigation of high frequency flame response on injector coupling in a perfectly premixed multi-jet combustor. J Eng Gas Turbines Power 2023; 145: 111022.