Thermo-acoustic Instabilities in the PRECCINSTA Combustor Investigated Using a Compressible LES-pdf Approach

Abstract

This work predicts the evolution of self-excited thermo-acoustic instabilities in a gas turbine model combustor using large eddy simulation. The applied flow solver is fully compressible and comprises a transported sub-grid probability density function approach in conjunction with the Eulerian stochastic fields method. An unstable operating condition in the PRECCINSTA test case—known to exhibit strong flame oscillations driven by thermo-acoustic instabilities—is the chosen target configuration. Good results are obtained in a comparison of time-averaged flow statistics against available measurement data. The flame’s self-excited oscillatory behaviour is successfully captured without any external forcing. Power spectral density analysis of the oscillation reveals a dominant thermo-acoustic mode at a frequency of 300 Hz; providing remarkable agreement with previous experimental observations. Moreover, the predicted limit-cycle amplitude is found to closely match its respective measured value obtained from experiments with rigid metal combustion chamber side walls. Finally, a phase-resolved study of the oscillation cycle is carried out leading to a detailed description of the physical mechanisms that sustain the closed feedback loop.

1 Introduction

Thermo-acoustic instabilities are a phenomenon often encountered during the development stages of modern low emission aero-engine or stationary gas turbine combustors designed to operate under lean, partially premixed conditions. Characterised by strong oscillations of the dynamic pressure and heat release rate, these instabilities can arise when a resonant feedback loop between the combustion chamber acoustics and the reactive flow dynamics is established (Lefebvre 1999; Lieuwen and Yang 2005). The time delay between an incident acoustic perturbation and the subsequent flame response must therefore be in phase to fulfil the so-called Rayleigh criterion (Rayleigh 1878). This will lead to an increase of the oscillation amplitude until a limit cycle is reached, where the amount of acoustic energy added to the system equals its energy losses. When operating in these limit cycle conditions, major problems such as higher sound and pollutant emissions, increased heat fluxes to the walls, flame extinction and flashback, or in severe cases even component failures may occur.

Advanced computational methods—essentially large eddy simulation (LES)—are potentially capable of facilitating the prevention of thermo-acoustic instabilities in the early design cycles of a combustion system. LES methods for reactive flows are based on the solution of spatially filtered transport equations governing mass, momentum, energy and species mass conservation. These equations must be closed through additional turbulence and combustion models accounting for the sub-grid scales (sgs). In inert flows away from solid boundaries, results are known to be relatively insensitive to the models used to represent the sgs stresses. Moreover, sgs contributions are generally expected to be small if a sufficiently fine domain discretisation is achieved. When dealing with self-excited thermo-acoustic phenomena, a fully compressible LES formulation is needed in order to allow for the propagation of acoustic waves represented by temporal fluctuations in the pressure field. Potential coupling mechanisms between the flame, flow and acoustics inside a combustion chamber can thereby be directly accounted for. However, such self-excited LES are often challenging to carry out and require careful consideration of all geometrical parts making up the domain, as well as all acoustic impedances at the inlets and outlets.

While initial modelling efforts of combustion instabilities have mainly focused on the simulation of rocket motors (Baum and Levine 1982) and ramjets (Kailasanath et al. 1991), more recent works have also targeted swirl-stabilised gas turbine combustors, see e.g. Huang and Yang (2009). Detailed reviews of numerical prediction and control methods for combustion instabilities can be found in Refs. Dowling and Morgans (2005) and Poinsot (2017). Only a very limited number of computational works in the open literature have utilised compressible LES to study the effects of self-excited thermo-acoustic instabilities in gas turbines, specifically. The majority of these studies (e.g. Refs. Roux et al. 2005; Schmitt et al. 2007; Staffelbach et al. 2009; Franzelli et al. 2012; Ghani et al. 2016; Schulz et al. 2019) were carried out with the widely-known AVBP solver, which couples LES with a thickened flame model (Colin et al. 2000) to represent sgs turbulence-chemistry interactions. Tachibana et al. (2015) and Chen et al. (2019) each employed different LES approaches with tabulated chemistry, where as Lourier et al. (2017) applied finite-rate chemistry in a scale adaptive simulation that transitions between LES and an unsteady Reynolds-averaged Navier–Stokes (URANS) solution depending on grid resolution.

The main advantage of the transported sgs probability density function (pdf) method utilised in the current work is its burning regime independent formulation; making the approach equally applicable to non-premixed, partially premixed and perfectly premixed turbulent flames (Jones and Prasad 2010; Jones et al. 2012; Bulat et al. 2013). This sgs-pdf method is implemented into an in-house LES code, which was extended for the application to fully compressible flows, and applied to examine the unstable behaviour of the lab-scale PRECCINSTA gas turbine model combustor. The combined LES-pdf methodology is presented below followed by an overview of the target test case and its numerical set-up. The method’s capabilities—especially in terms of predicting the combustor’s self-excited thermo-acoustic coupling—are assessed against available measurement data. A selection of relevant results are discussed including time-averaged flow statistics, dynamic pressure and heat release rate signals as well as a phase-resolved analysis of the observed oscillation cycle.

