Prediction of combustion instability limit cycle oscillations by combining ﬂame describing function simulations with a thermoacoustic network model

Accurate prediction of limit cycle oscillations resulting from combustion instability has been a long-standing challenge. The present work uses a coupled approach to predict the limit cycle characteristics of a combustor, developed at Cambridge University, for which experimental data are available ( Balachandran, Ph.D. thesis, 2005 ). The combustor ﬂame is bluff-body stabilised, turbulent and partially-premixed. The coupled approach combines Large Eddy Simulation (LES) in order to characterise the weakly non-linear response of the ﬂame to acoustic perturbations (the Flame Describing Function (FDF)), with a low order thermoacoustic network model for capturing the acoustic wave behaviour. The LES utilises the open source Computational Fluid Dy- namics (CFD) toolbox, OpenFOAM , with a low Mach number approximation for the ﬂow-ﬁeld and combustion modelled using the PaSR (Partially Stirred Reactor) model with a global one-step chemical reaction mecha- nism for ethylene/air. LES has not previously been applied to this partially-premixed ﬂame, to our knowledge. Code validation against experimental data for unreacting and partially-premixed reacting ﬂows without and with inlet velocity perturbations conﬁrmed that both the qualitative ﬂame dynamics and the quantitative response of the heat release rate were captured with very reasonable accuracy. The LES was then used to obtain the full FDF at conditions corresponding to combustion instability, using harmonic velocity forcing across six frequencies and four forcing amplitudes. The low order thermoacoustic network modelling tool used was the open source OSCILOS (http://www.oscilos.com). Validation of its use for limit cycle prediction was performed for a well-documented experimental conﬁguration, for which both experimental FDF data and limit cycle data were available. The FDF data from the LES for the present test case was then imported into the OSCILOS geometry network and limit cycle oscillations of frequency 342 Hz and normalised velocity amplitude of 0.26 were predicted. These were in good agreement with the experimental values of 348 Hz and 0.21 respectively. This work thus conﬁrms that a coupled numerical prediction of limit cycle behaviour is possible using an entirely open source numerical framework.


Introduction
The development of modern gas turbines requires high combustion performance as well as low emissions. In order to reduce NO x emissions, it is necessary to operate under lean combustion conditions. However, a serious issue related to lean combustors is susceptibility to damaging combustion instabilities [1]. Combustion instabilities generally refer to sustained pressure oscillations in the combustion chamber, resulting from the coupling of the system acoustics and the unsteady heat release [2]. Much of the early release calculations. The response of unsteady heat release to perturbations is modelled via a flame model [3,9], while the acoustic waves are captured by a low order combustor model [2,[10][11][12][13][14] or a Helmholtz solver [15], exploiting the fact that the acoustic wave behaviour, being linear, is well captured by analytical or simple numerical methods. The present study belongs to the second method, combining a low order network model for the combustor with a flame model obtained via high-fidelity CFD simulations.
The traditional flame model is defined via a (linear) Flame Transfer Function (FTF) [4]. However, this is restricted to small perturbations and thus cannot be used to predict limit cycle oscillations [10,16] nor other non-linear effects such as instability triggering and mode switching. More recently, it has been shown that the flame transfer function concept can be extended to the non-linear regime via a nonlinear FTF, also known as a Flame Describing Function (FDF) [9,17,18], in the form of: where Q /Q is the normalised heat release rate fluctuation and u /ū the normalised inlet velocity perturbation impinging on the flame. The FDF F(ω, |u |) is generally expressed in the frequency domain as gain (amplitude) G(ω, |u |) and phase ϕ(ω, |u |) which are functions of both forcing frequency ω and amplitude |u |. This approach makes the assumption of weak non-linearity, i.e. the flame response to harmonic forcing is assumed to be primarily at the same frequency as the forcing, but with a gain and phase shift which depend upon the forcing amplitude as well as the forcing frequency.
Several experimental studies [19][20][21][22][23][24][25][26][27] have been performed to determine non-linear flame models for the analysis of combustion instability, not only at lab scale, but also for real gas turbine combustors operated at high pressure [26]. The data reveal that the non-linearity of the FDF is of central importance as it governs the mechanisms leading to amplitude saturation into limit cycles [9,20,[23][24][25]28] -the behaviour of the acoustic waves remain linear. Different mechanisms causing saturation have been explored in experiments, such as the interactions of flame front with coherent structures [20,28], attachment point dynamics [29], and flame quenching [24]. Experimental data have also shown that the phase change of the FDF with forcing amplitude can cause saturation, due to the fact that it changes the Rayleigh source term which drives combustion instability [9].
Experimental flame model measurements are currently preferred for complex systems. However, computational simulations, if sufficiently accurate and fast, would offer the benefits of allowing stability predictions prior to experimental realisation. The approach of incorporating non-linear flame models provided by high-fidelity CFD into low order combustor models are only recently beginning to be exploited. Although some URANS (Unsteady Reynolds-Averaged Navier-Stokes) studies give reasonable results at specific conditions [30][31][32], large eddy simulation (LES) is capable of capturing unsteady flow and flame structures and is now widely used to investigate turbulent combustion problems [1,33]. However, it has recently drawn increasing concerns in the context of simulating combustion instability directly [8]. In order to resolve both the acoustic waves and the flame dynamics for low-Mach number flames, significant computation costs are required due to the small time step associated with the speed of sound appearing in the CFL time step limit, and the large computational domain needed to capture the whole acoustic system. The decoupled method is thus advantageous as computations are performed only for a small domain within the combustor to capture the flame dynamics, and no CFL time step limit based on the speed of sound is required. Consequently, low-Mach number or "incompressible LES" [34][35][36][37] can be used to determine the FDF as the flame response is well known to be unaffected by compressibility effects [20,29,38]. The use of low-Mach number or incompressible LES to identify the FDF is explained and justified in detail in other recent references [34][35][36][37]. Based on these ideas, the present study uses a similar low-Mach number LES solver to study an acoustically forced partially-premixed flame in order to identify the full FDF. The obtained FDF is then implemented in a low order network model for the combustor in order to investigate combustion instability, and in particular its non-linear features such as limit cycle frequency and amplitude.
Low-order thermoacoustic network models, which combine linear analytical models for the acoustic wave behaviour with an appropriate flame model, have been used fairly extensively for combustion instability studies [11,12,[39][40][41]. The basic idea is that, acoustically, the combustion system can be represented as a network of connected modules, each with simple geometry, which correspond to various components of the system. By combining with a well resolved flame model (a linear FTF or non-linear FDF), the frequencies and growth rates of the thermoacoustic modes, stability boundaries, and potentially limit cycle amplitudes, etc., can be determined. It is thus a useful and computationally efficient tool for combustion instability studies. Such a low-order network solver, OSCILOS [14,42], has been developed in the authors' group. It is written using Matlab © /Simulink © . The present study will combine it with a flame model from high fidelity CFD in order to predict limit cycle oscillations under unstable conditions.
The target case of the present study is the bluff-body stabilised flame investigated experimentally by Balachandran et al. [28] at Cambridge University, where series of experimental data are available for both partially-premixed and fully-premixed conditions. The case has an "acoustically short" flame within the combustor which simplifies determination of the FDF and its coupling with low order network models. For the fully-premixed case, LES has been recently performed [34] and the obtained FDF agrees well with experimental data. However, unstable thermoacoustic behaviour and hence limit cycle oscillations are only observed in the partially-premixed case (possible flame flash back has limited conditions considered under fullypremixed conditions). The present study thus considers the partiallypremixed case in the experiments by Balachandran et al. [28]. The objectives of the present paper are: (1) to perform LES studies of the turbulent partially-premixed flame for the first time and to compare these with experimental data, using the open-source CFD toolbox, OpenFOAM [43]; (2) to determine the full FDF of the flame response using LES for the first time; (3) to validate the OSCILOS network modelling tool by comparing limit cycle predictions with available experimental data; (4) to perform combustion instability analysis by predicting the limit-cycle amplitude for the unstable Cambridge configuration, and to compare with experimental measurements.
The target Cambridge experiment configuration will be described in Section 2, followed by the numerical details of the LES simulations in Section 3. The validation of the LES code for both cold flow and reactive flow results will be presented in Section 4. The determination of the full FDF will be given in Section 5. A brief description of the low-order network modelling tool and its validation is presented in Section 6. The limit cycle predictions for the Cambridge configuration are given in Section 7. Conclusions are presented in the last section.

Experimental test case
The burner considered in the present study has the simple construction illustrated in Fig. 1(a). This is described in detail by Balachandran et al. [28], and has featured in earlier studies [44][45][46][47].
In the experiments, the system can be operated in externally-forced or self-excited modes -the latter only under partially-premixed conditions. For the thermoacoustic network modelling, the whole system will be considered with the details given in Section 7. For the FDF determination, only part of the system (downstream of the plenum in Fig. 1(a)  The present study focusses exclusively on the lesser studied partially-premixed experiments, for which a combustion instability limit cycle was characterised. The burner geometry simulated consists of two concentric cylindrical ducts; one of inner diameter 35 mm carrying air and the other of outer diameter 8 mm carrying ethylene fuel. The latter leads to a conical bluff body of diameter d = 25 mm with a 45°cone angle. The fuel is injected radially through 6 injection holes along the circumference of the central pipe, situated 55 mm upstream of the exit of the bluff body. This setup results in partially-premixed conditions at the flame. The flame is enclosed using a quartz cylinder of inner diameter 70 mm which avoids equivalence ratio (φ) variations due to surrounding air entrainment. Note that the length of the enclosure is l = 350 mm under self-excited conditions, while the shorter length of l = 80 mm is used in the externally forced cases. Experimental observations [28] confirm that the length of the enclosure has neglectable impact on the flame dynamics, although it has a large effect on the acoustic characteristics. In the present LES, externally forced cases are considered and the length l = 80 mm is applied for the enclosure.
In the experiments, it should be noted that the global equivalence ratio is φ = 0.55 for the externally forced cases, and φ = 0.61 for the self-excited cases. Different experimental measurements are available at the different equivalence ratios (with not all measurements available at a single equivalence ratio). This leads to the present LES being performed for both equivalence ratios of φ = 0.55 and φ = 0.61.
For the externally forced cases, the acoustic forcing was generated experimentally by two loudspeakers mounted diametrically opposite one another on the circumference of the plenum chamber 100 mm downstream of the plenum inlet. This introduces a velocity oscillation of air along the chambers. As the mass flow rate of fuel is constant, the air mass flow rate oscillation leads to equivalence ratio oscillations prior to the combustor inlet. The combined oscillations of equivalence ratio and air velocity result in unsteady flame dynamics and heat release. The forcing amplitude (A) and frequency (f) were varied independently in the experiments. The forcing amplitude was as high as 70% of the mean velocity for certain forcing frequencies, and the forcing frequencies used ranged from 20 Hz to 400 Hz in the experiments.
Based on the definition of the FDF in Eq. (1), both the oscillating signals of heat release rate and reference velocity are needed to determine the FDF. In the experiments, the heat release rate was measured with OH * and CH * chemiluminescence. It was also possible to obtain a heat release rate predicting from the phase-averaged FSD (Flame Surface Density) images obtained from PLIF (Planar Laser-Induced Fluorescence), with images appropriately revolved around the burner central axis [28]. The values of OH * (f)/ OH * and CH * (f)/ CH * were used as estimates of Q (f)/ Q . The reference velocity in the experiments was taken at the combustor inlet (the position with the bulk velocity V b in Fig. 1(b)) and determined from acoustic pressure measurements using the two-microphone method.

Numerical method for flame simulations
In the present work, large eddy simulations are performed using the CFD toolbox, OpenFOAM. Specifically, a modified version of the reactingFOAM solver is used -this has been applied in previous LES studies of turbulent combustion [48][49][50]. The reactive flow equations are the Favre-filtered Navier-Stokes equations of mass, momentum, species mass fraction and energy. Following ref. [1], the gas mixture is presumed to be ideal, linearly viscous, with Fourier heat conduction and Fickian diffusion. The laminar viscosity is modelled by Sutherland's law. As the energy equation is solved, heat loss effects can be accounted for.
To close the governing equations, turbulence modelling is required. The popular Smagorinsky LES subgrid scale model [51] is applied, with the turbulent viscosity calculated by: where the model constant C s is equal to 0.167, | S| is the strain rate magnitude of the resolved velocity defined as | S| = 2 S i j S i j , and is the filter cutoff width, i.e. the characteristic length scale of the subgrid scale eddies. Note that the symbol¯denotes the spatial filtering used in the LES and the symbol˜denotes density-weighted filtering, defined asψ = ρψ/ρ for an arbitrary variable ψ.
A well-known problem for Smagorinsky based LES models is that the modelled turbulent viscosity, μ t , is too high in the near wall regions [52]. To improve the model performance near the wall, the (2)) should be damped by using the model for van Driest damping. In OpenFOAM, the damping is derived by changing the filter width, depending on the distance from the wall. Here, the filter width is calculated according to: where m is the cubic root of the cell volume, κ = 0.4187 is the von Karman constant, C = 0.158, A + = 26, y w represents the distance to the wall, and y + describes the dimensionless distance to the wall calculated from the wall shear stress. The van Driest damping has been used in previous studies [53][54][55] and improved results have been observed in calculations which are based on OpenFOAM.
For the target case, the air and fuel are not fully premixed prior to the combustor, which results in a partially-premixed flame. This is encountered quite often in gas turbines but is much less studied and understood than fully premixed flame problems. Recent developments on combustion modelling and relevant issues in the LES framework are reviewed in ref. [56] for gas turbine combustion. Some turbulent combustion models have been used in partially-premixed flame predictions based on LES, such as the PaSR (Partially Stirred Reactor) model [48][49][50], the transported Filtered Density Function or Probability Density Function model [57][58][59][60], the CMC (Conditional Moment Closure) model [61,62], etc. The present LES study applies the PaSR model to deal with the turbulence-combustion interactions.
The PaSR model accounts for finite rate chemistry. It solves the filtered LES equations using a model of the filtered combustion reaction rates,ω j for the j-th species. The reduction of the modelling ofω j is one of the core challenges in turbulent combustion modelling. The laminar Ahrrenius reaction rate is generally not valid for turbulent combustion, i.e.: where ρ is the filtered density, T is the Favre-filtered temperature and Y i is the Favre-filtered mass fraction of i-th species.
For the PaSR approach of modellingω j , the flow in a computational cell is split into two different parts; the fine structures in which mixing and reactions are assumed to take place, and the surroundings dominated by the large scale structures. Recent DNS (Direct Numerical Simulation) studies of combustion [63,64] confirmed the existence of fine structures and found that fine structure vortices at the flame are essentially parallel to the flame whereas those behind the flame are mostly perpendicular to the flame. In PaSR modelling, the fine structure component is treated like a perfectly stirred reactor, in which all present species are homogeneously mixed and reacted. After reactions have taken place, the species are assumed to be mixed with the surrounding component due to turbulence -this correspondingly takes place during the turbulent mixing time τ m . This then gives the final concentration in the entire, partially stirred, computational cell. The relative sizes of the two parts in the cell are governed by the combustion time and turbulent mixing time. The reaction rate for i-th species can then be scaled by the reactive volume fraction, κ, as in ref. [65]: where C 1 is the averaged concentration of the mixture leaving the computational cell, C 0 the initial averaged concentration of the mixture in the cell, and t the numerical time step. The term RR i is the laminar Arrhenius reaction source term, i.e. RR i =ω i (ρ, T , Y j ) (also refer to Eq. (4)). Correspondingly, the turbulence-combustion interaction in the PaSR model is reduced to the modelling of reactive volume fraction, κ. It is assumed that the final concentration C 1 is linearly related to the initial concentration C 0 and the unknown reacted concentration C, which results in the formulation of the reactive volume fraction (κ) as [65]: where τ m is the turbulent mixing time scale and τ c is the reaction time scale calculated by solving the fully coupled ODEs for the reaction system.
The reaction time scale τ c is determined by the chemical mechanisms which could be the global one-step reaction, reduced mechanism or comprehensive mechanisms. For the present ethylene/air reaction system, detailed mechanisms are available, such as the UCSD mechanism (46 species and 235 elementary reactions) [66] and Blanquart mechanism (149 species and 1651 elementary reactions) [67]. Such detailed mechanisms require significant computational cost and are mainly used for RANS simulations [68]. Based on these detailed mechanisms, reduced mechanisms have also been developed, such as the 19-species reduced mechanism [69], 22-species reduced mechanism [70], etc. These reduced mechanisms have been successfully used [71,72] in several LES studies of ethylene flames. However, the present study simulates many LES cases and those mechanisms remain too expensive. The global one-step (5 species) mechanism by Westbrook and Dryer [73] is applied in the present LES calculations. This global mechanism has been used for soot formation predictions by DNS [74]. Note that the present LES calculations mainly concern the global unsteady heat release of lean combustion induced from the unsteady flame dynamics due to acoustic forcing. The nonlinear flame dynamics is of key interest. Experiments on the partiallypremixed flame [28] that is the subject of the present study show that the heat release response as a function of forcing amplitude is similarly predicted from either the FSD, OH * or CH * measurements, for higher forcing frequencies (around f > 160 Hz). This suggests that, for this flame, the contribution of flame surface modulation to the total heat release response is significant and the flame area evolution is playing a key role. The variation in burning speeds due to timevarying equivalence ratio does not appear to be as significant. The present LES study mainly accounts for the flame dynamics at frequencies f > 150 Hz. Thus the experimental observations suggest that the reaction mechanism is not of major significance and the global onestep reaction mechanism should suffice for the present calculations.
Another important parameter in the PaSR model is the modelling of the turbulent mixing time, τ m , used to determine the reactive volume fraction (Eq. (6)). The turbulent mixing model is also a key component in the turbulent combustion models based on transported Probability Density Function methods [75,76]. In the original PaSR model within the RANS framework [77,78], the turbulent mixing time, τ m , is generally modelled using the integral time scale or Kolmogorov time scale. Recently, this modelling has been extended to the LES framework and some models have been developed. The default model of τ m in OpenFOAM uses the effective viscosity (μ e f f = μ + μ t ) and the dissipation rate, which has been used in previous studies of turbulent flames [50,79]. The present study applies a recently developed mixing time model which is used in the extended LES-PaSR model for high Reynolds, moderate Damköhler number turbulent flames [80,81]. The modelling is based on the subgrid velocity stretch time and Kolmogorov time scales, in the form of [80,81]: where the subgrid time scale τ and Kolmogorov time scale are calculated by: The turbulent mixing model shown in Eq. (7) was implemented in the OpenFOAM toolbox (version 2.3.0) and a modified low-Mach number reactingFOAM solver was developed for the present LES calculations. The C++ library OpenFOAM as a computational platform has previously been used for applications of varying complexity. The code employs an unstructured collocated Finite Volume Method (FVM) in which the discretisation is based on Gauss' theorem together with a semi-implicit time-integration scheme. The algorithm for pressurevelocity coupling is based on the PIMPLE method which results from combining the classic algorithms of SIMPLE and PISO (Pressure Implicit with Splitting of Operators) which is suitable for transient simulations. The convection divergence terms are discretised using a second order central difference scheme with the Sweby flux limiter to avoid unphysical oscillations. For temporal advancement, the second order implicit Crank-Nicolson scheme is used to discretise the unsteady terms, coupled with a fraction of first-order implicit Euler scheme to stabilise the calculations.
This section mainly concerns the numerical methods used to determine the FDF. As discussed previously, the decoupled method of calculating the acoustic waves and the unsteady heat release is used. Correspondingly, only a small part of the combustion system (see Fig. 1(a)) is employed in the LES calculations for FDF identification. The relatively small computational domain significantly reduces the computational cost of the LES. A z−cut of the computational domain is shown in Fig. 1(b) including the coordinate system used. It includes the whole combustor "enclosure" and the upper 80 mm of incoming circular duct. The six fuel injection holes are located 55 mm upstream of the combustor inlet. In the experiments, the diameter of the fuel injection holes is 0.25 mm which is small enough to lead to convergence problems in the simulations. To deal with this, the diameter of the six fuel injection holes is increased to 1.0 mm in the simulations, with the same total mass flow rate. Transverse fuel injection has complex interactions with the main air stream, and has been the topic of many studies both experimentally and numerically [82,83]. It still constitutes difficulties in numerical modelling due to the different flow features that are encountered in the complex instantaneous flow evolution. However, there exist some correlations describing the jet trajectory [84,85]. Based on these correlations [85], it is estimated that the jet penetration from the enlarged hole is around 79% of that from the original smaller hole. The effect on the fuel mixing is thus expected to be limited. The enlarged fuel injection holes reduce the complexity of the predicted flow structures, but the detailed effect is beyond the scope of the present study. An unstructured mesh is used for the present LES, with mesh independence checked and a final mesh containing about 3.48 million cells employed for all simulations presented. Meshes are clustered near the solid walls using a boundary-layer structured mesh with mean y + around 0.8. An unstructured mesh is used in all other regions. In the main combustion region, the mesh is nearly uniformly distributed.
For the base flow without acoustic forcing, the time-averaged bulk velocity at the combustor inlet is V b = 9.9 m/s, giving a Reynolds number of Re = dV b /ν = 17,000 [47]. In the experiments, the external forcing was introduced by two loudspeakers mounted upstream of the combustor. To emulate this forcing in the simulations, a single frequency harmonic velocity is superimposed on the mean flow at the computational inlet, with the form: (9) where A is the normalised velocity forcing amplitude and f the forcing frequency. A and f are varied independently in the simulations in order to obtain the FDF. Forcing of this form has been used to simulate harmonic loudspeaker forcing of a flame in previous numerical studies [30,[34][35][36][37]. This approach means that, at the computational inlet, the "acoustic" perturbations are mapped to "hydrodynamic" fluctuations for the purpose of the flame response as this is well known to be dominated by hydrodynamics. It should be noted that V b is the mean velocity entering the combustion chamber and V 0 the mean velocity at the inlet of computational domain -these are related by mass balance such that V 0 = 5.17 m/s (see Fig. 1

(b) for details).
In the simulations, all boundaries other than the inlet and outlet are treated as solid walls, where non-slip wall conditions are applied. In the experiments, it is observed that there is heat loss over the walls -the important walls are marked w1, w2 and w3 in Fig. 1(b). The heat loss effects on the unsteady flame dynamics and the final flame models have been studied previously [86][87][88]. It is found that heat loss tends to increase the heat release amplitude at low forcing frequencies and decrease it at high forcing frequencies. In the present simulations, the energy equation is solved and correspondingly heat loss effects on the unsteady heat release can be accounted for. A lower temperature than the adiabatic temperature is imposed on those walls with heat loss. As no experimental wall temperature measurements have been made, an estimated constant temperature is applied with T w1 = 1500 K, T w2 = 800 K and T w3 = 1000 K. Adiabatic conditions are applied for all other walls.
To determine the FDF defined in Eq. (1) from the present simulations, both the reference velocity signal and the heat release rate signal are required across different forcing amplitudes and frequencies. The forcing amplitude of the reference velocity is taken at the computational inlet (the position with bulk velocity V 0 in Fig. 1(b)), i.e. the value imposed of A in Eq. (9). For determining the phase of the FDF, a reference point P 0 is set at the combustor inlet (see Fig. 1(b)) in order to be consistent with the setup in experiments. The phase of the recorded velocity signal at point P 0 during the simulations is used as the phase of the reference velocity. For the heat release signal, the heat release rate Q calculated from the governing energy equation is integrated and recorded during the simulations. The signals are analysed spectrally using Fourier Transforms in order to determine the complex amplitude of the signals at the forcing frequency f. In the experiments, the values of OH * (f)/ OH * and CH * (f)/ CH * are used as estimates of Q (f)/ Q . The value of Q (f)/ Q is used to represent heat release rate fluctuations in the present simulations. Based on the heat release rate and reference velocity signals, the FDF defined in Eq. (1) can then be determined.

Validation of LES
No experimental data exists for the FDF at a global equivalence ratio of φ = 0.61 -the conditions under which combustion instability leads to limit cycle oscillations. It is thus important to validate the numerical method before it is used to calculate the FDF. Some experimental measurements for cold flow and partially-premixed reactive flow with a global equivalence ratio of φ = 0.55 are available, and are employed for validation of the present numerical method.

Case: cold flow
For the experimental configuration, the turbulent flow in the absence of fuel injection and a flame has been experimentally studied [89,90]. It is a pure turbulent flow past a bluff body, for which previous LES studies [34,61,90] have been performed. For comparisons, the previous LES studies are referred to as LES-1 and LES-2, corresponding to the results from refs. [61] and [34], respectively.
The main averaged flow structures can be seen from the timeaveraged axial (y-direction) velocity flow field, as shown in Fig. 2. The previous LES-2 result [34] is also included for comparison. The main flow structures captured by the present LES are similar to those from the previous LES-2 [34]. Both a central recirculation zone and a side recirculation zone can be observed, formed by the wake of the bluff body and the rearward-facing step, respectively. Shear layers are produced by these recirculation zones, which are very important for flame stabilisation. The agreement between the two sets of LES studies is qualitatively good. Differences can be observed, the main ones To further validate the present LES, direct comparison with experimental data is performed in terms of time-averaged and RMS axial (y-direction) velocities. Figure 3 shows velocity predictions from the present LES, two previous LES studies [34,61] and experimental measurements, for three distances from the bluff body of y/d = 0.22, 1.0 and 2.0. For the time-averaged velocity, the present LES predictions agree well with experimental data and previous LES calculations. For the RMS velocity, the previous LES studies generally underpredict compared with experimental data. The agreement of the present LES with experimental data is reasonably good, although the RMS velocity is a little overestimated in the far downstream region, e.g. y/d = 2.0. Note that the RMS velocity in Fig. 3(b) is

Case: partially-premixed reactive flow
The reactive flows for a partially-premixed flame with a global equivalence ratio of φ = 0.55, both without and with acoustic forcing are also used to validate the computational code. The forced case for which experimental data is available has a forcing frequency of f = 160 Hz. The response of the heat release rate and the non-linear flame dynamics are compared with available experimental measurements [28].
The natural flow in the absence of forcing is firstly studied. Figure 4 shows the time-averaged volumetric heat release rate from the present LES and the FSD image from experimental measurements [28], which represents a qualitative comparison of the heat release of the unforced reactive flow case. The agreement of the numerical prediction and the experimental data is reasonable. The present LES predicts a flame length of around 2.5d, slightly longer than is seen in  experiments, which implies that the speed of the combustion process is slightly under-estimated by the LES. Possible reasons for the difference are the global one-step reaction mechanism, the heat loss from the walls and the combustion modelling (especially the model constant c m ). The comparison is less good near the vertical walls, where the heat loss has large influence. Combustion in the present LES appears too strong close to the shear layers along the side recirculation zones even though low temperatures are imposed on the walls to partly account for the heat loss. Snapshots of the unforced reactive flow fields are shown in Figs. 5 and 6 for different flow quantities, including axial velocity (V in m/s), temperature (T in K), local equivalence ratio (φ) and turbulent/laminar viscosity ratio (μ t /μ), at different locations. It can be observed that the flame is anchored at the shear layers from the wake of the bluff body and the side recirculation zones, which is similar to the case of fully premixed flame [28,34]. Due to the intense heat release, the main central recirculation region behind the bluff body is enlarged compared to the cold flow condition (see Figs. 5 and 2). The two sets of shear layers are no longer joined together towards the top of the combustor chamber. Figures 5(c) and 6 show the flow fields of local equivalence ratio φ. The results in Fig. 5(c) demonstrate that the fuel is mixed with the incoming air after the injection holes. Prior to the combustor inlet, the fuel and the air are not fully premixed, which results in a partially-premixed mixture. This can be seen more clearly from Fig. 6. The fuel is spatially distributed and evolves with the main flow as it moves downstream. The fuel in regions directly influenced by the injection holes burns faster than in other regions, observed from the results in Figs. 6(c) and (d). Note that the results in Figs. 5(c) and 6 demonstrate that the local equivalence ratio in the main combustion region is mainly smaller than 0.8, within the lean combustion regime. This adds further weight to the suitability of using a global one-step reaction mechanism here.
The static Smagorinsky subgrid-scale LES model is used in the present study. Figure 5(d) shows the flow field of turbulent/laminar viscosity ratio (μ t /μ). It can be seen that the viscosity ratio is smaller than 0.3 in most regions, suggesting that the mesh resolution is quite fine and only a small part of the turbulence is modelled via the subgrid-scale model. The Smagorinsky model is thus suitable for the present study.
The forced reactive flow case is simulated to evaluate the performance of determining the heat release response to forcing. Velocity fluctuations are imposed on the mean velocity at the computational inlet (see Eq. (9)). A forcing frequency of f = 160 Hz is considered as experimental measurements are available at this frequency. Four forcing amplitudes are simulated, i.e. A = 0.05, 0.20, 0.40 and 0.65. For a harmonic velocity fluctuation at the computational inlet, an unsteady heat release rate result is recorded during the simulation. The simulation results presented in the following sections are based on at least 16 forcing cycles after transients have died away: phase averaging is thus carried out for at least 16 forcing cycles.
Fourier Transforms are used to process the time series of the heat release rate and the reference velocity, resulting in the gain and phase of the FDF in the frequency domain. The normalised amplitude of heat release rate fluctuation as a function of the forcing amplitude A is shown in Fig. 7, with simulation predictions compared to experimental measurements. The LES predictions agree reasonably well with the experimental measurements for both gain and phase. Note that the magnitude and phase of the heat release rate response measured by OH * and CH * chemiluminescence, and evaluated using the FSD based on OH PLIF are in good agreement in the experiments.  This implies that the variation in local equivalence ratio affecting the burning velocity probably results in increased flame area through increased flame wrinkling, which would be captured well by FSD measurements. Thus this validates our use of the global one-step reaction mechanism in the present LES calculations, although accuracy will always be slightly compromised compared to detailed mechanisms. It should be noted that the differences between the LES results and experimental data in Fig. 7 can be used to estimate the uncertainties of the predicted FDF gain/phase -these are estimated to be smaller than 11% and 9% for the gain and phase, respectively. The amplitude response from experiments in Fig. 7(a) demonstrates that the response is nearly linear up to a forcing amplitude of around A = 0.3 when non-linear effects start to develop. The LES predicts this linear to non-linear transition takes place earlier than in the experiments. It seems that the non-linear behaviour at relative high forcing amplitude is not accurately captured by the present LES. The observation here is similar to that in the fully-premixed flame at a forcing frequency of f = 310 Hz [34]. The differences between experiments and predictions may mainly come from the flame-wall interactions, specifically the heat loss effect from the walls and the global one-step chemical reaction in the near-wall regions. The phase results in Fig. 7(b) show a slight increase with forcing amplitude, which is captured reasonably well by the present LES.
The flame dynamics at a high forcing amplitude of A = 0.65 are visually shown in Fig. 8 at every 60°phase angle, for both the LES predictions and the experimental measurements. The image sequence shows clearly the deformation of the flame base, later resulting in radially inward rollup of the inner shear layer flame front and radially outward rollup of the outer shear layer flame front, which is similar to that in the fully premixed case [19,34]. This gives rise to the evolution of a mushroom-shaped flame contour. At phase angles of 280°and 340°, a new mushroom-shaped vortex starts to form at the base of the flame, and downstream vortex starts to weaken. This downstream mushroom-shaped vortex continues to weaken until it disappears, and only the new one remains in the flow field (see the results at the phase angles of 160°and 220°). It seems that the disappearance of the old mushroom-shaped vortex occurs slightly later in the present LES than that in the experiments, implying that the combustion process is slightly underestimated, consistent with the observation in the natural reactive flow case without forcing. It was observed in the experiments that the flame can impinge on the wall (see Fig. 8(d)) during the process. The wall-flame interactions are not captured well by the present LES.
To summarise, the present LES resolves the cold flow fields well compared with experimental measurements and previous LES. For the reactive flow cases without and with forcing, the present LES also captures the main flow structures very well, with the non-linear heat release response and the flame dynamics well predicted compared to available experimental data.

Full FDF determination by the present LES
The previous section confirms that the present LES code which uses the OpenFOAM toolbox can capture the partially-premixed reactive flow field well and can also predict the unsteady heat release response to acoustic forcing with good accuracy. The numerical method can thus be used for determination of the full FDF.
The target case exhibits combustion instability and hence limit cycle oscillations in the experiments [28]. The equivalence ratio is φ = 0.61 and no experimental data for the FDF is available. The same numerical method is applied to that in the previous section, except that the equivalence ratio is φ = 0.61 instead of φ = 0.55. LES calculations are carried out by varying the forcing frequency and forcing amplitude independently. The frequencies range from 150 Hz to 600 Hz. Frequencies below f = 150 Hz are not considered as the selfexcited oscillation observed has a frequency of around f = 348 Hz in the experiments [28]. For each frequency, four forcing amplitudes are performed, i.e. A = 0.1, 0.2, 0.3 and 0.4. Simulation data is based on at least 16 forcing cycles after the simulation transients have died away.
The dependence of the flame response on the forcing amplitude is shown in Fig. 9. The heat release saturates for all forcing frequencies except the lowest of f = 150 Hz, for which the response is approximately linear. The amplitude of the heat release response shows little variation with velocity amplitude, once above frequencies of f ≥ 250 Hz. At a forcing frequency of f = 350 Hz, very close to that of the self-excited oscillations (around f = 348 Hz), the amplitude of the heat release response is nearly constant across the four forcing amplitudes. For the phase results shown in Fig. 9(b), the phase shows little variation with velocity amplitude for the lower forcing frequencies of f = 150 Hz and f = 250 Hz. Above frequencies of f = 300 Hz, however, large phase variations with velocity amplitude are observed. For example, at a forcing frequency of f = 350 Hz, the phase jumps from around −0.67π at A = 0.1 to around −0.18π at A = 0.2, giving a phase change of around 0.49π . This large phase change may contribute to the limit cycle state of the oscillations.
The corresponding FDF from the present LES is given in Fig. 10 along with fitted lines for each forcing amplitude, obtained using a 15th order fit with the "fitfrd" command in Matlab. These fitted expressions will be used for analysis of the combustion instability and ensuing limit cycle. Note that fitted lines presented in Fig. 10 may not be valid below frequencies of f < 150 Hz. The gain is seen to fall off with increasing forcing frequency, other than for a small peak that appears at around f = 300 Hz for a forcing amplitude of A = 0.4 and around f = 350 Hz for the other three forcing amplitudes. With increasing forcing amplitude, the gain generally decreases. The nonlinearity of the gain is clearly evident -a linear response would not vary with forcing amplitude. For the phase response, a nearly constant time delay is evident up to a frequency of around f = 300 Hz for all forcing amplitudes, and even up to f = 400 Hz with the low forcing amplitude of A = 0.1. For higher forcing frequencies of f > 300 Hz, large phase changes are observed as already shown in Fig. 9(b). All the 24 LES runs use the same time step of t = 4.34 × 10 −6 s, with the total CPU time required to obtain the full FDF results around 12, 900 h.

Low order combustion instability network model, OSCILOS, and its validation
Combustion instability analysis is performed using the combustion instability low-order simulator (OSCILOS) [14,42], developed in the authors' group. It is written in Matlab © /Simulink © and has an incorporated Graphical User Interface (GUI). The solver is based on low-order network modelling. The thermoacoustic system is represented as a network of simple connected acoustic elements, where each element corresponds to a certain component of the system [91,92]. The acoustic wave behaviour is modelled analytically using linear wave-based methods, and a flame model is incorporated, capturing how the flame responds to acoustic waves. It provides predictions of the frequencies of resonant modes, their stability (positive and negative growth rates), mode shapes and the time evolution of disturbances. The fundamental basis of OSCILOS is similar to other thermoacoustic network models [2,12,26,35,93,94], which have been validated and used extensively in a variety of thermoacoustic problems [2,12,35,40,41].
The present study only concerns the 1-D plane (longitudinal) acoustic waves, as the ratio of the transverse dimensions of the elements to the acoustic wavelength is very small for the cases considered here. The flame is assumed "compact" compared to the acoustic wavelength. The acoustic analysis assumes that the combustor geometry can be represented as a network of connected modules, as schematically described in Fig. 11. Different modules have different sectional areas. For module k, the inlet and outlet are located at x = x k−1 and x = x k , respectively, where k = 1, 2, . . . N, with N the total number of elements. A compact flame located at x = x n (n is integer with 0 ≤ n ≤ N) is used as the heat source, and separates the unburned and burned gases. According to the linear acoustic theory, all flow and thermodynamic variables can be decomposed into a mean value and an acoustic perturbation [1,95,96], which is assumed to be small compared to the corresponding mean value. The acoustic field can be described by the superposition of forward and backward propagating plane waves. Considering acoustic waves propagating in both directions, the pressure, velocity and density in module k can be expressed as: where A + k and A − k mean the amplitude of the downstream and upstream propagating acoustic waves, respectively, E k = ρ k c 2 k /C p,k s k represents the amplitude of entropy waves, τ + Substituting Eqs. (10)-(12) into standard reactive flow balance equations, assuming low Mach number and weak (linear) disturbances, high order terms can be neglected and it is possible to get the steady and first order term balance equations to relate the upstream and downstream acoustic waves. Applying the process to all the elements in the system, a linear system can be obtained to describe the acoustic waves within the system. It is then possible to derive a global matrix G 1→N (s) to link the waves in the first (k = 1) and the last modules (k = N), which can be written as: (13) where the superscript · represents the Laplace transform, s = σ + i2π f indicates the Laplace variable, σ is the growth rate and f represents the frequency. The attenuation of entropy waves due to shear dispersion can be accounted for [97]. Entropy waves are considered to disappear when they reach the end of combustor, and the indirect noise due to entropy waves is neglected in the present study, due to the presence of open downstream boundaries. Pressure reflection coefficients R 1 and R 2 are employed to characterise the inlet and outlet acoustic boundary conditions. Readers can refer to ref. [42] for more details. At the inlet of the combustor, A − 1 (s) = 1 and E 1 (s) = 0 are set, along with A +  (13). The error at the outlet boundary (a measure of how well the outlet acoustic boundary condition is met for this choice of frequency and growth rate (i.e. s)) can be mathematically expressed as: which then makes it possible to plot the contour map of 20log 10 |δe(s)| with the growth rate σ and the frequency f, e.g. the results shown in Fig. 15. The eigenvalues are those frequencies and growth rates which satisfy the outlet boundary condition and are hence located at minima.
To close the thermoacoustic system, it is necessary to link the acoustic velocity perturbations upstream of the flame to the unsteady heat release rate, Q , integrated over the flame volume. This is provided in a flame model. In OSCILOS, different flame models can be selected, ranging through from simple linear n − τ models to nonlinear flame describing functions, either prescribed analytically or loaded from CFD data/experimental measurements. More details on OSCILOS can be found in ref. [42].
Before OSCILOS is applied for predicting the thermoacoustic characteristics of the present unstable combustion system, it is firstly validated using a well documented case; the experimental configuration developed at Laboratory EM2C, which benefits from a variety of combustion instability studies [12,[98][99][100]]. The combustor system is shown in Fig. 12, and includes a plenum, an injection unit and a combustion chamber with an open end. A compact flame is stabilised at the beginning of the combustion chamber. Experiments were carried out with different lengths of the plenum and chamber in order to vary the eigenvalues of the combustor system. Herein, we only take one unstable case, designated as L − 400 [12], for the comparisons between the combustion instability predictions from OSCILOS and the experimental measurements. The plenum comprises a straight cylindrical container with a length of 224 mm and diameter of 65 mm, and a smoothly convergent cylindrical unit with a length of 60 mm (this is divided into 50 subelements in the OSCILOS calculations). The diameters of the inlet and outlet are 60 mm and 35 mm, respectively. The injection unit has a length of 56 mm and a diameter of 22 mm. The length and the diameter of the combustion chamber are 400 mm and 70 mm, respectively. More details can be found from refs. [12,99]. The mean velocity at the outlet of the injection unit isū 2 = 4.13 m/s, the mean pressure is p 1 = 1 bar and the mean temperature is T 1 = 300 K. Methane is used as the fuel and the equivalence ratio is φ = 0.7. The measured mean temperature of the burned gases is around 1600 K, which is consistent with the predicted mean temperature in Fig. 13(b) by setting the combustion efficiency η as 0.825 in OSCILOS. It is thus possible to calculate the mean thermal properties and mean flow velocity in different sections. Figure 13 shows the plots of the mean velocityū and the mean temperature T in different sections calculated by OSCILOS, where T 3 = 1601 K and T 1 = T 2 = 300 K. In OSCILOS, acoustic losses occur at boundaries and at area increase interfaces between modules. The presence of mean flow also makes contributions to acoustic losses, which is for example discussed in ref. [101]. At the boundaries, acoustic losses increase when the magnitude of pressure reflection coefficient |R| decreases [14]. For the acoustic boundary conditions of the EM2C combustor, the inlet can be considered as a rigid wall and the pressure reflection coefficient is adjusted to account for the slight acoustic loss within the plenum. The outlet of the combustion chamber is open to atmosphere. Analyses of the sensitivity of predicted limit cycle to dissipations at the boundaries were carried out by successively varying the magnitudes of pressure reflection coefficients R 1 and R 2 . When R 2 is changed from −1 to −0.9, the predicted normalised velocity perturbationû u /ū u of the limit cycle varies from 0.69 to 0.74, and the corresponding eigenfrequency varies from 126.3 Hz to 124.2 Hz. The inlet pressure reflection coefficient plays a more important role on the dissipation of the system. When R 1 decreases from 0.95 to 0.85, the predictedû u /ū u varies from 0.69 to 0.36, and the corresponding eigenfrequency varies from 126.3 Hz to 131.7 Hz. The inlet and outlet pressure reflection coefficients thus should be carefully selected to account for acoustic losses of the system. The damping rate of the system, α, can be evaluated with the stable flame using the shortest chamber, which is similar to those used in the experiment [12]. The predicted growth rate of the first mode σ 1 is negative and equals to σ 1 = −α = −45 rad s −1 when the inlet pressure reflection coefficient R 1 is set to 0.95. The outlet pressure reflection coefficient is R 2 = −1. An end correction of 0.4 times of the diameter of the combustion chamber diameter is taken into account to consider the sound radiation at the outlet of the chamber [15,102,103].
The flame describing functions are provided by experimental measurements which are imported into OSCILOS. A fitting procedure to obtain the FDF in mathematical form is then performed within OS-CILOS. Figure 14 shows the experimentally measured and fitted flame describing functions. The fitted FDF has order 16 and captures the shape of FDF for the most "dangerous" frequency range (0−400 Hz) where combustion instability is known to occur. Note that the ratio of forcing amplitude to its mean valueû u /ū u is used as the normalised velocity perturbation in the present study, while the ratio of RMS value to its mean value u u,rms /ū u is used in the experiments [12,99].
Substituting the fitted FDF results into the thermoacoustic network model, we can obtain the evolution of the eigenvalues with changing velocity perturbation levels. Figure 15 shows the contour maps of 20 lg |δe(s)| in the s-plane for the four normalised velocity perturbations. For weak velocity perturbations (such as in Fig. 15(a)), the growth rate of the main mode equals 21.5 rad s −1 , meaning that the system is unstable with disturbances oscillating at the corresponding eigenfrequency of f = 130.8 Hz. With increasing velocity perturbation, the growth rate decreases (see the evolution of the main mode's growth rate in Figs. 15(b)-(d)  finally established when the growth rate is equal to zero, between the normalised velocity perturbations of 0.679 and 0.848. Figure 16 shows the evolution of the eigenfrequency and its corresponding growth rate with the normalised velocity perturbations. With increasing velocity perturbations, both the eigenfrequency and the growth rate decrease. We can then predict the normalised velocity perturbation for the limit cycle based on a linear interpolation for which the growth rate is zero, giving 0. 69

Thermoacoustic analysis of the present Cambridge experimental configuration
The discussions thus far have addressed the two separate aspects needed for a coupled approach to combustion instability analysis: the FDF calculations and the low order thermoacoustic network model. Both have been validated against experimental measurements, and we now combine them for the target test case.
The target Cambridge configuration can be represented as the network of connected modules shown in Fig. 17 where the geometry details are adopted from ref. [41]. The combustor system constitutes four sections: the upstream settling chamber, the plenum, the inflow pipe and the combustor chamber, from upstream to downstream (also see Fig. 1(a) for reference). The upstream settling chamber has a length of 80 mm and a diameter of 34 mm. The plenum is represented as having three parts with lengths of 50 mm, 200 mm and 50 mm, respectively. The part with constant cross sectional area has a diameter of 100 mm. The other two parts are assumed to be linearly connected to the upstream settling chamber and the inflow pipe, respectively. In OSCILOS, these two parts are divided into 10 sub-elements for the calculations. The inflow pipe has a length of 450 mm and a diameter of 17 mm. Different combustor chamber lengths were used in the experiments [28]; self-excited oscillations were observed for a length of 350 mm (the diameter is 70 mm) and this will hence be the length used here.
The experiments were carried out at atmosphere pressure. The mean velocity in the upstream settling chamber is the same as in the inflow pipe, i.e. 5.17 m/s, accompanied by a mean temperature of 300 K. For the present partially-premixed flame with a global equivalence ratio of φ = 0.61, assuming a combustion efficiency of η = 0.95, giving the mean temperature of the burned gas to be around 1766 K in OSCILOS. The mean velocities and temperatures along the four modules as calculated by OSCILOS are shown in Fig. 18. The acoustic boundary at the inlet is considered as a rigid wall, and the outlet of the combustor chamber is treated as open to atmosphere. Analyses of the sensitivity of predicted limit cycle to dissipations at the boundaries were carried out as well. When R 1 and R 2 are changed by 10%, variations in the predicted normalised velocity perturbation  u u /ū u of the limit cycle and corresponding eigenfrequency are within 1%. The stability of the system thus does not depend on the dissipations at the two boundaries. We thus set R 1 = 0.95 and R 2 = −1, which are the same as those in Section 6. The unsteady heat release response to perturbations is modelled via the flame model obtained via the present LES, as discussed in Section 5. The full FDF is shown in Fig. 10, which will be used for the present thermoacoustic analysis. Note that for the present target case of Cambridge configuration under partially-premixed conditions, a kinematic flame model has been developed based on the G-equation [41,104]. Although the main features of the flame response were reasonably captured, the accuracy was not sufficient for accurate limit cycle prediction.
The obtained FDF data are imported into OSCILOS and the 15th order fits for each forcing amplitude obtained (see Fig. 10). The thermoacoustic modes of the Cambridge configuration are then calculated. Figure 19 presents contour maps of 20 lg |δe(s)| in the s-plane for the four normalised velocity perturbations, i.e.û u /ū u = 0.1, 0.  Fig. 19 implying that the mode is always stable. For the mode around f = 342 Hz, the growth rate is positive for the lower two perturbation levels ofû u /ū u = 0.1 and 0.2, and Table 1 Comparisons between experimental measurements and predictions at limit cycle for the Cambridge configuration.
To determine the frequency and the amplitude at the limit cycle of the present thermoacoustic system, the evolution of the eigenfrequency and its corresponding growth rate with increasing velocity perturbation are presented in Fig. 20. Using linear interpolation, both the eigenfrequency and the amplitude of the limit cycle can be predicted. The predictions are shown in Table 1 along with the experimental measurements. It can be seen that both the frequency and the amplitude of the limit cycle resulting from combustion instability are well predicted. An analysis of the sensitivity of the predicted limit cycle to variations in the gain and time delay of the flame describing function [105,106] were performed by changing these values by ±10% (this being a reasonable estimate for the uncertainty in the FDF). The results are summarised in Table 2. The change in the predicted limit cycleû u /ū u varies from 5.8% to −3.1% and that of corresponding eigenfrequency varies from 0.9% to 1.3%. The change in the time delay of the FDF therefore plays an more important role, but   20. Evolution of the eigenfrequency (marked •) and the corresponding growth rate (marked ) with normalised velocity perturbationsûu/ūu predicted by OSCILOS for the Cambridge configuration. The blue marker + represents the predicted normalised velocity perturbation for the limit cycle with the marker × for the corresponding frequency. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
the sensitivity to uncertainty in both variables is small. These analyses indicate that the saturation of the combustion process has been successfully captured, and combining determination of the FDF from the present LES with low order thermoacoustic modelling via OSCILOS accurately predicts the unstable behaviour of the present Cambridge configuration.

Conclusions
Predicting the limit cycle behaviour resulting from combustion instability has been a long-standing challenge in the combustion instability community. This is primarily because it depends on accurate characterisation of both the non-linear flame response to acoustic waves, and the acoustic wave behaviour and losses within the combustion system. The present study has successfully performed a limit cycle analysis of combustion instability involving a partially-premixed flame combustor. This was achieved by combining a weakly non-linear flame model in the form of a Flame Describing Function (FDF), obtained using large eddy simulations (LES), with a low order thermoacoustic network modelling tool. The target case was a bluff body stabilised, lean, partially-premixed flame combustor developed at Cambridge University, for which previous experimental data are available. Low Mach number LES was used to determine the non-linear heat release rate response to acoustic forcing, i.e. the FDF, using the CFD toolbox OpenFOAM. The low order thermoacoustic network modeling tool, OSCILOS, developed in the authors' group, then captured the acoustic wave behaviour within the combustion system. This is the first work, to the authors' knowledge, which studies this particular partially-premixed flame combustor using LES.
The LES method which used the CFD toolbox OpenFOAM was firstly validated. Turbulent combustion was modelled using the Partial Stirred Reactor (PaSR) model with a global one-step reaction mechanism. Both the unforced cold flow and the unforced/forced reactive flows were simulated and compared with available experimental data. The results demonstrated that both the flow and flame dynamics, as well as the unsteady heat release, were captured well. Simulations were then performed with varying inlet velocity in order to determine the full FDF. Both the forcing frequency and the forcing amplitude were varied independently, with six frequencies and four normalised forcing amplitudes considered. Non-linearity of the obtained FDF was clearly visible. Heat release saturation and phase changes were also observed at higher forcing frequencies, which could contribute to saturation into the limit cycle.
The solver OSCILOS is a combustion instability low-order simulator based on Matlab © /Simulink © . It models a combustor as a network of connected simple modules and combines wave-based linear analytical models for the acoustic waves with a more complex flame model, which for the present study was based on importing data from the LES flame study. Its validation for limit cycle studies was conducted for a well-documented unstable combustor system. The eigenfrequencies and growth rates were predicted for different levels of velocity perturbation upstream of the flame. Based on interpolation to the zero growth rate condition, the frequency and forcing amplitude of the limit cycle oscillations were predicted and both agreed well with the experimental measurements. This confirmed that OS-CILOS is capable of capturing limit cycle behaviour.
For the target Cambridge experimental configuration, the obtained FDF and validated thermoacoustic network modelling tool, OSCILOS, were then combined to study combustion instability. Two modes were predicted, the first being stable with a frequency of around f = 49 Hz and the second one unstable with a frequency of around f = 342 Hz. The limit cycle frequency and amplitude were predicted to be 342 Hz and 0.26, respectively, which agreed well with the experimental measurements of around 348 Hz and 0.21, respectively. This is the first work, to the authors' knowledge, which predicts the limit cycle frequency and amplitude under unstable conditions by combining high-fidelity CFD methods with a low order network modelling tool for the target case. This confirms that an open-source software framework combining OpenFOAM and OSCILOS, can be used to study combustion instability problems numerically, with good accuracy. The study also suggests that a sufficiently accurate flame model can be deduced from high-fidelity LES using the open source CFD toolbox OpenFOAM.