<figure><img src="image_1.png" /></figure>

Dedication

A dedicated paper for Pr Sébastien Candel EM2C laboratory, CentraleSupélec

Christophe Bailly<sup>a</sup>, Daniel Durox<sup>b</sup>, Matthew Juniper<sup>c</sup>, Nicolas Noiray<sup>d</sup>, Thierry Poinsot<sup>e,*</sup>,
Thierry Schuller<sup>e,f</sup>, Denis Veynante<sup>b</sup>

<figure><img src="image_2.png" /></figure>

ARTICLE INFO

Keywords:

Turbulent combustion

Combustion noise

Combustion theory

Instabilities

Gas turbine

Annular chambers

Nonlinear modes

A B S T R A C T

This dedicated paper presents an overview of six fields of research in combustion in which Pr Sébastien Candel has made remarkable contributions: cryogenic combustion and liquid rocket engine dynamics, turbulent combustion modeling, combustion noise, flame dynamics and combustion instability, flame describing functions in combustion instability analysis, combustion dynamics and instabilities in annular systems. The paper is designed to (1) provide a review of research in each of the six fields, which can be used by all members of our community and (2) show how Pr Candel's work has oriented and shaped research in these fields over the last few decades.

1. Introduction

Pr Sébastien Candel has been a leader in the field of combustion for decades as well as a remarkable member of the Combustion Institute. This dedicated paper, written by some of his students and colleagues, presents a summary of his achievements in six fields of research, starting from his initial ideas and showing how these ideas have influenced multiple generations of researchers in combustion. While Pr Candel has also worked on steady flames, his main interest has always been unsteady combustion as found in turbulent flames, combustion noise and combustion instabilities. He has addressed all types of flames, from laminar Bunsen type flames (in the kW range) to cryogenic hydrogen-oxygen rocket engines (in the MW range):

* Cryogenic combustion and liquid rocket engine dynamics

* Turbulent combustion modeling

* Combustion noise

* Flame dynamics and combustion instability

* Flame describing functions in combustion instability analysis

* Combustion dynamics and instabilities of annular systems

For each topic, there are two objectives: (1) to provide a historical view of combustion research in each field so that readers will have an overview of recent progress, and (2) to show the impressive impact of Pr Candel in all of these fields. Each section is designed to provide a

summary of the field over the last forty years and to be used as an entry point into the domain for all interested researchers. Details can be found in the corresponding papers.

2. Cryogenic combustion and liquid rocket engine dynamics

Rocket engines are the most extreme application of combustion. They combine high pressures, huge powers, unusually low temperatures, and extremely high temperatures. The design of fuel and oxidizer injection systems is a crucial component in rocket engines. Early research into combustion in liquid rocket motors relied heavily on experiments [1,2]. By the 1970's, tests on many injector designs had led to adoption of the coaxial injector with a slow central liquid oxygen (LOx) jet surrounded by a fast annular fuel jet. This was used in the Space Shuttle main engine (SSME) with swirling flow at 200 bar, and in the Ariane family with non-swirling flow at 100 bar. The more complex design of the SSME was justified by the fact that each engine was expected to be re-used up to 50 times, although in practice each was used around 5 times.

By the standards of today, numerical analysis before 1990 was rudimentary [3]. Furthermore, experiments were limited to measurements of global outcomes such as combustion efficiency, thrust, and mixing rates because rigs did not have good optical access. Until the

<footnote><sup>a</sup> Ecole Centrale de Lyon, CNRS, Univ Claude Bernard Lyon 1, INSA Lyon, LMFA, UMR5509, 69130 Ecully, France</footnote>

<footnote><sup>b</sup> Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 91192 Gif-sur-Yvette cedex, France</footnote>

<footnote><sup>c</sup> Engineering Department, University of Cambridge, Trumpington Street, Cambridge, CB2 1PZ, United Kingdom</footnote>

<footnote><sup>d</sup> CAPS Laboratory, Department of Mechanical and Process Engineering, ETH Zürich, 8092, Zürich, Switzerland</footnote>

<footnote><sup>e</sup> Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS, Toulouse, France</footnote>

<footnote><sup>f</sup> Institut Universitaire de France (IUF), France</footnote>

<footnote>* Corresponding author.
E-mail address: thierry-poinsot@imft.fr (T. Poinsot).</footnote>
1990's, rocket engine design therefore relied on empirical correlations and simple axisymmetric models of the flows from fuel injectors [4]. These models assumed that shear layers grow linearly with downstream distance and considered the streams as quasi-1D incompressible steady flows. In this way, the shear experienced by the LOx jet could be estimated from the injector geometry, which went some way to explaining the experimental observations. Nevertheless, many details remained unknown. For example, scientists did not know the shape of the flame in a coaxial injector, how it was anchored, how the LOx jet was atomized, or how trans-criticality influenced jet break-up.

With advances in diagnostics, instrumentation, data acquisition, and numerical simulations available in the 1990's, it became possible to start to answer these questions systematically. In 1993 Sébastien Candel and others convened the Groupement de Recherche 'Combustion dans les Moteurs Fusées' to (i) investigate the physical processes involved in mixing and atomization, (ii) investigate the physical processes involved in combustion of sprays, (iii) develop experimental diagnostics applicable to cryogenic combustion, (iv) perform experimental measurements on a trans-critical cryogenic combustion facility and (v) develop numerical simulations of cryogenic combustion. One of the most instrumental tools was ONERA's MASCOTTE test rig: a high pressure (> 70 bar) coaxial test rig for cryogenic propellants with good optical access and the ability to be acoustically forced [5]. The information available from experiments on MASCOTTE enabled advances in models and simulations that would otherwise have been impossible. This section contains a description of the work that Sébastien Candel led from the 1990's, which has changed our understanding of the structure and stability of cryogenic spray flames.

## 2.1. Unanswered questions about cryogenic combustion in the early 1990's

Cryogenic combustion occurs between trans-critical species over such a wide range of spatial and temporal scales that it is too numerically expensive to simulate directly. Instead, reduced order models must be used for many aspects of the physics, meaning that experiments are required to test and improve models. Estimates of characteristic times [6] showed that chemical reaction is faster than mixing, which is faster than vaporization, which in turn is governed by atomization of the LOx jet, meaning that atomization is the rate limiting process. For atomization, cold flow experiments showed that the most influential dimensionless groups are the momentum flux ratio $J \equiv (\rho_g u_g^2)/(\rho_l u_l^2)$ (also known as the dynamic head ratio) and the Weber number $We \equiv (\rho_g(u_g - u_l)^2 d_l)/\sigma$, where g denotes gas, l denotes liquid, u is speed, ρ is density, d is jet diameter, and σ is surface tension. For supercritical flows, the Weber number is undefined because there is no surface tension and it was not clear what, if anything, should replace it. Despite the fast chemical time, the position of the base of the flame was also not known before the 1990's. Four possible flame anchoring configurations were thought possible [7], shown in Fig. 1. The main questions were therefore (i) what controls the rate of atomization or mixing, which impacts combustion efficiency and (ii) what controls the flame anchoring, which impacts reliability.

## 2.2. Flame structure and stabilization at subcritical pressures

The MASCOTTE test rig took up three rooms at ONERA. Students would spend one week setting up and testing optical diagnostics in the laser room adjacent to the outdoor test rig. Each test would begin by launching the optical diagnostics, confirming to Lucien Vingert that they were running, hearing Lucien mutter “c’est parti”, listening to the start-up sequence, then a click as the LOx and H₂ valves opened, a swoosh as the reactants poured into the chamber, and a deafening roar as they fired a 5 MW cloud of oxygen and water vapor into the Palaiseau forest.

The first published experimental images from MASCOTTE were of a single coaxial injector with a central liquid oxygen (LOx) jet

and a surrounding gaseous hydrogen ($GH_2$) annulus at 5 bar and 10 bar with no recess [6]. The flame was imaged with (i) planar laser-induced fluorescence (PLIF) of OH radicals, (ii) elastic scattering of the LOx jet, and (iii) spontaneous emission of OH*. Diagnostics (i) and (ii) were acquired simultaneously, showing that the flame front was attached to the LOx injector lip, which came as a surprise. It had been assumed previously that the flame was lifted and that some non-reacting mixing was taking place in the space between the injector and the flame, as would be the case for a classic triple flame. Furthermore, these experiments showed that the momentum flux ratio, $J$, had no influence on flame anchoring over the range examined and that the flame was a continuous sheet containing no discernible holes. This qualitative information became useful later for modeling and experimental diagnostics.

A systematic sweep of the two most influential parameters was performed [7]. These parameters were pressure (1, 5, 10 bar) and momentum flux ratio $J$ (~14, ~10, ~6.5). This necessitated varying the diameter ratio, the mass flux ratio, the velocity ratio, and the Weber number. As in [6] the flame was imaged with (i) OH PLIF, (ii) scattering from the LOx jet, and (iii) spontaneous OH* emission. This revealed that the flame remained anchored at all conditions, as expected. More importantly, it showed that the flame structure was similar at all pressures (which were all subcritical), that the flame shortened as the pressure increased, and that the flame diameter and flame intensity increased as the momentum flux ratio increased. This is consistent with the LOx jet breaking up and evaporating faster as the momentum flux ratio increases, and with the reaction rate being dominated by this process. This combination of theoretical reasoning combined with sparse but informative experiments greatly increased physical understanding of cryogenic rocket engine flames compared with previous tests. Qualitative visualization available through MASCOTTE was key. Further, it was possible to estimate the mean heat release rate per unit volume in the rig to be 4.5 GW/m³, which is an order of magnitude higher than typical hydrocarbon flames at atmospheric pressure. As expected, this value is similar to those in industrial rocket engines [8]. This large volumetric heat release rate explains why thermoacoustic instability is such a problem in rocket engines: the heat release rate is so high that, even if only a small portion is converted to work each cycle through the thermoacoustic mechanism, the acoustic energy can build up quickly.

Through the tests up to 10 bar, it became clear that PLIF and laser sheet scattering from the LOx jets were struggling to be effective at 10 bar and would not work well at higher pressures. Firstly, LIF imaging suffers from more rapid quenching of the excited molecules at higher pressures, leading to a lower signal. Secondly, LIF suffers from dispersion of the laser sheet due to refractive index gradients in the turbulent flow around the flame and jet. This dispersion increases at higher pressures. On the other hand, the natural OH* emission images remained informative at all pressures. The qualitative information that the flame was known to be a thin hollow shell now proved useful because it showed that it would be acceptable to perform computerized tomography (an Abel transform) of the OH* natural emission images in order to deduce the mean position of the flame.

If the radiating region of a flame is axisymmetric and transparent then a single side-on line-of-sight image contains sufficient information to deduce the radial emission profile [9,10] (Fig. 2). This is done with the Abel transform, which can be implemented numerically in several different ways, but is easiest to describe as ‘onion-peeling’: the contribution of the outer layer of emission is deduced from the outer pixels of the side-on image; the outer layer’s contribution to the adjacent (inner) layer of pixels is then subtracted from these pixel values to reveal the contribution of the adjacent inner layer of emission; this process is repeated to the centerline of the flame. Background emission needs to be eliminated, which is straightforward. Measurement noise accumulates during the onion-peeling process, meaning that the Abel-transformed images become less reliable towards the centerline.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Possible flame stabilization configurations: (a) the flame is anchored on the injector rim; (b) the flame is lifted and stabilized at a finite lift-off distance; (c) the flame is stabilized as a premixed front anchored in the liquid core break-up region; (d) a weak pilot reaction develops near the injector rim initiating reactions which after a certain distance run away and produce flame stabilization.</figcaption></figure>

Source: Reproduced from [7].

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Time average and Abel transform of OH\* emission at 10 bar, *J* = 6.5, no recess from Figs. 3b and 5b of [10].</figcaption></figure>

For several reasons, the Abel transform of natural OH\* emission is an ideal diagnostic given the cryogenic flame structure. Firstly, there is no laser sheet to be scattered by the high refractive index gradients in the flow and the LOx jet. Secondly, the accumulation of measurement noise towards the center is unimportant because most emission takes place in the outer shells and there is no emission at all from the LOx jet. Thirdly, any scattering of natural emission from the LOx jet does not influence the emission measurement from the outer shells. The only significant drawback is the axisymmetric assumption. The images in [7], which were then Abel transformed in [9], were instantaneous snapshots. Between 20 and 140 of these were averaged to obtain an approximate time average, but this image had to be treated with a low pass spatial filter in order to obtain an image that could be Abel-transformed without measurement noise exceeding the signal. Because of this limitation, later experiments on Mascotte used a longer exposure time such that the OH\* natural emission images were naturally time-averaged during acquisition.

### 2.3. Flame structure and anchoring at transcritical pressures

Experiments at subcritical pressures had shown that the flame sheet was continuous and was stabilized at the LOx post lip (Fig. 2). A series of careful 1D and 2D calculations on counterflow $GH_2/LOx$ diffusion

flames [11–14], showed that the extinction strain rate of this flame remains an order of magnitude higher than flow-induced strain rates, even when the flame is close to the liquid oxygen surface. This reduces the question of flame anchoring to a study of the zone behind the LOx post lip because the subsequent diffusion flame cannot be extinguished.

The next series of experiments on MASCOTTE [15] was at 70 bar pressure, which is greater than the critical pressure of $O_2$ at 50.5 bar. Natural OH\* emission processed with the Abel transform was used to identify the mean flame position. Backlighting (shadowgraphy) was used to identify the LOx jet position, quantified by the proportion of time that the LOx jet was in the line of sight of each pixel. The flame was visualized from 0 to 22 LOx jet diameters downstream of the injection plane to obtain the overall spray flame shape. A zoom of the injection region from 0 to 4 LOx jet diameters was visualized to examine flame holding on the LOx post in more detail (Fig. 3).

