## Abstract

Predictive modelling of turbulent combustion is important for the development of air-breathing engines, internal combustion engines, furnaces and for power generation. Significant advances in modelling non-reactive turbulent flows are now possible with the development of large eddy simulation (LES), in which the large energetic scales of the flow are resolved on the grid while modelling the effects of the small scales. Here, we discuss the use of combustion LES in predictive modelling of propulsion applications such as gas turbine, ramjet and scramjet engines. The LES models used are described in some detail and are validated against laboratory data—of which results from two cases are presented. These validated LES models are then applied to an annular multi-burner gas turbine combustor and a simplified scramjet combustor, for which some additional experimental data are available. For these cases, good agreement with the available reference data is obtained, and the LES predictions are used to elucidate the flow physics in such devices to further enhance our knowledge of these propulsion systems. Particular attention is focused on the influence of the combustion chemistry, turbulence–chemistry interaction, self-ignition, flame holding burner-to-burner interactions and combustion oscillations.

## 1. Introduction

Predictive modelling of turbulent combustion is becoming increasingly important for the development of air-breathing engines, internal combustion engines, furnaces and for power generation. The increase in computational power in the past decade has made many of these flow configurations numerically accessible. Nevertheless, the interaction of turbulence with physical processes, such as chemical kinetics and thermal radiation, is a great challenge. In the gas turbine industry, Reynolds-averaged Navier–Stokes (RANS) models, together with flamelet or eddy break-up combustion models, are the primary means of analysing combusting flows (Hsu *et al*. 1997), mainly due to their fast turnaround time and success in providing design guidance for meeting exit temperature profile requirements. RANS has also been used to predict emissions, but such predictions have met with mixed success. For high-speed ramjet and scramjet engines, similar methods are currently used, often together with more detailed chemical reaction mechanisms (Oevermann 2000), in order to better predict the ignition delay time of the mixture. A key issue for the successful operation of such engines is rapid mixing between fuel and air prior to self-ignition, which however is difficult to predict with RANS. Significant advances in modelling non-reactive turbulent flows are now possible with the development of large eddy simulation (LES; Sagaut 2001; Grinstein *et al*. 2007) and similar methods. The philosophy behind LES is to explicitly solve for the large- (energetic) scale flow, directly affected by boundary conditions, while modelling the small- (less energetic) scale flow. The development of LES has so far primarily been based on ordinary turbulence theory and on RANS (Wilcox 1993), resulting in improved predictions due mainly to that, in LES, the large scales are resolved. Considering the use of LES in modelling combustion for propulsion applications, we recognize the direct benefits offered by resolving the large-scale flow, but need to develop a better understanding for what is required in terms of modelling the turbulence–chemistry interactions, in particular for reduced reaction mechanisms, and how well such simulations can predict mixing, self-ignition, flame stabilization and combustion oscillations in complex geometries.

The objectives of this study is to present a short overview of the current status of finite-rate chemistry LES and to examine how such LES models can handle propulsion applications such as gas turbine engines and scramjet engines. The cases selected are simplified, yet realistic, and associated with laboratory studies providing experimental data for validation.

## 2. Large eddy simulation combustion modelling

The reactive flow equations are the balance equations of mass, momentum and energy describing convection, diffusion and chemical reactions (Poinsot & Veynante 2001). In LES, all variables are decomposed into resolved and unresolved components by a spatial filter such that , where is the Favre-filtered variable (Sagaut 2001). Filtering yields(2.1)Here, , , , and , respectively, are the resolved density, pressure, velocity, temperature and species mass fractions, , , , , and are the constitutive equations for a linear viscous mixture with Fourier heat conduction and Fickian diffusion. In addition, are the species' reaction rates, *P*_{ij} the stoichiometric coefficients and the Arrhenius reaction rates, with and being the rate constants, and , and are the subgrid stress tensor and flux vectors, representing the small-scale unresolved flow physics that needs to be modelled. The viscosity, *μ*, is given by Sutherland's law, and the species and thermal diffusivities are modelled as *D*_{i}=*μ*/*Sc*_{i} and *κ*=*μ*/*Pr*, respectively, in which *Sc*_{i} and *Pr* are the Schmidt and Prandtl numbers. Here, the mixed model (MM; Sagaut 2001; Bensow & Fureby 2007) is used in which , and , with *μ*_{k} being provided by the one-equation eddy viscosity model (Schumann 1975). To reduce the computational cost, we use *wall-modelled* LES, in which a subgrid wall model (Fureby 2007*a*,*b*) is used to model the near-wall flow physics in ** B**,

