Analysis of the Influence of the Mean Flow Field Mach Number on Thermo-Acoustic Combustion Instability

This work deals with numerical simulations of an annular burner with cross-section variations by using the FEM method. The main novelty is the introduction of the axial component of the mean flow in 3D simulations, which is usually neglected in the literature. The geometry investigated consists of an annular plenum chamber and an annular combustion chamber connected by 12 swirl burners. New insights are provided for a wide range of mean flow values, thus of Mach numbers. Firstly, a numerical assessment has been performed. Indeed, a simple annular geometry, which accurately approximates a real LP combustor, has been simulated and the results compared with the ones obtained with the software OSCILOS. The results achieved agree with the ones presented by Dowling and Stow [1], which demonstrated by applying a loworder one-dimensional model that the higher the mean flow velocity, the lower the frequencies and the growth rates gained.


Introduction
Modern gas turbines for power generation equipped with lean premixed (LP) dry low NOx combustion systems are susceptible to thermo-acoustic instability. These phenomena are characterized by self-sustained pressure oscillations, which arise from the coupling mechanism of unsteady heat source produced by the flame and the acoustic waves of the system. Over the years, annular configuration have been largely analyses in literature [1][2][3][4][5]. Despite the large numbers of theoretical, numerical and experimental studies presented on these configurations, very few are the works in which numerical simulations are developed with the introduction of the mean flow for real experimental test rigs [1,3,[6][7][8][9]18]. The present paper aims to fill this gap by developing a numerical procedure relying on the combination of a Helmholtz solver and a linear flame law with the introduction of the mean flow. In this framework a sensitivity study of the thermo-acoustic instability by varying the Mach number (M), which represents the velocity set at the burner inlet, has been carried out. In general, the mean flow is not considered as the Mach number of the flow entering the combustion chamber is usually low. Usually the mean flow is small compared to that of the sound that is generated in the combustion chamber for gas turbines (in general the Mach number does not exceed 0.2 [1]), but at the same time it is not fully correct to neglect it. Mean flow has two main consequences. First of all, it affects the propagation speed of acoustic waves, with the traveling downstream wave at a speed c + u and upstream at a speed c − u. Furthermore, it affects the existence of entropy and vorticity waves. Moreover, the contribution of the mean flow is of fundamental importance to describe the behavior of compact elements such as a burner. The introduction of the velocity in a computational model, able to solve the convective wave equation allows to evaluate with greater accuracy the effect of the fluid dynamics and acoustic losses on the propagation of the waves within a compact element. At this point, it is worth briefly reviewing some paramount investigations of the influence of mean flow in longitudinal devices. Dowling and Stow in [1] used a low-order model to demonstrate the damping of frequencies as the number of Mach increases. Besides this, Dowling and Stow, demonstrated that the introduction of mean flow on simple geometries, is not significant for Mach number less than about 0.2, without mentioning its influence for real test rig applications. The mean flow would introduce the possibility of a new mode of oscillation at much lower frequency in which the time taken for the convection of entropic waves or hot spots sets the period of fluctuation [10]. Campa in [6] introduced the mean flow in a methodology based on FEM in order to analyze the instability of thermo-acoustic combustion of gas turbines by solving a convective Helmholtz equation on simple longitudinal geometries with constant cross section. In this work, a simple geometry, which consists of an annular combustion plenum connected by 12 swirl burners and an annular combustion chamber has been analyzed by means of FEM technique. This geometry has been modeled in COMSOL Multiphysics ® . Simulations have been carried out with the introduction of a flame sheet in the first part of combustion chamber. The main novelty lies in the application of different levels of mean flow magnitude. To do this, the velocities deriving from 7 different Mach numbers have been imposed in the plenum according to the continuity equation. The results show a good stability of the analysis method despite of the complexity of the numerical model, showing a decrease of both the frequency and growth rate, with the increase of the mean flow magnitude.

GOVERNING EQUATION
The acoustic analysis is based on the resolution of the wave equation. To derive the wave equation, a compressible and viscous fluid is considered; in absence of external force, the continuity and momentum equations can be expressed, respectively, as: where ρ is the density, p is the pressure, u is the velocity vector, σ i,j is the tensor of viscous stress and e i is the versor in direction of coordinate i. D/Dt represents the material derivative operator defined as ∂/∂t + u · ∇. Considering the fluid as a perfect gas, the equation of state p = RρT can be applied, where T represents the temperature and R = c p − c v the perfect gas constant, equal to the difference between specific heat at constant pressure and volume, respectively. Calling the specific energy as e = c v T and the enthalpy h = c p T = e + p/ρ, the energy equation can be written as: whereq represents the heat supplied to the fluid per volume unit and k is the thermal conductivity. Exploiting Eq.2, the Eq.3 becomes: Replacing, the entropy expression S , defined by dh = T dS + (1/ρ)d p, Eq.4 becomes: Under the assumption of non-viscous flow σ i,j = 0, and fluid as ideal gas γ = it is possible to decompose the flow in one component of mean flow, constant and uniform, identified by overbar, and a small perturbation, identified by apex. Therefore, it is possible to write the flow variables as: Combining the continuity, the momentum and the energy equations for the fluctuation fields, the non-homogeneous wave equation reads as follows (equation derived from equation 8 in [1]): For analysis in the frequency domain, the fluctuation variables are expressed in the form: By dividing equation Eq.7 by density and using Eq.8, writing the eigenvalue −iω as λ, the wave equation is expressed as: where Q CM represents the heat monopole source (that represents a domain heat source causing pressure variation) that the user has to introduce in the Helmholtz solver. In this study, the mean flow is characterized by the solely axial velocity vector component, which corresponds to the x direction (u x ). Thus, in absence of the heat release fluctuations, the monopole source added on the right side of the wave Eq.9 is equal to: where Q u accounts for the mean flux and Q ht for the heat release.