The ligaments and droplets that were visible in the experiments below 10 bar were replaced by packets of high density oxygen with indistinct edges at 70 bar. The 70 bar flames were several times more intense than the 1 to 10 bar flames but expanded more slowly with downstream distance and were less influenced by the momentum flux ratio, *J*. A comparison of typical timescales showed that, at 70 bar, large scale turbulent mixing becomes the slowest and most influential parameter. Below the critical pressure of oxygen, on the other hand,
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Transcritical cryogenic LOx/GH2 flame at 70 bar and *J* ≈ 11 from [15]. Color levels: Abel transform of the time-averaged OH* emissions. Gray levels: proportion of time that the LOx jet was present along the line of sight from 0% (transparent) to 100% (black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** OH PLIF images of the flame holding region, from [17]. LOx is injected below the lip and GH₂ is injected above the lip.</figcaption></figure>

atomization is the slowest and most influential parameter. The main conclusion is that, in combustion chambers operating above the critical pressure, one should attempt to increase mixing rates rather than improve atomization [16]. This detail aside, the supercritical flame structure was found to be similar to the subcritical flame structure.

Attention now turned to the flame anchoring mechanism because of its importance for reliability and, potentially, resistance to thermo-acoustic instability. By considering the canonical problem of two types of cross-flow flame, Juniper & Candel [18] showed that the stand-off distance of cross-flow flames controlled by the strain rate increases in proportion to the strain rate, while the stand-off distance of cross-flow flames controlled by convection increases in proportion to the cube of the flow velocity. In other words, a flame head exposed to a high velocity stream will blow off and should be avoided in a rocket engine.

The influence of heat transfer to, and evaporation from, a liquid oxygen stream was then considered [19], in which the evaporation velocity of the fuel is coupled to the heat release rate. A flame edge above a flat pool was found, like the cross-flow flames, to be sensitive to the incoming flow velocity. The addition of a step, modeling the LOx post, was found to anchor the flame, as long as the flame thickness was smaller than the LOx post thickness. The ratio of LOx post thickness to flame thickness was found to be the most important parameter for the stabilization of this type of cryogenic flame.

Although no systematic experiments were performed to test the above conclusion, two experimental campaigns were conducted on cryogenic flame stabilization. The first [15,20], which was published in journal form in [16], showed that the flame was indeed tucked behind the LOx post lip, adjacent to the LOx jet. A later study [17] showed that, contrary to prior expectations, OH PLIF could be used successfully up to 63 bar in GH₂/LOx flames, although not in high pressure CH₄/LOx flames [21]. PLIF worked even in the flame holding region a few millimeters from the LOx post lip, which was around 0.3 mm wide. Ref. [17] also showed, in higher resolution, that OH radicals from the H₂/LOx flame are tucked in behind the LOx post lip (Fig. 4). On the other hand, if H₂ is replaced with CH₄, the flame (as measured by OH radical location) is thicker and the anchoring point is blown downstream and oscillates back and forth, depending on local

flow conditions. The main conclusion is that the width of the LOx post lip needs to be greater than the flame thickness for the flame to be anchored.

When the hydrogen stream is replaced with methane, it is possible for both the fuel and the oxidizer to be transcritical. A systematic study at various pressures was performed by [22]. This revealed the influence of the criticality of the outer (fuel) stream. When the outer stream is subcritical, whether the fuel be GH₂ or GCH₄, the flame structure is the same as that shown previously. When the outer stream is transcritical, however, the flame structure is different. A flame surrounds the LOx core, as before, but a second flame surrounds the outer methane annulus, reacting with residual oxygen in the chamber.

## 2.4. Effect of transverse oscillations

Thermo-acoustic instability remains one of the greatest challenges facing liquid rocket engine designers [23]. It occurs when local heat release rate fluctuations in a combustion chamber are sufficiently in phase with local pressure fluctuations. This mechanism converts heat into work and, due to the vast heat release rate and low damping in rocket engines, it needs only to be slightly efficient in order to induce large amplitude acoustic vibrations. Great effort is made to avoid it.

One influential element of the thermo-acoustic feedback loop is the mechanism through which the acoustic waves modulate the heat release rate. In rocket engines, in which radial and azimuthal (rather than longitudinal) acoustic chamber modes tend to be most dangerous, this mechanism is through transverse oscillations of the flames. With this in mind, a set of experiments were performed on MASCOTTE in order to examine the influence of transverse oscillations on the break-up and atomization of LOx jets and their subsequent combustion with the fuel. In rocket engines, the fuel injectors sit within a showerhead of other injectors, so the interaction between jets could be influential. For this reason, three or five injectors in a line were studied. Backlighting, OH* natural emission, and photomultipliers were used.

The transverse oscillations were forced with a siren placed behind a nozzle on the side of MASCOTTE. This was designed to excite an acoustic resonance in the chamber, thereby causing high amplitude
<figure><img src="image_7.png" /><figcaption>Fig. 5. Instantaneous OH* images without forcing (left) and with transverse forcing (right) from [26].</figcaption></figure>

oscillations from relatively low amplitude excitation. This worked well in cold flow. In [24], particle image velocimetry (PIV) measurements of the flow during transverse forcing were compared with the acoustic mode shapes found with a Helmholtz solver (AVSP developed by Cerfacs). The acoustic amplitude was high when the system was forced close to its acoustic resonant frequency. Furthermore, modulation close to the most amplified hydrodynamic frequency in the jet augmented the jet spreading angle.

Transverse forcing around the resonant frequency was, however, less effective in hot flow. This was studied carefully in [25] for 3 flame and 5 flame configurations both experimentally and with low order models. This paper showed that the chamber resonance is sharp in cold flow but distributed in hot flow. This is because temperature fluctuations inside the chamber cause speed of sound fluctuations, and hence resonant frequency fluctuations within the chamber. Resonant bursts then appear during a firing, at moments when the chamber resonant frequency happens to coincide with the forcing frequency. These resonant bursts make quantitative interpretation difficult (see Fig. 5).

Nevertheless, valuable qualitative information could be obtained from hot fire tests. In [27], 3 coaxial fuel injectors placed side by side in MASCOTTE were forced transversally at 9 bar. The flame spread and heat release rate were enhanced by the forcing, and the wall temperature increased. This study was not able to discover the fundamental mechanisms leading to the strong coupling between acoustic forcing and heat release rate but a later experiment [26] on a 5 flame configuration provided more information. This experiment forced at 20% of the mean chamber pressure. At moderate amplitude forcing, secondary atomization is enhanced, augmenting the vaporization rate. At higher amplitude forcing, the primary liquid jet break-up is also enhanced. Both mechanisms increase the evaporation rate of liquid oxygen, which is the rate controlling factor below the critical pressure of oxygen, so the heat release rate is enhanced.

## 2.5. Large eddy simulations in non-reacting transcritical conditions

The complexity of cryogenic combustion, with steep density and temperature gradients, high velocities, and the need for detailed equations of state, made numerical simulation difficult. An early large eddy simulation (LES) [28] was remarkable at the time, and was extended further over the next decade [29,30]. A Reynolds-average Navier-Stokes (RANS) approach [31] with a $k-\epsilon$ model for turbulence, models for either transcritical mass transfer or primary and secondary atomization, and a flame surface density transport equation, yielded qualitatively correct mean flame structures, quantitatively correct heat release rates, and the correct dependence on parameters. Nevertheless, this approach was not able to model unsteady effects such as transverse forcing. Progress with LES during that decade, combined with

experimental data from MASCOTTE provided the motivation for further attempts to use LES as an investigative tool.

In a collaboration between EM2C and Edwards Air Force Base, a supercritical pressure nitrogen coaxial jet was visualized with backlighting and compared with LES of this jet from the AVBP solver of CERFACS equipped with the Peng-Robinson equation of state [32]. Supercritical nitrogen was fed into the inner and outer flows of the coaxial injector, with different temperatures and velocities in the inner and outer flows. The momentum flux ratio varied from 1 to 10. The flow patterns were qualitatively the same, with the central jet shortening as the momentum flux ratio increases. The flow was then forced transversally, both experimentally and numerically, and the flow patterns were found to be qualitatively the same. Transverse forcing shortened the jet. The experimental diagnostic, backlighting, is integrated over the line of sight and is therefore impossible to process quantitatively. The LES on the other hand can be probed quantitatively. From the LES it was observed that the outer jet responds strongly to external forcing but that the inner jet does not, probably because the chosen forcing frequency was far from the natural hydrodynamic frequency of the inner jet. A further study [33] used LES to examine the influence of transverse forcing of a supercritical jet injected into a hot environment. Perhaps the most important conclusion of these studies, however, was that AVBP could simulate supercritical jets, paving the way for the reacting simulations that followed.

In parallel to the above studies, LES were performed of LOx/GH4 flames under transcritical conditions with an infinitely-fast combustion model, as suggested by experimental results showing that the flame was attached to the injector lip [35]. A single operating point was considered, taking 70,000 CPU hours. At that time, the computational expense precluded an investigation of parameters. In [34] a single LOX/GH4 coaxial injector was considered at $Re > 10^5$, $J = 7.7$ and 114 bar. The LES showed the flame and flow structure in far greater detail than could be obtained from experiments (Fig. 6). This flame was then subjected to transverse acoustic forcing, which flattened the dense LOx jet into a plane whose normal was in the direction of the forcing. At low forcing frequencies, the jet flapped from side to side and the flame followed. At high forcing frequencies, the jet remained stationary while the flame flapped from side to side around it. In all cases, the contribution of each region to the thermoacoustic growth rate (the Rayleigh index) was calculated. This revealed extraordinary detail about the flame dynamics and what their contribution to thermoacoustic instability would have been if the acoustic oscillations had been self-excited. This was repeated for five aligned coaxial jets in [36], using 8.5 million CPU hours, and for two operating conditions of a showerhead of 42 coaxial injectors modeling DLR's BKD combustor in [37] (Fig. 7). In the latter study, the system was linearly stable but could be triggered to a stable limit cycle involving the first transverse mode of the combustion chamber. This mode was coupled between the $O_2$ feed system and the combustion chamber and its structure was the same as that found with a Helmholtz
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Top: Longitudinal slice of instantaneous axial velocity (blue to red) with a super-posed slice of density (deep blue to light blue) and the contour at which the O₂ mass fraction is 0.5. Bottom: Instantaneous iso contour of temperature (1500 K) colored by axial velocity (blue 0 ms⁻¹ to red 90 ms⁻¹), from [34]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Top: Longitudinal slice of instantaneous temperature, from LES of the BKD experiment at DLR Lampoldshausen [37].</figcaption></figure>

solver. The main point of these studies is that, with sufficient computing power, transcritical cryogenic combustion with real gas effects can be simulated with LES. Parametric studies are possible, limited only by computational expense.

## 2.6. Effect of recessing the LOx tube inside the hydrogen injector

During tests in the 1960s and 1970s it was discovered that combustion efficiency of rocket engines was enhanced if the LOx tube was recessed inside the surrounding fuel tube: the efficiency increases because a larger fraction of the oxygen burns before reaching the throat. The mechanism for this was not known but further tests in cold flow [2, p45], showed that recess of the LOx tube enhances mixing of the LOx jet with the surrounding flow. This mixing is maximal when the LOx post is recessed by 1 LOx jet diameter but reduces when the recessed length is greater or less than this.

An early study [4] proposed that, when the inner tube is recessed, the outer tube prevents the annular flow from expanding and thereby retains a higher speed for a longer distance. This increases the distance over which high shear can act on the inner flow, thereby increasing jet breakup and atomization. Although plausible, this mechanism predicts that mixing would continue to be enhanced as the recessed distance is increased, contrary to observations. An alternative explanation was proposed in [38] for 2D planar flows and in [39] for axisymmetric flows. They showed that the flow in the recessed portion of the injector is strongly absolutely unstable [40] to long wavelength helical (azimuthal

wavenumber 1) oscillations. When the recessed length is sufficiently long, this absolute instability triggers a global helical oscillation in the recessed portion, which enhances primary jet break up. This helical oscillation is visible in shadowgraphy images of MASCOTTE [38]. This mechanism explains the observation that the mixing efficiency starts abruptly as the LOx tube is recessed and does not increase when the LOx tube is recessed further. The same mechanism is found in gas turbine combustion chambers, causing the precessing vortex core [41,42].

# 3. Turbulent combustion modeling

## 3.1. Introduction

The interaction of turbulence and combustion is probably one of the most discussed issue in our community in the last sixty years, especially with the development of Computational Fluid Dynamics (CFD). Sébastien Candel introduced CFD for turbulent reacting flows and related modeling activities at Laboratoire EM2C in the end of the seventies, starting from an original approach devised by his professor at Caltech, Franck Marble, for non-premixed turbulent flames, the so-called Coherent Flame Model (CFM) deriving a balance equation for the flame surface density [43].

S. Candel and his coworkers extended the applicability range of this model and continuously improved its ingredients, intensively performing RANS simulations of turbulent combustion. In the 1990s, taking advantage of spectacular progress in computational power, modeling
<figure><img src="image_10.png" /><figcaption>Fig. 8. Schematic drawing of flame stretching and shortening as described in the Coherent Flame Model. From Marble and Broadwell [43].</figcaption></figure>

turbulent reactive flows naturally turned to LES. Once again, S. Candel suggested an original approach based on the flame thickening procedure proposed by Butler and O'Rourke [44,45] in another context and whose closure, based on flame surface wrinkling factor, is closely related to CFM concepts.

This section summarizes the main contributions in turbulent combustion modeling conducted or inspired by S. Candel. Sections 3.2–3.4 are devoted to the initial version and main developments of the CFM while Section 3.5 focuses on LES. Most of the work presented here is due to Sébastien Candel, his disciples and his students but many references are also given regarding the use and improvements of these models by other teams.

## 3.2. The original Coherent Flame Model (CFM) and primary extensions

The CFM was developed under the auspices of the Project SQUID, a cooperative program of basic research on Jet Propulsion supported by the US Office of Naval Research and administered by Purdue University [46,47] and, to the best knowledge of the authors, no further investigated by its creators, Frank E. Marble and James E. Broadwell. The model considers turbulent flames between coflowing fuel and oxidizer streams as sketched in Fig. 8. Its originality is to explicitly take into account the coherent structures evidenced few years sooner by Brown and Roshko [48] at a time where combustion models relied on reaction rates to turbulent mixing speeds [49–51]. Reaction layers are then identified to thin sheets separating fuel and oxidizer. These sheets are stretched by the flow, increasing the total flame area, and destroyed by “mutual annihilation”, i.e. the consumption of the reactant, fuel or oxidizer, between two adjacent flame sheets. “A formulation of the coherent flame model for the turbulent diffusion flame then requires (i) a model for inhomogeneous turbulence, including closure conditions; (ii) a model for flame surface distribution over the turbulent region, which leads to a corresponding reactant consumption and heat release; and (iii) calculation or measurement of the reactant consumption rate for a laminar diffusion flame that is undergoing strain in its own plane” [43].

The key ingredient of CFM is the balance equation for the flame surface density, i.e. the flame area per unit volume, $\Sigma$, scaling as the inverse of a length. This balance equation is built from the formal calculation by Batchelor [52] who examined the behavior of line and surface segments in homogeneous and isotropic turbulence, adding a consumption term describing the mutual annihilation. This term is written as proportional to the volume consumption rate of reactant i, $V_{D_i}\Sigma$, where $V_{D_i}$ is the mean volume consumption rate per unit of flame surface area of the ith species and inversely proportional to the distance between adjacent flame sheets estimated as proportional to $v_i/\Sigma$ where $v_i$ is the volume fraction of species i. For a single step irreversible $F + sO \rightarrow (1 + s)P$ reaction, where F, O and P denote

fuel, oxidizer and products, respectively, and s the mass stoichiometric coefficient, the $\Sigma$-balance equation is postulated as:

$$ \frac{\partial \Sigma}{\partial t} + \tilde{u}_k \frac{\partial \Sigma}{\partial x_k} = \frac{\partial}{\partial x_k} \left( \frac{v_t}{\sigma_\Sigma} \frac{\partial \Sigma}{\partial x_k} \right) + \epsilon_s \Sigma - \beta \left( \frac{V_{D_F}}{v_F} + \frac{V_{D_O}}{v_O} \right) \Sigma^2 \quad (1) $$

where $\tilde{u}_k$ denotes the mass-weighted mean component of the flow velocity in the direction $x_k$, $v_t$ the turbulent kinematic viscosity, $\sigma_\Sigma$ the turbulent Schmidt number for the flame surface density, $\epsilon_s$ the strain rate induced by turbulent motions and acting to the flame surface, while $\beta$ is the model parameter. The three right hand side terms correspond to the turbulent transport of the flame surface, modeled through a turbulent viscosity, the flame surface increase due to strain rate and its consumption by mutual annihilation.

The balance equation for mass-weighted mean fuel and oxidizer mass fractions, $\tilde{Y}_F$ and $\tilde{Y}_O$, respectively, are then, in conservative form (written in terms of volume fractions in [43]):

$$ \frac{\partial \bar{\rho} \tilde{Y}_F}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_k \tilde{Y}_F}{\partial x_k} = \frac{\partial}{\partial x_k} \left( \bar{\rho} \frac{v_t}{\sigma_F} \frac{\partial \tilde{Y}_F}{\partial x_k} \right) - \bar{\rho} V_{D_F} \Sigma \quad (2) $$

$$ \frac{\partial \bar{\rho} \tilde{Y}_O}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_k \tilde{Y}_O}{\partial x_k} = \frac{\partial}{\partial x_k} \left( \bar{\rho} \frac{v_t}{\sigma_O} \frac{\partial \tilde{Y}_O}{\partial x_k} \right) - \bar{\rho} V_{D_O} \Sigma \quad (3) $$

where $\bar{\rho}$ is the mean flow density, $\sigma_F$ and $\sigma_O$ fuel and oxidizer turbulent Schmidt numbers. The mean mass reaction rate of species i is written:

$$ \bar{\omega}_i = \bar{\rho} V_{D_i} \Sigma \quad (4) $$

and split into two contributions separating chemical aspects and turbulence: $\bar{\rho}V_{D_i}$, the reaction rate per unit of flame surface area, estimated from one-dimensional laminar flame and possibly including complex chemistry features, and the flame surface density $\Sigma$.

The model is then complemented by (i) a turbulence model. Marble and Broadwell chose the Saffman model [53] based on turbulent kinetic energy and vorticity, this latter quantity being linked to coherent flow motions, and (ii) a model for the reaction rates per unit of flame surface area, $V_{D_i}$. For one-dimensional steady-state counter flow laminar diffusion flames assuming one-step irreversible and infinitely fast chemical reaction and unit species Lewis numbers, these reaction rates are given by analytical expressions [43,54]:

$$ V_{D_F} = Y_F^\infty \sqrt{\frac{D\epsilon_s}{2\pi}} \left(\frac{\phi+1}{\phi}\right) \exp\left(-\left[\text{erf}^{-1}\left(\frac{\phi-1}{\phi+1}\right)\right]^2\right) \quad (5) $$

$$ V_{D_O} = Y_O^\infty \sqrt{\frac{D\epsilon_s}{2\pi}} (\phi + 1) \exp\left(-\left[\text{erf}^{-1}\left(\frac{\phi-1}{\phi+1}\right)\right]^2\right) \quad (6) $$

where $Y_F^\infty$ and $Y_O^\infty$ are the fuel and oxidizer mass fractions in the pure fuel and oxidizer streams, respectively, D the species mass diffusivity, $\epsilon_s$ the strain rate and $\phi = sY_F^\infty/Y_O^\infty$ the equivalence ratio. As expected, $V_{D_F}/V_{D_O} = 1/s$. The last model ingredient is the strain
rate $\epsilon_s$. Marble and Broadwell [43], looking for self-similar solutions in two-dimensional mixing layers, related it to the absolute value of the gradient in the transverse direction $x_2$ of the mean downstream velocity (direction $x_1$):

$$ \epsilon_s = \left| \frac{\partial \tilde{u}_1}{\partial x_2} \right| \qquad (7) $$

Sébastien Candel and his coworkers extended the CFM to turbulent premixed combustion to perform numerical simulations of a two-dimensional dump combustor [55]. Even recognizing that a detailed analysis of premixed laminar flames under strain is essential to the determination of the volume rates of consumption per unit of flame surface area, the authors, facing the lack of information about this point, decided to use constant fuel ($V_{D_F}$) and oxidizer ($V_{D_O}$) consumption rates of the same order as the corresponding laminar flame speed. The strain rate $\epsilon_s$ is assumed proportional to the rate of strain of the mean flow:

$$ \epsilon_s = \alpha (\tilde{S}_{ij} \tilde{S}_{ij})^{1/2} \qquad (8) $$

where $\alpha$ is a model constant, set here to unity, and:

$$ \tilde{S}_{ij} = \frac{1}{2} \left( \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial \tilde{u}_j}{\partial x_i} \right) \qquad (9) $$

Transient ignition is mimicked setting a Gaussian shaped initial flame surface density.

Two alternatives are proposed in [56,57]. The first is to rely on the strain rate $\epsilon_s$ to the inverse of an integral turbulence time scale $\tau_t = k/\epsilon$, where $k$ and $\epsilon$ denote the turbulent kinetic energy and its dissipation rate, respectively:

$$ \epsilon_s = \alpha \epsilon / k \qquad (10) $$

The second evolution recasted the mutual annihilation of flame surface as:

$$ \beta \frac{\dot{q}}{(-\Delta H_F^0) \tilde{Y}_F} \Sigma^2 \qquad (11) $$

where $\dot{q}$ is the heat release rate per unit of flame surface area, estimated from laminar flame calculations, and $(-\Delta H_F^0)$ the heat release per unit mass of fuel.

When flame elements are identified to one-dimensional unstrained premixed laminar flames propagating at the laminar flame speed $S_L^0$, the fuel consumption term may be recast as:

$$ \bar{\omega}_F = -\rho_0 S_L^0 Y_F^0 \Sigma \qquad (12) $$

where $\rho_0$ and $Y_F^0$ denote the density and fuel mass fraction in fresh gases, respectively. The flame surface density consumption term is reformulated as:

$$ \beta \frac{S_L^0 Y_F^0}{\tilde{Y}_F} \Sigma^2 = \beta \frac{S_L^0}{1 - \tilde{\theta}} \Sigma^2 \quad \text{or alternatively} \quad \beta \frac{S_L^0}{\tilde{\theta}(1 - \tilde{\theta})} \Sigma^2 \qquad (13) $$

where $\tilde{\theta} = 1 - \tilde{Y}_F/Y_F^0$ is the mean progress variable, defined such as $\tilde{\theta} = 0$ is fresh gases and $\tilde{\theta} = 1$ in fully burnt products. The right hand side proposal in Eq. (13) ensures that $\Sigma$ vanishes in pure reactants and products. Retaining source and sink terms given by Eqs. (10) and (13) leads to the $\Sigma$-balance equation:

$$ \frac{\partial \Sigma}{\partial t} + \tilde{u}_k \frac{\partial \Sigma}{\partial x_k} = \frac{\partial}{\partial x_k} \left( \frac{v_t}{\sigma_\Sigma} \frac{\partial \Sigma}{\partial x_k} \right) + \alpha \frac{\epsilon}{k} \Sigma - \beta \frac{S_L^0}{\tilde{\theta}(1 - \tilde{\theta})} \Sigma^2 \qquad (14) $$

while the progress variable source term is $\bar{\omega}_\theta = \rho_0 S_L^0 \Sigma$. This simple formulation leads to two comments:

* Assuming an equilibrium between source and sink terms in Eq. (14) gives:

$$ \bar{\omega}_{\theta} = \rho_0 S_L^0 \Sigma = \rho_0 \frac{\alpha \epsilon}{\beta k} \tilde{\theta} (1 - \tilde{\theta}) = \rho_0 \frac{\alpha \epsilon}{\beta k} \frac{\tilde{Y}_F^0}{Y_F^0} \left( 1 - \frac{\tilde{Y}_F^0}{Y_F^0} \right) \qquad (15) $$

recovering the well-known Eddy-Break-Up (EBU) model [50,54, 58]. Accordingly, the CFM is coherent with the EBU model corresponding to an equilibrium limit but is able to handle transient and non-equilibrium situations through its flame surface density balance equation.

* Multiplying each term in Eq. (14) by $\rho_0 S_L^0$ assumed constant gives a balance equation for the mean progress variable reaction rate $\bar{\omega}_\theta$:

$$ \frac{\partial \bar{\omega}_\theta}{\partial t} + \tilde{u}_k \frac{\partial \bar{\omega}_\theta}{\partial x_k} = \frac{\partial}{\partial x_k} \left( \frac{v_t}{\sigma_\Sigma} \frac{\partial \bar{\omega}_\theta}{\partial x_k} \right) + \alpha \frac{\epsilon}{k} \bar{\omega}_\theta - \beta \frac{\bar{\omega}_\theta^2}{\rho_0 \tilde{\theta}(1 - \tilde{\theta})} \qquad (16) $$

This CFM version is not sensitive to the consumption rate per unit of flame surface area, $\rho_0 S_L^0$. Increasing $S_L^0$ will decrease $\Sigma$ in the same proportion and vice versa, preserving $\bar{\omega}_\theta$ values. This sensitivity can be achieved only through initial conditions, non-uniform reaction rates per unit of flame surface area or more sophisticated expressions for $\Sigma$-balance equation source and sink terms.

In parallel, the CFM was also used in the Candel's group for non-premixed turbulent combustion under infinitely fast chemistry assumption in mixing layers following the original Marble and Broadwell formulation [59] or modeling the turbulent strain rate with Eq. (10) [60]. A phenomenological multi-surface model to describe mixing between reactants, premixed and non-premixed combustion in the same flow was proposed to extend the CFM to situations such as the transient ignition of non-premixed flames [61,62].

### 3.3. Investigating key ingredients of the Coherent Flame Model

S. Candel and his coworkers then focused on two key ingredients of the CFM: (i) its ability to incorporate detailed chemistry features, separating chemistry details ($V_D$) from the flame/turbulence interactions ($\Sigma$), thanks to Eq. (4); (ii) the mutual annihilation of adjacent flame fronts in a vortex, described by the consumption term in the flame surface density balance equation.

#### 3.3.1. Reaction rate per unit of flame surface area

Simple expressions of reaction rates per unit of flame surface area are analytically given in non-premixed flames for one-step chemistry and unity Lewis numbers (Eqs. (5) and (6)) or can be linked to the laminar flame speed in premixed combustion ($-\rho_0 S_L^0 Y_F^0$ in Eq. (12)), but complex chemistry details may be easily taken into account. Further investigations require the selection of a model configuration for laminar flame elements. One of the most convenient is counterflow flames (Fig. 9). Two flows, fuel and oxidizer for non-premixed flames, fresh gases against fresh or burnt gases in premixed combustion, are injected in the opposite directions. This configuration is rather simple and the flow is steady, allowing experiments [63], the flame strain rate is directly adjusted through flow rates, while one-dimensional calculations, along the flow injection axis, are possible. For these reasons, such configurations have been widely studied for a long time [64,65] and still are today [66]. A complete review of these studies is out of the scope of this contribution and the reader is referred to some review papers [66–68].

To improve CFM and incorporate complex chemistry features in simulations, S. Candel and his coworkers investigated extensively counterflow laminar flames with detailed chemical schemes, both in premixed [69–74] and non-premixed [75–79] configurations, including spray combustion [80], soot precursor formation [81] or hydrogen/oxygen transcritical conditions for rocket engines [12,13,82].

Fig. 10 displays evolution of the heat release as a function of strain rate in counterflow premixed flames between propane/air fresh mixture and burnt gas streams (Fig. 9 middle) [72]. A chemical scheme based on 33 species and 126 reactions proposed by Warnatz (private communication) is considered. The curves are monotonic when the burnt gas temperature $T_b$ exceeds 1530 K and exhibit the so-called
<figure><img src="image_11.png" /><figcaption>Fig. 9. Counterflow flames generally selected to model laminar flame elements in the Coherent Flame Model. Left: fuel and oxidizer flows for diffusion flame; center: counterflow of fresh and burned gases (issued from a previous combustion); right: double flame in a counterflow of fresh gases.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Heat release rate versus strain rate in counterflow premixed flames between propane/air fresh mixture and burnt gas streams (see Fig. 9, middle). Left: equivalence ratio $\phi = 0.75$ and burnt gas temperature $T_b = 1450$ K; Right: $\phi = 0.75$, different values of the burnt gas temperature, from $T_b = 1400$ to $T_b = 1965$ K. From [72].</figcaption></figure>

S-shape below this value. In this later situation (Fig. 10, left), the upper branch corresponds to standard flames, the lower branch to extinguished flames while the middle branch is unstable. Extinction when increasing strain takes place at B while ignition occurs at D when decreasing strain rate. The occurrence of extinction under excessive strain rates suggested to add a related consumption term, $D_{ext}$, in the flame surface density balance equation Eq. (1) [56,57]:

$$D_{ext} = -\gamma_e \epsilon_s h (\epsilon_s - \epsilon_s^{ext}) \Sigma \quad \text{or} \quad D_{ext} = -\gamma_e (\epsilon_s - \epsilon_s^{ext}) h (\epsilon_s - \epsilon_s^{ext}) \qquad (17)$$

where $\gamma_e$ is a model constant and h the Heaviside function ($h(x) = 0$ when $x < 0$; $h(x) = 1$ for $x \ge 0$). However, these formulations were not used in practice because strain rates predicted in turbulent flames are often larger than extinction strain rates $\epsilon_s^{ext}$ of laminar flames. Three types of reasons can be invoked. First, the model to estimate strain rates $\epsilon_s$ in turbulent flows, generally related to mean velocity gradients (Eqs. (7) and (8)) or to the inverse of a turbulent time scale (Eq. (10)) might be deficient. A second point lies in the orientation factors (see Eqs. (20) and (24) in Section 3.4.1 below): $\epsilon_s$ corresponds to the strain rate acting in the plane tangential to the flame surface, not to the overall strain rate induced by the flow field. Moreover, detailed analysis shows that laminar flames respond to strain rate changes with a delay of the order of $1/\epsilon_s$ and act as low pass filters (see [54,75,83] for diffusion flames and [84] for premixed flames) so that flames can sustain strain rates exceeding extinction strain rates as long as they are not applied during a too long time. Haworth et al. [83] incorporated this response time through a balance equation for the effective strain rate seen by the flame to be used to estimate the reaction rate per unit of flame surface area from expressions such as Eqs. (5) and (6), established for steady flames.