*b*_{i}and

*b*_{E}by adjusting

*μ*

_{k}in the first near-wall cell.

The reaction-rate model is required to handle non-unity Lewis number1 mixing, extinction and reignition, fuel modulation and different modes of combustion, including accurate predictions of transient species. Mixing-based (,Bray & Moss 1977; Magnussen 1981), or flamelet (Weller *et al*. 1998; Pitsch 2006), models are unsuitable for this devoid of modifications, and we thus focus on finite-rate chemistry models (Magnussen 1981; Karlsson 1995; Colin *et al*. 2000; Givi 2006; Porumbel & Menon 2006). Here, we use the eddy dissipation concept (EDC) model (Magnussen 1981) adapted for LES, in which , with and , respectively, being the conditions in the fine structures and the surroundings and the (internal) intermittency factor. The conditions in the fine structures and surroundings are related by the subgrid balance equations(2.2)in which is the subgrid residence time. In order to provide an independent evaluation of the sensitivity of the results to the reaction-rate model, the filtered reaction rates are also modelled by the thickened flame model (TFM; Colin *et al*. 2000), in which(2.3)where *F* is the flame thickening factor and *E* is the ratio of the true flame wrinkling to the resolved flame wrinkling, here modelled by the fractal model (Fureby 2004).

## 3. Numerical methods

The reactive LES equations are discretized using Reynolds transport theorem in a finite-volume framework, designed to handle arbitrary cell shapes, using the OpenFOAM C++ library (Weller *et al*. 1997). Depending on the flow speed, different application software will be used: for low Mach number flows, a semi-implicit code based on a second-order accurate Crank–Nicolson time-integration scheme and a pressure implicit with splitting of operators procedure based on a Rhie & Chow interpolation for the cell-centred data storage structure is used (Issa 1986), and for high Mach number flows, a fully explicit code based on a second-order accurate total variational diminishing Runge–Kutta scheme (Gottlieb & Shu 1998) is employed. In both codes, the reconstruction of the convective flux function is performed by linear interpolation between neighbouring control volumes, essentially resulting in a second-order accurate algorithm, with unphysical high-frequency high-wavenumber oscillations damped by a sixth-order hyperviscosity term. For the remaining viscous and subgrid fluxes, linear interpolation is used to reconstruct the inner gradient, essentially resulting in a second-order accurate scheme. Both codes use collocated cell-centred variable arrangements, and the equations are solved sequentially with a fixed time step corresponding to a Courant number of approximately 0.3.

## 4. Verification and validation

The computational methodology used, involving the OpenFOAM platform (Weller *et al*. 1997), has been verified for a range of cases using systematic grid refinement studies, the method of manufactured solutions and modified equations analysis to estimate the truncation errors. The application codes are extensively validated against direct numerical simulation or experimental data, and many such cases have been reported earlier, most performed with LES, but also some with RANS to provide a reference to current engineering capabilities. Some of the flows considered includes the Taylor–Green flow (Drikakis *et al*. 2007), homogeneous isotropic turbulence (Fureby *et al*. 1997), fully developed turbulent channel flow (Fureby *et al*. 2004), flow past various blunt or streamlined bodies (Fureby *et al*. 1999; Wikström *et al*. 2004; Persson *et al*. 2006), supersonic flows (Fureby *et al*. 2007, 2008) and reacting flow in different laboratory combustors. Figure 1 shows some representative results from the low-swirl burner of Shepherd *et al*. (2002), studied experimentally by Petersson *et al*. (2007) and computationally by Nogenmyr *et al*. (2008), and the triple annular research swirler (TARS) studied experimentally by Li & Gutmark (2006) and computationally by Fureby *et al*. (2006). These and other similar verification and validation studies give us confidence in the computational models and methods, as well as the procedures used to investigate and evaluate the numerical results.

## 5. Gas turbine combustor application