2 LES-pdf Approach

2.1 Filtering Operation

The evolution of any given compressible, turbulent reacting flow can be described via the governing equations for mass, momentum, energy and species mass conservation. In LES, a spatial filter is applied to these equations with the aim of filtering out fine-scale sub-filter motions. The spatial filter of a function $f=f(\mathbf {x},t)$ is defined as its convolution with a filter function G, according to:

$$\begin{aligned} \overline{f}(\mathbf {x},t) = \int _\Omega G(\mathbf {x} -{\mathbf {x'}} ; \Delta (\mathbf {x})) f({\mathbf {x'}},t) d {\mathbf {x'}} \end{aligned}$$

(1)

where ${\mathbf {x}}$ and t are the spatial coordinate and time, respectively. The filter function must be positive definite in order to maintain filtered values of scalars such as mass fraction within bound values and to preserve the nature of the source terms associated with chemical reaction (a filter that changes sign may change consumption terms to formation terms). The integration is defined over the entire flow domain $\Omega$ and the condition that the filter kernel G should be positive definite implies that it has the properties of a pdf. The filter function has a characteristic width $\Delta$ which, in general, may vary with position—if it varies smoothly then the filtering operation (1) commutes with spatial differentiation. In the present work, a box or ‘top hat’ filter defined by:

$$\begin{aligned} G(\mathbf {x} -\mathbf {x'})&= \frac{1}{\Delta ^3} \quad \text{ for }\quad |\mathbf {x} -\mathbf {x'}|<\frac{\Delta }{2} \\&= 0 \quad \text{ otherwise } \end{aligned}$$

is used. The density variations in the unresolved scales that arise in combusting flows can be treated through the use of density-weighted, or Favre, filtering, defined by $\widetilde{f} (\mathbf {x},t) = \overline{\rho f}/ \overline{\rho }$, where ${\rho }$ denotes the density.

2.2 Filtered Equations of Fluid Motion

The governing equations of fluid motion are obtained by applying a spatial, density weighted filter to the continuity and Navier–Stokes equations:

$$\begin{aligned}&\frac{\partial \overline{\rho }}{\partial t } + \frac{\partial \overline{\rho }\widetilde{u}_i }{\partial x_i } = 0 \end{aligned}$$

(2)

$$\begin{aligned} \frac{\partial \overline{\rho }\widetilde{u}_i}{\partial t } + \frac{\partial \overline{\rho }\widetilde{u}_i\widetilde{u}_j }{\partial x_j } = -\frac{\partial \overline{p} }{\partial x_i } + \frac{\partial }{\partial x_j } \left[ \mu \left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} \right) - \frac{2}{3} \mu \frac{\partial \overline{u}_k}{\partial x_k} \delta _{ij} \right] - \frac{\partial \tau _{ij}^{sgs}}{\partial x_j } \end{aligned}$$

(3)

where ${u_i}$ represents the velocity in the i-direction, p is the pressure, ${\mu }$ is the viscosity and ${\delta _{ij}}$ is the Kronecker delta. The sgs stress tensor ${\tau _{ij}^{sgs}}$ is closed through a dynamic version of the well-known Smagorinsky turbulence model (Piomelli and Liu 1995), introducing an sgs viscosity ${\mu _{sgs}}$ to mimic the diffusion process on the dissipative small-scales. The isotropic parts of the viscous and sgs stresses are absorbed into the pressure.

2.3 Scalar Transport Equations

A system’s instantaneous chemical state can be expressed in terms of transport equations for, respectively, the specific mole number ${n_{\alpha }}$ of each of the $N_s$ species ${\alpha }$ present in the mixture and the enthalpy h. The species diffusive fluxes are represented by Fick’s law together with the assumption of equal diffusivities and in the enthalpy equation the heat flux is represented by Fourier’s law with a unity Lewis number assumption (Poinsot and Veynante 2005). The resulting equations are:

$$\begin{aligned}&\frac{\partial \rho n_{\alpha }}{\partial t} + \frac{\partial \rho u_i n_{\alpha }}{\partial x_i} = \frac{\partial }{\partial x_i} \left( \frac{\mu }{\sigma } \frac{\partial n_\alpha }{\partial x_i} \right) + \rho \dot{\omega }_\alpha ; \quad 1 \le \alpha \le N_s \end{aligned}$$

(4)

$$\begin{aligned}&\frac{\partial \rho h}{\partial t} + \frac{\partial \rho u_i h}{\partial x_i} = \frac{\partial }{\partial x_i} \left( \frac{\mu }{\sigma } \frac{\partial h}{\partial x_i} \right) +\frac{\partial p}{\partial t} + u_i \frac{\partial p}{\partial x_i} \end{aligned}$$

(5)

The above species and enthalpy equations can be conveniently combined and written in the following generic form:

$$\begin{aligned} \frac{\partial \rho \phi _\alpha }{\partial t} + \frac{\partial \rho u_i \phi _{\alpha }}{\partial x_i} = \frac{\partial }{\partial x_i} \left[ \left( \frac{\mu }{\sigma } \right) \frac{\partial \phi _{\alpha }}{\partial x_i} \right] + \rho \dot{\omega }_{\alpha }; \quad 1 \le \alpha \le N_s +1 \end{aligned}$$

(6)

For chemical species, the last term on the right-hand side of Eq. (6) represents the chemical source term, while for enthalpy, it includes the pressure gradient terms and radiative heat loss (neglected here). However, the filtered forms of Eqs. (4) and (5) are not used in the present work because of the difficulties encountered in evaluating the filtered values of the chemical source terms in Eq. (4).

2.4 Transported pdf Approach

A joint sgs pdf, ${\widetilde{P}_{sgs}}$, for the $N_s$ species and enthalpy is introduced to account for the influence of the sub-grid scale contributions on the species formation rates. It represents the filtered product of the fine-grained one-point pdf of each scalar as follows:

$$\begin{aligned} \widetilde{P}_{sgs}(\underline{\psi }; \mathbf{x} , t) = \frac{1}{\overline{\rho }} \int _{\Omega } \rho G(\mathbf{x} - \mathbf{x} ^{\prime }, \Delta (\mathbf {x})) \times \prod _{\alpha =1}^{N_s+1} P_\alpha = \prod _{\alpha =1}^{N_s+1} \delta \left[ \psi _{\alpha } - \phi _{\alpha }(\mathbf{x} ^{\prime },t) \right] \text {d}{} \mathbf{x} ^{\prime } \end{aligned}$$

(7)

where ${\psi _{\alpha }}$ is the sample space of the scalar ${\phi _{\alpha }}$ and ${\delta }$ is the Dirac-function. Following Gao and O’Brien (1993) and taking into account the approach discussed in Brauner et al. (2016), an exact equation describing the evolution of ${\widetilde{P}_{sgs}}$ can be derived:

$$\begin{aligned} \begin{aligned}&\overline{\rho } \frac{\partial \widetilde{P}_{sgs}(\underline{\psi })}{\partial t} + \overline{\rho } \widetilde{u}_i \frac{\partial \widetilde{P}_{sgs}(\underline{\psi })}{\partial x_i} + \sum _{\alpha =1}^{N_s+1} \sum _{\beta =1}^{N_s+1} \frac{\mu }{\sigma } \frac{\partial \widetilde{\phi }_\alpha }{\partial x_i} \frac{\partial \widetilde{\phi }_\beta }{\partial x_i} \frac{\partial ^2\widetilde{P}_{sgs}(\underline{\psi })}{\partial \psi _\alpha \partial \psi _\beta } \\&\quad = \,- \underbrace{ \sum _{\alpha =1}^{N_s+1} \frac{\partial }{\partial \psi _\alpha } \left[ \overline{\rho }\dot{\omega }_\alpha (\underline{\psi })\widetilde{P}_{sgs}(\underline{\psi })\right] }_\mathbf{I: \text {source term (closed)}} + \underbrace{ \frac{\partial }{\partial x_i}\left[ \left( \frac{\mu }{\sigma }+\frac{\mu _{sgs}}{\sigma _{sgs}}\right) \frac{\partial \widetilde{P}_{sgs}(\underline{\psi })}{\partial x_i}\right] }_\mathbf{II: sgs \; pdf \; \text {transport (gradient closure)}} \\&\qquad - \underbrace{ \frac{C_d}{2 \tau _{sgs}} \sum _{\alpha =1}^{N_s+1} \frac{\partial }{\partial \psi _\alpha } \left[ \overline{\rho }(\psi _\alpha - \widetilde{\phi }_\alpha (\mathbf{x} ,t))\widetilde{P}_{sgs}(\underline{\psi })\right] }_\mathbf{III: sgs \; \text {micro-mixing (LMSE closure)}} \end{aligned} \end{aligned}$$

(8)

where the sgs micro-mixing constant is ${C_d = 2.0}$ and the sgs micro-mixing time scale is given by ${\tau _{sgs} = \overline{\rho }\Delta ^2 /\mu _{sgs}}$. In this formulation, the source term ${\dot{\omega }_\alpha }$ (I) appears in closed form. The second and third term on the right-hand side of the equation represent sgs transport of the pdf (II) and sgs micro-mixing (III). These terms are approximated by a gradient closure analogous to the Smagorinsky model and the linear mean square estimation (LMSE) closure Dopazo and O’Brien (1974), respectively.

2.5 Stochastic Fields Method

The fully Eulerian stochastic fields method is employed as a solution method to the sub-grid pdf evolution equation (8) as proposed in Valiño (1998) and Sabel’nikov and Soulard (2005). For this purpose, ${\widetilde{P}_{sgs}(\underline{\psi })}$ is represented by an ensemble of N stochastic fields ${\xi _{\alpha }^n(\mathbf{x} ,t)}$ assuming that ${1 \leqslant n \leqslant N}$ and ${1 \leqslant \alpha \leqslant N_s+1}$:

$$\begin{aligned} \begin{aligned} \bar{\rho } \text {d}\xi _{\alpha }^n + \bar{\rho }\widetilde{u}_i \frac{\partial \xi _{\alpha }^n}{\partial x_i} \text {d}t =\,&\overline{\rho }\dot{\omega }_{\alpha }^n(\underline{\xi }^n) \text {d}t + \frac{\partial }{\partial x_i} \left[ \left( \frac{\mu }{\sigma }+\frac{\mu _{sgs}}{\sigma _{sgs}} \right) \frac{\partial \xi _{\alpha }^n}{\partial x_i}\right] \text {d}t \\&- \frac{C_d }{2 }\frac{\bar{\rho } }{\tau _{sgs}} \left( \xi _{\alpha }^n - \widetilde{\phi }_{\alpha } \right) \text {d}t + \left( 2 \bar{\rho } \frac{\mu _{sgs}}{\sigma _{sgs}} \right) ^{1/2} \frac{\partial \xi _{\alpha }^n}{\partial x_i} \text {d}W_i^n \end{aligned} \end{aligned}$$

(9)

where ${\text {d}W_i^n}$ denotes increments of a Wiener process, different for each stochastic field but independent of the spatial location $\mathbf{x }$. It is approximated by time-step increments defined as ${\eta _i^n (\text {d}t)^{1/2}}$ where ${\eta _i^n}$ is a ${\{-1,1\}}$ dichotomic random vector. Averaging over the stochastic fields results in the filtered value of each scalar.

3 PRECCINSTA Combustor

3.1 Test Case Review

The swirl-stabilised PRECCINSTA gas turbine model combustor was studied in an experimental test campaign (Meier et al. 2007) at the German Aerospace Center (DLR), where it was shown to exhibit flame oscillations driven by thermo-acoustic instabilities under certain operating conditions. A significant number of experimental investigations have since been carried out attempting to identify the underlying mechanisms involved in the flame oscillation. Among the different phenomena studied in these works are the dynamics of coherent vortex structures (Steinberg et al. 2010; Caux-Brisebois et al. 2014; Oberleithner et al. 2015; An et al. 2016; Zhang et al. 2019), the interaction between velocity and equivalence ratio fluctuations (Stöhr et al. 2017) and the transition between bistable flame shapes (Yin et al. 2017; Stöhr et al. 2018).

Most computational studies on the other hand, have focused on validating LES combustion models based on the reference case involving ‘stable’ operating conditions with a global equivalence ratio of $\phi _{glob}$ = 0.83 and thermal power of $P_{th}$ = 30 kW. Only very few investigations dealt with the oscillating flame behaviour of the PRECCINSTA combustor and its underlying mechanisms at play. Roux et al. (2005) utilised a compressible LES method in combination with acoustic analysis to compute both a non-reacting and reacting flow $\phi _{glob}$ = 0.75, $P_{th}$ = 27 kW. Comparison of their results to experimental data showed reasonable agreement for both the predicted acoustic and hydrodynamic mode found in the cold flow as well as the thermo-acoustic mode obtained with combustion. Their study also underlined the importance of resolving the radial swirler instead of starting all computations from the combustion chamber inlet plane. Franzelli et al. (2012) were the first to simulate the oscillating flame case $\phi _{glob}$ = 0.7, $P_{th}$ = 25.1 kW measured in Meier et al. (2007) and considered in the current work. Their compressible LES reproduced a self-excited frequency of 390 Hz compared to the experimental value of 290 Hz. The overestimation was argued to be the result of a falsely modelled impedance at the fuel injection inlets. It was furthermore concluded that resolving the fuel-air mixing—as opposed to prescribing a fully premixed mixture at the domain inlet—is crucial for the prediction of thermo-acoustic instabilities in this test case.

Finally, Lourier et al. (2017) performed a scale adaptive simulation (SAS) and additional experimental measurements of the oscillating flame case with adjusted acoustic boundary conditions. In these experiments, the inflow nozzle was equipped with a chocked orifice plate 15 cm upstream of the air plenum to create a well-defined acoustic inlet boundary. This however, led to a marginal decrease of the thermo-acoustic frequency by about 15 Hz. Furthermore, rigid metal combustion chamber side walls were securely mounted to mitigate acoustic damping effects observed in the initial measurement campaign due to the use of loosely fitted quartz glass windows. As a result, power spectral density (PSD) analysis of the air plenum and combustion chamber pressure signals yielded an increase in the oscillation amplitude of around 10 dB. The employed SAS showed very good agreement for both the thermo-acoustic frequency and its amplitude with a relatively small over-prediction of about 15–25 Hz and 5 dB, respectively. The simulation included the fuel supply plenum and channels according to the findings of Franzelli et al. (2012).

Fig. 1

Schematic of the PRECCINSTA combustor showing the two pressure probe locations

3.2 Computational Set-up

The developed flow solver (see also Appendix) is based on the in-house, pressure-based, second-order accurate finite volume method BOFFIN (BOundary Fitted Flow INtegrator). It utilises a fully compressible formulation in order to account for acoustic wave propagation and eight stochastic fields to describe the influence of the sub-grid scale contributions (Jones and Navarro-Martinez 2007; Noh et al. 2019). The solver’s default model parameters remain unchanged (Fredrich et al. 2019) and all chemical reactions are represented by an augmented reduced mechanism (ARM) (Sung et al. 2001) involving 15 steps and 19 species. A non-reflective outflow boundary condition based on the work of Rudy and Strikwerda (1980) is implemented to mitigate wave reflection at the combustor outlet:

$$\begin{aligned} \begin{aligned} \left( \frac{\partial p}{\partial t} + \rho c \frac{\partial u_i}{\partial t} \right) + \eta _p \frac{c \left( 1 - \mathcal {M}^2 \right) }{L_x} \left( p - p_\infty \right) + \beta _t \mathbb {T} = 0 \end{aligned} \end{aligned}$$

(10)

where c is the sound speed, $\mathcal {M}$ is the average Mach number in the outlet plane, $L_x$ is the characteristic size of the computational domain and $p_\infty$ represents the target far-field pressure equal to 1 bar. The pressure relaxation factor $\eta _p$ is a constant of 0.28 (Rudy and Strikwerda 1980), chosen to manipulate (i.e. minimise) the amplitudes of incoming waves. The transverse damping parameter $\beta _t$ is set equal to $\mathcal {M}$, as suggested by Yoo et al. (2005), and combined with the transverse term $\mathbb {T}$ defined by:

$$\begin{aligned} \begin{aligned} \mathbb {T} = u_t \frac{\partial p}{\partial x_t} + \gamma p \frac{\partial u_t}{\partial x_t} - \rho c u_t \frac{\partial u_i}{\partial x_t} \, , \ \ t = [j,k] \end{aligned} \end{aligned}$$

(11)

where $\gamma$ denotes the isentropic expansion factor. In the current formulation, all three velocity components plus the density are obtained on the boundary using zero-th order extrapolation from interior points. The pressure can then be calculated from the characteristic boundary condition and the temperature is obtained from the ideal gas law. A more detailed description of non-reflective boundary conditions for the simulation of compressible flows can also be found in Ref. Poinsot and Lele (1992).

The investigated operating condition $\phi _{glob}$ = 0.7, $P_{th}$ = 25.1 kW and detailed geometry (Fig. 1) of the PRECCINSTA combustor are identical to case (1) of the initial experiments of Meier et al. (2007). The velocities and density of all incoming air and fuel streams are fixed to prescribe the correct mass flow rates at a temperature and pressure of 320 K and 1 bar, respectively. The computational mesh consists of 2.7 million grid points with refinements in the reaction zone and fuel-air mixing regions. All wall boundaries use a logarithmic law-of-the-wall formulation (Hoffmann and Benocci 1995). Previous studies (Benard et al. 2019; Gövert et al. 2018; Fredrich et al. 2019) of the ‘stable’ flame case have shown that an application of non-adiabatic combustion chamber walls noticeably improves the near-wall scalar field and can induce flame lift-off in line with experimental observations. Hence, wall heat transfer is included in the current work via isothermal temperatures of 1400 K and 700 K prescribed at the combustion chamber side and base plate walls. No forcing of the inlet velocity is applied in this work to ensure that any potential oscillatory behaviour observed in the simulation is fully self-excited.

Fig. 2

Radial profiles of the mean axial and radial velocity from the LES (—) and experiments ($\blacklozenge$) at different downstream positions

Fig. 3

Radial profiles of the mean and RMS temperature from the LES (—) and experiments ($\blacklozenge$) at different downstream positions

4 Results and Discussion

4.1 Time-Averaged Flow Statistics

Radial profiles of the mean axial and radial velocity at different downstream positions are shown in Fig. 2. Overall, good agreement with the initial experiments of Meier et al. (2007) can be reported. Both the location and strength of the inner recirculation zone (IRZ) and outer recirculation zone (ORZ) are accurately reproduced. The incoming jet of fresh gases, represented by two pronounced velocity peaks, experiences a minor inwards shift with downstream position. This shift in combination with underestimated peak values in the radial velocity indicates a smaller mean spreading angle and expansion of the jet into the combustion chamber.

The accompanying under-prediction of the mean length and spreading angle of the flame can be determined from the temperature and methane (CH₄) mass fraction profiles displayed in Figs. 3 and 4. Nevertheless, the general temperature and CH₄ distributions in the IRZ and ORZ are retrieved, suggesting the applied wall heat transfer treatment is realistic, and the mean flame position and shape are qualitatively captured. A marginally increased mean flame lift-off height is observed close to the centreline, implying a longer duration of flame detachment (Hermeth et al. 2014) from the central cone bluff body compared to the measurements. Temperature fluctuations of up to 500 K, visualised as root mean square (RMS) profiles in Fig. 3, are expectedly high in both simulation and experiment due to the oscillating nature of the flame. However, LES results show an almost constant level of fluctuation inside of the inner shear layer (ISL), $-16$ mm $\leqslant$ r $\leqslant$ 16 mm, whereas experimental data reveals a dip of RMS values in between the inner and outer shear layer (OSL), 10 mm $\leqslant$ r $\leqslant$ 16 mm. The computed location of the incoming jet is therefore less steady over the full oscillation cycle. A stronger flame angle oscillation—or flapping of the flame—can be concluded, which likely contributes to the under-predicted mean spreading angle observed for the jet and flame. The increased flapping motion is potentially induced by swirl number fluctuations due to stronger acoustic perturbations (Candel et al. 2014). These are caused by a higher amplitude limit-cycle oscillation obtained in the simulation, as discussed in the following section.

Fig. 4

Radial profiles of the mean and RMS methane mass fraction from the LES (—) and experiments ($\blacklozenge$) at different downstream positions

4.2 Dynamic Pressure and Heat Release Rate Signals

Fig. 5

Dynamic pressure and integrated heat release rate (HRR) fluctuations and their respective power spectral density (PSD) signals in the air plenum and combustion chamber. Dominant peaks measured in the experiments are plotted as vertical dotted lines (frequencies only) for quartz glass walls (Meier et al. 2007) and as ${\varvec{\times }}$-markers for rigid metal walls (Lourier et al. 2017)

Time series of the dynamic pressure in the air plenum and combustion chamber as well as the integrated heat release rate fluctuations are displayed in Fig. 5. Both the heat release rate and combustion chamber pressure are in phase satisfying the Rayleigh criterion. Constant mean amplitudes in the pressure signal translate into a growth rate equal to zero, thus indicating a limit-cycle oscillation. The existing modes are revealed via PSD analysis providing a dominant peak for the air plenum, combustion chamber and heat release rate at approximately 300 Hz. This peak represents the thermo-acoustic frequency and is in remarkable agreement with experimental measurements of about 290 Hz (Meier et al. 2007). A second peak is clearly visible at around 600 Hz representing the first harmonic, also reported experimentally. However, the computed PSD amplitudes of the thermo-acoustic mode are overestimated by about 13 dB in both air plenum and combustion chamber. This overestimation is a result of acoustic damping (Lourier et al. 2017) occurring at the combustion chamber side walls in the initial experiments. Excellent agreement is obtained when comparing the present PSD amplitudes to the additional measurements carried out by Lourier et al. (2017) who used tightly fitted rigid metal walls to mitigate the effects of acoustic damping. Overall, the dynamic behaviour of the combustor is shown to be well captured as self-excited thermo-acoustic instabilities are successfully reproduced at the correct frequency and amplitude.

4.3 Phase-Resolved Analysis of the Oscillation Cycle

Fig. 6

Definition of eight different phase angles in relation to one full oscillation cycle of the integrated heat release rate (HRR) and dynamic pressure in the air plenum and combustion chamber

Fig. 7

Phase-resolved temporal evolution of the heat release rate (HRR), velocity magnitude and CH₄ mass fraction for one full oscillation cycle

While the quantitative velocity and scalar fields are influenced by the acoustic pressure amplitude in the combustor, the physical feedback mechanisms driving the thermo-acoustic coupling remain largely unaffected by this and can still be investigated. A phase-resolved analysis of the oscillation cycle is therefore carried out and eight different phase angles (Ph 1–8) are introduced as illustrated in Fig. 6. Starting from Ph 1, the evolution of the main flow properties at each phase is described below, based on Fig. 7.

At Ph 1, fresh gases enter the combustion chamber with a high mass flow rate and equivalence ratio, as characterised by the high velocities and CH₄ mass fraction concentrations in this region. Large amounts of unburnt reactants are thereby convected into the main flame zone enhancing combustion. The increasing heat release rate approaches its maximum between Ph 2 and 3 elevating the pressure inside the combustion chamber and generating a pressure wave that propagates up- and downstream within the domain. A negative pressure drop between the combustion chamber and air plenum develops decreasing the mixture mass flow rate and hence supply of fresh gases into the reaction zone. As a result, the bulk flow inside the swirler becomes almost stagnant. CH₄, on the other hand, is still injected at an almost constant rate since the fuel jets are only weakly affected by the pressure oscillation due to their high acoustic impedance. This leads to fuel accumulation inside the swirler during phases 4 and 5, raising the equivalence ratio of the local mixture. Minor flame flashback occurs as the heat release rapidly decreases followed by a reduction of the local pressure. A phase difference of about ${80}^{\circ }$ between the recorded pressure signals in the combustion chamber and air plenum is determined at this point, similar to the shift reported in the experiments Meier et al. (2007). This pressure drop in turn generates an increased mass flow rate within the swirler from Ph 6 and 7 while the integrated heat release rate reaches its minimum. The high equivalence ratio mixture finally exits the burner nozzle at Ph 8 pushing the flame front downstream and closing the feedback loop as the heat release rate begins to increase again. In summary, fluctuations of the mixture mass flow rate due to a periodically varying pressure drop between the air plenum and combustion chamber can be identified as the dominant driver of the oscillation. This mechanism is superimposed by equivalence ratio fluctuations (Hermeth et al. 2013; Lourier et al. 2017; Stöhr et al. 2017 )in the local mixture entering the reaction zone as a result of fuel accumulation inside the swirler at one specific period of the oscillation cycle.

5 Conclusions

Self-excited thermo-acoustic instabilities in the PRECCINSTA gas turbine model combustor were investigated using a novel, fully compressible LES-pdf approach comprising the Eulerian stochastic fields method. Comparisons with available measurement data revealed good agreement for the predicted velocity, temperature and CH₄ mass fraction radial profiles. Small deviations in the mean and RMS results were attributed to an increased oscillation of the flame spreading angle resulting from higher swirl number fluctuations. These are likely caused by an overestimated amplitude of the acoustic pressure oscillation compared to the original experiments, which incurred damping effects at the vibrating combustion chamber side walls.

Excellent agreement was achieved in reproducing the dynamic behaviour of the combustor. Recorded signals of the fluctuating heat release rate and pressure in the combustion chamber were in phase and showed growth rates equal to zero indicating a thermo-acoustic limit-cycle oscillation. PSD analysis of the pressure fluctuations revealed a dominant frequency of 300 Hz, identified as the thermo-acoustic mode, as well as its first harmonic at 600 Hz; closely matching the respective experimental values. The amplitude of these two modes was further found to be in very good agreement with measurements obtained using stiff combustion chamber side walls to suppress acoustic damping. Finally, a phase-resolved analysis of the oscillation cycle was carried out to identify its underlying coupling mechanisms. A periodically changing pressure drop between the air plenum and combustion chamber due to acoustic wave propagation led to a fluctuating mixture mass flow rate entering the reaction zone. The mixture further showed local equivalence ratio variations caused by fuel accumulation inside the swirler at one specific period of the cycle. The superposition of these two effects enhanced the unsteadiness of the heat release rate, which in turn generated acoustic wave propagation closing the feedback loop.

Appendix: Solution Algorithm

1.1 Preliminaries

If the spatial derivatives in the continuity, Navier–Stokes and stochastic fields equations are replaced by discrete approximations, then Eqs. (2), (3) and (9) can be written in the form:

$$\begin{aligned} \frac{\underline{\rho }^{n+1} - \underline{\rho }^{n}}{\delta t} + \nabla _d (\underline{\rho }\underline{u}) = 0 \end{aligned}$$

(12)

$$\begin{aligned} \frac{\underline{u}^{n+1}-\underline{u}^n}{\delta t}=\underline{A}_u(\underline{\rho },\underline{u}) \underline{u} + \underline{b}_u(\underline{u}) - \frac{1}{\rho } \nabla _d (\underline{p}) \end{aligned}$$

(13)

$$\begin{aligned} \frac{\underline{\phi }^{n+1}-\underline{\phi }^{n}}{\delta t}=\underline{A}_\phi (\underline{\rho },\underline{u}) \underline{\phi } + \underline{b}_\phi (\underline{\phi }); \ \ \ 1 \le \alpha \le N_s+1 \end{aligned}$$

(14)

where $\nabla _d(\ )$ is the discrete form of the divergence or gradient operator. If transporting the total enthalpy, i.e. $\underline{\phi }= h_t = h + 0.5u_iu_i$, then $\underline{b}_\phi = \frac{\partial p}{\partial t}$.

All the spatial terms in the continuity and Navier–Stokes equations are approximated by standard central differences. In the case of the stochastic fields equations, a total variation diminishing (TVD) scheme is used for the advection terms and the equations are solved using an approximately factored Euler–Maruyama scheme (Kloeden and Platen 1999), details of which are given in Jones et al. (2012). Hence, the Navier–Stokes equations and stochastic fields species equations are approximately ‘factored’ to yield:

$$\begin{aligned}&\underline{u}^{n+1}(\underline{I}-\underline{A}_u \delta t)= \underline{u}^{n}+\underline{b}_u \delta t - \frac{1}{\rho } \nabla _d(p^n) \delta t - (\underline{I}-\underline{A}_u \delta t) \frac{1}{\rho } \nabla _d(p^{n+1}-p^{n}) \delta t \end{aligned}$$

(15)

$$\begin{aligned}&\underline{\phi }^{n+1}(\underline{I}-\underline{A}_\phi \delta t)= \underline{\phi }^{n}+(\underline{I}-\underline{A}_\phi \delta t) \underline{b}_\phi \delta t; \ \ \ 1 \le \alpha \le N_{s} \end{aligned}$$

(16)

The solution to Eq. (15) yields:

$$\begin{aligned} \underline{u}^{n+1}+\mathcal {O}(\delta t^2) \times \nabla _d(p^{n+1}-p^{n}) \end{aligned}$$

(17)

and thus $\underline{u}^{n+1} + \mathcal {O}(\delta t^3)$, while Eq. (16) gives $\underline{\phi }^{n+1} + \mathcal {O}(\delta t^2)$.

1.2 Algorithm

Step (I): Stochastic fields equations.

Solve the discrete form of the stochastic fields species Eq. (16), using $\underline{\rho }^n$ and $\underline{u}^n$, as follows:

$$\begin{aligned} \begin{aligned} \underline{\hat{\phi }}^{n+1}=&\ \underline{\phi }^{n}+\underline{b}_\phi (\underline{\hat{\phi }}^{n+1})\delta t \\ (\underline{I} - \underline{A}_\phi \delta t)\underline{\phi }^{n+1} =&\ \underline{\hat{\phi }}^{n+1}; \ \ \ \ 1 \le \alpha \le N_{s} \end{aligned} \end{aligned}$$

(18)

Solve the stochastic fields equation for total enthalpy:

$$\begin{aligned} (\underline{I} - \underline{A}_\phi \delta t)\underline{h}_t^{n+1} = \underline{h}_t^n + \frac{\partial \underline{p}^n}{\partial t} \delta t \end{aligned}$$

(19)

The solution to Eqs. (18) and (19) can be used to obtain ${\hat{{\rm T}}}^{n+1}$ and $\underline{\rho }^{*} = \underline{\rho } ({\hat{{\rm T}}}^{n+1} ,\hat{ n}^{n+1}_\alpha , p^n)$ and $\underline{u}^{*}$ is set equal to $\underline{u}^n$.

Step (II): Navier–Stokes equations.

The approximate factorisation Eq. (15) is complemented as follows:

$$\begin{aligned} \frac{\underline{\hat{u}}^{n+1} - \underline{{u}}^{n}}{\delta t} = \underline{A}_u (\underline{\rho }^{*}, \underline{u}^{*}) \underline{\hat{u}}^{n+1/2} + \underline{b}_u (\underline{u}^{*}) - \frac{1}{\rho ^*}\nabla _d(p^n) \end{aligned}$$

(20)

The solution of Eq. (20) can then be used to obtain:

$$\begin{aligned} \begin{aligned} \frac{\underline{{u}}^{n+1} - \underline{\hat{u}}^{n+1}}{\delta t} = - \frac{1}{\rho ^*}\nabla _d( p^{n+1}-p^n) \end{aligned} \end{aligned}$$

(21)

and

$$\begin{aligned} \begin{aligned} \underline{u}^{n+1/2} =\underline{\hat{u}}^{n+1/2} - \frac{1}{2 \rho ^{*}}\nabla _d (p^{n+1} - p^n) \delta t \end{aligned} \end{aligned}$$

(22)

Step (III): continuity equation.

The continuity equation may be written as:

$$\begin{aligned} \frac{\rho ^{n+1} - \rho ^n}{\delta t} + \frac{\partial }{\partial x_i} \left( \rho ^{n+1/2} \underline{u}^{n+1/2} \right) = 0 \end{aligned}$$

(23)

If the pressure dependence of the density is taken into account, the densities in Eq. (23) can be obtained from:

$$\begin{aligned} \begin{aligned} \rho ^{n+1}&= \rho ^{*} + \frac{\partial \rho }{\partial p} (p^{n+1} - p^n) \\ \rho ^{n+1/2}&= \rho ^{*} + \frac{1}{2} \frac{\partial \rho }{\partial p} (p^{n+1} - p^n) \end{aligned} \end{aligned}$$

(24)

On substitution of Eqs. (22) and (24) into the continuity equation (23) and with the neglect of the second order terms $(p^{n+1}-p^n)\times (u^{n+1/2}- \hat{u}^{n+1/2})$ the following may be obtained:

$$\begin{aligned} \begin{aligned} \frac{\partial \rho }{\partial p} \frac{\underline{p}'}{\delta t} + \frac{1}{2} \frac{\partial }{\partial x_i} \left\{ \underline{\hat{u}}^{n+1/2} \frac{\partial \rho }{\partial p} \underline{p}'- \nabla _d (\underline{p}') \delta t \right\} = - \left[ \frac{\underline{\rho }^{*} - \underline{ \rho }^n}{\delta t} + \frac{\partial }{\partial x_i} \left( \rho ^{*} \hat{u}^{n+1/2}_i \right) \right] \end{aligned} \end{aligned}$$

(25)

where $p' \equiv p^{n+1} - p^n$.

On solution of Eq. (25) the pressure and density can be updated, viz: $p^{n+1} = p^n + p'$ and $\rho ^{n+1} = \rho ^{*} + \frac{\partial \rho }{\partial p}p'$. Finally set:

$$\begin{aligned} \begin{aligned} \underline{\rho }^{*}&= 1/2 (\rho ^{n+1} + \rho ^n) \\ \underline{u}^{*}&= 1/2 (\underline{u}^{n+1} + \underline{u}^n) \end{aligned} \end{aligned}$$

(26)

Repeat Steps (II) to (III) once.

The accuracy of $\underline{\rho }^{n+1}, \underline{u}^{n+1} \sim \mathcal {O} (\delta t^3)$, i.e. second order accuracy for the partial differential equations.