To compute laminar flame structures in view of integration in turbulent flame models such as CFM leaded later naturally to tabulation methods (see review paper [85]) where reaction rates are

estimated as functions of progress variable, mixture fraction (or equivalence ratio), strain rates, heat losses, etc. Simultaneously, Gicquel et al. [86] (FPI: Flame prolongation of Intrinsic Low Dimensional Manifold) and Oijen et al. [87] (FGM: Flamelet generated manifold) tabulated one-dimensional laminar premixed flames. This approach was then extended to partially premixed and non-adiabatic flames. [88–93]. The size of databases might become prohibitive in terms of computational resources, but can be reduced using flame self-similarity properties [94–97].

### 3.3.2. Flame in a vortex

According to Fig. 8, a configuration of interest is the rolling of a flame in a vortex. This configuration is well-suited to study the straining of flame in mixing layer, the increase in flame surface area due to strain as well as mutual annihilation of adjacent flame elements. S. Candel and his coworkers [98–101] investigated this situation using numerical simulations to assess the scaling laws theoretically predicted by Marble and coworkers [102–105] and compare results to prediction of the CFM [59]. The agreements were excellent and confirmed the relevance of the CFM concept, even if simulations were limited by strong assumptions such as infinitely fast chemistry for diffusion flame or thermodiffusive approximation (no action of the flame to the flow field, due to a constant density assumption), because of the very limited computational resources available at that time.

### 3.4. Theoretical investigations

The CFM was initially developed intuitively and the $\Sigma$-equation (Eq. (1) or Eq. (14)) was essentially postulated without rigorous derivation. In the 1990s, S. Candel revisited the whole framework of flame surface model using theoretical analysis, combined with direct numerical simulation (DNS) data, briefly summarized here.
### 3.4.1. Balance equation for the flame surface density

Candel and Poinsot [106] derived an exact expression for the flame stretch:

$$ \phi_s = \frac{1}{\delta A} \frac{d(\delta A)}{dt} = \nabla \cdot \mathbf{u} - \mathbf{n}\mathbf{n} : \nabla \mathbf{u} + S_L \nabla \cdot \mathbf{n}. \quad (18) $$

where $\delta A$ denotes a flame surface element propagating at the speed $S_L$ in the normal direction $\mathbf{n}$. The two first RHS terms, $\nabla \cdot \mathbf{u} - \mathbf{n}\mathbf{n} : \nabla \mathbf{u}$, correspond to the strain rate, i.e. the velocity gradient in the flame tangential plane, while the last term is related to the flame curvature effects and may be expressed in terms of the principal radii of curvature of the surface, $R_1$ and $R_2$:

$$ S_L \nabla \cdot \mathbf{n} = S_L \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \quad (19) $$

and is positive (negative) when the flame is convex (concave) towards fresh gases. Then, an exact balance equation for the instantaneous local flame area per unit volume, $\Sigma_l$ is deduced:

$$ \frac{\partial \Sigma_l}{\partial t} + \nabla \cdot \mathbf{u} \Sigma_l = (\nabla \cdot \mathbf{u} - \mathbf{n}\mathbf{n} : \nabla \mathbf{u}) \Sigma_l - \mathbf{n} \cdot \nabla (S_L \Sigma_l) \quad (20) $$

The authors suggested that considering the flame surface area per unit mass, $a_s$, is more appropriate for simulations. As $\Sigma_l = \rho a_s$, the corresponding balance equation is:

$$ \frac{\partial (\rho a_s)}{\partial t} + \nabla \cdot (\rho \mathbf{u} a_s) = \rho (\nabla \cdot \mathbf{u} - \mathbf{n}\mathbf{n} : \nabla \mathbf{u}) a_s - \mathbf{n} \cdot \nabla (\rho S_L a_s) \quad (21) $$

Pope [107] showed that the mean flame surface density can only be defined statistically. For a premixed flame described by a progress variable $\theta$ (reduced mass fraction or temperature such as $\theta = 0$ in fresh gases and $\theta = 1$ in fully burnt ones), the exact definition of $\Sigma^*$ is:

$$ \Sigma^* = \overline{|\nabla\theta| \delta(\theta - \theta^*)} = \left(\overline{|\nabla\theta| |\theta = \theta^*}\right) \rho(\theta^*) \quad (22) $$

where $\theta = \theta^*$ identifies the flame surface. The Dirac delta function $\delta(\theta - \theta^*)$ gives a local measure of the probability to be on the $\theta = \theta^*$ isosurface, $(|\nabla\theta| |\theta = \theta^*|$) is the conditional average of $|\nabla\theta|$ for $\theta = \theta^*$ and $\rho(\theta^*)$ is the probability to find $\theta = \theta^*$ at the given location. Recasting the $\theta$-balance equation under a propagative form:

$$ \frac{\partial \theta}{\partial t} + \mathbf{u} \cdot \nabla \theta = \frac{1}{\rho} \left[ \nabla \cdot (\rho D_{\theta} \nabla \theta) + \dot{\omega} \right] = S_d |\nabla \theta| \quad (23) $$

where $S_d$ denotes the displacement speed of the $\theta$-isosurface in its normal direction relatively to the flow field and $D_{\theta}$ the diffusivity of $\theta$, a balance equation for the mean (or filtered in LES) flame surface density can be derived [107–112]:

$$ \frac{\partial \Sigma^*}{\partial t} + \nabla \cdot (\langle \mathbf{u} \rangle_s^* \Sigma^*) + \nabla \cdot \left[ \langle S_d^* \mathbf{n} \rangle_s \Sigma^* \right] = \langle \nabla \cdot \mathbf{u} - \mathbf{n}\mathbf{n} : \nabla \mathbf{u} \rangle_s^* \Sigma^* + \langle S_d^* \nabla \cdot \mathbf{n} \rangle_s^* \Sigma^* \quad (24) $$

where $S_d^*$ is the $\theta^*$-surface displacement speed. The surface averaged $Q$-quantity, $\langle Q \rangle_s^*$, is:

$$ \langle Q \rangle_s^* = \frac{\overline{Q |\nabla c| \delta(\theta = \theta^*)}}{|\nabla c| \delta(\theta = \theta^*)} = \left( \frac{\overline{Q |\nabla c| |\theta = \theta^*}}{|\nabla c| |\theta = \theta^*|} \right) \quad (25) $$

The notation $\langle Q \rangle_s^*$ indicates that the average is taken for $\theta = \theta^*$ and is a function of $\theta^*$. The three terms in LHS of Eq. (24) are, respectively, the unsteady term, the convection of the flame surface by the flow field, that can be split into mean ($\tilde{\mathbf{u}}$) and unresolved ($\mathbf{u}''$) contributions using $\langle \mathbf{u} \rangle_s^* = \tilde{\mathbf{u}} + \langle \mathbf{u}'' \rangle_s^*$, and the flame propagation at the displacement speed $S_d^*$. The two RHS terms correspond to the action of stretch, decomposed in strain rate and curvature effects (Eq. (18)). The strain rate is often split into a contribution due to the mean flow and a contribution due to turbulent velocity fluctuations:

$$ \langle \nabla \cdot \mathbf{u} - \mathbf{n}\mathbf{n} : \nabla \mathbf{u} \rangle_s^* = \langle \nabla \cdot \tilde{\mathbf{u}} - \langle \mathbf{n}\mathbf{n} \rangle_s^* : \nabla \tilde{\mathbf{u}} \rangle + \langle \nabla \cdot \mathbf{u}'' - \mathbf{n}\mathbf{n} : \nabla \mathbf{u}'' \rangle_s^* \quad (26) $$

The flame surface density concept was built for infinitely thin flame fronts but can be extended to non-infinitely thin flames, integrating

Eq. (22) across isosurface levels [113,114]:<fnref n="1" />

$$ \Sigma = \int_{0}^{1} \Sigma^{*} d\theta^{*} = \int_{0}^{1} \left( \overline{|\nabla\theta| |\theta = \theta^{*}|} \right) p(\theta^{*}) d\theta^{*} = \overline{|\nabla\theta|} \quad (27) $$

that may be viewed as a “generalized flame surface density”, following a balance equation similar to Eq. (24), replacing surface averages $\langle Q \rangle_s^*$ by “generalized surface averages” $\langle Q \rangle_s$ which do not depend on $\theta^*$:

$$ \langle Q \rangle_s = \frac{1}{|\nabla\theta|} \int_0^1 \langle Q \rangle_s^* \Sigma^* d\theta^* \quad (28) $$

Transport by unresolved flow motions, flame propagation, strain rate and curvature contributions in Eq. (24) require modeling [54,58, 108,111,112,114,115]:

* The convection of flame surface by unresolved flow motions is generally closed assuming a usual gradient turbulent transport even if Bidaux and Bray (unpublished contribution, 1994) showed that this transport may be of counter-gradient type [116].

* The strain rate due to the mean flow, first RHS term in Eq. (26), requires only the modeling of the orientation factors [108,115].

* The strain rate due to unresolved flow motions, last term in Eq. (26), is generally closed as the inverse of a turbulent time scale, either the integral Eq. (10) or the Kolmogorov [108] time scales.

* The curvature term is usually modeled as a sink term but there is evidence that this term can be positive in the fresh gas side of the turbulent flame brush [111,112,115].

Eqs. (20) and (24) have provided a theoretical support initiated by S. Candel to a nice idea primarily based on intuitive arguments. DNS of reacting flows were developed to investigate various aspects of this concept [109,112,113,115,117–127]. Few studies discussed the extraction of flame surface densities and related quantities from two-dimensional experimental data [128–132].

### 3.4.2. Efficiency functions

Following flame/vortex studies by S. Candel and coworkers discussed in Section 3.3.2, Poinsot et al. [117] investigated laminar premixed flame / vortex interactions using realistic DNS. An example of results is plotted in Fig. 11 for a case where the vortex pair size and speed are high enough to induce quenching of the flame front. Results for a set of simulations varying vortex characteristics are displayed in Fig. 12 in terms of vortex intensity and size, compared to laminar flame speed and thickness, respectively. Four situations are identified, depending on velocity and length ratios: no effects, wrinkled flame front, formation of fresh gas pockets in burnt gases and flame quenching. This figure shows clearly that small and/or low intensity vortices have no effects on the flame while large and/or intense vortices induce flame quenching. These results were primarily used to revisit turbulent premixed combustion diagrams [54,117], showing that the flamelet regime domain, where the turbulent flame may be identified as a collection of laminar flame elements, sustaining CFM concept, is larger than expected from a simple analysis comparing orders of magnitude in Borghi-type diagrams.

Later, Meneveau and Poinsot [133] used these DNS results to estimate the strain rate entering the flame surface density balance equation Eq. (24). The first step is to extract from DNS the vortex efficiency, C,

<footnote n="1">¹ This generalized flame surface density may be of interest when extracting information from direct numerical simulations resolving the inner flame structure. To identify the flame to a given $\theta^*$-isosurface might be misleading when the behavior of this surface is not representative of the overall behavior of the flame, for example when the local flame thickness is evolving, possibly leading to unexpected negative displacement speeds $S_d^*$ [54].</footnote>
<figure><img src="image_13.png" /><figcaption>Fig. 11. Temperature (left) and reaction rate (right) fields at four instants during flame vortex interaction [54,117].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Spectral diagram summarizing DNS results of flame vortex interaction [54, 117].</figcaption></figure>

defined as:

$$ \frac{1}{A} \frac{dA}{dt} = C \left( \frac{r}{\delta_L^0} \right) \frac{u'(r)}{r} \qquad (29) $$

where A is the flame surface, r and u'(r) the vortex size and velocity, δ_L^0 the flame thickness and u'(r)/r the theoretical stretch acting on the flame. The RHS term in the above equation is the effective stretch due to the (r, u'(r)) vortex. The idea is then to integrate these stretches, assumed cumulative, over all the turbulence eddies from the Kolmogorov to the integral (I_t) length scales using multifractal theories [134,135], leading to the Intermittent turbulent net flame stretch (ITNFS) model:

$$ \epsilon_s = \alpha \Gamma_s \left( \frac{u'}{S_L^0}, \frac{l_t}{\delta_L^0} \right) \frac{\epsilon}{k} \qquad (30) $$

where u' is the turbulent RMS velocity in fresh gases, α a model constant. Γ_s is the global efficiency function fitted from DNS and displayed in Fig. 13, evidencing the reduced ability of small vortices to wrinkle the flame front. This model provides significant improvements to the CFM [58]. It can also be combined with Eddy-Break-Up (EBU) and Bray-Moss-Libby (BML) models [136,137] to introduce a sensitivity to chemistry and decrease the mean reaction rate in highly strained regions when turbulent scales are small.

This efficiency function was later revisited and improved for LES by Colin et al. [138] and Charlette et al. [139] while Bougrine et al. [140]

<figure><img src="image_15.png" /><figcaption>Fig. 13. ITNFS efficiency function Γ, Eq. (30) as a function of the length scale ratio l_t/δ_L^0 for several values of the velocity ratio u'/s_L^0 (0.1; 1.0; 10; 100) [54].</figcaption></figure>

introduced Lewis number effects. Thiesset et al. [141] proposed a refined analysis combining numerical simulations and experimental data. Suillaud et al. [127] derived a simple analytical expression assuming an Heaviside function for the vortex efficiency function C in Eq. (29).<fnref n="2" />

### 3.4.3. Kolmogorov-Petrovski-Piskunov (KPP) analysis

The Kolmogorov-Petrovski-Piskunov (KPP) analysis is a simple, powerful theoretical tool providing analytical expressions for turbulent flame speed. It is conducted under some restrictive assumptions such as homogeneous and frozen turbulence (turbulent flow field is not affected by combustion) and can only be derived for some combustion models. However, it is an efficient way to investigate basic model trends. The idea of the KPP analysis is to search an exponential solution for θ and Σ at the flame leading edge, keeping first order terms (see details in [54,58,142,143]). Fichot et al. [143] used this approach to analyze the behavior of the CFM and showed that the analytical expression for the flame speed derived by the KKP analysis agrees very well with numerical simulations conducted under the same assumptions. Duclos et al. [58] compared various formulations of the flame surface density balance equation issued from the literature.

<footnote n="2">² In LES, the integration over vortex size is limited to the larger subgrid scale contribution, i.e. to the filter size Δ. Of course, the corresponding efficiency function can be used in RANS, replacing Δ by the integral length scale l_t.</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Modeling approaches for turbulent combustion: geometrical analysis (flame front is viewed as a surface), mixing description (reaction rate is assumed to be controlled by the mixing rate, i.e. by scalar dissipation rate) or pure statistical approach (probability density function formalism) [114]. From [54].</figcaption></figure>

### 3.4.4. Relations between models

Beginners are often surprised by the very different nature of turbulent combustion models. Despite their differences, models can be classified in three categories (Fig. 14) [114]:

*   **Turbulent mixing:** assuming that chemical time scales are shorter than turbulent time scales (large Damköhler numbers), reaction rates are controlled by turbulent mixing rate, usually expressed in terms of scalar dissipation rates $\tilde{\chi}$, defined as:

$$ \bar{\rho} \tilde{\chi} = \bar{\rho} D \nabla \theta \cdot \nabla \theta \quad (31) $$

where $\theta$ is a progress variable or a mixture fraction in premixed and non-premixed combustion, respectively. Simple and popular pioneering models such as Eddy-Break-Up [50] and Eddy-Dissipation Concept [51] are derived from this analysis.

*   **Geometrical description:** the flame front is identified as a surface evolving in the turbulent flow field. This surface may be related to the global flame brush as done in the level set G-equation framework [144–146] but is often linked to the instantaneous isosurface of a mass or mixture fraction as in the CFM described in this paper.

*   **One-point statistics:** mean reaction rates are expressed combining instantaneous reaction rates given from Arrhenius laws, $\dot{\omega}_k(\Psi_1, \dots, \Psi_N)$ with the joint probability density function $p(\Psi_1, \Psi_2, \dots, \Psi_N)$ to have given values of the thermochemical variables (i.e. species mass fractions, temperature...) $\Psi_1, \Psi_2, \dots, \Psi_N$:

$$ \bar{\omega}_k = \int_{\Psi_1, \dots, \Psi_N} \dot{\omega}_k(\Psi_1, \dots, \Psi_N) p(\Psi_1, \Psi_2, \dots, \Psi_N) d\Psi_1 d\Psi_2 \dots d\Psi_N \quad (32) $$

This expression may be recast in terms of conditional means:

$$ \bar{\omega}_k = \int_0^1 \left( \overline{\dot{\omega}_k(\Psi_1, \dots, \Psi_N) | \theta = \theta^*} \right) p(\theta = \theta^*) d\theta^* \quad (33) $$

where the conditional quantity, $\theta$, is generally either a progress variable or a mixture fraction.

The Conditional Moment Closure (CMC) model [147] may be viewed either as a probability density function or a multi-surface model.

During the 1980s, vigorous debates took place in the combustion community regarding these concepts and modeling approaches. The flame surface density concept, initially built on intuitive ideas, took more time to find its place. Later, it was shown that previous approaches, even based on different physical concepts, are closely related [114,148].

<caption>**Table 1**</caption>

Exact expressions for the averaged (or filtered) reaction rate as a function of the mean scalar dissipation rate, $\tilde{\chi}$ (Eq. (31)), probability density function, $p(\theta^*)$ or surface density, $\Sigma^*$ (Eq. (22)) of iso-surface $\theta = \theta^*$, $\theta$ being a progress variable (premixed combustion) or a mixture fraction (non-premixed combustion). $\dot{\omega}$ is the local instantaneous reaction rate, $\langle Q \rangle_s^*$ a surface average (Eq. (25)), $\langle Q \rangle_s$ corresponds to generalized surface average and is defined in Eq. (28).

Source: From [114].

<table><thead><tr><th>Tool</th><th>Averaged reaction rate</th></tr></thead><tbody><tr><td>Scalar dissipation rate</td><td>$$ \bar{\rho} \frac{\int_0^1 \left\langle \frac{\dot{\omega}}{|\nabla\theta|} \right\rangle_s^* d\theta^*}{\int_0^1 \langle (\rho D |\nabla\theta|) \rangle_s^* d\theta^*} = \tilde{\chi} = \bar{\rho} \frac{\left\langle \frac{\dot{\omega}}{|\nabla\theta|} \right\rangle_s}{\langle (\rho D |\nabla\theta|) \rangle_s} \tilde{\chi} $$</td></tr><tr><td>Probability density function</td><td>$$ \int_0^1 \dot{\omega}(\theta^*) p(\theta^*) d\theta^* $$</td></tr><tr><td>Flame surface density</td><td>$$ \int_0^1 \left\langle \frac{\dot{\omega}}{|\nabla\theta|} \right\rangle_s^* d\theta^* = \left\langle \frac{\dot{\omega}}{|\nabla\theta|} \right\rangle_s \frac{|\nabla\theta|}{|\nabla\theta|} $$</td></tr></tbody></table>

Table 1 gives the exact expressions of mean (RANS) or filtered (LES) reaction rates in terms of scalar dissipation rate, probability density function and flame surface density [114]. These relations are useful to understand the implications of the physical hypothesis underlying the model developments and can be combined to devise hybrid models [149,150].

No approach displays a decisive advantage over others and results mainly depend on closure models: when the flame behavior is mainly governed by flame propagation, the flame surface density concept is easier to close, while scalar dissipation rate or probability density functions are more suited to describe mixing controlled combustion. The CFM and its variants have been widely used outside Candel's team, mainly for industrial applications such as internal combustion engines or gas turbines [108,151–172].

### 3.5. Large eddy simulations

#### 3.5.1. Introduction

At the end of the 1990s, thanks to the continuous increase of computational power, LES emerged as a very efficient tool for combustion for three reasons. First, LES, describing the behavior of the largest turbulent motions, gives access to the flow dynamics, of primary importance to predict and study combustion instabilities in burners or cycle-to-cycle variations in internal combustion engines. The second point is that LES models only the effects of smallest (unresolved) scales, expected to behave in a more universal way, reducing the modeling effort (at the expense of computational costs). Moreover, LES identifies fresh and burnt gas zones at the resolved scale level. These zones behave differently in terms of turbulence, radiative heat transfer or pollutant emissions and can be better described.<fnref n="3" /> For these reasons, very good results are achieved even with rather simple sub-grid scale models.

There is no difficulty to extend the flame surface density concept of Section 3.4 to LES and derive a balance equation for the filtered flame surface density. Most balance equations hold, replacing mean values used in RANS by filtered quantities in LES [54,111–113,123,173]. However, the flame surface density equation, needs closures and is rarely used in practice [112,174–182]: (i) good results are achieved with algebraic models where the filtered flame surface density or

<footnote n="3">³ The main drawback of Reynolds averaged Navier-Stokes simulations (RANS) lies in the turbulence model. Most combustion models are too sophisticated compared to their connection with the turbulent model through a single integral time scale.</footnote>
the sub-grid scale flame wrinkling factor are expressed from known quantities; (ii) solving this equation increases costs significantly; (iii) LES give access to resolved scales that can be used to model sub-grid scale contributions using scale-similarity models [183,184] or dynamic formalism as described below.<fnref n="4" />

Flame fronts are generally much thinner than usual computational mesh sizes, making combustion mainly a sub-grid scale combustion phenomenon [54,186].<fnref n="5" /> This difficulty can be circumvented using a flame front tracking technique [144,146,187] or filtering with a Gaussian filter larger than the mesh size [113,188]. S. Candel suggested to propagate artificially thickened premixed flame fronts on a coarse grid following an original idea developed by O'Rourke and coworkers [44, 45] in another context.

### 3.5.2. The Thickened Flame model for LES (TFLES)

Following simple theories of laminar premixed flames [189,190], the flame speed $S_L^0$ and thickness $\delta_L^0$ scale as:

$$S_L^0 \propto \sqrt{D_{th} W} \quad \text{and} \quad \delta_L^0 \propto \sqrt{\frac{D_{th}}{W}} \qquad (34)$$

where $D_{th}$ is the thermal diffusivity while $W$ measures the reaction rate intensity, for example the pre-exponential factor when using a single Arrhenius law. If the thermal diffusivity is multiplied by a factor $F$ while $W$ is divided by the same factor $F$, the flame thickness $\delta_L^0$ increases by $F$ while the flame speed $S_L^0$ is conserved. For sufficiently large $F$ values, the thickened flame front can be resolved on any LES computational mesh.<fnref n="6" /> Since reaction rates are still expressed using Arrhenius laws, various phenomena such as ignition, flame stabilization or flame/wall interactions are expected to be accounted for without requiring dedicated submodels. Using this original approach, Thibaut and Candel [191] performed LES of flame flashback due to instabilities in a turbulent premixed dump combustor, reproducing behavior experimentally observed by Keller et al. [192].

However, there is a prize to pay when the flame front is thickened: the Damköhler number $D_a$ comparing turbulent, $\tau_t$, and chemical, $\tau_c$, time scales:

$$D_a = \frac{\tau_t}{\tau_c} = \frac{l_t}{u'} \frac{s_l}{\delta_L^0} \qquad (35)$$

decreases by a factor $F$, becoming $D_a/F$, meaning that interactions between turbulence and combustion are modified. As discussed in Section 3.4.2 and displayed in Figs. 12 and 13, when the ratio between the turbulence integral length scale $l_t$ and the laminar flame thickness $\delta_L^0$ is decreased by $F$, the flame becomes less sensitive to turbulent motions and, being thicker, can be less wrinkled. Angelberger et al. [193] and Colin et al. [138] investigated this point using DNS and, following the ITNFS formalism (see Section 3.4.2), derived an efficiency function to compensate the loss in flame surface due to the thickening process, i.e. to model sub-grid scale combustion. Charlette et al. [139] refined later this function, starting from the same DNS.<fnref n="7" />

In practice, the thickened flame approach is implemented by a “mapping” where the diffusivity and the reaction rate are changed according to:

$$ \begin{array}{l} \text{Diffusivity: } D_{th} \rightarrow F D_{th} \rightarrow \Xi F D_{th} \\ \text{Reaction rate: } W \rightarrow W/F \rightarrow \Xi W/F \\ \text{thickening} \quad \rightarrow \quad \text{wrinkling} \end{array} $$

where $\Xi$ is the efficiency function measuring the sub-grid scale surface, i.e. flame wrinkles due to turbulent motions smaller than the flame thickness $F\delta_L^0$, lost in the thickening procedure.<fnref n="8" />

According to Eq. (34), the flame speed $S_L^0$ and thickness $\delta_L^0$ become respectively:

$$S_T^0 = \Xi S_L^0 \quad \text{and} \quad \delta_T = F \delta_L^0 \qquad (36)$$

where $S_T^0$ is the sub-grid scale turbulent flame speed. The lost surface is then replaced by an increased flame speed. Regarding the classification of models proposed in Section 3.4.4, the thickened flame model is related to the geometrical approach as the wrinkling factor $\Xi$ measures the sub-grid scale flame surface. Note also that the effective diffusivity $\Xi F D_{th}$ may be split as:

$$\Xi F D_{th} = F D_{th} + (\Xi - 1) F D_{th} \qquad (37)$$

where the last RHS term models the sub-grid scale turbulent contribution.

This model has been derived and is theoretically supported for premixed flames. However, it has also been successfully used in partially premixed and non-premixed combustion, introducing a dynamic version where thickening is only applied in the reaction zones, identified with a dedicated sensor [194–196] to preserve diffusion in non-reacting regions. The thickening factor may also depend on the local grid resolution [197].

The TFLES, also called artificially thickened flame (ATF) model by other authors, has been successfully used in a large number of works [138,197–232] among others.

### 3.5.3. Flame wrinkling factor models

Following Charlette et al. [233], the filtered reaction rate, $\overline{\omega}$, entering species and energy filtered balance equations may be written under the generic form:<fnref n="9" />

$$\overline{\omega}(\theta) = \Xi_{\Delta} W_{\Delta}(\tilde{\theta}) / \Delta \qquad (38)$$

where $\Delta$ is the LES filter size, $W_{\Delta}(\tilde{\theta})$, depending only on filtered quantities, is the reaction rate per unit of flame surface area, while $W_{\Delta}(\tilde{\theta})/\Delta$ corresponds to the resolved reaction rate (i.e. when $\Xi_{\Delta} = 1$). $\theta$ and $\tilde{\theta}$ stand here for any quantity entering the reaction rate. The wrinkling factor $\Xi_{\Delta}$ measures the ratio of total to resolved flame surface area in the filtering volume.

Table 2 summarizes $\overline{\omega}(\theta)$ and $W_{\Delta}(\tilde{\theta})$ expressions for several models: the Boger et al. [113,235] algebraic model, extending the Bray–Moss–Libby formalism [236] to LES, the level-set, or “G-equation”,

$$\overline{\omega}(\theta) = \Xi_{\Delta} W_{\Delta}(\tilde{\theta}) / \sqrt{\Delta^2 + \delta_L^0} \qquad (39)$$

<footnote n="4">⁴ Integrating the balance equation for the filtered flame surface density over the combustion chamber volume, Richard and Veynante [185] derived a 0-D flame wrinkling balance equation for simulations of internal combustion engines.</footnote>

<footnote n="5">⁵ Note that turbulent flows and flames behave differently: most of the flow turbulence kinetic energy lies in largest scales and is expected to be resolved in LES [186], while combustion takes place at a sub-grid scale level. This difference must be kept in mind when developing models, especially for scale-similarity or dynamical approaches.</footnote>

<footnote n="6">⁶ This result is a property of convection/diffusion/reaction balance equations and is easily proved by replacing the spatial coordinate $x$ by $x/F$ in the balance equation for a 1D steady propagating flame.</footnote>

<footnote n="7">⁷ Charlette et al. [139] performed the integration over turbulence scales using the energy spectrum and not a multi-fractal description of turbulence [133,138], corrected some limiting behaviors of the Colin et al. function,</footnote>

<footnote>introduced the model parameter as an exponent and suppressed the a priori prescription of the Reynolds turbulence number.</footnote>

<footnote n="8">⁸ Colin et al. [138] defined the efficiency function, comparing the wrinkling factors of the actual flame of thickness $\delta_L^0$, $\Xi(\delta_L^0)$, and the artificially thickened flame of thickness $F\delta_L^0$, $\Xi(F\delta_L^0)$, as $E = \Xi(\delta_L^0)/\Xi(F\delta_L^0)$. Because the thickened flame is not expected to be wrinkled at scales lower than its thickness $F\delta_L^0$ (see Figs. 12 and 13), $\Xi(F\delta_L^0) \approx 1$ and $E \approx \Xi(\delta_L^0)$ [139].</footnote>

<footnote n="9">⁹ The original expression Eq. (38) has been slightly modified by Veynante and Moureau [234], introducing the laminar flame thickness $\delta_L^0$ to ensure a correct behavior and suitable comparisons with DNS when $\Delta \to 0$:</footnote>

<footnote>This correction is not considered here as some models in Table 2 are derived for infinitely thin flame fronts.</footnote>
<caption>**Table 2**</caption>

Filtered reaction rate $\dot{\omega}(\bar{\theta})$ and $W_{\Delta}(\bar{\theta})$ entering the generic expression Eq. (38) for four turbulent premixed combustion LES models. $\rho_a$ is the fresh gases density, $S_L^0$ and $\delta_L^0$ the laminar flame speed and thickness, respectively. G is the level-set field, usually defined as a signed distance to the flame front, $S_T = \Xi_{\Delta}S_L^0$ the sub-grid scale turbulent flame speed, F the TFLES thickening factor and $\dot{\omega}$ the progress variable reaction rate (laminar flame). $\dot{\Omega}_L(\bar{\theta}, \Delta)$ is given by filtering one-dimensional laminar premixed flames [237]. Units are given in mass (M), length (L) and time (T).

Source: Table adapted from [234].

<table><thead><tr><th>Model</th><th>ω̃(θ̃) [ML<sup>-3</sup>T<sup>-1</sup>]</th><th>W<sub>Δ</sub>(θ̃) [ML<sup>-2</sup>T<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Boger et al.</td><td>4ρ<sub>a</sub>S<sub>L</sub><sup>0</sup>Ξ<sub>Δ</sub>√(6/π)(θ̃(1-θ̃)/Δ)</td><td>4ρ<sub>a</sub>S<sub>L</sub><sup>0</sup>√(6/π)(θ̃(1-θ̃)/Δ)</td></tr><tr><td>Level set</td><td>ρ<sub>a</sub>Ξ<sub>Δ</sub>S<sub>L</sub><sup>0</sup>|V̇G|</td><td>ρ<sub>a</sub>S<sub>L</sub><sup>0</sup>Δ|V̇G|</td></tr><tr><td>TFLES</td><td>Ξ<sub>Δ</sub>(ω̃(θ̃)/F)</td><td>δ<sub>L</sub><sup>0</sup>ω̃(θ̃)</td></tr><tr><td>F-TACLES</td><td>Ξ<sub>Δ</sub>δ<sub>L</sub>(θ̃, Δ)</td><td>ΔΩ<sub>L</sub>(θ̃, Δ)</td></tr></tbody></table>

approach [187], the thickened flame model (Section 3.5.2) where the thickening factor is here expressed in terms of flame thickness and filter size as $F = \Delta/\delta_L^0$ or the F-TACLES model where resolved reaction rates, as well as transport and diffusive terms in the filtered progress variable balance equation, are modeled from filtered one-dimensional laminar premixed flames [92,93,217,237,238].<fnref n="10" />

Expression Eq. (38) holds as long as flame/turbulence interactions are described in terms of flame surface wrinkling factor or sub-grid scale turbulent flame speed $S_T = \Xi_{\Delta}S_L^0$ (flamelet assumption). Note that if, strictly speaking, the thickening process is not equivalent to a filtering operation, models based on filtering laminar flame elements, such as F-TACLES or Duwig model [239,240], and TFLES can be closed using the flame surface wrinkling factor $\Xi_{\Delta}$.

The wrinkling factor $\Xi_{\Delta}$ is then the only quantity to be closed [234]. This can be done using the efficiency function proposed by Charlette et al. [139] and discussed in Section 3.5.2:

$$ \Xi_{\Delta} = \left( 1 + \min \left[ \max \left( \frac{\Delta}{\delta_L^0} - 1, 0 \right), \Gamma \left( \frac{\Delta}{\delta_L^0}, \frac{u'_\Delta}{S_L^0}, Re_\Delta \right) \frac{u'_\Delta}{S_L^0} \right] \right)^\beta \quad (39) $$

where $\Gamma$ is an efficiency function [139], $\beta$ the model parameter<fnref n="11" /> and $Re_\Delta = u'_\Delta\Delta/\nu$ the sub-grid scale turbulence Reynolds number, $\nu$ being the fresh gas kinematic viscosity.<fnref n="12" /> A model is needed for the sub-grid scale intensity $u'_\Delta$ [138].

However, expressions such as Eq. (39) assume an equilibrium between turbulence motions and flame surface which is clearly not reached at least in two situations: (i) close to the injector lips in jet flames, where the flame is initially laminar and wrinkled by turbulence motions while convected downstream [242]; (ii) during the early stages of the development of a flame kernel following ignition, as in spark-ignited combustion engines: the initially laminar flame kernel is progressively wrinkled by turbulence motions as it grows [112]. These situations can be handled solving a balance equation for the flame surface density [111–113,173] or the wrinkling factor [243,244], relaxing the equilibrium assumption.<fnref n="13" />

An alternative is to automatically adjust the model parameter $\beta$ in Eq. (39), both in time and space, from the known resolved flame motions, following the dynamic formalism proposed by Germano

et al. [245] for unresolved momentum fluxes. The idea is to estimate the flame surface density, averaged over a given volume, referred here as $\langle \cdot \rangle$, from two routes: (i) filtering the LES flame surface density field at the test-filter scale $\hat{\Delta} > \Delta$; (ii) estimating directly the flame surface density at the test-filter scale from the model [234]:

$$ \langle \Xi_{\Delta} | \nabla \bar{\theta} | \rangle = \langle \Xi_{\gamma \Delta} | \hat{\nabla} \bar{\theta} | \rangle \quad (40) $$

where $|\nabla\bar{\theta}|$, $\Xi_{\Delta}|\nabla\bar{\theta}|$, $|\hat{\nabla}\bar{\theta}|$ and $\Xi_{\gamma\Delta}|\hat{\nabla}\bar{\theta}|$ measure resolved and total generalized flame surface densities at LES and effective filter scales, respectively [54,114]. The length $\gamma\Delta$ corresponds to the effective filter width when combining LES and test filter operators $(\gamma\Delta = (\Delta^2 + \hat{\Delta}^2)^{1/2}$ for Gaussian filters). Modeling the wrinkling factor can be made using a simple fractal-like expression:<fnref n="14" />

$$ \Xi_{\Delta} = (\Delta/\delta_c)^{\beta} \quad (41) $$

where $\delta_c$ is the inner cut-off scale of flame wrinkles, of the order of $\delta_L^0$. Reporting Eq. (41) into Eq. (40) provides $\beta$, uniform over the averaging volume and independent on the filtering scale:

$$ \beta = \frac{\ln\left(\left\langle\left|\nabla\bar{\theta}\right|\right\rangle/\left\langle\left|\hat{\nabla}\bar{\theta}\right|\right\rangle\right)}{\ln(\gamma)} \quad (42) $$

By construction, this expression recovers unity wrinkling factors ($\Xi_{\Delta} = \Xi_{\gamma\Delta} = 1$) when the flame wrinkles are fully resolved in simulations ($\beta = 0$), a constraint required to recover the DNS limit when the filter size tends towards 0. Direct numerical simulations show that the filtered progress variable $\bar{\theta}$ can be approximated using mass-weighted filtered progress variable $\tilde{\theta}$ [234]. This approach has also been combined with a scale similarity flame surface density model [184].

These dynamic formulations based on flame surface description provided very good results, especially for transient situations that cannot be handled assuming equilibrium between turbulence motions and flame wrinkles or when operating conditions are changed [220, 222,234,241,242,246–257].

## 4. Combustion noise

Combustion is not the only field of expertise of Pr. Sébastien Candel. His first interest was actually acoustics, explaining his impressive contribution in fields where combustion and acoustics are closely linked such as combustion noise and thermoacoustics. Combustion noise is an issue not to be underestimated for aircraft engine, helicopter engine, auxiliary power unit or gas turbine, to name but a few examples where his work has been essential.

This section is organized as follows. The concepts of direct noise combustion and indirect noise combustion are introduced in Section 4.1 and in Section 4.2 respectively. The seminal study by Candel [258] and Marble and Candel [259] on the response of a nozzle flow to entropy disturbances is discussed in Section 4.3. The contribution of combustion noise in core noise is finally presented in Section 4.4. Fig. 15 shows a diagram of the combustion chamber, turbines stages and nozzle exhaust of an aircraft engine to introduce the two main components of combustion noise.

### 4.1. Direct combustion noise

The first noise source leading to the acoustic field ($p_1^+$) is called direct combustion noise and is generated by the unsteady heat released by the flame (F) in the combustion chamber. In free space for an open

<footnote n="10">¹⁰ The Boger et al. model [113,235] may be viewed as the limit of F-TACLES for infinitely thin flame fronts.</footnote>

<footnote n="11">¹¹ $\beta$ can be related to the fractal dimension D of the flame surface as $\beta = D - 2$ when $\beta$ is spatially uniform and lower than unity, but Eq. (39) is more general [139].</footnote>

<footnote n="12">¹² The original expression [139] is slightly modified to both ensure a maximum value $(\Delta/\delta_L^0)^\beta$ for large turbulence intensities and a unity wrinkling factor when $\Delta \le \delta_L^0$ [234,241].</footnote>

<footnote n="13">¹³ Filtered flame surface density, $\Sigma_{\Delta}$, and wrinkling factor, $\Xi_{\Delta}$, are linked by $\Sigma_{\Delta} = \Xi_{\Delta}|\nabla\bar{\theta}|$ [114,243].</footnote>

<footnote n="14">¹⁴ This expression corresponds to the saturated value of Eq. (39), when the min value is given by $\Delta\delta_L^0 - 1$. However, Eq. (39) can be recast under the form of Eq. (41) [234].</footnote>
<figure><img src="image_17.png" /><figcaption>Fig. 15. Noise in a jet engine: direct combustion noise in solid line generated by the flame ($p_1^+$) to ($p_2^+$); indirect combustion noise ($p_{w2}^+$) generated by the sudden acceleration and deceleration of entropy spots ($w_1^+$) through the turbine stages and ($w_2^+$) possibly through the exhaust nozzle.</figcaption></figure>

flame, it can be shown from Lighthill's theory and by assuming isobaric combustion that the acoustic far field is given by [260,261]:

$$p'(\mathbf{x}, t) \approx \frac{\gamma - 1}{4\pi x c_{\infty}^2} \int_{V} \frac{\partial \dot{\omega}_T}{\partial t} \left( \mathbf{y}, t - \frac{|\mathbf{x} - \mathbf{y}|}{c_{\infty}} \right) d\mathbf{y} \quad (43)$$

where $p' = p - p_{\infty}$ is the fluctuating pressure perceived at the observer position $\mathbf{x}$ at time $t$, $\gamma$ is the constant ratio of specific heats, $c$ is the speed of sound and $\dot{\omega}_T$ the heat release per unit volume of the flame. The subscript $\cdot_{\infty}$ indicates a suited constant reference value, the ambient medium at rest here. This source term largely dominates all molecular processes involved in an exact formulation [262], this point is discussed in the next paragraph. In Eq. (43), the integration is performed over the combustion volume $V$ and $|\mathbf{x} - \mathbf{y}|/c_{\infty}$ corresponds to the time necessary for acoustic waves to propagate from a current source point $\mathbf{y}$ to the observer $\mathbf{x}$. This integral formulation and its many subtleties are well known to the aeroacoustic community [263].

Direct combustion noise is monopolar in terms of acoustic efficiency, broadband and concentrated on low frequencies. In other words, this noise mechanism can often be modeled with compact monopole sources: the typical acoustic wavelength $\lambda_u$ is large compared to the size of the source ($F$) and the retarded time in Eq. (43) is then reduced to $|\mathbf{x}|/c_{\infty}$. Experimental demonstration of the link between acoustic pressure and the heat release fluctuations predicted by Eq. (43) is generally achieved by measuring the intensity of light emission from free radicals excited during combustion. It is now possible to make highly detailed comparisons with the improvement of optical techniques and acquisition rate. For an acoustically excited, laminar premixed flame impinging on a water-cooled plate, the spontaneous emission of CH$^+$ radicals is used in Schuller et al. [264] to compute the radiated pressure, which matches the measured signal by a microphone. A recent review as well as an experimental characterization of the direct combustion noise of a low swirl turbulent flame can be found in Shoji et al. [265]. Numerical simulations combining direct or large eddy simulation and an acoustic analogy can also provide insightful analyses [260,261,266,267].

Direct combustion noise ($p_1^+$) in Fig. 15 may be reflected at the combustion chamber outlet ($p_1^-$) and propagate to the flame ($F$), and may also excite resonance modes of the combustion chamber. These acoustic feedbacks lead to undesirable self-sustained oscillations and thermoacoustic instabilities when the Rayleigh criterion is satisfied [268], representing a major risk of destroying the engine itself. The topic is presented in a general framework in Section 5, and for a specific geometrical arrangement of gas turbine combustion chambers in Section 7.

## 4.2. Indirect combustion noise

Unsteady combustion also produces entropy spots ($w_1^+$), regions of non-uniform temperature, which are convected by the low Mach number flow through the combustor towards the turbine stage (Fig. 15). The successive acceleration and deceleration of entropy perturbations in the turbine stages will also produce reflected ($p_{w1}^-$) and transmitted acoustic waves ($p_{w2}^+$) contributing to a secondary combustion noise. This indirect noise was identified early on in the literature by Candel [258] during his PhD thesis as a possible additional source of combustion noise, associated with the passage of these disturbances through the turbine blade channels. This was all the more remarkable as no experimental evidence were available.

Indirect combustion noise can be seen as noise generation by the non-uniform convection of entropy, as well as changes in the mixture composition and possible variation of thermodynamic and transport properties of species [261,262,266]. The contribution of entropic disturbances is however first-order compared to all the other mentioned sources, when breakdown of combustion noise sources is performed on realistic cases. Very large temperature difference can indeed exist at the outlet of a combustion chamber.

A elegant manner of introducing the subject is to consider Lighthill's analogy. Although this formulation is not always appropriate for numerical predictions [269], it provides a good understanding of the physical mechanisms involved in aeroacoustics and indirect combustion noise [260,261]. Lighthill's equation is obtained from the fundamental governing equations of fluid dynamics:

$$\frac{1}{c_{\infty}^2} \frac{\partial^2 p}{\partial t^2} - \nabla^2 p = \nabla \cdot \nabla \cdot \left( \rho \mathbf{u} \otimes \mathbf{u} - \bar{\bar{\tau}} \right) - \frac{\partial^2}{\partial t^2} \left( \rho - \frac{1}{c_{\infty}^2} p \right) \quad (44)$$

where $\rho$ denotes the density, $\mathbf{u}$ the velocity vector and $\bar{\bar{\tau}}$ the viscous tensor. This exact equation needs to be interpreted in two ways. Source terms are identified in the right-hand side of Eq. (44), forcing an acoustic wave equation, but the source domain is localized in space since sound wave propagation takes place in a medium at rest with a constant speed of sound $c_{\infty}$. Since the reformulation is exact, one can also say that for propagation in a more complex medium, the linear complementary terms needed to form the associated operator on the left-hand side are embedded in the source terms [262,263,270]. Let us return to the source terms now. The first term is the Lighthill stress tensor $\rho \mathbf{u} \otimes \mathbf{u}$ leading to the well-known eighth power law of the velocity for subsonic aerodynamic noise [263,271]. The viscous term $\bar{\bar{\tau}}$ plays no role in noise generation, but is responsible of sound attenuation during propagation. The generation of direct and indirect combustion noise is contained in the third term on the right-hand side and its expression is
usually recast by introducing the so-called excess density $\rho_e$ [272]:

$$ \rho_e = (\rho - \rho_\infty) - \frac{1}{c_\infty^2}(\rho - \rho_\infty) \quad (45) $$

which leads [262] for this third source term to:

$$ -\frac{\partial^2 \rho_e}{\partial t^2} = \frac{\partial}{\partial t} \left[ \frac{\rho_\infty}{\rho} \frac{\gamma - 1}{c^2} \left( \dot{\omega}_T + \nabla \cdot (\lambda \nabla T) + \nabla u : \bar{\tau} - \rho \sum_k Y_k c_{p,k} \mathbf{v}_k \cdot \nabla T \right) + \rho_\infty \frac{D}{Dt} (\ln r) \right. \\ \left. + \frac{1}{c_\infty^2} \frac{\partial}{\partial t} \left[ \left( 1 - \frac{\rho_\infty c_\infty^2}{\rho c^2} \right) \frac{Dp}{Dt} - \frac{p - p_\infty}{\rho} \frac{D\rho}{Dt} \right] + \nabla \cdot \frac{\partial}{\partial t} (\rho_e u) \right] \quad (46) $$

for an ideal gas. In this expression, $D/Dt = \partial_t + u \cdot \nabla$ is the material derivative, $\lambda$ is the thermal conductivity, $T$ the temperature; $Y_k$, $c_{p,k}$ and $\mathbf{v}_k$ are the mass fraction, the specific heat and the diffusion velocity of species $k$, and $r$ is the ratio between the perfect gas constant and the molar mass of the mixture. The first term associated with unsteady heat release corresponds to the integral solution given by Eq. (43) for combustion at constant pressure. This monopole source term is found to be dominant in experiments [260,264,265] and numerical simulations [273] for low Mach number flows. On the second line, the three first source terms are the contributions of the heat flux, viscous dissipation and diffusion of species whereas the last term corresponds to noise induced by changes in molar composition of the mixture. The contribution of all these noise components to indirect combustion noise is expected to remain small, as already indicated above. Finally, the third line contains the source terms associated with entropy noise.

The reader can find valuable details and references in Candel et al. [260] and Dowling and Mahmoudi [261] regarding Eq. (46) and alternatives expressions of the source terms for low Mach number flows. Equivalent expressions to Eqs. (45) and (46) can also be derived for Phillips and Lilley's wave equation [262] to better identify and separate propagation terms from noise source terms [269]. These scientific topics go beyond the scope of this section and fall within the field of aeroacoustics, an area of research in which Sébastien Candel's contribution is also outstanding.

### 4.3. Response of a nozzle flow to entropy disturbances

Three families of waves can be identified from the linearization of the Euler equations around an arbitrary steady base flow, with no gravitational force. These waves only remain decoupled in two particular cases: either the base flow is homogeneous, or the waves are very high-frequency compared with the length scales of this base flow [258,274] leading to geometrical acoustics [275] for instance. To investigate the main source of indirect combustion noise linked to entropy fluctuations, Marble and Candel [259] consider the propagation of plane acoustic and entropic perturbations through a one-dimensional inviscid nozzle flow. The base flow is assumed to be steady and the problem is analytically solved for the low frequency case, when the length of the nozzle $l_N$ is small with respect to the characteristic wavelength of perturbations, or equivalently when the dimensionless angular frequency $\Omega = \omega l_N/c_1$ goes to zero where $M_1 = u_1/c_1$ is the upstream Mach number, the subscript 2 denoting the downstream part of the nozzle, as shown in Fig. 16.

For an incoming convected entropy wave $w_1^+$ of the form $s'/c_p = w_1^+ \exp(i\omega(t - x_1/u_1))$, the propagating acoustic wave $p_{w2}^+$ associated with the indirect combustion noise is given by:

$$ \frac{p_{w2}^+}{w_1^+} = \frac{M_2 - M_1}{1 + M_2} \frac{(1/2)M_2}{1 + [(\gamma - 1)/2]M_1 M_2} \quad (47) $$

for a subcritical flow. The coupling or conversion of an entropy wave to an acoustic wave is induced by the flow acceleration in the nozzle. Matching conditions for the wave amplitudes, without phase shift for a compact nozzle, can be established. Leyko et al. [276] have extended

<figure><img src="image_18.png" /><figcaption>Fig. 16. Sketch of a nozzle flow to investigate indirect noise combustion.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Ratio in logarithmic scale between the indirect and the direct combustion noise for a compact nozzle, using the relations of Tab. 3 for the subcritical case and the supersonic case without shock wave, adapted from Leyko et al. [276]. The solid line corresponds to $\eta = 1$.</figcaption></figure>

this approach by introducing a modeling of the acoustic response to a heat release under the same assumptions, $w_1^+/p_1^+ = M_1/(1 + M_1)$, leading to the following expression for the ratio between the indirect and direct combustion noise:

$$ \eta = \frac{p_{w2}^+}{p_2^+} = \frac{p_{w2}^+}{w_1^+} \frac{w_1^+}{p_1^+} \frac{p_1^+}{p_2^+} = \frac{1}{4M_1} \frac{M_2^2 - M_1^2}{1 + [(\gamma - 1)/2]M_2^2} \quad (48) $$

where $p_2^+/p_1^+$ is the transfer function for an acoustic wave through the compact nozzle. Similar expressions can be derived for a supersonic nozzle flow and also the case including a normal shock in the divergent [259]. The results for the two first cases are gathered in Fig. 17.

The ratio $\eta$ is plotted in Fig. 17 as a function of the two Mach numbers $M_1$ and $M_2$ using a logarithmic scale. The solid line is provided by the condition $\eta = 1$. For a given upstream Mach number $M_1$, the indirect combustion noise dominates for downstream Mach number $M_2$ above this solid line. For the sake of completeness, numerical simulations for $0 \le \Omega \le 1$ indicate that the threshold $\eta = 1$ decreases with the reduced frequency $\Omega$ for a non compact nozzle [276,277].

The experimental counterpart to Marble and Candel's work is provided by Bake et al. [278] with the entropy wave generator. As always in acoustics, many precautions must be taken to perfectly control or characterize in terms of impedance the boundary conditions. The electrical heat source must also be carefully designed. A typical signal for temperature measured upstream of the nozzle ($w_{1s}^+$) is shown in Fig. 18, and is compared to the pressure signal measured downstream, ($p_{w2}^+$), but including installation effects, that is a reflection coefficient imposed by the duct outlet. The duration of the wavepacket is about 0.1 s, the nozzle length is about $l_N = 263$ mm which leads to a dimensionless angular frequency $\Omega \approx 0.05$. Numerical simulations are able to reproduce the measured pressure signals, taking into account the effective boundary conditions of the channel [277,279], or using an acoustic analogy [278,280] based on the main source term $\nabla \cdot \partial(\rho_e u)\partial t$ of the indirect combustion noise source term in Eq. (46).

Returning to the original problem, it is physically useful to relate the entropy wave to the expression of the indirect combustion noise source term Eq. (46) in the previous section. By differentiation of the
<table><caption>**Table 3**
Determination of the ratio η between indirect and direct combustion noise [276] (including the correction of numerical factors) for a compact nozzle using Marble and Candel relations [259]. Same notations are used as in Fig. 15 and Fig. 16, the sign +/− indicates a downstream/upstream propagation along x₁.</caption><thead><tr><th colspan="3">Subsonic compact nozzle flow</th><th colspan="3">Supersonic compact nozzle flow</th></tr></thead><tbody><tr><td><sup>p</sup><sub>w2</sub><br>ω<sub>1</sub><sup>+</sup></td><td>M<sub>2</sub> - M<sub>1</sub><br>1 + M<sub>2</sub></td><td>(1/2)M<sub>2</sub><br>1 + [(γ - 1)/2]M<sub>1</sub>M<sub>2</sub></td><td><sup>p</sup><sub>w2</sub><br>ω<sub>1</sub><sup>+</sup></td><td>M<sub>2</sub> - M<sub>1</sub><br>2</td><td>(1/2)<br>1 + [(γ - 1)/2]M<sub>1</sub></td><td>≠ Eq.(29) in [276]</td></tr><tr><td><sup>p</sup><sub>1</sub><sup>+</sup></td><td>2M<sub>2</sub><br>1 + M<sub>1</sub></td><td>1 + [(γ - 1)/2]M<sub>2</sub><br>1 + [(γ - 1)/2]M<sub>1</sub>M<sub>2</sub></td><td><sup>p</sup><sub>1</sub><sup>+</sup></td><td>1 + [(γ - 1)/2]M<sub>2</sub><br>1 + [(γ - 1)/2]M<sub>1</sub></td><td></td><td></td></tr><tr><td>ω<sub>1</sub><sup>+</sup></td><td>M<sub>1</sub><br>1 + M<sub>1</sub></td><td></td><td>ω<sub>1</sub><sup>+</sup></td><td>M<sub>1</sub><br>1 + M<sub>1</sub></td><td></td><td></td></tr><tr><td>η<br>≠ Eq. (34) in [276]</td><td><sup>p</sup><sub>w2</sub><br>ω<sub>1</sub><sup>+</sup></td><td><sup>p</sup><sub>w2</sub><br>ω<sub>1</sub><sup>+</sup></td><td><sup>p</sup><sub>1</sub><sup>+</sup></td><td>M<sub>2</sub> - M<sub>1</sub><br>1 + [(γ - 1)/2]M<sub>2</sub></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Reproduction of Fig. 4 in Bake et al. [278]. The temperature signal measured by a fast thermocouple upstream of the nozzle is plotted in red dashed line and the pressure signal measured by a microphone downstream of the nozzle is plotted in solid black line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

thermodynamic relationship $p = p(\rho, s)$ and using the ideal gas law $p = \rho rT$, one gets:

$$ \frac{p'}{\gamma p_{\infty}} = \frac{\rho'}{\rho_{\infty}} + \frac{s'}{c_p} \qquad \frac{T'}{T_{\infty}} = \frac{s'}{c_p} + (\gamma - 1) \frac{p'}{\gamma p_{\infty}} \qquad (49) $$

The second relation reminds us that entropy fluctuations are not temperature fluctuations, and the first identity corresponds to the excess density Eq. (45) introduced earlier in the text. The excess density, in other words, the amount of density fluctuations not balanced by acoustic pressure fluctuations, corresponds to entropy fluctuations.

The Candel-Marble problem, as described in Ffowcs-Williams and Howe [281], has led to a large number of developments still aimed at characterizing indirect combustion noise. A few references are given here for this second golden age of entropy noise, among many others that can be found in these studies. Huet & Giauque [282] have extended the analytical work by Marble and Candel [259] to the nonlinear regime, and Giauque et al. [283] has investigated numerically the effect of a finite nozzle length. Duran & Moreau [284] have proposed an asymptotic expansion to take into account a finite frequency and the detailed nozzle geometry. Azimuthal modes including vorticity waves are considered in Duran & Morgans [285]. The contribution of the compositional noise is investigated in Jain & Magri [286]. The response of a nozzle flow to entropic disturbances has now become a classic problem to study numerically [287–290].

New insights have also been gained experimentally. In Weilenmann et al. [291], measurements of entropy transfer function relating the acoustic forcing to the production of entropy wave by technically premixed flames have been obtained with laser-induced fluorescence thermometry of OH radical between 40 Hz and 90 Hz for two types of burners. This work demonstrates the highly nonlinear character of

the entropy wave production. The entropy and flame transfer functions have been measured in [292] for a turbulent swirled technically premixed flame with a new technique using tunable diode laser absorption spectroscopy with wavelength modulation spectroscopy and multi-microphone method. The transfer function relating backward acoustic waves to forward entropy waves in a highly turbulent channel flow terminated by a choked nozzle has been measured between 60 and 180 Hz in [293]. It is the first time that sound production from entropy waves accelerating in a choked nozzle is quantified in this frequency range and at such large Reynolds number in the upstream duct. The indirect sound production was found approximately 60% lower than the result predicted by the classic compact model [259] for supercritical compact nozzles. The difference is caused by the three dimensional deformation of the coherent entropy spots during their acceleration in the nozzle.

### 4.4. Contribution to core noise

Fan noise, jet noise, core noise and airframe noise make up the classical breakdown for commercial aircraft noise. Core noise is the additional noise component to jet noise which can be identified from engine or gas turbine exhaust [294] and combustion noise is a component of core noise. There is no obvious dimensional law for simply separating subsonic jet noise from direct and indirect combustion noise, even though the latter increases less rapidly than jet noise as a function of the exhaust velocity.

Entropy waves propagate more slowly than acoustic waves in the engine so that a cross-correlation procedure can be used to identify direct and indirect combustion noise [295]. The time delay between probes in the combustor and far field microphones increases as the frequency decreases, indicating the relative domination of direct combustion noise. As an illustration, the measured noise spectrum of a turboshaft engine is plotted in Fig. 19. The different components of core noise are reported, with the domination of combustion noise in the lower part of the spectrum.

Cumpsty and Marble [294] have generalized in one sense the one-dimensional compact nozzle solution presented in Section 4.3 by replacing the nozzle by a turbine blade row. Wavelengths associated with acoustic, entropy and vorticity disturbances are again large by comparison with the blade chord. The blade row is then replaced by an infinitesimally thin actuator disk and two-dimensional waves are now considered with the main direction of the engine ($x_1$ previously) and a propagation angle $\theta$ with respect to this main axis. There is no radial dependence for waves and the mean flow in this approach.

Leyko [296] has first proposed a comprehensive procedure for dealing with direct and indirect combustion noise in realistic configurations, see Fig. 20. The first step is to carry out a simulation of the combustion chamber, taking into account all the specifics of the case under consideration. Entropy, acoustic and vorticity waves are then identified in a second step, by searching for a consistent expression with the assumptions of the Cumpsty and Marble model [294]. An average in the radial direction is for instance required to comply with
<figure><img src="image_21.png" /><figcaption>Fig. 19. Sound pressure level spectrum (SPL, relative dB) of a turboshaft engine in one-third octave bands: combustion noise, compressor; turbines; total noise (courtesy of Éric Bouty from Safran Helicopter Engine). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Overview of the workflow to compute direct and indirect combustion noise. The analytical model of Cumpsty and Marble [294] is implemented in CHORUS [296] and the global methodology from the large eddy simulation of the combustor to the acoustic far field is called CONOCHAIN [297].</figcaption></figure>

the analytical model. Furthermore, the calculations are performed in Fourier space. A matrix formulation is then used for each blade row of the engine to finally obtain the direct and indirect acoustic field associated with combustion. The first application to the engine demonstrator DEM21 has shown the importance of indirect combustion noise compared with direct noise. The full methodology has been improved by Liveardon [297] with a real breakthrough in combustion noise prediction. Liveardon et al. [298] have applied the CONOCHAIN of Fig. 20 to an helicopter engine. For the low frequency part of spectra, a good agreement is found between numerical predictions and measurements. This work also confirms the importance of indirect combustion noise for helicopter turboshaft engine.

As in the case of the compact nozzle model, many different teams have worked on the limits of the analytical model [299]. The topic has been reviewed by Ihme [266] and discussed by Brouzet et al. [300].

## 5. Flame dynamics and combustion instability

Among all fascinating phenomena found in flames, combustion instabilities coupled by pressure waves are probably the most spectacular: they are not only a very wide field of studies for fundamental research but also a crucial issue in many real systems where they lead to considerable problems [301,302]. Pr Candel has brought contributions to this field since the 80's.

The next three sections present how the field of combustion instabilities has evolved in the last 50 years and shows the contribution of Pr Candel and his team. It begins with a general description of combustion

instability physics and approaches (Section 5) mainly considering the growth of linear perturbations on a baseline flame. Section 6 then tackles the issue of nonlinearity and of limit cycles observed in many combustors. Finally Section 7 presents studies in annular chambers, typical of gas turbines, where instabilities couple to azimuthal acoustic modes.

### 5.1. Combustion instability coupled by pressure waves

Thermo-acoustic instabilities, which manifest as undesired large self-sustained combustion oscillations, arise from resonant feedback between the unsteady heat released by a flame and synchronized pressure waves, which, in turn, alter the flow upstream of the flame as described by the feedback loop in Fig. 21(a) taken from Candel's Hottel lecture in 2002 [303].

In the 1970s, the issue of combustion instabilities gained prominence due to dramatic failures in rocket engines, as well as the push to develop new generations of afterburners, ramjets, and, to a lesser extent, large-scale boilers. Sébastien Candel, who had just joined ONERA at the time, began addressing these problems using similar analytical tools he had developed during his PhD to analyze sound waves in jet engines [258]. At the core of thermoacoustic analysis back then was the foundational two-parameter $n-\tau$ time-lag model, originally introduced by Crocco for rocket engines [304] and adapted to industrial burners by Putnam [305]. This heuristic model captured the key insight that fluctuations in the heat release rate $\dot{Q}_1$ arise from variations in fuel or air injection velocity $u$, after a certain time delay:

$$ \frac{\dot{Q}_1}{\dot{Q}_0} = n \frac{u_1(t - \tau)}{u_0} \qquad (50) $$

where the interaction index $n$ quantifies the strength of this coupling between relative velocity $u_1/u_0$ and relative heat release rate $\dot{Q}_1/\dot{Q}_0$ fluctuations around their mean values, while $\tau$ represents the time lag for a flow perturbation at the injector to produce a change in the flame heat release rate. Despite its simplicity, the $n-\tau$ model played a crucial role in enabling stability predictions [8] that matched many experimental observations. These approaches made it possible to predict the stability characteristics and dynamic behavior of the acoustic modes as functions of the selected values of the $n-\tau$ parameters.

In this context, Pr Candel's contributions were twofold. First, his work laid the foundation for active control strategies aimed at modifying the time delay $\tau$ or reducing the interaction index $n$ through feedback mechanisms that integrate sensing and actuation (Fig. 21(b)). He demonstrated how to realize these strategies by combining acoustic and optical sensors, designing new actuators, and developing customized control methodologies (Fig. 22). His reviews [306,307] remain seminal references in the field and have served as a springboard for subsequent advances [308–310].

On the other hand, a significant effort by Pr Candel and his team at Ecole Centrale Paris was dedicated to uncovering the variety of physical mechanisms that govern the global delay $\tau$ and define the strength of the coupling (represented by the index $n$) in flames of increasing technological complexity. This research effort ranged from studies of the response of laminar premixed flames to low-frequency acoustic waves in canonical, few-kilowatt setups operating at atmospheric conditions, to investigations of high-pressure cryogenic flames subjected to high-frequency waves in a multi-megawatt facility. Work before 2013 is summarized in [311], while later developments discussed in Section 7 focus mainly on the dynamics of annular combustion chambers [312]. The scientific approach integrates analytical modeling, carefully designed experiments, the development of novel actuators for acoustic forcing (see for example the Very High Actuator Modulator used to perturb flames transversely in liquid fueled rocket combustion chambers [313,314]), and high-fidelity LES tailored to capture acoustic wave propagation in reactive flows (see for example [37] for transverse instabilities in cryogenic conditions).
<figure><img src="image_23.png" /><figcaption>Fig. 21. (a) Combustion/acoustics feedback loop. (b) Active control loop to stabilize the system. From [303].</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22. Active control methods. (a) Open loop control with fixed controller parameters. (b) Closed loop control with fixed controller parameters. (c) Closed loop adaptive control. Source: Adapted from [306].</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 23. (a) Ramjet configuration studied in [315,316] with (b) reduced frequencies $\omega L/u_0$ associated to large acoustic amplifications as a function of normalized flame speed $S_d/u_0$ of the ducted flame [315] and (c) acoustic pressure, acoustic velocity and resulting velocity perturbation along the TFS simulated in [316] at normalized time $t = 0.18$. (d) Experimental setup studied in [317], with phase averaged Schlieren (e) and heat release rate (f) images.</figcaption></figure>

The following section retraces the major developments that led to the identification of key driving mechanisms and culminated in the introduction of the flame transfer function (FTF) as tool for instability predictions. The final section highlights how this framework has since been further refined and adopted by research teams and industries across the globe.

## 5.2. Elementary driving mechanisms

Similar to nozzles or turbine blades, the flame reaction layer can often be considered compact relative to the acoustic wavelength of the thermo-acoustic instability and modeled as a thin interface separating fresh reactants from burnt gases. This low-frequency approximation forms the foundation of the Thin Flame Sheet (TFS) model, first introduced by Marble and Candel [315] for studying a model ramjet combustor with a ducted flame stabilized on a bluff body that could be influenced by the resonant acoustic field (Fig. 23(a)). They showed that, within a linear framework, strong flame responses occur at specific reduced frequencies, defined as $\omega L/u_0$, where $u_0$ is the incoming flow velocity, $L$ the flame length, and $\omega$ the pulsation frequency of the acoustic wave. Furthermore, the amplitude of these flame responses was found to depend on the ratio of the normal flame speed to the local flow velocity $S_d/u_0$, which determines the shape of the flame (Fig. 23(b)). In their analytical study, the observed peak responses were attributed to disturbances convected along the flame front, which could

later be directly visualized through direct numerical simulations [316] (Fig. 23(c)).

When S. Candel joined Ecole Centrale Paris and the EM2C laboratory, he set out to develop new experimental techniques and innovative simulation tools, specifically designed to study combustion instabilities and flame responses to acoustically externally forced conditions. One of the breakthroughs in the late 1980s was the successful use of phase-averaged analysis to investigate vortex-driven, acoustically coupled combustion instabilities in a multi-slit turbulent combustor [317] (Fig. 23(d)). Photomultiplier tubes and microphones were used to map the local unsteady heat release rate, allowing the distribution of the Rayleigh thermoacoustic source term to be deduced (Figs. 23(e,f)). Simultaneously, one of the earliest direct compressible numerical simulations of unsteady reactive flow enabled the visualization of convective disturbances along the flame front (Fig. 23(c)). In this study, the flame was not resolved directly but modeled using the TFS approach [316]. These developments demonstrated that the dynamics of combustion instabilities could be reasonably well predicted using reduced-order models based on a linear acoustic framework provided some information about the flame response is known. The central challenge lies in accurately identifying and characterizing the mechanisms that govern the flame response, specifically how to model the interaction index $n$ and the global time lag $\tau$.

This then motivated systematic investigations of canonical flame configurations to understand how incoming flow perturbations from
<figure><img src="image_26.png" /><figcaption>Fig. 24. Left: elementary processes leading to a delayed heat release rate response [303]. Right: Interfering mechanisms governing heat release rate fluctuations in swirling premixed flames [318].</figcaption></figure>

the injector lead to a delayed flame response [303,318,319]. The main elementary processes are summarized in Fig. 24(a). Among them, vortical disturbances generated in acoustically perturbed shear layers [320, 321], as well as mixture composition fluctuations [322] arising from the differential response of fuel and air supply lines to pressure disturbances, were shown to play a central role. Furthermore, interactions between the flame and solid boundaries [323], as well as the mutual interactions between neighboring flame fronts [324], were found to significantly impact the flame surface area. These interactions lead to substantial surface area destruction rates, which are closely linked to intense noise emission as described in Section 4 and play a crucial role in triggering instability [260].

Finally, the acoustic response of the injector itself was shown to generate flow perturbations that must be taken into account to understand the overall delayed flame response [318]. This mechanism is often called 'mode conversion' [325] and is particularly important when the injector includes a swirler, which converts incident acoustic waves into azimuthal vortical disturbances caused by flow separation at the swirler trailing edge (Fig. 24(b)). These vortical modes can subsequently interfere with incident acoustic waves, potentially leading to constructive or destructive interferences in the flame response [326, 327]. Later investigations demonstrated that the bluff body responsible for anchoring the flame plays a pivotal role by limiting the motion of the leading flame edge [328]. More recently, this mechanism of acoustic-to-vorticity mode conversion has been deliberately exploited to tailor the flame response to incident flow perturbations to specific ranges of forcing frequencies [329].

### 5.3. Flame Transfer Function

In many burners, the combustion reaction layer is inclined with respect to the main flow direction. Flames subjected to acoustic waves have been shown to be significantly more sensitive to velocity disturbances than to pressure fluctuations. Pressure effects become relevant primarily for normal flames at very high frequencies, i.e. when the acoustic wavelength approaches the flame thickness and the compact flame approximation breaks down [330].

At low Mach numbers, and under the assumption of compact flames, the dynamics of a TFS is governed by the displacement and deformation of the flame interface in response to periodic flow perturbations. This interface can be represented by a level-set function $G(x, t) = 0$, which separates fresh reactants ($G < 0$) from burnt gases ($G > 0$), and moves at an absolute velocity $w = v + S_d n$. In other words, the flame is convected by the flow at the velocity $v$ and propagates toward the unburnt gas along the normal vector $n = -\nabla G / |\nabla G|$, with a flame speed $S_d$, as described in Section 3. The flame speed $S_d$ depends on the fuel-air mixture and is influenced by thermodynamic conditions, local mixture composition, and local flame and flow properties—including effects of flame stretch and turbulence.

For small perturbations of the flow velocity, $v_1 = v - v_0$, and small disturbances in flame speed, $S_{d1} = S_d - S_{d0}$, the resulting perturbation in flame position, defined as $G_1(x, t) = G(x, t) - G_0(x)$, satisfies the linearized transport equation [303]:

$$ \frac{\partial G_1}{\partial t} + v_0 \cdot \nabla G_1 = w_1^{\perp} |\nabla G_0| \quad (51) $$

This expression highlights two key features of flame sheet dynamics: (i) perturbations generated at a point on the flame front are convected along the front with velocity $v_0 = v_0 - (v_0 \cdot n_0)n_0$, which corresponds to the component of the mean flow velocity tangential to the mean flame front; (ii) only absolute velocity perturbations normal to the flame front, given by $w_1^{\perp} = w_1 \cdot n_0 = S_{d1} + v_1 \cdot n_0$, directly influence flame displacement.

For a single-valued function $G_1(x, t) = \xi(X, t)$, integration of Eq. (51) along an inclined TFS shows that the flame displacement $\xi$ at a curvilinear position $X$ along the flame and at time $t$, relative to the unperturbed steady flame, results from the cumulative effect of absolute velocity disturbances normal to the mean flame, $w_1^{\perp}$, and the displacement at the flame root, $\xi_0$. These contributions are delayed by the time required for the mean flow to advect the perturbations to the current flame location [331]:

$$ \xi(X, t) = \frac{1}{v_0^{\parallel}} \int_{0}^{X} w_1^{\perp} \left( X', t - \frac{X - X'}{v_0^{\parallel}} \right) dX' + \xi_0 \left( t - \frac{X}{v_0^{\parallel}} \right) \quad (52) $$

Integration of this displacement along the flame front yields the resulting perturbed flame surface area $\int dA_1$, from which resulting heat release rate perturbations can be inferred at globally lean conditions [332,333]:

$$ \frac{\dot{Q}_1}{\dot{Q}_0} = \frac{1}{\rho_f} \frac{\int \rho_f dA_0}{\int dA_0} + \frac{1}{S_{d0}} \frac{\int S_{d1} dA_0}{\int dA_0} + \frac{\int dA_1}{\int dA_0} \quad (53) $$

where the first term, associated to perturbation of the fuel density $\rho_f$ can be neglected in low Mach flows. Fig. 25 exemplifies how four color Schlieren experiments, direct numerical flow simulations and predictions by integration of the nonlinear G-equation compare for the case of a premixed conical flame submitted to an harmonic flow modulation.

Within the linear framework (Eqs. (52) and (53)), the various contributions to heat release rate disturbances $\dot{Q}_1$ on the right-hand side of Eq. (53) can all be related to flow velocity perturbations $u_1$ at the burner outlet. It has proven to be highly effective, enabling, for example, the determination of the impulse response of flames in the time domain [334]. But generally it is used to analyze the forced flame response to harmonic perturbations of the form $a_1 = Re[\tilde{a}_1 \exp(-i\omega t)]$ at angular frequency $\omega$ in the frequency domain:

$$ \frac{\dot{Q}_1}{\dot{Q}_0} = G(\omega) \exp(i\varphi) \frac{\tilde{u}_1}{u_0} \quad (54) $$
<figure><img src="image_27.png" /><figcaption>Fig. 25. Response of a premixed conical flame submitted to harmonic modulation at reduced frequency $\omega L/u_0 \approx 15$ and modulation level $u_1/u_0 \approx 0.2$ [303]. (a) Schlieren images. (b) Direct numerical simulations. (c) Predictions from Eq. (51) for velocity perturbations advected by the mean flow [323].</figcaption></figure>

where G and $\varphi$ denote, respectively, the gain and phase lag of the complex Flame Transfer Function (FTF). This linear approach enables the derivation of expressions for the $n-\tau$ parameters under harmonic flow perturbations and shows that these parameters are inherently frequency-dependent.

Marble and Candel [315] were the first to propose using the G-equation to describe flame response to acoustic perturbations, based on rigorous flame dynamics calculations. While their approach was analytically tractable, obtaining a fully closed-form solution was not possible due to the need to solve simultaneously for both velocity and flame position. This study was followed up a few years later by Yang and Culick [335] using similar analytical developments. However, the topic saw little further development until the influential study by Fleifil et al. [336]. In contrast to earlier efforts, they assumed an isothermal flame by directly imposing the velocity disturbances, rather than solving for velocity and flame position concurrently. This simplification removed the need for the transverse integration required in [315, 335], enabling analytical progress with pen-and-paper calculations and sparking hundreds of subsequent studies.

Fleifil et al. [336] derived an expression for the FTF of elongated conical premixed flames followed up by Dowling [337] for a ducted V-flame submitted to low frequency acoustic disturbances. Theoretical predictions for conical flames were experimentally validated at EM2C by Ducruix et al. [338]. The relative ease with which FTF could be determined [323], combined with the development of flow forcing and FTF measurement techniques [339], and also improved signal post-processing methods for extracting FTFs from time-domain simulations [340], led to the rapid adoption of these tools for analyzing flame responses in increasingly complex flow configurations and under more sophisticated flow perturbation scenarios.

The flame response was confirmed to be highly sensitive to flame geometry. Certain configurations act as filters that attenuate incoming perturbations, while others, such as M and V-shaped flames, can amplify disturbances over specific frequencies [339,341]. Systematic efforts have been made to derive analytical expressions for the FTF

in canonical flame geometries, including dihedral, V-shaped, and M-shaped flames [342]. Parallel experimental investigations, primarily conducted at EM2C, have focused on measuring FTFs for these configurations [339].

Concurrently, research has expanded to consider various types of flow perturbations, such as advected coherent structures [343], fluctuations in mixture composition [344,345], and the effects of flame anchoring dynamics and unsteady heat transfer [346–348]. On the flame dynamics description side, attention has been given to incorporating the influence of flame stretch [349], turbulence [350].

Early experimental studies revealed that nonlinear effects were dominated the flame response as the amplitude of flame motion increased and the instability evolved toward a limit cycle [317]. This recognition prompted efforts to develop models that could capture this essential behavior [351]. As a result, researchers began investigating the impact of the forcing amplitude in externally forced flames [320, 339], which, in the context of principal harmonic analysis of the flame response, led to the development of the Flame Describing Function (FDF) framework [337,352], as will be discussed in Section 6.

The case of swirling V- and M-shaped flames has received particular attention, as these configurations are representative of flames aerodynamically stabilized in gas turbine combustors, which are especially susceptible to thermoacoustic instabilities [327,353–357]. These studies confirmed that the reduced frequency $\omega L/u_0$ could be used to scale FTF and FDF data for different operating conditions and fuels. It was also shown that the FTF of swirled V-flames anchored on a bluff body can exhibit both constructive and destructive interference patterns at specific frequencies [358]. Further studies investigated the influence of helical structures mimicking the Precessing Vortex Core (PVC) [359], the impact of confinement [360], slight anchoring flame asymmetries [361,362], and the effects of transverse acoustic forcing [363, 364]. The latter is particularly relevant for high frequency combustion instabilities where coupling with transverse acoustic modes becomes significant, such as in liquid fueled rocket engines (see Section 2) and, to a lesser extent, in annular combustion chambers (see Section 7).
More recently, these approaches have been extended to address the frequency response of diffusion-controlled flames as well [365].

These efforts laid the groundwork for a robust theoretical framework. They enabled identification of the key dimensionless parameters that govern the FTF, providing a basis for scaling the flame response and extrapolating the results to other flow conditions. This research also informed the development of heuristic models, such as the distributed time-lag model [366], helping to minimize the number of modeling parameters required. The work of S. Candel made it possible to develop these tools, validate them, and extend their capabilities to configurations increasingly closer to real injector setups. New challenges address extension of FTF concepts to high frequency instabilities (see for example [367]) and more generally to the problem of the acoustic response of non-compact flames.

# 6. Flame describing functions in combustion instability analysis

## 6.1. Introduction

The emergence of thermoacoustic instabilities in combustion chambers is a nonlinear process. Their onset can be predicted by linearizing the mass, momentum and energy equations governing chemically reacting flows, and using a suitable description of the two-way interaction between the acoustic field and the unsteady heat release rate of the flames. The key element of such linear stability analysis is the flame transfer function (FTF), which describes the frequency response of a flame when it is subject to small amplitude acoustic perturbations [333, 366]. The use of FTF addresses the crucial need of predicting the linear stability of a combustor after it has been equipped with a new burner, and the eigenfrequencies at which the perturbations are exponentially amplified. However, for linearly unstable conditions, using FTFs in a thermoacoustic network model cannot provide information on the amplitude and frequency of the self-oscillations. It also cannot be used to predict the co-existence of stable limit cycles and linearly stable non-oscillating state. For the development of practical combustors, this knowledge can however be crucial. It is therefore not surprising that the flame describing function (FDF) framework, which was presented in 2008 by Noiray, Schuller, Durox and Candel [352] and which enables predictions of limit cycle amplitudes and nonlinear phenomena such as instability triggering, received significant attention from the research community working on thermoacoustic instabilities.

The research efforts on nonlinear aspects of combustion instabilities preceding the above-mentioned article are now briefly presented. Early studies on the nonlinearities of thermoacoustic oscillations, such as the ones from Culick and co-workers, e.g. [368], naturally aimed at finding approximate analytical solutions of the nonlinear dynamics using classic averaging methods [369]. In rocket engines, the sound pressure level reached during combustion instabilities is so high that the propagation of acoustic waves becomes a nonlinear process. The description of the corresponding nonlinear gas dynamics, which can lead to phenomena such as wave steepening and shock formation must be incorporated in the models of rocket engines instabilities. However, it was shown in [370] that these purely acoustic nonlinearities are not sufficient to explain the observed pulse-triggered thermoacoustic limit cycles in rocket engines, and that they can only be reproduced when a nonlinear model of the heat release rate response of the flame to velocity perturbation is used. In gas turbines combustors, where the sound pressure level associated to thermoacoustic instabilities is lower than in rocket engines, this nonlinear flame response becomes the main nonlinearity of the thermoacoustic system, while in most cases, linear acoustic propagation can be considered as a good approximation.

The crucial importance of the nonlinear nature of the coupling between acoustic velocity and heat release rate had been already identified by Dowling [351] in 1997 in the context of premixed flames. In 1999 [337], she showed, by combining a kinematic representation of an archetypal premixed flame and a simple acoustic model of the

duct containing this flame, that the delayed response of the flame to finite amplitude acoustic perturbations is nonlinear. She performed time-domain simulations of the onset of oscillations for a set of model parameters corresponding to a linearly unstable situation, and of the saturation of these oscillations to a limit cycle. Furthermore, she introduced for the first time the describing function analysis [371], a technique from control engineering, in the field of combustor thermoacoustics. The FDF is a quasi-linear approach, which enables prediction not only of linear stability borders, but also of amplitudes and frequency of limit cycles. In [337], she showed that the decrease of the flame response gain when the amplitude increases was responsible of the saturation to a stable limit cycle. About ten years later, another formulation of the FDF was proposed in [352]. In contrast with [337], where the nonlinear frequency response is expressed as the product of the linear response and an effectively real amplitude dependent gain, it was demonstrated experimentally in [352] that not only the FDF gain can nonlinearly depend on the forcing amplitude but also the FDF phase. This observation motivated the following general formulation of the FDF:

$$F(\omega_r, |u'|) = \frac{Q'/\bar{Q}}{u'/\bar{u}} = G(\omega_r, |u'|) \exp [i \varphi(\omega_r, |u'|)], \quad (55)$$

where $\omega_r$ is the real part of the angular frequency, $|u'|$ is the modulus of the acoustic velocity fluctuations at the base of the flame, $\bar{u}$ is the mean flow velocity, $Q'$ and $\bar{Q}$ respectively denote the fluctuations and the mean of the heat release rate, and $G$ and $\varphi$ are the FDF gain and phase. It was then used for predicting the amplitude and frequency of limit cycles emerging from sub- and super-critical Hopf bifurcations. Remarkably, it enabled the prediction of experimentally observed phenomena associated to the nonlinear FDF phase such as frequency shift during transient to limit cycle, mode switching when triggering amplitude of a low frequency mode is reached from the growth of a linearly unstable high frequency mode, or hysteresis of stable thermoacoustic state when the combustor geometry is varied (see Fig. 26a). Furthermore, this FDF framework enabled the identification of different types of trajectories taken by the thermoacoustic system in the amplitude-frequency space. Indeed, it was shown that stable limit cycles can be reached when the nonlinear gain becomes too low, or when the nonlinear phase approach the stability criterion, and that triggering is caused, for these flames, by the nonlinear phase (see Fig. 26b). The FDF framework presented in [352] therefore constituted a significant advancement for the modeling and prediction of nonlinear phenomena in thermoacoustic instabilities.

## 6.2. Contributions from EM2C laboratory

A selected subset of the key achievements of Sébastien Candel and his coworkers with regard to the FDF framework is now discussed. Following the findings presented in [352], the FDFs of four distinct fully premixed flame geometries were examined in [339]: a single conical flame, a V-flame, an M-flame, and a collection of small conical flames stabilized on a perforated plate, all with the same fuel and air mass flows. The study highlights that the nonlinear characteristics of the flame response to finite amplitude acoustic forcing significantly depend on their geometry. For the conical flame, the gain is relatively insensitive to modulation levels, while the phase changes linearly with frequency at low amplitudes and saturates at higher amplitudes. In contrast, the FDF gain of the V-flame decreases as the amplitude increases while the phase evolves quasi-linearly with frequency and is largely amplitude-independent. For the collection of small conical flames, the gain decreases with increasing amplitude and the slope of the phase increases with amplitude, indicating an augmented effective time lag. These different FDF characteristics can lead to frequency shifts during thermoacoustic transients, mode switching, and hysteresis observed in e.g. [352,372], and these phenomena can be effectively predicted with the FDF framework.
<figure><img src="image_28.png" /><figcaption>**Fig. 26.** Comparison between FDF-based predictions of limit cycle amplitudes and experimental measurements from [352]. (a) Experimentally measured limit cycle amplitudes $u'_{\text{max}}/\bar{u}$ as function of the combustor plenum length $L$ superimposed on the predicted positive growth rates for the first three eigenmodes. (b) Collection of calculated positive growth rate trajectories superimposed on the FDF phase contours. It highlights the key importance of the nonlinear phase for the saturation amplitude of type III trajectories, which is obtained for vanishing growth rate (open symbols), and for the triggering threshold of type I trajectories (dark symbols). The stable limit cycles (open symbols) in the case of types I, II and IV trajectories are caused by the nonlinear saturating FDF gain.

Source: Figure adapted from [352].</figcaption></figure>

In [373], the FDF of fully premixed swirled flames is investigated. The findings of this study have been essential for the understanding of thermoacoustic instabilities involving swirled flames which are found in most of the practical gas turbine combustors. In parallel, the FDF framework was also successfully applied to predict combustion instabilities in turbulent premixed swirled flames [374] and the one of the corresponding FDF is shown in Fig. 27a.

While in [352,372,374] the FDF were used with combustors exhibiting a sufficiently simple geometry for accurately representing their acoustic field in 1D, practical combustors exhibit complex 3D geometries whose acoustic field cannot be described with analytical functions. In [375], the authors extended the FDF framework by combining a Helmholtz solver for the 3D acoustic field and FDFs (see Fig. 27a) in order to predict the thermoacoustic limit cycles of a swirled combustor. For the validation of this extended approach that can be used for 3D geometries, they considered with a combustor having a simple 1D geometry, which allowed them to compare the growth rate trajectories of the combustor's fundamental mode to analytical results. They also show that the prediction of the linearly unstable modes and of the corresponding stable limit cycles crucially depends on the modeling of the system's acoustic dissipation. More recently, the FDFs obtained from single and multiple injector setups [376] were investigated. This research was based on a new experimental setup called TICCA-Spray, which features a linear arrangement of three injectors, to study the flame response. The central flame is surrounded by two identical side flames, mimicking the conditions of an annular combustor. It is shown that the FDFs from single-flame experiments can be used as a first step for predicting thermoacoustic oscillations of annular combustors featuring a collection of flames, but that they may not fully capture the differences in boundary conditions between the single-flame system and the annular one. In particular, the study compares FDFs from two configurations with different swirler types and swirl directions (co- and counter-swirl). Differences in FDFs between single and multi-injector setups (see Fig. 27b) depend on the interaction between flames and

on their spacing. One of the most recent examples of the contribution from Candel and his co-workers on the subject of FDF framework is the paper [377], in which they investigate how explicit rotational symmetry breaking of the flames in the annular combustor MICCA affects the thermoacoustic instabilities. They use FDFs to predict the linear stability and the limit cycle amplitudes in this annular combustor. It is also worth to mention the recent review paper from Schuller, Poinsot and Candel [333], which provides a comprehensive overview of the modeling strategies of thermoacoustic instabilities in combustion chambers, with emphasis on the prediction capabilities enabled by the use of FDFs.

### 6.3. Contributions from other research groups

The early work on FDF at the EM2C laboratory has been followed by contributions from other research groups. The study presented in [357] is an excellent example of experimental research on FDF, where the behavior of a swirl-stabilized flame subject to high-amplitude simultaneous perturbations of equivalence ratio and velocity was investigated. The authors used a decomposition ansatz for the FDF that separates contributions from the velocity perturbations and from the equivalence ratio ones. In particular, they show that the forcing amplitude influences the phase between the responses of the flame to these two types of perturbations, which ultimately influence the overall response of the flame to the acoustic excitation. Recently, in [378], the effect of hydrogen enrichment on the FDF of V-shape swirled lean premixed methane/air flames was studied experimentally. This study examines two different hydrogen injection schemes: premixed with the main methane/air flow and injected pure as a pilot jet directly at the flame base. It is shown that the corresponding FDFs significantly differ, despite similar flame shapes. In the fully premixed strategy, the frequency range where the flame responds strongly increases as the flame shortens due to the higher reactivity of the hydrogen-enriched mixture. Conversely, the pilot injection strategy leads to a rebalancing of the heat
<figure><img src="image_29.png" /><figcaption>Fig. 27. Examples of measured FDF. (a) FDF of the swirled premixed flame investigated in [374] and in [375]. (b) FDFs measured in [376] showing the significant differences of nonlinear response of a spray flame caused by the presence (TICCA) or not (SICCA) of neighboring flames.</figcaption></figure>

release rate distribution along the flame brush, with higher reaction rates near the flame base, which reduces the FDF gain. These findings highlight the complex interactions between hydrogen content, injection schemes, and flame dynamics, and show that localized low-power hydrogen injection may be considered to control the instabilities.

Numerical research work on FDF has also been extensive in the last decade. For instance, in [379], the authors investigated the linear and nonlinear flame dynamics in the second stage of a sequential combustor. In contrast with previous FDF studies, which were dealing with relatively cold mixtures corresponding to propagating flames, this numerical work investigated the dynamics of hot mixtures prone to autoignite. This specificity makes the flame heat release rate not only sensitive to velocity and equivalence ratio perturbations, but also to temperature ones. For such technically premixed sequential burner, it was shown that the FDF gain governing the flame response to mixture temperature oscillations can increase with the amplitude of forcing due to autoignition kernels impinging on the flame (see Fig. 28a). This highly nonlinear gain increase can lead to instability triggering and sub-critical Hopf bifurcations. In [380], the authors integrate FDF simulations into a thermoacoustic network model. The study quantitatively predicts the onset and the amplitude of limit cycle oscillations.

Another example of numerical study on FDF is [383]. The use of the FDF framework was limited to studies on combustors operated at atmospheric pressure until this article was published. This study demonstrated the prediction of thermoacoustic limit cycles of a turbulent swirled flame combustor at pressures of 3 and 6 bar. The FDF was obtained with incompressible LES and the instability prediction results were validated with experimental data.

Important extensions of the classic FDF framework were also proposed in several studies. For instance, Ref. [384] presents a methodology for achieving a state-space realization of a given describing function. This approach allows for the conversion of the frequency domain information of the FDF obtained numerically or experimentally into a nonlinear time-domain description of the thermoacoustic system. This time-domain translation of a given FDF facilitates the study and simulation of the system dynamics using state-space methods, which are well-established in control theory. In [385], the authors propose to combine FDFs with three dimensional Green's functions in order to

perform parametric studies for situations where one dimensional description of the thermoacoustic problem would fail. Another extension of the FDF approach is the Flame Double Input Describing Function (FDIDF) proposed in [381]. The authors of this work address situations where thermoacoustic oscillations are quasiperiodic, involving multiple incommensurate frequencies. This extension allows for a more comprehensive analysis of the nonlinear response of flames by using two amplitudes and two frequencies, albeit at a higher computational cost than the traditional FDF. The FDIDF is embedded into a thermoacoustic network to predict the nature and amplitude of these oscillations using the harmonic balance method. The study outlines a criterion for the stability of these oscillations and compares the results with classical FDF analysis and self-excited time domain simulations. The FDIDF improves stability predictions over the FDF by identifying the onset of Neimark–Sacker bifurcations and determining the frequency of oscillations around unstable limit cycles (see Fig. 28b).

In [382], the role of higher harmonics in thermoacoustic limit cycles is studied. The work focuses on improving the predictive capability of the FDF framework by including additional transfer functions relating higher harmonics of the heat release rate to the forcing velocity. No additional data are needed to compute the extended FDF (xFDF), which provides accurate predictions in situations where higher harmonics play an important role in the thermoacoustic limit cycles and where the classic FDF approach would underestimate the overall sound pressure level. Indeed, the implementation of a xFDF in a thermoacoustic network enables the prediction of the amplitude and phase at the fundamental frequency and also of those of the harmonics, allowing to reconstruct the full coherent oscillations (see Fig. 28c). In [386], essential ingredients for low order modeling of flame response delay and nonlinearities are presented and applied to a swirled turbulent premixed flame, and allow prediction of intermittency of the thermoacoustic dynamics. Regarding the nonlinearities, it is shown that both of the FDF gain and phase are nonlinear, and that the acoustic spectrum of the limit cycle can be reproduced only if the higher harmonics of the delayed flame response are correctly accounted for.

Recently, machine learning algorithms were considered in combination to the FDF approach [387]. The authors investigated how their algorithms performed in interpolating and extrapolating xFDFs using
<figure><img src="image_30.png" /><figcaption>Relative forcing amplitude 8.5% leads to FDF gain of 2</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>b) Bifurcation diagram from FDIDF</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>c) Reconstructed time traces from xFDF</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Fig. 28.** (a) Nonlinear flame response to 125 Hz temperature oscillations at the inlet of a sequential combustor from [379]. Left: Normalized time signals of volume integrated heat release rate fluctuations $\dot{Q}'$ and temperature fluctuations at domain inlet $T'$. Right: Instantaneous heat release rate fields at maximum positive heat release rate oscillations. (b) Bifurcation diagram obtained from the FDIDF in [381] superimposed to the probability density function of velocity fluctuations obtained from time-domain simulations of the thermoacoustic dynamics. (c) Simulated time series for a combustor plenum length $L = 0.18$ m (solid gray line) compared to prediction by xFDF (solid black line) and prediction by standard FDF (dashed line) from [382].</figcaption></figure>

both experimental and simulation data. They showed that cubic spline interpolation and Gaussian processes (GPs), which also provide confidence estimates, are the preferred option. Their study also highlights that extrapolation performance deteriorates rapidly as the distance from the known domain increases.

Finally, in [388], a novel neural network architecture is proposed in the context of thermoacoustic instability analysis based on FDF. This multilayer perceptron-convolutional autoencoder (MLP-CAE) combines the strengths of multilayer perceptrons and convolutional autoencoders, and it enables to quantitatively capture spatial and temporal features of a forced laminar flame, as well as its FDF. The fact that the local flame dynamics is also accessible with the MLP-CAE opens the way for the prediction of instabilities in situations where the flame is not compact with respect to the acoustic wavelength of the instability peak. This research highlights the potential of advanced neural network architectures for FDF-based thermoacoustic instability analysis.

## 7. Combustion dynamics and instabilities of annular systems

Most studies on combustion dynamics focus on single flames in simplified geometries, as discussed in Sections 5 and 6. Fifteen years ago, it became clear that time had come to examine more complex configurations and in particular annular geometries that are most often found in gas turbines. In 2008, Sébastien Candel's team at EM2C Laboratory secured initial funding to develop a laboratory-scale annular system with large optical access. The setup enabled detailed analysis of the pressure field inside the plenum and combustion chamber together with direct visualization of the multiple injection units. These injection units were simplified models based on real-world devices. The facility, named MICCA, was fully designed by 2010 and became operational in 2012 [389,390].

Notably, during the same period, similar advancements were achieved at the University of Cambridge by Nicolas Worth and James Dawson [391]. A few years earlier, a separate facility was established at the Technical University of Munich under the leadership of Thomas Sattelmayer [392,393], albeit with limited optical access. As is often the case, theoretical investigations in this field preceded experimental developments, beginning around the year 2000 and subsequently

pursued across several institutions [394–396]. On the numerical front, Thierry Poinsot and his team at Cerfacs were the first to demonstrate that azimuthal instabilities could be captured using LES. Their pioneering work in 2009 [397] reveal a spinning instability in a small combustion chamber representative of those used in helicopter engines. Complementing these numerical advances, the MICCA facility provided the combustion community with invaluable experimental insights, as well as high-quality data essential for the development of reduced-order models and high-fidelity simulation tools.

This brief overview highlights the early milestones of over a decade of research on annular systems, which has significantly deepened the understanding of their instability mechanisms. Central to this effort has been the EM2C team, working collaboratively around the MICCA facility. The research was led by S. Candel's group, with D. Durox heading the experimental activities, alongside key contributions from T. Schuller, R. Vicquelin and more recently from A. Renaud, and several doctoral students. The following section summarizes major achievements and the collaborative efforts that made them possible.

Annular configurations are of particular importance, as they are widely used in aircraft engines and heavy-duty gas turbines. Ground-based gas turbines often favor can-type or hybrid can-annular designs, primarily due to ease of maintenance. However, in aerospace propulsion, the annular geometry is especially prevalent due to its compactness and performance advantages. As new engine architectures and advanced combustion concepts are developed to meet environmental regulations and reduce pollutant emissions, a range of flame dynamics challenges must be carefully addressed. In particular, achieving precise ignition control, ensuring reliable relight under high-altitude conditions, and preventing the onset of combustion instabilities across the entire operating envelope are critical considerations.

Ignition must be achieved safely while minimizing energy consumption. Due to the lack of laboratory-scale experiments, the light-round mechanism governing ignition in annular configurations remained poorly understood. As a result, the design and placement of ignitors relied heavily on trial-and-error testing and empirical knowledge. A significant breakthrough came with the first LES performed by CERFACS, which successfully demonstrated that light-round ignition could be numerically predicted in a helicopter combustion chamber [203].
<figure><img src="image_34.png" /><figcaption>**Fig. 29.** Comparison between experiment and numerical simulation during an ignition sequence in MICCA. Propane and air mixture at equivalence ratio $\phi = 0.74$, bulk velocity in each injector $U_b = 17.1 \text{ m/s}$ corresponding to a total power of $P = 52 \text{ kW}$. Images of flames are set in false colors to improve the visualization of the flame front. Yellow and dark red correspond to the highest and lowest values of light intensity. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Source: Adapted from [389,398].

However, it became evident that experimental validation was essential to support these simulations, necessitating a well-controlled laboratory-scale setup. This need motivated the development of the MICCA annular combustion chamber to provide a suitable platform for detailed experimental investigation and validate LES as described in Section 7.1.

A second distinctive feature of annular combustion chambers is their susceptibility to trigger combustion instabilities coupled to azimuthal acoustic modes. The MICCA annular chamber enabled in-depth investigation of these phenomena using combined acoustic and advanced optical diagnostic techniques. In particular, Section 7.2.1 introduces the spin ratio, developed to characterize the nature of azimuthal acoustic modes observed in such configurations. Key results are then synthesized in Sections 7.2.2 to 7.2.6 for different fuel injectors and fuel injection schemes. The key role played by Sébastien Candel in advancing the understanding of these complex phenomena and contributing to the broader knowledge in the field will be highlighted.

## 7.1. Ignition of annular combustion chambers

Until the early 2010s, there was no academic experimental study on ignition in annular combustors. The first experimental paper on the ignition in an annular combustion chamber was published in early 2013 by Bach et al. [399], and extended a few months later by the EM2C team using the MICCA annular chamber [389]. In this study, the flames were premixed and swirled, and they were stabilized aerodynamically by recirculation of the swirled flow generated by the injection units. Experiments reported in [389] provided new insight on the light round process: (i) ignition spreads in a smooth fashion from injector to injector (Fig. 29) without flame jumps and back propagation to one or more injectors of the kind observed in experiments on linear arrays of injectors [400]; (ii) expansion pushes gases residing in the chamber; (iii) the flame propagation velocity increases with the flow rate in each injector, reducing the total light round time. Despite different geometries and fuels, simulations in [203] already resembled the propagation observed in the MICCA chamber. Later simulations conducted at EM2C Laboratory of the exact geometrical configuration with the AVBP code from CERFACS led to a remarkable comparison (Fig. 29) between experiments and simulations [398].

In 2014, MICCA was modified to allow liquid fuel injection through ‘hollow-cone’ atomizer nozzles, while retaining the principle of swirling

the air flow to stabilize the flame aerodynamically. This new design that also allowed injection of air and propane mixtures, was used to study the influence of the fuel type on the light round process with two different liquid fuels of high and low volatility (n-heptane and dodecane) and compare results with those measured for propane-air mixture injection [401]. In particular, ignition was shown to be faster in the premixed case and that the total light-round time increases with the least volatile fuel.

LES were further pursued, enabling comparisons with experimental data and offering deeper insight into the phenomena occurring within the annular chamber. In particular, two key processes were elucidated. First, the flushing of fuel droplets on the upstream side of the flame front — triggered by the ejection of fresh gases — was shown to create local fuel concentration heterogeneities, resulting from the velocity mismatch between the fresh gases and the droplet spray. Second, a brief flashback phase was identified, during which the flame propagates into the injection unit as combustion of the air-fuel spray is initiated [402]. High-speed visualizations showed that the flame entered the injector before being ejected from this unit and taking its final shape above it. This effect is due to the back pressure induced by a transient increase in pressure in the chamber each time a burner is ignited. These two phenomena could be well retrieved by LES [403,404]. Effects of temperature of the backplane and of the lateral walls on light round ignition could also be reproduced by LES [405,406]. The data, theoretical interpretation, and simulations were primarily guided by the ideas of Sébastien Candel and Ronan Vicquelin. This foundation enabled the development of a simplified model describing the variation in flow rate through a burner, driven by the back pressure peak that is induced by ignition. The model also accounts for the influence of heat transfer on the overall process.

## 7.2. Thermo-acoustic instabilities in annular combustion chambers

In many high performance systems, the most dangerous instabilities are coupled by transverse acoustic waves. Combustion instabilities of this nature first emerged during the early development of rocket engines, with several failures attributed to the coupling between combustion dynamics and the azimuthal or radial acoustic modes of cylindrical thrust chambers. A notable incident occurred during the second launch of the Ariane rocket on May 23, 1980, when a failure
<figure><img src="image_35.png" /><figcaption>Fig. 30. Annular combustion instabilities calculated in a helicopter combustion chamber with the solver AVBP at Cerfacs [397,420,421]. The radius is approximately 20 cm, and the chamber is equipped with 15 burners. The two images present two different instants of the pressure field. The azimuthal pressure oscillations induce axial flow fluctuation through each burner.</figcaption></figure>

caused by combustion instability led to the mission's loss. Sébastien Candel was part of the investigative panel tasked with determining the root cause of the accident. The panel concluded that the failure resulted from an azimuthal resonant mode within the thrust chamber, which caused the propellant injection system to melt down. Following this diagnosis, modifications were made to the Viking engines, and the Ariane launcher successfully returned to flight in June 1981.

In aeroengines and gas turbines equipped with annular combustors, azimuthal modes tend to be less effectively damped. Because purely azimuthal modes have eigenfrequencies that are inversely proportional to the mean perimeter of the chamber, their resonances often fall within the frequency range where flames are most sensitive to perturbations. These challenges have become more pronounced with increasing performance demands and stricter environmental regulations. Analytical studies of thermoacoustic instabilities in annular systems began around the turn of the century [394,407], while the first experimental evidence of standing and spinning azimuthal modes in heavy-duty gas turbines was reported in [408].

Around a decade later, a comprehensive study based on measurements in a heavy-duty gas turbine revealed that the self-sustained azimuthal modes are a mixture of both spinning and standing modes, rather than purely one or the other. The study also highlighted that the modal dynamics are significantly influenced by stochastic forcing from turbulent combustion noise [409]. Numerous theoretical studies have since addressed azimuthal acoustic modes and related combustion instabilities [395,396,410–412]. Many of these investigations have focused on predicting the nature of the modes (standing, spinning, or mixed) as a function of instability amplitude and injector parameters [409,413–415]. Concurrently, Helmholtz solvers combined with FTF models have been developed to ease these predictions [416–419].

The first LES was reported in 2009 [397,420,421]. While the results were impressive (see Fig. 30), experimental validation was still lacking. This simulation revealed that the azimuthal pressure wave in the combustion chamber induces axial flow oscillations at the injector outlets. At that time, the only well-instrumented laboratory-scale annular combustor was the TUM ring, developed by Professor Thomas Sattelmayer [422,423]. This slightly pressurized annular chamber featured small portholes, which provided limited access to the combustion region.

Two atmospheric pressure facilities were independently designed at Cambridge [424] and at EM2C [390] to enable pressure field measurements along with optical access to the flame region. The Cambridge configuration was further developed after Professor Dawson and Dr. Worth moved to NTNU in Norway. They then successfully develop a high-pressure facility at NTNU which is operational since 2021 [425], and numerous other annular combustors have been commissioned in China [426–428] and India [429] after 2018.

### 7.2.1. Definition of the spin ratio

One aspect that has received much attention has been that of the modal identification in annular systems. With a flow crossing the combustion chamber in the axial direction and low velocity components in the azimuthal direction, the pressure field is formed by the sum of four waves traveling in the upstream or downstream direction, and having a direction of rotation clockwise (CW) or counter-clockwise (CCW) in the annular chamber [395]. These waves combine to give different modes of oscillation in the chamber. Limit cases are the spinning CW mode and its opposite, the spinning CCW mode. When CW waves have exactly the same frequency and amplitude as CCW waves, they combine and the mode is standing, with an oscillation of the pressure on both sides of a nodal line. If the amplitudes of CW and CCW waves are not identical, the mode is 'mixed'.

Classification of these modes is facilitated by a characteristic parameter known as the spin ratio, as introduced by Evesque et al. [395]. Defined based on energy considerations, the spin ratio takes a value of 1 for purely spinning modes, regardless of their direction (clockwise or counterclockwise), and a value of 0 for purely standing modes.

The ambiguity in the rotation direction was resolved later on by Sébastien Candel with a revised definition of the spin ratio that provided the missing information. It was first used to characterize the azimuthal instabilities arising in the annular MICCA combustion chamber [390]. The pressure field analysis was simplified by assuming that the mean flow in the azimuthal direction could be neglected in view of the large difference between flow velocity and speed of sound. For the first azimuthal mode and assuming time variation for the pressure the pressure field can be written [430]:

$$p(\theta, t) = \text{Re}\{A_{+}e^{i\theta - i\omega t + i\phi_{+}} + A_{-}e^{-i\theta - i\omega t + i\phi_{-}}\} \quad (56)$$
where $A_+$ and $A_-$ are positive real amplitudes, the spin ratio defined as follows:

$$s_r = \frac{A_+ - A_-}{A_+ + A_-} \qquad (57)$$

characterizes the relative levels of CW and CCW waves. When $s_r = 1$, the wave is spinning in the CCW direction, when $s_r = -1$, the spinning is in the CW direction, and when $s_r = 0$, the mode is purely standing.<fnref n="15" />

The $s_r$ index was subsequently adopted by other research teams [431,432]. In 2018, Ghirardo and Bothien [433] proposed a more general formulation, based on the Bloch sphere of quantum mechanics in combination with the quaternion formalism, to fully characterize azimuthal instabilities. The quaternion takes into account the amplitude of the instability, its nature (spinning, standing or mixed), the orientation of the nodal line and phase. The nature angle $\chi$, used to define the latitude on the Bloch sphere, is directly related to the spin ratio, so that the antinodal line location $\theta_0$ and phase $\varphi$ may be linked to the phases of the two waves appearing in the pressure field expression Eq. (56) through:

$$\tan \chi = s_r, \quad \theta_0 = \frac{1}{2}(\phi_- - \phi_+), \quad \text{and} \quad \varphi = \frac{1}{2}(\phi_- + \phi_+) \qquad (58)$$

Recently, this formalism was combined with the wave equation by Faure-Beaulieu and Noiray [434] to describe the combustion dynamics in an ideal annular combustor, unifying previous modeling approaches with the effects of time delayed flame response, asymmetries, mean azimuthal flow, and additive stochastic forcing from the turbulent combustion noise that is inherent to practical combustors. In [435, 436], the theoretical model presented in [434] was used to explain the predicted spontaneous reflectional symmetry breaking, and the unexpected explicit rotational symmetry breaking that were observed experimentally. Furthermore, a new type of modal dynamics characterized by periodic switching of the spinning direction was observed and presented in [437]. The origin of this so-called beating mode is theoretically explained by the inescapable presence of tiny asymmetries of the geometry, the mean temperature field, or the thermoacoustic response of the flames.

### 7.2.2. Instabilities with premixed swirling flames

Laboratory-scale experiments on instabilities in annular combustors, including flame visualizations, were first conducted at Cambridge by Worth and Dawson [424,438]. In these experiments, premixed swirled injectors were fitted with central bluff bodies to stabilize the flames. Azimuthal instabilities in the MICCA annular chamber were reported shortly thereafter at EM2C [390]. In this study, the premixed flames were stabilized aerodynamically by the swirling flow at the outlet of each injector. Initial tests with MICCA indicated stable operation, but Sébastien Candel proposed modifying the configuration by terminating each burner with a cup (bowl) to alter the flame shapes and enhance interactions between the flame fronts established by the different injectors. Invited to the EM2C lab by Thierry Schuller, James Dawson suggested varying the lengths of the quartz tubes forming the lateral boundaries of the chamber, a geometry that had been implemented in the Cambridge annular chamber [424]. These two modifications were introduced, leading to the onset of azimuthal instabilities of spinning, mixed, or standing types in MICCA. In [390], it was experimentally demonstrated that pressure oscillations at the injector outlet caused fluctuations in the flow rate, which, in turn, led to the roll-up of the flame front by a large vortical structures released from the injector.

### 7.2.3. Instabilities with laminar flame matrix burners

Daniel Durox proposed simplifying the configuration by replacing the swirling injection units with matrix injection plates to investigate

instabilities in a chamber with multiple laminar flames, established by a regular matrix of small orifices connecting the upstream plenum to the combustion chamber. This modification leveraged the expertise at EM2C on the dynamics of premixed laminar flames produced by matrix burners [352,439]. Previous studies on single matrix systems had provided valuable insights into nonlinear effects such as frequency shifting during transient phases, leading to limit-cycle oscillations, mode switching, instability triggering, and hysteresis phenomena that could be predicted within the FDF framework, as described in Section 6. Notably, it was found that these phenomena were linked to changes in the slope of the FDF phase as a function of perturbation amplitude. A set of 16 matrix burners was installed in the MICCA system, which was operated with premixed propane and air. Unstable operation was easily achieved, resulting in robust oscillations even with equal-length chamber walls. Depending on the flame settings, either pure standing or pure spinning modes could be observed for extended durations (Fig. 31) [440,441], enabling detailed investigations of flame motion. This setup allowed direct comparison of flame behavior at the nodal and antinodal locations when the unstable system exhibits a standing mode.

It was also shown that the final behavior of the system also depended on the path followed to get to the operating point [442]. An unusual mode of instability was unraveled corresponding to a combination of longitudinal and standing modes, both having the same frequency [443]. Sébastien Candel suggested to designate this as the 'slanted mode'. A close examination of the spinning and standing modes obtained with matrix burners, which are acoustically transparent, showed that for these modes, the strongest pressure disturbances were found in the plenum. When the spinning mode was present, it was found that the pressure waves in the plenum and in the combustion chamber were not in phase indicating, according to Sébastien Candel, that the time delay between the chamber and the plenum could be linked to the phase of the FDF, and this was included in a comprehensive model of the plenum, injectors and combustion chamber [441].

### 7.2.4. Instabilities with swirling spray flames

In 2014, to investigate aero-engine applications, the EM2C team opted to return to swirling injectors with liquid fuel injection, while maintaining the capability to operate under premixed air/propane conditions. The new version of MICCA, designated MICCA-Spray, enabled extensive studies on the combustion dynamics of swirl spray flames. A unique aspect of this approach was the establishment of a productive dialogue between measurements taken in a sector featuring a single injection unit (in a facility known as SICCA) and experiments conducted on MICCA.

Determining FDFs is simpler in the single-sector configuration than in the annular system. The single-sector geometry also facilitates a more straightforward analysis of the flame structure and is more suitable for detailed characterizations of the droplet spray. Additionally, it was hoped that measurements from the single-sector setup would provide insights into how to adapt the injection unit to promote unstable operation in the annular system, while also providing guidance on modifications required to achieve stable operation.

Recent investigations have benefited from the integrated operation of these facilities, involving former Ph.D. students: Kevin Prieur, Guillaume Vignat, Preethi Rajendram Soundararajan, and Véranika Latour, as well as current PhD candidate Nicolas Vaysse. Ongoing collaborations continue with Laurent Gicquel's team at CERFACS and Françoise Baillot's team at CORIA.

Current work aims at identifying the role of the injector itself on combustion dynamics. The injector not only governs the combustion quality in the chamber, largely determining combustion efficiency, operability, and pollutant emissions of the unit, but it also has a substantial impact on the combustion dynamics of the system, as revealed by FDF data. These data were used to derive stability predictions in SICCA, analyze the effects of different fuels, study the influence of head loss and swirl number, examine the impact of the atomizer position

<footnote n="15">¹⁵ The analysis that led to the definition of a new spin ratio was carried out on a beach in Juan-les-Pins.</footnote>
<figure><img src="image_36.png" /><figcaption>Fig. 31. Spinning mode in MICCA equipped with matrix burners. Eight images obtained by phase averaging showing a full cycle. $\phi = 0.96$ and $u_b = 1.49$ ms$^{-1}$. The direction of reading is from left to right and from top to bottom. From the first to the last image, the maximum light intensity evolves circularly counterclockwise.</figcaption></figure>

Source: From [441].

<figure><img src="image_37.png" /><figcaption>Fig. 32. Blow-off of several flames in MICCA-Spray due to a strong azimuthal instability when the acoustic mode is standing. For a few moments the flames are extinguished in the zone of the antinode acoustic velocity. The pressure nodal line passes between burners 1 and 16, and between burners 8 and 9. P = 117 kW; $\phi = 1.09$.</figcaption></figure>

Source: Courtesy of [447].

<figure><img src="image_38.png" /><figcaption>Fig. 33. Flame blow out boundary for three regimes of operation. Liquid n-heptane fuel. The blue symbols correspond to $U_b = 17.8$ m/s, f = 745–760 Hz, the green symbols to $U_b = 22.4$ m/s, f = 755 Hz and the red symbols correspond to $U_b = 26.7$ m/s, f = 750–762 Hz. Filled symbols indicate the maximum achievable level of $P'/\rho_0 c_0 U_b$.</figcaption></figure>

Source: Courtesy of [449].

relative to the injector outlet, and document changes in FDFs and instabilities when the fuel is a blend of two pure substances with varying mixture compositions [444–446].

### 7.2.5. Strong combustion instabilities and dynamical blow-off

High amplitude oscillations have long been thought to induce flame blow-off but this phenomenon is not well documented in annular combustors. When MICCA was equipped with a long outer wall (500 mm) and a much shorter inner wall (200 mm), an instability coupled by the 1A-1L combustion chamber mode led to violent oscillations with peak amplitudes of more than 5000 Pa [447,448]. High-speed video records indicated that, when the oscillation corresponded to a standing mode for a sufficient amount of time, the flames established near the pressure nodal line were blown-off temporarily. At the maximum of instability, nearly half of the flames were extinguished (Fig. 32) for a few periods leading to a reduction of the instability level accompanied by a re-ignition of the extinguished flames. This phenomenon could be repeated a few instants later, when the oscillation had again reached large amplitudes.

In Sébastien Candel’s view of this dynamical blow-off, the controlling parameter was the transverse acoustic velocity $v'_x$ acting on the flame in the lateral direction. Using dimensional analysis, the process could be described in terms of four dimensionless groups: $v'_x/U_b$, the transverse relative velocity perturbation amplitude, $b = (\tau_c U_b)/d$, the ratio of the chemical conversion time to a typical mechanical time, $St = \omega d/U_b$, the Strouhal number and the swirl number S. The initial version of the dynamical blow-off criterion formulated in [448] was later cast in a more practical form [449]:

$$ \frac{v'_x}{U_b} > \psi_3(S, St)(b_* - b)^n \qquad (59) $$

where n is an exponent determined from experiments, $\psi_3(S, St)$ is a function of the swirl S and Strouhal St numbers, and $b_*$ is a critical value corresponding to static blow-off and has a value of order of unity. Note that b is the inverse of a Damköhler number $b \sim 1/Da$. When b is equal to $b_*$ the flame is blown-off. Thus the difference $b_* - b$ represents a margin with respect to static blow-off.

Experiments carried out in MICCA were not used to test this model. It was done in collaboration with the CORIA laboratory in Rouen and
<figure><img src="image_39.png" /><figcaption>Fig. 34. (Left) Growth rate ($\omega_i$) maps obtained from FDF data. (Right) Amplitude stability map obtained in the annular combustor MICCA-Spray.</figcaption></figure>

their linear array facility, called TACC, equipped with three aligned injection units, that were similar to those used in MICCA. TACC may be submitted to high amplitude transverse oscillations using lateral driver units allowing variations of the relative transverse velocity $v'_x/U_b$ and a direct verification of the model (Fig. 33). It was found that the criterion in Eq. (59) suitably reproduces the flame blow-off at the velocity antinode for different values of the equivalence ratio [449].

### 7.2.6. Prediction of instabilities in annular combustors

Prediction of combustion instabilities has been a permanent objective of research in combustion dynamics, especially with reduced order models that could be used for design at engineering level. Many of these models rely on FTF or FDF concepts derived from control systems theory and the pioneering work of Luigi Crocco [304]. The EM2C team has tackled this problem by making predictions and comparing them with a large set of data from various experimental setups. A family of swirlers was designed allowing to vary the head loss and swirl number of the injection units. These injectors mounted in the single injector SICCA were then used to test theoretical predictions with respect to experimental data [450]. A good agreement with experiments on SICCA was obtained but it was found that this required a knowledge of the FDF, of the impedance at the flame and a reliable estimate of the damping rate. By shifting the location of the instability band, impedances of feedlines and exhaust pipes play a key role. In the low frequency range, the upstream impedance was found to be controlled by that of the injector unit [450]. The transfer matrices of injection units were measured and checked by making use of acoustic energy conservation principles [451]. At higher frequencies that are typical of the 1A1L mode in MICCA, it was found that the effective impedance at the flame $\zeta$ is also controlled by the presence of the rigid backplane and that it was important to be able to represent the impact of the injector and backplane [446]. Physically, the flame is submitted to the joint effects of the injector through its impedance $\zeta_i$ and of the rigid backplane which imposes an impedance $(i \tan(ka))^{-1}$ at the location a of the flame with respect to the backplane. S. Candel proposed to model these two influences by defining a bridging function g depending on the ratio between injector outlet surface $S_o$ and backplane surface area of a single injector $S_{bp}$, so that the effective impedance $\zeta$ in front of the flame becomes:

$$ \zeta = \frac{\zeta_i}{g + (1 - g) i \tan(ka) \zeta_i} \qquad (60) $$

where a is the distance between the chamber backplane and the flame barycenter. If $S_o/S_{bp}$ is small, this model yields an impedance $\zeta = (i \tan(ka))^{-1}$ which is that of a fully rigid backplane. When $S_o/S_{bp} = 1$, the injector outlet occupies all the surface area available in the backplane, $\zeta = \zeta_i$, and the injector imposes its impedance to the flame. Using this unified impedance model in combination with measured FDFs, experimental gathered in MICCA can be interpreted by changing

the fuel composition from pure heptane to pure dodecane to explain observations made when the combustor is fed with these various blends. This model was also used to derive instability predictions for MICCA relying on a full mapping of the FDF over a broad range of operation (for $0.75 < \phi < 1.05$ and $92 \text{ kW} < P < 118 \text{ kW}$) [452]. The maps retrieve the location and general layout of the unstable region determined experimentally in MICCA (Fig. 34).

This agreement is partly due to the fact that there flame-to-flame interactions between different injectors are limited in MICCA, and that the overall response of the burners can be considered as the sum of the responses of the individual burners. These latest results also show that a complete prediction of the stability maps of an annular combustion chamber becomes possible. The point that remains to be improved is the modeling of the impedances upstream of the flame fronts.

## 8. Conclusions

Through six fields of research, this paper demonstrates the impact of Pr Candel on the combustion community over the last fifty years. Many of his intuitions have become standard methods for our community. These developments have changed our understanding of multiple crucial mechanisms in combustion, the design of many real combustion systems and the numerical simulations of gas turbines, rocket engines, internal combustion engines and industrial furnaces.

Sébastien Candel’s contribution was considerable in all these studies. He always made sure to closely follow experimental tests, suggest ideas to resolve technical problems, propose theoretical concepts to model physical phenomena, provide mathematical methods for processing the results, developing the theory or performing simulations.

This section omits the strong interest and investment of Sébastien in numerical methods, high performance computing and data computing as well as his service to the French and international scientific and industrial communities.

Throughout these years, Pr Candel has also been an eminent member of the Combustion Institute, helping to organize and structure it, acting as Editor in Chief for Combustion and Flame for many years. He was also vice-president and then president of the French Academy of Sciences (2015–2018). We hope that the present paper will be useful for young researchers not only to identify these contributions but also to have a proper view of the field and guide them in their own work. We also hope that this paper will convey the admiration of all his colleagues for his remarkable contribution.

## References

[1] M. Barrere, A. Jaumotte, B.F. de Veubeke, J. Vandenkerckhove, Rocket Propulsion, Elsevier, 1960.

[2] R. Burick, Space storable propellant performance program coaxial injector characterization, NASA Rep. (1972).

[3] P. Liang, Liquid rocket combustor computer code development, NASA NAS8-34928 (1986) 696–716.
[4] G. Gill, A qualitative technique for concentric tube element optimization, utilizing the factor (dynamic head ratio - 1), American Institute of Aeronautics and Astronautics (AIAA), 1978, pp. 1-16,

[5] L. Vingert, M. Habiballah, J. Traineau, Mascotte, a research test facility for high pressure combustion of cryogenic propellants, in: Technical Report, Office National d'Etudes et de Recherches Aérospatiales (ONERA), 92 - Chatillon (France), 2000.

[6] G. Herding, R. Snyder, P. Scouflaire, C. Rolon, S. Candel, Flame stabilization in cryogenic propellant combustion, Proc. Comb. Inst. 26 (1996) 2041-2047.

[7] R. Snyder, G. Herding, J. Rolon, S. Candel, Analysis of flame patterns in cryogenic propellant combustion, Combust. Sci. Tech. 124 (1997) 331-370.

[8] F.C.T. Culick, Unsteady motions in combustion chambers for propulsion systems, AGARD AG-AVT-039, NATO, 2006.

[9] G. Herding, R. Snyder, C. Rolon, S. Candel, Investigation of cryogenic propellant flames using computerized tomography of emission images, J. Prop. Power 14 (1998) 146-151.

[10] D. Kendrick, G. Herding, P. Scouflaire, C. Rolon, S. Candel, Effects of a recess on cryogenic flame stabilization, Combust. Flame 118 (1999) 327-339.

[11] M. Juniper, N. Darabiha, S. Candel, The extinction limits of a hydrogen counterflow diffusion flame above liquid oxygen, Combust. Flame 135 (2003) 87-96.

[12] G. Ribert, N. Zong, Y. V., L. Pons, N. Darabiha, S. Candel, Counterflow diffusion flames of general fluids: Oxygen/hydrogen mixtures, Combust. Flame 154 (3) (2008) 319-330.

[13] L. Pons, N. Darabiha, S. Candel, T. Schmitt, B. Cuenot, The structure of multidimensional strained flames under transcritical conditions, Comptes Rendus Mécanique 337 (6-7) (2009) 517-527.

[14] L. Pons, N. Darabiha, S. Candel, G. Ribert, V. Yang, Mass transfer and combustion in transcritical non-premixed counterflows, Combust. Theor. Model. 13 (2009) 57-81.

[15] M.P. Juniper, A. Tripathi, P. Scouflaire, J.C. Rolon, S. Candel, Structure of cryogenic flames at elevated pressures, Proc. Comb. Inst. 28 (2000) 1103-1109.

[16] S. Candel, M.P. Juniper, G. Singla, P. Scouflaire, C. Rolon, Structure and dynamics of cryogenic flames at supercritical pressure, Combust. Sci. Tech. 178 (2006) 161-192.

[17] G. Singla, P. Scouflaire, C. Rolon, S. Candel, Planar laser-induced fluorescence of OH in high-pressure cryogenic LOx/GH2 jet flames, Combust. Flame 144 (2006) 151-169.

[18] Combust. Theor. Model. 7 (2003) 563-577.

[19] M.P. Juniper, S. Candel, Edge diffusion flame stabilization behind a step over a liquid reactant, J. Prop. Power 19 (2003) 332-341.

[20] M. Juniper, Ecole Centrale Paris, Structure and Stabilization of Cryogenic Spray Flames (Ph.D. thesis), Ecole Centrale Paris, 2001.

[21] G. Singla, P. Scouflaire, J. Rolon, S. Candel, L. Vingert, OH planar laser-induced fluorescence and emission imaging in high-pressure LOx/methane flames, J. Prop. Power 23 (2007) 593-602.

[22] G. Singla, P. Scouflaire, C. Rolon, S. Candel, Transcritical oxygen/transcritical or supercritical methane combustion, Proc. Comb. Inst. 30 II (2005) 2921-2928.

[23] J.C. Oefelein, V. Yang, Comprehensive review of liquid-propellant combustion instabilities in F-1 engines, J. Prop. Power 9 (1993) 657-677.

[24] Y. Méry, S. Ducruix, P. Scouflaire, S. Candel, Injection coupling with high amplitude transverse modes: Experimentation and simulation, Comptes Rendus Mécanique 337 (2009) 426-437.

[25] F. Richecour, S. Ducruix, P. Scouflaire, S. Candel, Effect of temperature fluctuations on high frequency acoustic coupling, Proc. Comb. Inst. 32 (2009) 1663-1670.

[26] Y. Méry, L. Hakim, P. Scoufalre, L. Vingert, S. Ducruix, S. Candel, Experimental investigation of cryogenic flame dynamics under transverse acoustic modulations, Comptes Rendus Mécanique 341 (2013) 100-109.

[27] F. Richecoeur, S. Ducruix, P. Scouflaire, S. Candel, Experimental investigation of high-frequency combustion instabilities in liquid rocket engine, Acta Astronaut. 62 (2008) 18-27.

[28] J. Oefelein, V. Yang, Modeling high-pressure mixing and combustion processes in liquid rocket engines, J. Prop. Power 14 (1998) 843-857.

[29] J. Oefelein, Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure, Combust. Sci. Tech. 178 (2006) 229-252.

[30] N. Zong, V. Yang, Near-field flow and flame dynamics of LOX/methane shear-coaxial injector under supercritical conditions, Proc. Comb. Inst. 31 II (2007) 2309-2317.

[31] S. Jay, F. Lacas, S. Candel, Combined surface density concepts for dense spray combustion, Combust. Flame 144 (2006) 558-577.

[32] T. Schmitt, J. Rodriguez, I.A. Leyva, S. Candel, Experiments and numerical simulation of mixing under supercritical conditions, Phys. Fluids 24 (2012) 055104.

[33] M. Gonzalez-Flesca, T. Schmitt, S. Ducruix, S. Candel, Large eddy simulations of a transcritical round jet submitted to transverse acoustic modulation, Phys. Fluids 28 (2016) 055106.

[34] Y. Méry, L. Hakim, P. Scouflaire, L. Vingert, S. Ducruix, S. Candel, Dynamics of a transcritical coaxial flame under a high-frequency transverse acoustic forcing: Influence of the modulation frequency on the flame response, Combust. Flame 162 (2015) 3482-3502.

[35] T. Schmitt, Y. Méry, M. Bileau, S. Candel, Large-Eddy simulation of oxygen/methane flames under transcritical conditions, Proc. Comb. Inst. 33 (2011) 1383-1390.

[36] L. Hakim, A. Ruiz, T. Schmitt, M. Boileau, G. Staffelbach, S. Ducruix, B. Cuenot, S. Candel, Large eddy simulations of multiple transcritical coaxial flames submitted to a high-frequency transverse acoustic modulation, Proc. Comb. Inst. 35 (2015) 1461-1468.

[37] A. Urbano, L. Selle, G. Staffelbach, B. Cuenot, T. Schmitt, S. Ducruix, S. Candel, Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine, Combust. Flame 169 (2016) 129-140.

[38] M. Juniper, S. Candel, The stability of ducted compound flows and sequences for the geometry of coaxial injectors, J. Fluid Mech. 482 (2003).

[39] M. Juniper, The effect of confinement on the stability of non-swirling round jet/wake flows, J. Fluid Mech. 605 (2008) 227-252.

[40] P. Huerre, P.A. Monkewitz, Open shear flow instabilities, Cambridge University Press, 2000.

[41] M. Juniper, Absolute and convective instability in gas turbine fuel injectors, in: Proceedings of the ASME Turbo Expo, Copenhagen, Denmark, June 11-15, 2012, (GT2012-68253) 2012.

[42] K. Oberleithner, M. Sieber, C.N. Nayeri, C.O. Paschereit, C. Petz, H.C. Hege, B.R. Noack, I. Wygnanski, Three-dimensional coherent structures in a swirling jet undergoing vortex breakdown: Stability analysis and empirical mode construction, J. Fluid Mech. 679 (2011) 383-414.

[43] F. Marble, J. Broadwell, The Coherent Flame Model for Turbulent Chemical Reactions, Technical Report, (TRW-9-PU) Project Squid, 1977.

[44] T. Butler, P. O'Rourke, A numerical method for two-dimensional unsteady reacting flows, Proc. Comb. Inst. 16 (1977) 1503-1515.

[45] P. O'Rourke, F. Bracco, Two scaling transformations for the numerical computation of multidimensional unsteady laminar flames, J. Comput. Phys. 33 (2) (1979) 185-203.

[46] Project squid: Introductory report, 1946, https://apps.dtic.mil/sti/tr/pdf/ADA952695.pdf.

[47] Project squid: Field survey report, 1947, https://apps.dtic.mil/sti/tr/pdf/ADA952998.pdf.

[48] G. Brown, A. Roshko, On density effects and large structure in turbulent mixing layers, J. Fluid Mech. 64 (4) (1974) 775-816.

[49] D. Spalding, Mixing and chemical reaction in steady confined turbulent flames, in: Symp. (Int.) Combust., vol. 13, 1971, pp. 649-657.

[50] D.B. Spalding, Development of the Eddy-break-up model of turbulent combustion, in: Symp. (Int.) Combust., vol. 16, The Combustion Institute, 1976, pp. 1657-1663.

[51] B. Magnussen, B. Hjertager, On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion, in: Symp. (Int.) Combust., vol. 16, 1976, pp. 719-727.

[52] G. Batchelor, The effect of homogeneous turbulence on material lines and surfaces, Proc. Royal Soc. London. Series A, Math. and Phys. Sci. 213 (1114) (1952) 349-366.

[53] P. Saffman, A model for inhomogeneous turbulent flow, Proc. Royal Soc. London. Series A, Math. and Phys. Sci. 317 (1530) (1970) 417-433.

[54] T. Poinsot, D. Veynante, Theoretical and Numerical Combustion, fourth ed., 2022, https://www.amazon.fr/Theoretical-Numerical-Combustion-Thierry-POINSOT/dp/2746639904, ISBN 978-2746639904.

[55] S. Candel, N. Darabiha, E. Esposito, Models for a turbulent, premixed dump combustor, in: AIAA/SAE/ASME 18th Joint Propulsion Conference, 1982, AIAA Paper 82-1261.

[56] N. Darabiha, V. Giovangigli, A. Trouvé, S. Candel, E. Esposito, Coherent flame description of turbulent premixed ducted flame, in: France - USA Workshop on Turbulent Combustion, Rouen, 1987.

[57] N. Darabiha, V. Giovangigli, A. Trouvé, S. Candel, E. Esposito, Flamelet modeling of turbulent premixed flames, in: AGARD (Ed.), Combustion and Fuels in Gas Turbine Engines, in: Conference proceedings No 422, 1987.

[58] J. Duclos, D. Veynante, T. Poinsot, A comparison of flamelet models for premixed turbulent combustion, Combust. Flame 95 (1/2) (1993) 101-118.

[59] D. Veynante, S. Candel, J. Martin, Coherent flame modeling of chemical reactions in a turbulent mixing layer, in: 2nd Workshop on Modelling of Chemical Reaction System, Heidelberg, Springer-Verlag, 1986, pp. 386-398.

[60] F. Lacas, S. Zikikout, S. Candel, Comparison between calculated and experimental mean source terms in non-premixed turbulent combustion, in: AIAA/SAE/ASME/ASEE 23rd Joint Propulsion Conference, 1987, pp. AIAA-87-1782.

[61] D. Veynante, F. Lacas, S. Candel, A new flamelet combustion model combining premixed and non-premixed turbulent flames, in: 27th Aerospace Sciences Meeting, AIAA 89-0487, 1989.

[62] D. Veynante, F. Lacas, S. Candel, Numerical simulation of the transient ignition regime of a turbulent diffusion flame, AIAA J. 29 (5) (1991) 848-851.