Flame model
Considering as reference the flame model proposed by Dowling and Stow [1], the release of non-stationary heat, coming from the flame, can be assumed concentrated on an axial plane x = b called flame sheet, Figure 1. Based on this assumption, the heat fluctuation q (x) is related to the air velocity u i , calculated before the cross section of the flame, coming from the burner, with a time delay τ: where Q (t) is the heat fluctuation for unit area of cross section of duct and the subscript i refers to the conditions before the area of heat release. The Eq.11 relates the heat fluctuations per unit of volume, q (x), to the heat fluctuations per unit area of cross section, Q (t), by means of the Dirac delta δ(x − b). The factor β provides a measure of the intensity of the heat release. The sensitivity study of the model to the factor β underlined that it has little influence on the Growth Rate and it causes a decrease in the value of frequency on vibrating axial modes, but much less on vibrating circumferential modes [1]. τ is defined as the convective time from the injection of fuel to its combustion, i.e., the time taken by the fuel to reach the flame front. Time delay greatly influences the stability of the system, therefore it is possible to identify the value of τ corresponding to a stable condition for the individual vibrating modes of the system. In the FEM eigenvalue analysis, it is assumed that the heat fluctuations occur in a thin volume of thickness s, and the Dirac delta δ(x − b) can be approximated as follows: Combining Eq.11 and Eq.12, one can write: Introducing the eigenvalue λ = −iω, Eq.14 can be written in the frequency domain as follows: Its simplified version is:q

HELMHOLTZ SOLVER: COMSOL Multiphysics ®
The software COMSOL Multiphysics ® is one of the most advanced and complete finite element software on the market. It is a software environment for modeling and simulating several physical systems. Indeed, the program features a series of modules capable of simulating physical, magnetic, gravitational, fluid dynamic phenomena, etc.. In the present work, the frequency domain module has been used. This allows to detect complex self-frequencies that characterize the system and, therefore, to perform a stability study. It consists in solving a problem with differential equations, converted into a problem with eigenvalues in the frequency domain.

LOW ORDER CODE: OSCILOS
OSCILOS is a low-order simulator developed at the Imperial College London by Prof Aimee Morgans and Dr. Dong Yang. It is written in Matlab ® /Simulink ® . OSCILOS permits to draw a combustor as a network of connected modules and to model the acoustic waves whether in 1-D or 2-D plane. OSCILOS is a user friendly low-order code and it returns reliable results by considering simple geometries. The mean flow is calculated by assuming 1-D flow condition with changes only across flames or module interface. It is necessary to connect the thermal properties and the mean flow between the neighboring combustor sections. In OSCILOS, between neighboring combustor sections, the mass and the energy flux are unchanged while the momentum flux is increased by the axial force on the walls [1,11]. Therefore, the firstorder energy equation can be written in function of ratio of sectional surface areas [11].

NUMERICAL ASSESSMENT
First of all, a simple cylinder with a constant cross section has been investigated in order to validate the numerical set up in COMSOL Multiphysics ® , see Figure 2. The geometry reproduces a well-known case study presented in the literature [1,6]. The geometric parameters are shown in Table 1.  Table 2. The flame sheet has been positioned at 1/10 of the cylinder length; moreover, the monopole source has been set only in the domain of the flame sheet. The velocity has been imposed in order to set the mean flow in the respective domains (upstream and downstream of the flame sheet) considering them constant along the streamwise direction. To discretize the model an unstructured grid of 18218 elements has been used. For this model, the inlet section has been considered as an open wall (zero acoustic pressure), while the exit section has been considered as a closed wall (u = 0). The  Figure 3.
The results obtained for the frequency and the growth rate are shown in the Figure 4, where the results of the FEM analysis (COMSOL) are compared with those obtained by means of OSCILOS. As it can be seen, the results obtained with the two solvers are in good agreement. The main difference between the Helmholtz solver (COMSOL Multiphysics ® ) and the low order solver (OSCILOS) is that COMSOL allows the user to discretize much   On the other hand, the main discrepancies are found for the GR. Indeed, the gap between the curves is 540% at M = 0.2 and GR changes its sing at M = 0.2, where the mode becomes stable.
Also, in the case of stable mode (GR>0), in accordance with the results obtained by low order approach (OSCILOS) the frequency (normalized by f max =294.713 Hz) and the growth rate (normalized by |GR max |=1.0588) decrease as the number of Mach in the plenum increases, Figure 5. Differently from the previous case, the frequency and the GR curves show roughly the same gap (=78%). In this case both COMSOL and OSCILOS show a frequency reduction of 7.6% from M = 0 to M = 0.2. The discrepancies between the frequencies and the growth rate calculate by using COMSOL and OSCILOS can be due to the fact that in COMSOL a constant density has been considered in the flame sheet.