Annular gas turbine combustors are mainly used in aeropropulsion, but also for power generation and for marine applications. The current trend is to operate fuel lean and premixed, close to the lean blow-out limit. Under such conditions, low emissions and improved performance can be obtained, but as they operate close to the lean stability limit, they are vulnerable to combustion oscillations that may lead to flashback, blow-out and acoustical and mechanical vibrations. Combustion oscillations have been studied theoretically, experimentally and computationally for can combustors and single burner arrangements (Candel 1992), and are reasonably well understood. However, little is known about their behaviour in multi-burner annular combustor configurations, representative of real engines. Fureby (submitted) made an attempt to examine whether LES could be used to shed some light on annular multi-burner combustor dynamics using a model annular combustor derived from the General Electric (GE) LM600 laboratory combustor (Hura *et al*. 1998), used in several swirl combustor studies (Kim *et al*. 1999; Grinstein & Fureby 2004). This model annular combustor consists of 18 burners characterized by velocity and fuel mass-fraction profiles provided by GE (Hura *et al*. 1998), and representative of the DACRS burner, used in the single burner laboratory experiments (Hura *et al*. 1998), which will be used as reference. According to table 1, we present results for both CH_{4} and *n*-C_{10}H_{22} combustion at *ϕ*=0.56, a combustor pressure of *p*_{0}=6.18×10^{5} Pa, a reactant temperature of *T*_{0}=644 K and a centreline axial jet inflow velocity of *v*_{0}=110 m s^{−1}. By changing the velocity profiles, the effects of changing swirlers are examined, as characterized by different swirl and radial numbers, *S* and *R*, respectively, and by modifying the impedance of the outflow boundary conditions (characterized by the parameter *K*), the sensitivity of the outflow turbine boundary conditions is examined. The chemistry consists of a two-step reduced mechanism CH_{4}+1.5O_{2}←CO+2H_{2}O or C_{10}H_{22}+10.5O_{2}←10CO+11H_{2}O and CO+0.5O_{2}+CO_{2}, with rate parameters and Strouhal numbers determined from the Gri Mech 3.0 mechanism.

As an illustration of the combustor flow fields, being fairly similar at a first glance, the results of case 4 are presented in figure 2*a*,*b*. Between the dump plane and the swirling jets, toroidal-shaped vortices are formed that merge in the region between the jets. Adjacent helicoidal vortex systems form in the combustor when the swirling jets break-up, whereby the recirculation zones merge to one in the annular combustor. As reactants enter the combustor through the swirling jets, to mix with the hot combustion products, reactions occur along the wrinkled flame whereby chemical energy is released, causing volumetric expansion and flow acceleration. Each flame takes the shape of a short wrenched expanding tube, along which intermediate species and NO_{X} are formed, which folds back on itself to collapse in a tongue twisted by the helicoidal vortex. The flames are stabilized by the combined effect of swirl and axial velocity reduction enforced by the dump plane as the streamlines fan out. Vorticity is produced within the swirling jet shear layers, along the annular cooling air slits and through the baroclinic torque (∇*ρ*×∇*p*) due to exothermicity. As the flames are located between the clockwise swirling jets and the adjacent rotating vortex sheets, they experience both shear-layer instabilities and high-strain effects, affecting the turbulence–chemistry interactions. Viewed at one instant, the flames are very different and the outer (semi-transparent) combustor wall, coloured by pressure, reveals intertwined high- and low-pressure regions, emphasizing significant pressure variations.

As an example of how the mean flow is affected by modifying the fuel, outflow impedance conditions, swirl number and radial number, figure 2*c* shows cross-sectional profiles of the axial velocity, , from the cases in table 1. For case 4 (with *K*=1), the shallow recirculation region is located around *x*/*D*≈1.7; for case 5 (with *K*=10), it is repositioned upstream to *x*/*D*≈1.3 and weakened; for case 6 (burning *n*-C_{10}H_{22}) and case 7 (with *S*≈0.49), it is strengthened and located at *x*/*D*≈1.4 and 2.0, respectively; whereas for case 8 (with *R*≈0.04), it is eliminated. The shape of the cross-sectional -profiles are similar, with the exception of case 7 (with *S*≈0.49), which results in a higher peak velocity at *x*/*D*=0.18 due to the lower swirl number and is 5 per cent higher for the 18 burner annular combustor cases 4–8 when compared to the single burner cases. The lower swirl number of case 7 also results in a narrower distribution at *x*/*D*≈0.72. Compared with the single burner laboratory data, the differences between all computed cases are most evident at the flame tip, at *x*/*D*=0.72, where the present LES shows somewhat lower core velocities than the data and the LES results of Kim *et al*. (1999).