NUMERICAL SETUP OF ANNULAR COMBUSTOR
Once the numerical assessment, a simplified model of an annular combustor has been considered. This model is characterized by a diffusion chamber (plenum) and a combustion chamber, both of them being annular connected to each other thanks to 12 swirled burners. A mixture of air and fuel is injected in the plenum from the outside. The numerical model does not take into account what is present upstream and downstream of the combustor. The dimensions of the various components of the system have used having as reference the model used by Pankiewitz and Sattelmayer in their study [2]. In the numerical model used in COM-SOL Multiphysics ® , see Figure 6a, the burners have been reproduced in a simplified way (ducts with constant cross section area) and the nozzles, located at the end of the combustion chamber, have been omitted. This leads to a simplified geometry characterized by an annular plenum, 12 thin ducts, representing the burners, and an annular combustion chamber.   The plenum inlet and the combustion chamber outlet behave like walls acoustically closed (u = 0). In the plenum, and therefore in the burners, there is the mixture of air and fuel at a controlled temperature of T p = 774 K. The absolute pressure is assumed constant and equal to p 0 = 101325 Pa. The combustion chamber is modeled as an annular duct at the temperature of T CC = 2350 K. In the first part of the combustion chamber, a sheet of 0.01 mm has been introduced to model the flame as proposed by Dowling and Stow [1]. The density of the flow inside both the plenum and the burners is equal to 0.45 kg/m 3 , whereas, in the combustion chamber the density is equal to 0.15 kg/m 3 . The speed of sound inside both the plenum and the burners is equal to 556 m/s, whereas in the combustion chamber c 2 = 945.2 m/s. A grid independence study has been performed. Three unstructured grids have been generated, see Figures 6b,7a and 7b. Table 3 sums up their number of cells and the values of frequency and growth rate obtained for the third vibrating mode, in the case of Mach number equal to zero. We note the independence of the frequencies from the number of elements of the grid. This is explained by the fact that the wavelengths involved are sufficiently large compared to the cell size used, thus making its influence negligible. It is worth to highlight that a grid with less than 68601 numbers of elements, is not able to discretize correctly the flame sheet zone. Finally, the unstructured grid, made of 68601 tetrahedron elements, has been chosen for the further simulations.
The mean flow in the plenum is calculated by the Mach number (Mach number which varies from 0 to 0.125 with a range equal to 0.025) and the speed of sound imposed in the plenum. Pankiewits and Sattelmayer in [2] consider a model with a Mach number equal to 0.017, in this preliminary work a wide range of values has been investigated. In the combustion chamber, the mean flow is derived by applying the continuity equation, u 2 = (u p ρ p A p )/(ρ 2 A 2 ). The velocity in the 12 burners is considered equal to the velocity in the plenum; its variation in the passage between the plenum and the swirler is neglected in this preliminary no zero Mach number study. In the future, corrected values of the flow velocity will be taken into account for each domain.

RESULTS OF THE ANNULAR COMBUSTOR
The influence of the Mach number on the third stable mode is shown in Figure 8. The Mach number in the plenum has been changed from 0 to 0.125. This range of Mach number and the thermodynamic properties of the flow are representative of the flow field at the inlet of real gas burner applications (i.e., u=80 m/s). The results of the sensitivity analysis on the frequency and the growth rate of the third stable mode, by varying the Mach number in the plenum, are shown in Figure 8 and summed up in Table 4.  Table 4, the frequency is not affected by great variations (at least 1.4%), on the contrary the GR undergoes a remarkable variation. As shown in Figure 8, the frequency (Blue line) and the growth rate (Red line) decrease while the Mach number increases, according to Dowling and Stow [3]. This means that in this case the third stable mode becomes more stable by introducing the mean flow.

CONCLUSIONS
In this work, the FEM method has been applied with the aim to investigate the influence of different mean flow magnitudes on the frequency and GR of the system of an annular combustor. For the first time the mean flow is introduced in a 3D model. Firstly, a numerical assessment has been performed by studying a simple cylindrical burner by means of the low order code OSCILOS and the Helmholtz solver in COMSOL. In these models the flame model is defined by applying the flame sheet in the first part of the combustion chamber. The results of both a stable and unstable modes have been compared. Given the good results a real annular combustor has been investigated. Increasing the Mach number in the plenum, and therefore the flow velocity at the inlet of the domain, both the frequency and the relative Growth Rate, decrease. This means that the mean flow leads to a decrease of the burner instability. The mean flow introduces energy dissipation in the burner and, for this reason, stabilizes the oscillation mode of the combustor. In the future, numerical simulations will be carried out to perform comparison between results on the annular combustor obtained with COMSOL Multipysics ® and OSCILOS, to take into account likely discontinuities in the continuity, momentum and energy equations in the passage from one section to another.