Combustion oscillations can be due to several reasons and the LES predictions will next be used to study the influence of outflow boundary condition and fuel and premixer characteristics on the combustor dynamics. The wave-transmissive boundary condition (Poinsot & Lele 1992) is based on linear relaxation, so that the incoming characteristic is , where *p*_{∞} is the pressure at infinity and *K* a stiffness parameter. For small values of *K*, *p* becomes close to *p*_{∞} and the outlet becomes non-reflecting, whereas for large *K*, *p* is equal to *p*_{0} and the outlet becomes reflecting. In case 4 (with *K*=1), the outlet is non-reflecting and the acoustic waves exit with small reflection levels, and in case 5 (with *K*=10), the outlet is less reflective. For *K*=1, the acoustic feedback is small with r.m.s. pressure fluctuations of *p*_{r.m.s.}≈2.4 kPa, whereas for *K*=10, the acoustic feedback is larger, with *p*_{r.m.s.}=6.4 kPa. The flames look similar in both cases, but the flame in case 5 is shorter and more compact than in case 4. Comparing the pressure frequency spectra evaluated at (*x*/*D*, *y*/*D*, *z*/*D*)=(1, 0, 0) and the pressure distributions on the combustor wall in figure 2*d*, we find that case 4 (with *K*=1) is dominated by *azimuthal* pressure fluctuations at 2240 Hz, with additional peaks at 686 and 1963 Hz, with magnitudes of approximately 20 per cent of the main frequency. Case 5 (with *K*=10) is dominated by *radial* pressure fluctuations at 2410 Hz with additional peaks at 1010, 1920 and 2600 Hz. For case 6 (burning *n*-C_{10}H_{22}), we find this case to be dominated by *radial* pressure fluctuations at 2170 Hz, with additional peaks at 520, 837 and 1250 Hz, and r.m.s. pressure fluctuations of *p*_{r.m.s.}≈5.1 kPa. For case 7 (with *S*≈0.49), we find this case to be subject to *mixed radial and azimuthal* pressure fluctuations at 2187 Hz, having additional peaks at 490, 850 and 1440 Hz, and r.m.s. pressure fluctuations of *p*_{r.m.s.}≈2.6 kPa. For case 8 (with *R*≈0.04), we find this case to be dominated by *longitudinal* pressure fluctuations at 724 Hz and *azimuthal* pressure fluctuations at 2124 Hz, with additional peaks at 1245 Hz, and r.m.s. pressure fluctuations of *p*_{r.m.s.}≈4.6 kPa. The laboratory combustor (cases 1–3) has r.m.s. pressure fluctuations of *p*_{r.m.s.}≈1.7 kPa that are dominated by *longitudinal* pressure fluctuations at 645 Hz.

## 6. High-speed ramjet and scramjet combustion

High-speed flight has, for a long time, been of interest to man—both for terrestrial travel and space exploration and transport. Key issues for these types of travel are the vehicle design and the propulsion system, and how to integrate them. Two types of propulsion systems, ramjets (Fry 2004) and scramjets (Curran 2001), are suitable for such vehicles. In a ramjet, the flow is decelerated to subsonic levels before it enters the combustor, allowing an efficient operational regime of 3<*Ma*<5 (where *Ma* is the Mach number), above which the deceleration leads to excessive thermal losses, whereas in a scramjet, the flow through the engine remains supersonic, thus allowing an operational regime of 6<*Ma*<15. Seamless integration of ramjet and scramjet operations is possible in the same engine (Andreadis 2007) by allowing it to sequentially operate as a ramjet, dual-mode ramjet, dual-mode scramjet and scramjet. Below *Ma*≈3, a turbine-based combined cycle or a rocket-based combined cycle, with externally or internally integrated rockets, can be used.

Part of the reason for the slow pace in high-speed combustion applications is the complex nature of the high-temperature reacting flow in the engine. Most research has been performed in ground-based research facilities, but recently, test engines have flown on conventionally powered vehicles (Andreadis 2007). Modelling high-speed combustion is a challenge due to the complex aerothermodynamics, but is currently a viable alternative to laboratory experiments due to the short test times in shock-tube facilities and the inaccessibility of appropriate operating conditions in blowdown facilities. LES has proven successful in this (Génin *et al*. 2003; Berglund & Fureby 2006), and recently the influence of the chemical kinetics on the self-ignition in vitiated hot supersonic air was also studied using one-, two- and seven-step mechanisms (Davidenko *et al*. 2006). The ONERA/JAXA supersonic combustion rig (Sunami *et al*. 2005) was used by Berglund *et al*. (2008) and in the present study, and consists of an air-vitiation heater (providing an N_{2}/H_{2}O mixture with a mass-fraction distribution of 0.23 : 0.70 : 0.07 at a temperature of 830 K, a velocity of 1449 m s^{−1} and a pressure of 34.2 kPa) and a constant-area duct connected to the diverging area combustor. The flame is stabilized by a wedge-shaped flameholder of height *h*, from which H_{2} is injected at the top and bottom walls and at the base into the vitiated hot supersonic airstream, figure 3*a*.

Here, we report more recent results focusing on the combustion physics and the influence of the turbulence–chemistry interactions using EDC LES and TFM LES and the seven-step mechanism of Davidenko *et al*. (2006) In figure 3*a*, the scramjet combustor flow field is illustrated. The flow can be divided into three regions: *the induction zone*, where turbulence controls mixing and combustion, in which the flame essentially follows the shear layers; *the transitional zone*, in which the flame escapes the shear layers to develop three-dimensional undulating structures, affecting the convective mixing and exothermicity; and *the turbulent combustion zone*, in which the three-dimensional coherent structures dominate the convective mixing and exothermicity. Longitudinal vortex structures develop primarily in the boundary layers, just upstream of the flameholder, vortex rings are generated around the H_{2} jets emanating from the three strut-based injectors, horseshoe vortices develop around the injectors situated at the top and bottom walls of the flameholder, and transverse vortical structures are shed off the upper and lower trailing edges of the flameholder. The leading-edge flameholder shock reflects at the upper and lower combustor walls before it impacts on the flame, creating further vortical structures through shock–boundary-layer interactions and shock–flame interactions, respectively. Similarly, the weaker shock waves emanate from the kink in the flameholder geometry, whereas the rarefaction waves from the strut base contribute to the generation of vorticity. In figure 3*b*, we show measured and predicted spontaneous flame images from the EDC LES, suggesting good qualitative agreement. Intense and stable flame emission appears downstream of the flameholder, with the flame gradually growing in thickness until it suddenly widens at *x*_{te}≈9 h. The blunt flame front oscillates back and forth with an amplitude of 3 h, and separate into an upper and lower branch at *x*_{te}≈14 h. George *et al*. (2007) suggested that this behaviour may be a consequence of thermal blockage.

Figure 3*c* presents the time-averaged bottom wall pressure distributions from the laboratory data and from the seven-step EDC LES and TFM LES. The first peak in the wall pressure is due to the incidence of the shock wave from the leading edge of the flameholder, the second peak is due to the weaker shock wave emanating from the kink in flameholder geometry, whereas the third peak is due to the volumetric expansion due to exothermicity. Included in figure 3*c* are also experimental data and predictions without combustion to demonstrate the difference in wall pressure due to exothermicity. The significant increase in wall pressure due to combustion indicates that the mixing is sufficient to burn most of the fuel, which is a prerequisite for successful scramjet operation. We have also evaluated the combustion efficiency at *x*_{te}=60 h, which is found to be 96 per cent for the EDC LES and 81 per cent for the TFM LES, quantifying some of the differences between the two turbulence–chemistry interaction models.

## 7. Concluding remarks

Here, we present LES predictions of both multi-burner annular gas turbine combustor and high-speed vitiated hot confined supersonic air in a simplified scramjet engine. The LES models used are based on finite-rate chemistry models using the eddy dissipation concept and the TFM. Good agreement with the limited reference data available for these two cases is observed. In order to further validate the LES models, comparisons are also made with experimental data in various laboratory combustors, such as an open low-swirl burner, resulting in a stratified mixture, and a modern gas turbine premixer connected to a round can combustor. The good agreement obtained in the validation studies gives confidence to the LES methodology and the ability to predict combustion under realistic, yet simplified, conditions. The predictions obtained for the multi-burner annular combustor and scramjet engine are used to further elucidate the complex flow fields in these devices, which in turn may be used to direct the design of new combustor devices, to identify critical issues and to further improve the general understanding of the processes involved in gas turbine and scramjet combustion.

## Acknowledgments

The author wishes to acknowledge scientific support from Suresh Menon, Fernando Grinstein, Hukam Mongia, Vladimir Sabelnikov, Jon Tegnér and Magnus Berglund, and the financial support of the Swedish armed forces.

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Applied large eddy simulation’.

↵The Lewis number compares the thermal and species diffusion speeds.

- © 2009 The Royal Society