Emilien Courtine

Laurent Selle

Thierry Poinsot

DNS of Intrinsic ThermoAcoustic modes in laminar premixed flames

Keywords:

come

Introduction

Thermoacoustic instabilities have been a topic of strong interest in the aerospace industry as well as many other engineering applications for the past decades. Rayleigh [START_REF] Rayleigh | The Theory of Sound, Mac Millan[END_REF] first described the underlying mechanism of thermoacoustic instabilities as a constructive build up of acoustic energy by the product of acoustic pressure and unsteady heat released by the flame. Such accumulation is usually found in systems that can store acoustic energy, thereby involving acoustic eigen modes. Resonance has become a mainstay of combustion instability analysis, leading to the idea that there is no such thing as a thermo-acoustically unstable flame in the absence of acoustic modes in the combustor, and that there has to be some sort of acoustic coupling with the surrounding system for an instability to develop.

However, recent theoretical studies [START_REF] Hoeijmakers | [END_REF]3] and experiments [START_REF] Hoeijmakers | Flame-acoustic coupling in combustion instabilities[END_REF] suggest that, even in an anechoic environment, thermoacoustic instabilities may exist. These modes are called here Intrinsic Thermo-Acoustic (ITA) modes because they do not require any acoustic feedback from the boundaries of the burner to exist 1 . They correspond to a feedback mechanism inherent to the flame and its anchoring station, completely independent of the acoustic behavior of the surrounding systems (combustion chamber, injectors, nozzle, compressors, turbines, * Corresponding author.

E-mail address: emilien.courtine@imft.fr (E. Courtine). 1 The intrinsic terminology used here does not refer to flame front instabilities leading to cellular flames such as the Darrieus-Landau instability or thermodiffusive effects [START_REF] Williams | Combustion Theory[END_REF][START_REF] Matalon | [END_REF], which do not require acoustics to develop. It is neither related to parametric instabilities [7] which apply to planar flames. etc.). Even though ITA modes have been evidenced only recently by studying anechoic combustors, they may play a major role in most combustors where they can interact with other feedback mechanisms encountered in standard thermoacoustic modes.

ITA modes have been predicted theoretically and observed experimentally [START_REF] Hoeijmakers | [END_REF][3][START_REF] Hoeijmakers | Flame-acoustic coupling in combustion instabilities[END_REF]. Analyzing them using Direct Numerical Simulation (DNS) is an obvious next step [START_REF] Courtine | Proceedings of the Summer Program[END_REF]. This is the subject of this paper. In Section 2, the theoretical derivation of ITA modes is recalled, based on the work of Hoeijmakers et al. [START_REF] Hoeijmakers | [END_REF] and Emmert et al. [3], and further emphasis is given on the key parameters controlling this instability. Section 3 presents the numerical strategy to capture intrinsic thermoacoustic instabilities, and details are given on the laminar premixed flame setup used in this study. Since theory suggests that ITA modes are controlled by the cross-section ratio S 2 /S 1 between chamber and injection duct (Fig. 1), the Flame Transfer Function (FTF) is computed for four cross-section ratios ranging from 1.5 to 6, and used in the theoretical model for stability and frequency predictions in Section 4. The corresponding DNS are performed in Section 5 and compared to theory in terms of stability and frequency. Good agreement is found with theory, but practical limitations for real configurations are also discussed. Section 6 focuses on the particular case S 2 /S 1 = 2, which is found to be unstable and used for comparison with theory in terms of frequency and mode structure. The acoustic properties of the burner play no role on ITA modes. This is confirmed by performing the simulation of the same flame in two different burner configurations where the lengths of the injection duct and combustion chamber are changed. Finally, in Section 7, the results from DNS are used to capture and evidence the underlying physical phenomena responsible for intrinsic thermoacoustic instabilities.

Theoretical background

The generic setup studied here corresponds to the classical problem of a laminar premixed flame stabilized at a sudden expansion plane separating the injection duct (length l 1 , section S 1 ) and the combustion chamber (length l 2 , section S 2 ). The configuration is twodimensional (dihedral flame). The acoustic properties of the inlet and outlet boundaries are represented by reflection coefficients R j , where the index j ∈ {1, 2} respectively refers to the inlet and outlet ducts.

The linearized Euler equations lead to a simple 1D acoustic wave equation in each duct. The solutions of this equation are forward and backward propagating planar waves A + and A -(see Fig. 1) that couple the acoustic pressure p ′ and velocity u ′ fields as follows (assuming harmonic fluctuations):

p ′ j (x, t) = ℜ( p j (x)e -iωt ) with p j (x) = A + j e ik j x + A - j e -ik j x (1) u ′ j (x, t) = ℜ( û j (x)e -iωt ) with û j (x) = 1 ρ j c j A + j e ik j x -A - j e -ik j x ( 2
)
where ω is the considered angular frequency, k j = ω/c j the associated wave number, and ρ j and c j the density and sound speed that change between cold and hot gases. The prime ( ′ ) and hat (^) symbols respectively denote temporal harmonic fluctuations and the associated complex amplitude at the corresponding angular frequency2 . In the wave formulations of Eqs. ( 1) and ( 2) the Mach number is assumed to be zero so that acoustics is not affected by the mean flow. The amplitudes of these planar waves in the injection and combustion ducts are coupled via jump relations to account for the crosssection area change and the combustion source term at x = 0, using the assumption that the flame is compact with regards to the acoustic wave lengths. These are known as the acoustical Rankine-Hugoniot jump relations, expressed here in the limit of zero Mach number [START_REF] Chu | [END_REF][10][11][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Bauerheim | [END_REF]:

p ′ 2 (x = 0 + , t) -p ′ 1 (x = 0 -, t) = 0 (3) S 2 u ′ 2 (x = 0 + , t) -S 1 u ′ 1 (x = 0 -, t) = γ -1 γ p 0 ˙ ′ T (t) (4)
where γ is the specific heat capacity ratio, p 0 is the reference pressure, assumed to be constant across the flame, and ˙ ′ T is the fluctuating component of heat release rate integrated over the flame domain. In the velocity jump relation (Eq. ( 4)), the acoustic emission is due to dilatation induced by the unsteady reaction rate. This relation is expressed in terms of volumetric flow rate, as effects of mean flow on the acoustics are ignored3 [10,[START_REF] Bauerheim | [END_REF].

To close the set of Eqs. ( 3) and ( 4), ˙ ′

T must be expressed as a function of the acoustic field. The simplest model for the unsteady reaction rate is to link ˙ ′ T with the upstream acoustic velocity u ′ 1 [START_REF] Crocco | Proceedings of the 10th International Symposium on Combustion[END_REF].

This model combines a time delay τ and a dimensionless interaction index n:

˙ ′ T (t) = ρ 1 c 2 1 γ -1 S 1 nu ′ 1 (t -τ ) or ˆ T = ρ 1 c 2 1 γ -1 S 1 n û1 e iωτ (5)
The term of velocity sensitive has been coined for the description of flames whose behavior matches Eq. [START_REF] Williams | Combustion Theory[END_REF]. Experimental data suggest that the heat release of non-planar premixed laminar or turbulent flames is indeed mostly driven by velocity fluctuations rather than pressure fluctuations [START_REF] Langhorne | [END_REF]16]. The parameters n and τ contain all the convective mechanisms controlling the flame response (e.g. vortex formation caused by the acoustic velocity surge, vortex convection by the mean flow followed by vortex breakdown and combustion in turbulent flows [17]). Thus, the model of Crocco (Eq. ( 5)) links the mechanisms in the "convective world" to those in the "acoustic world", u ′ 1 being supposed to contain only acoustical components [START_REF] Crocco | Proceedings of the 10th International Symposium on Combustion[END_REF]. One may note the peculiar property of this model as it seems to be acoustically non causal 4 : the input quantity being the reference velocity fluctuations in the inlet duct u ′

1 = (A + 1 -A - 1 )/(ρ 1 c 1 )
, it is the difference of upstream and downstream propagating waves. At first sight it seems implausible that an upstream oriented wave A - 1 , propagating away from the flame front, can affect the combustion process. This causality concern has been addressed by [START_REF] Courtine | Proceedings of the Summer Program[END_REF] who showed that Crocco's model holds for a laminar premixed flame in this configuration, proving that this flame is indeed velocity sensitive. In fact, A - 1 in the inlet duct is the trace of an acoustic wave, generated by the flame or coming from the outlet of the combustor, propagating upstream and that produced a convective wave when passing through the backward step, leading to mode conversion as discussed more thoroughly in Section 7.

Using Crocco's model (Eq. ( 5)) and the wave decomposition (Eqs. ( 1) and ( 2)), the acoustic jump relations (Eqs. ( 3) and ( 4)) become:

A + 2 + A - 2 = A + 1 + A - 1 (6) A + 2 -A - 2 = (A + 1 -A - 1 )(1 + ne iωτ ) (7)
where = (S 1 ρ 2 c 2 )/(S 2 ρ 1 c 1 ) is a coupling parameter that reduces to = S 1 T 1 /(S 2 T 2 ) for isobaric flames using the perfect gas law [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Palies | Dynamique et instabilités de combustion de flammes swirlées[END_REF][START_REF] Schuller | [END_REF], T 1 and T 2 being respectively the fresh and burnt gases temperatures.

Eqs. ( 6) and ( 7) are the basic elements of all network models for thermoacoustics. ITA modes can be constructed by considering the specific case where both terminations are anechoic: R 1 = R 2 = 0. In this case, the amplitudes of the waves propagating towards the flame vanish:

A + 1 = R 1 A - 1 = 0 and A - 2 = R 2 A + 2 = 0
, and the jump relations (6) and (7) reduce to:

A + 2 = A - 1 and A + 2 [1 + (1 + ne iωτ )] = 0
which leads to the dispersion equation derived by [START_REF] Hoeijmakers | [END_REF] and [3]:

1 + (1 + ne iωτ ) = 0 (8)
The solutions of Eq. ( 8) correspond to the ITA modes and have analytical expressions for the real part ω r (frequency) and the imaginary part ω i (growth rate):

ω q r = (2q -1)π τ (q ∈ N * ) (9)
ω i = 1 τ ln n 1 + (10)
The angular frequency of the first mode (Eq. ( 9) with q = 1) simply writes:

ω 1 r = π τ so that T = 2τ (11)
where T = 2π /ω 1 r is the mode period. As noted by Emmert et al. [3] this mode is quite unusual in thermoacoustics since it is obtained in a situation where acoustic losses at the burner's ends are maximum (anechoic terminations). The period of the first ITA mode is twice the flame delay and is not linked to any acoustic mode of the duct. Note that no acoustic mode can exist here because both terminations are anechoic: therefore the only time scale in this problem is the flame delay τ .

Another interesting characteristic of ITA modes is that their stability, i.e. the sign of their growth rate, does not depend on the flame delay τ but only on the gain n. In classical thermoacoustic modes, stability is primarily determined by the delay τ , n playing only a secondary role: in the limit where small interaction indices are considered (weak flames) and acoustic losses are not accounted for, the sign of the growth rate is even independent of n [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. For ITA modes however, a critical gain n c appears when writing the stability criterion:

ω i > 0 ⇐⇒ n > n c = + 1 (12)
This critical value may also be found using Rayleigh's criterion on the compact flame element, as shown in Appendix A.

Crocco's model with constant parameters n and τ (Eq. ( 5)) constitutes a rather limited representation of real flames. Their dynamics can be modeled more accurately by using frequency dependent parameters n(ω) and τ (ω) or, equivalently, normalized Flame Transfer Functions, F(ω), that have been introduced for this purpose and are now used by most authors:

ˆ T ¯ T = F (ω) û1 ū1 (13)
where ¯ T is the total, time-averaged heat release rate. By taking ū1 as the bulk velocity of the inlet flow, one can see that the n -τ model is equivalent to the FTF, using the relation between steady heat release rate and bulk velocity ¯ T = ρ 1 S 1 ū1 c p (T 2 -T 1 ), where c p is the specific heat capacity at constant pressure of the fresh gas mixture 5 :

θ F = ne iωτ where θ = T 2 T 1 -1 (14)
Eq. [START_REF] Crocco | Proceedings of the 10th International Symposium on Combustion[END_REF] shows that F(ω) carries the same information as the parameters n(ω) and τ (ω). In the general case, n and τ being functions of the frequency, Eq. ( 8) can be recast in terms of FTF as:

F (ω) = - + 1 θ ( 15
)
5 cp is assumed to be constant through the flame. or, equivalently:

φ = arg (F(ω)) = (2q -1)π( q ∈ N * ) (16) |F(ω)| = + 1 θ (17)
Eqs. ( 16) and ( 17) are implicit equations that require further knowledge on the transfer function of the flame, i.e. the FTF, to solve for ω r and ω i . A useful particular case corresponds to flames for which τ is supposed to be independent of ω (an assumption which is reasonable for many premixed laminar flames at low frequencies [20]).

In this case, φ = ω r τ so that Eq. ( 16) reduces to the simple expressions of Eqs. ( 9) and (11). In other cases, Eqs. ( 16) and ( 17) require a numerical resolution, which in turn requires the Flame Transfer Function F(ω) to be known on the entire complex plane (i.e. at ω i ̸ = 0) [21].

However, in practical situations, only the frequency response F(ω r ) is known. This frequency response, although not sufficient to obtain the exact solutions of Eq. ( 15), is sufficient to derive a stability criterion for ITA modes, according to Nyquist's criterion 6 . Henceforth, F will refer to the frequency response of the flame and will indistinctly also be called FTF. The stability criterion is twofold: a necessary condition on the gain, and a criterion for the phase.

A necessary condition for the instability to start is that the modulus of the frequency response |F(ω r )| exceeds a critical threshold F c , as can be inferred from Eqs. ( 12) and [START_REF] Crocco | Proceedings of the 10th International Symposium on Combustion[END_REF].

|F| > F c with F c = + 1 θ (18)
For a perfect gas and an isobaric flame, = S 1 T 1 / (S 2 T 2 ) so that the threshold F c can be written:

F c = 1 T 2 /T 1 -1 1 + S 2 S 1 T 2 T 1 (19)
This is only a necessary condition for an unstable ITA mode to appear: for real flames, the gain of the FTF, |F(ω r )|, is a function of the frequency. Therefore the criterion on gain given in Eq. ( 18) may be satisfied only in a restricted frequency range, within which the phase criterion arg (F(ω)) = π (Eq. ( 16)) must also be satisfied to trigger an unstable ITA mode 7 . Eq. [START_REF] Schuller | [END_REF] shows that the ITA stability threshold goes down linearly when the section ratio between inlet duct and combustion chamber S 2 /S 1 decreases, or when the temperature ratio T 2 /T 1 (or θ ) increases, as shown in Fig. 2. Intense flames in chambers with small section changes (strong confinement) should therefore be more prone to intrinsic thermoacoustic instabilities.

Figure 2 displays a stability map of ITA modes for one particular parameter: confinement S 2 /S 1 . The stable and unstable regions are separated by the line |F| = F c (S 2 /S 1 ). A similar map can be drawn for the temperature jump across the flame, the stability boundary becoming, in that case, the curve |F| = F c (θ) .

The final part of this section on theory is devoted to the derivation of the spatial structure of ITA modes. Using the definition of the FTF (Eq. ( 13)), and recalling that arg (F) = π [2π ] for ITA modes so that F = -|F |, it is found that ˆ T and û1 are completely out of phase. Moreover, using the null pressure jump condition p1 (x = 0) = p2 (x = 0) = p (Eq. ( 3)) and recalling that anechoic boundaries are considered, it can be shown that the velocity fluctuations at the reference point û1 and the acoustic pressure at the flame location p are out of phase too, meaning that p and ˆ T are in phase. This is consistent with Rayleigh's criterion that states that production of acoustic energy is only possible when pressure and heat release oscillations are in phase (Appendix A). The relationships between û1 , p and ˆ T for an ITA mode are the following:

ˆ T = -|F | ¯ T ū1 û1 and û1 = - A - 1 ρ 1 c 1 = - p ρ 1 c 1 ( 20
)
With the compact flame assumption, and anechoic terminations, the acoustic velocity and pressure in each duct are proportional to a single planar harmonic wave, e.g. A - 1 , thereby their phase k j |x| = ω|x|/c j unwraps linearly with the distance to the flame |x| (Eqs. ( 1) and ( 2)). The pressure jump condition (Eq. ( 3)) shows continuity of the pressure phase and amplitude across the flame. Regarding acoustic velocity, the jump condition (Eq. ( 7)) combined with the definition of the FTF (Eq. ( 13)), together with the fact that arg (F) = π [2π ] for ITA modes, imply a phase difference of π across the flame and an amplitude ratio S 1 /S 2 (θ |F| -1) between hot and cold gases 8 . Within each duct acoustic losses are not taken into account so that the amplitudes of velocity and pressure fluctuations remain constant with the distance to the flame. This spatial structure of the first ITA mode is summarized in Eq. ( 21), and shown in Fig. 3.

⎧ ⎪ ⎨ ⎪ ⎩ | p2 | | p1 | = 1 | û2 | | û1 | = S 1 S 2 (θ |F| -1) ⎧ ⎪ ⎨ ⎪ ⎩ arg[ p1 ] = - π c 1 τ x arg[ û1 ] = - π c 1 τ x ⎧ ⎪ ⎨ ⎪ ⎩ arg[ p2 ] = π c 2 τ x arg[ û2 ] = π c 2 τ x + π (21)
8 The π phase shift for the velocity fluctuations between hot and cold gases comes from the fact that 1 + θF = 1 -θ|F| with θ|F| -1 > 1/ > 0 for an unstable ITA mode (Eq. ( 18)). This property is also true for stable ITA modes as long as |F| > 1/θ , otherwise the phase of u ′ is continuous through the flame.

Numerical setup

The DNS configuration corresponds to a laminar 2D dihedral flame previously studied in [START_REF] Mejia | [END_REF]24], on which a combustion chamber has been added for the study of intrinsic thermoacoustic modes (Fig. 1). The cross-section of the inlet duct S 1 = 1 cm is fixed, while the crosssection of the combustion chamber can be adjusted in the DNS. The lengths of the inlet duct l 1 = 5 cm and combustion chamber l 2 = 10 cm are fixed, except in one case where they are doubled to ensure that they have no impact on the mode. An homogeneous methane-air mixture at equivalence ratio φ = 0.95 and temperature T 1 = 300 K is fed through the inlet duct at ū1 = 1.8 m/s, leading to a Reynolds number Re 1 ≃ 2000, typical of a laminar flow. The adiabatic temperature jump induced by the flame at this equivalence ratio is θ adiab (φ = 0.95) = T 2 /T 1 -1 = 6.3, and the laminar flame speed S L = 0.41 m/s will be constant throughout the study. The resulting flame tip halfangle of the flame is α = sin -1 (S L / ū) = 13.2 • , as shown in Fig. 1.

The numerical solver used for the computation of the Flame Transfer Functions and the unsteady simulations is AVBP. A highorder fully explicit scheme is used to advance the compressible reacting Navier-Stokes equations [25]. Chemistry is modeled using a two-step chemical scheme for methane-air flames, yielding a laminar flame speed S L = 0.41 m/s, and the thickened flame model [26] with a flame thickening factor of 2 for computational cost issues. Numerical meshes are composed of approximately 500 000 cells (depending on the geometry considered) all with identical smallest cell volume V cell min ≃ 5.5 × 10 -10 m 2 (dimensioned to resolve temperature gradients at the flame anchor station), leading to time steps 0.015 dt 0.017 µs 9 , with a CFL condition of 0.7 and a Fourier number F o = 0.1. It was checked on an unstable case that using finer grids, with different flame thickness factors, had very limited to no impact on the results both in terms of flame shape and quantitative results for the mode. Acoustic boundary conditions at the inlet and outlet of the combustor are imposed using the NSCBC method [27,28]. For all present simulations, fully non-reflecting conditions are used to allow intrinsic thermoacoustic modes to develop. Walls are modeled with no-slip isothermal boundary conditions at 300 K.

In order to capture an unstable ITA mode numerically, first the FTF of the flame has to be computed so that it can be used in the theoretical stability criterion (Eq. ( 18)) to evaluate the range of the confinement parameter (S 2 /S 1 ) required for triggering an unstable ITA mode 10 , θ being fixed by the chemistry and the equivalence ratio that are constant throughout the study. The FTF will also serve to predict, based on theory, the frequency of the instability (Eq. ( 16)), or 9 A typical physical time of interest of 100 ms requires about 65 h on 192 Intel Harpertown 3.00 GHz cores (totaling 12 500 CPU hours) on the SGI Altix Ice 8200 scalar supercomputer Jade of CINES. 10 Throughout this study, triggering will refer to the onset of a linear instability starting from a steady flame configuration, and not to high-amplitude external excitations required to set off non-linear instabilities as can be found in studies dedicated to nonlinear systems.

- ) and experimental results for the unconfined flame ( ). From left to right: modulus and argument of F against forcing frequency. The reference point P is located 1 cm upstream of the combustion chamber in the inlet duct.

l 1 0 l 2 0 1 2 3 S1 S2 (θ|F | -1) x G -l 1 0 l 2 0 π 50 π c2τ π c1τ x arg(p) [rad] 0 π arg(û) [rad]

Table 1

Summary of the various configurations used in numerical simulations. Corresponding case names refer to standard duct length l j (REF) or enlarged duct length l ′ j = 2l j (DOUBLE), and confinement expressed in terms of cross-section ratio S 2 /S 1 . Then a set of simulations is performed for different cross-section ratios (S 2 /S 1 ) (see Table 1) and intrinsic thermoacoustic stability will be assessed using temporal signals of heat release rate, pressure and velocity fluctuations. These are compared with theoretical predictions in terms of stability (Fig. 2), frequency (Eq. ( 11)), and mode structure (Fig. 3).

Configuration name Case

REF

Measurement of the FTF in an intrinsically unstable system.

The gain of Flame Transfer Functions of premixed laminar flames is typically of order unity at low frequencies and goes down to zero at high frequencies (low-pass filter behavior). It was shown analytically from the linear solution of the G-equation [29] that |F| stays below unity at every frequency for conical flames, and may reach values up to |F| = 2 for very shallow V-flames. The case of 2D dihedral flames, such as the one investigated in this study, is not discussed in [29] but is expected to behave qualitatively more like an axisymmetric conical flame. Therefore, even though slight discrepancies may be found between experimental or numerical results and the linearized G-equation theory, unstable ITA modes should be found for values of S 2 /S 1 such that the critical FTF gain is lower than 1, i.e. F c < 1. This corresponds to S 2 /S 1 < 2, according to Fig. 2 and using θ adiab = 6.3.

Although being just a rough estimate, this order of magnitude of critical cross-section ratio S c 2 /S 1 ≃ 2 sets off the procedure for the entire quantitative analysis of ITA modes.

First, the FTF was computed at very low confinement (S 2 /S 1 = 6) and compared to the experimental measurements performed by [START_REF] Mejia | [END_REF] on the same flame but in an unconfined configuration. For S 2 /S 1 = 6 (case REF-6.0 on Table 1), the critical threshold for ITA modes is F c = 2.72 (Eq. ( 19)) while FTF values never exceed 1.2. Therefore, as suggested in the previous paragraphs, the flame is stable for this case: no ITA mode is observed. The FTF measured by DNS (Fig. 4) shows reasonable agreement with the unconfined experimental FTF. The computation of the FTF was performed by enforcing harmonic acoustic waves through the inlet boundary condition using the NSCBC formalism [27]. Results do not depend on the origin of the acoustic excitation (upstream or downstream) as long as the reference point is located inside the inlet duct [START_REF] Courtine | Proceedings of the Summer Program[END_REF].

The comparison of Fig. 4 between numerical and experimental FTF results may seem unsuitable considering they were performed with different confinements (S 2 /S 1 = 6 for DNS, and S 2 /S 1 = ∞ for the experiment). This is a small limitation as shown by recent work on the effect of confinement on laminar conical flames [30]. When the flame is unconfined, the plume of hot gases expands in the transversal direction only up to a given surface value, noted S b and called the unconfined plume surface [31]. Confinement effects may be quantified using the dimensionless parameter C b = S b /S 2 (where S 2 is the surface of the hot gas jet in a confined burner) that can be related to the geometry of the burner and thermodynamic quantities:

C b = S b S 2 = C r 1 - θ θ + 1 cos α -1 2 ( 22
)
where C r = S 2 /S 1 and α = sin -1 (S L / ū) is the flame tip half-angle.

For the FTF shown in Fig. 4, C b (S 2 /S 1 = 6) = 1.02 ≃ 1 so that confinement is expected to have negligible impact on the FTF in this case. However, when S 2 /S 1 decreases, the FTF is expected to change. This was accounted for by additional FTF computations, which can become increasingly difficult to perform when S 2 /S 1 reaches values of the order of 2, i.e. when ITA modes become unstable and can spontaneously be excited (cases REF-2.0 and REF-1.5). Contrary to typical thermoacoustic studies where the computation of the FTF at the frequency of the instability is usually possible by isolating the flame in an anechoic combustor (therefore removing the coupling with the acoustics of the surrounding elements and making the flame stable), unstable ITA modes can appear even for cases where R 1 = R 2 = 0. However, Hoeijmakers et al. [START_REF] Hoeijmakers | [END_REF] suggest that ITA modes may be damped by enforcing partially reflecting boundary conditions with a well-chosen set of reflection coefficients and duct lengths. In other words, unstable ITA modes can be stabilized by using reflecting boundary conditions if R 1 and R 2 are properly chosen. This procedure allows the computation of the FTF for virtually any confinement, even those that would trigger unstable ITA modes in anechoic environments.

0 150 0.5 1 F c f 1 IT A ∆f f [Hz] |F | 0 150 0 π 2π f 1 IT A f [Hz] arg [ F ]
The flame was stabilized using this technique for cases REF-1.5 and REF-2.0 (the less confined cases REF-3.0, REF-4.0, and REF-6.0 being intrinsically stable). The appropriate reflection coefficients for stabilization R j were obtained from the relaxation coefficients K j enforced through Navier-Stokes Characteristics Boundary Conditions at the inlet and outlet 11 [32]. By choosing K 1 = K 2 = 800 s -1 , the associated complex-valued frequency-dependent reflection coefficients R j = ±1/(1 -2iω/K j ), that have phases arg (R j ) ≃ π /2 and modulus values |R j | ranging from 0.1 and 0.2 over the frequency range of interest (50 ࣠ f ࣠ 100 Hz), are sufficient to stabilize the flame for cases REF-1.5 and REF-2.0 12 . FTFs were then computed using acoustic forcing from the inlet boundary condition, following the same procedure as for case REF-6.0 (Fig. 4). It was also checked that the reflection coefficients R j , that behave like low-pass filters of cut-off frequency f c = K/(4π ) = 32 Hz [32], let the acoustic fluxes of the harmonic forcing leave and not accumulate in the cavity. This method yielded constant amplitude harmonic signals suitable for the computations of the FTFs, represented in Fig. 5.

Figure 5 shows that, as S 2 /S 1 decreases, the slope of the FTF phase is reduced, while the overall gain curve is also reduced. The flame delay τ = arg (F)/ω is reduced because of the deformation of the unperturbed flame front and the increased velocity of the fresh stream of reactants along the central axis due to confinement. While theory predicts a shift of the gain peak towards higher frequencies [30], DNS 11 At the inlet boundary the NSCBC formalism enforces a characteristic wave entering the domain L 5 = K 1 ρ 1 c 1 (uu T ), whereas a characteristic wave L 5 = K(pp T ) is imposed at the outlet boundary condition, the "T" subscript referring to target values. 12 Complex-valued reflection coefficients are obtained via this procedure, whereas only real-valued reflection coefficients were considered in the parametric study of [START_REF] Hoeijmakers | [END_REF]. Nevertheless, similar results hold with complex-valued reflection coefficients in terms of stabilization of ITA modes. results only show a gain reduction with increasing confinement. The FTFs presented in Fig. 5 are computed only on the frequency range relevant for the study of ITA modes, i.e. close to arg (F( f )) = π , as explained hereafter.

Once the FTF has been computed at a given confinement, it is straightforward to obtain the frequency and the stability criterion of any ITA mode. The critical gain F c is given by Eq. ( 18) as a function of confinement S 2 /S 1 (see Fig. 2). The theoretical frequency of the first ITA mode corresponds to arg (F( f 1 ITA )) = π (Eq. ( 16)), as shown in Fig. 6. Let f be the range of frequencies for which |F(f)| ≥ F c . If f 1 ITA ∈ f, as it is the case on the example in Fig. 6, the first mode should be unstable. The corresponding flame delay given by the FTF at this frequency is τ = 1/(2 f 1 ITA ). The frequencies of higher order ITA modes satisfy arg (F( f q ITA )) = π + 2(q -1)π with q > 1 so that they are generally out of the frequency range where the FTF is above the critical gain. Thus, on the example of Fig. 6, only the first ITA mode should be unstable. This method of prediction based on the FTF and the critical gain is equivalent to the Nyquist criterion applied to a Bode diagram.

Stability limits

The stability of an ITA mode can be predicted by using the procedure of Fig. 6 on the FTFs of Fig. 5. The corresponding theoretical stability map is represented in Fig. 2 as a function of confinement. The purpose of this section is to check whether the DNS results confirm this stability map in terms of critical gain for various confinements. To do so, DNS are repeated for several combustion chambers of different cross-sectional area ratios S 2 /S 1 = 1.5, 2.0, 2.5, 3.0, 4.0 and 6.0 with fixed lengths l 1 and l 2 (see Table 1). The procedure described in Section 6 is used: first the flame is stabilized using appropriate reflection coefficients at the inlet and outlet boundaries, then the acoustic feedback is removed and the flame slightly perturbed from equilibrium by a low amplitude acoustic impulse. Stability is assessed by considering the ensuing growth or decay of oscillations.

The results of these DNS are in qualitative agreement with theoretical predictions, as summarized in Fig. 7. Stable simulations are represented with squares, while unstable simulations are marked by triangles. Their location on the stability map is determined by their respective cross-section ratio S 2 /S 1 and the corresponding FTF gain computed at the frequency of the instability as shown in Fig. 5 (see computation of the FTFs in Section 4). Some cases (S 2 /S 1 = 2.5 and ). Unstable ( ) and stable ( ) results from DNS. Ordinates are FTF gains at the frequency of the first ITA mode. [START_REF] Hoeijmakers | Flame-acoustic coupling in combustion instabilities[END_REF].0) are displayed in Fig. 7 even though their FTFs were not computed. Their FTF gains are interpolated from neighboring points.

The trend is correctly reproduced: for confinements at which the FTF gain exceeds the critical gain F c , i.e. for S 2 /S 1 ࣠ 2, ITA modes are found unstable in the DNS. Conversely, for S 2 /S 1 ࣡ 3, the gain of the FTF is lower than the critical value and the corresponding simulations are stable. Some DNS, however, do not comply with theory. While the case S 2 /S 1 = 2.0 is predicted to be marginally stable by theory, it is found to be unstable in the DNS (see details in Section 6). Similarly, case S 2 /S 1 = 2.5 is predicted to be stable by theory and yet is unstable in the DNS. Possible explanations for these discrepancies are twofold, and directly linked to each constitutive block of the ITA theory. The first one is the determination of the critical FTF gain F c (Eq. ( 19)), based on the Rankine-Hugoniot acoustic jump relations (Eqs. ( 3) and ( 4)). While the classical assumption of acoustical compactness of the flame is reasonable at the frequencies of interest (the associated wavelengths, λ 3.5 m, being much larger than the average flame length H = 2 cm), it may be invalidated for the intrinsic thermoacoustic feedback loop which is entirely comprized within the flame scale, where the acoustic velocity field is highly discontinuous (see Section 7). This acoustic field is also clearly two-dimensional (see Fig. 15), questioning the 1D assumption used in the acoustic jump relations. Moreover, the normalized adiabatic temperature jump across the flame, θ adiab = T 2 /T 1 -1 = 6.3, present in the expression of F c (Eq. ( 19)), should be modified to account for heat losses at the walls (even though reducing θ would increase F c and therefore widen the gap between theory and DNS in Fig. 7). The second aspect of the ITA theory that may explain discrepancies is the computation of the FTF gain. Even though FTF measurements have been performed with upstream and downstream acoustic excitation, giving very consistent results and proving the velocity sensitive assumption [START_REF] Courtine | Proceedings of the Summer Program[END_REF], the FTF gain is sensitive to the location of the reference point used for velocity measurements [33,[START_REF] Mejia | Wall-temperature effects on flame response to acoustic oscillations[END_REF]. Bringing the reference point closer to the dump plane would increase the FTF gain, thereby reducing the gap between theoretical predictions and DNS results.

Despite discrepancies in terms of stability limits, Fig. 7 confirms that ITA modes are expected to be amplified only in chambers with small cross-section ratios (strong confinement). In terms of frequency of the instability, it is interesting to compare the theoretical prediction of the frequency of ITA modes based on FTFs (Eq. ( 16)), and the frequency of the flame oscillations actually observed in DNS (Fig. 8). As explained in Section 4, FTFs were computed for four confinements: S 2 /S 1 = 1.5, 2.0, 3.0 and 6.0, so theoretical predictions are available only for these cases. For the unstable anechoic simulations, the frequency of the ITA mode was evaluated in the early cycles of instability in the DNS, i.e. at low oscillation amplitudes. By doing so, the effects of non-linearities are discarded (see Section 6). A very good agreement is obtained between theoretical predictions and the frequency of ITA modes in the DNS, for both stable and unstable configurations. The results presented in Fig. 8 also show that it is necessary to con- ). Unstable ( ) and stable ( ) results from DNS. Ordinates are oscillation frequencies in the linear regime. sider the influence of confinement on the phase of the FTF in order to have a correct prediction of the instability.

Control, growth and spatial structure of an ITA mode for

S 2 /S 1 = 2
This section focuses on a particular case, S 2 /S 1 = 2. For this case, theory predicts that the frequency of the first ITA mode is f 1 ITA = 71 Hz (Eq. ( 16) and Fig. 5), and that it is marginally stable, using the stability criterion (Eq. ( 18), Figs. 5 and6):

|F( f 1 ITA )| = 1.0 ≃ F c (S 2 /S 1 ) = 1.01.
However, the corresponding anechoic DNS exhibits a strong unstable ITA mode, which is used for more thorough investigation of ITA mode characteristics and comparison with theory.

As mentioned in Section 4, with a well-chosen set of relaxation coefficients K j for the inlet and outlet boundaries of the DNS, it is possible to stabilize an intrinsically unstable flame. Once the stable solution is obtained, relaxation coefficients can be set to zero, corresponding to a perfectly anechoic environment, and a low amplitude Gaussian acoustic impulse can be injected at the inlet boundary condition to trigger the instability 13 . The impulse is a velocity perturbation of maximum amplitude u ′ pulse = ū1 × 10% and characteristic time 1.25 ms τ × 10%. Immediately after the flame is moved away from equilibrium by the acoustic disturbance, it starts pulsating at a distinct frequency as shown by the evolution of heat release rate in Fig. 9 for case REF-2.0. These self-sustained oscillations grow rapidly in time and reach very high amplitudes within a few periods: the fluctuations of global heat release rate are of the same order of magnitude as the mean value after only 5 periods. Following the growth phase, a limit-cycle establishes in the combustion chamber at a very high amplitude. The frequency of the oscillation in the limit-cycle (t ≥ 100 ms), f 1 DNS,l-c = 50 Hz, differs from the one in the phase of linear growth (t ≤ 60 ms), f 1 DNS, lin = 68 Hz. The latter is close to the frequency of the first ITA mode predicted theoretically: f good agreement between theory and DNS in this linear domain (as seen in Fig. 8). At 71 Hz, the time delay τ , measured in the FTF of Fig. 5 for this confinement

S 2 /S 1 = 2, is τ = π /ω r = 1/(2 f 1 ITA ) = 7.0
ms. This is consistent with the theoretical evaluation of Eq. ( 11) that states that the period T = 14.7 ms must be equal to twice the flame delay τ .

In order to verify that the mode observed here is indeed an ITA mode, i.e. that it is intrinsic to the flame and does not depend on the acoustics of the burner, the lengths of the feeding duct l 1 and combustion chamber l 2 were doubled for another simulation: l ′ j = 2l j (case DOUBLE-2.0 on Table 1), and the exact same simulation was conducted on this enlarged setup. The acoustic eigenmodes of REF-2.0 and DOUBLE-2.0 configurations do not overlap, making this test relevant to disprove the acoustic origin of the modes observed in the DNS. Moreover, their frequency is much larger than the frequency range of interest in the DNS (the acoustic mode with lowest frequency in the DOUBLE-2.0 configuration is at 723 Hz). All parameters except the lengths were kept unchanged, as well as the mesh in the flame region, and the onset procedure for triggering the instability was identical. Results (Fig. 10) show that, indeed, the frequency and growth rate of the oscillations were kept unchanged, proving that ITA modes do not depend on the acoustics of the burner.

The DNS mode structure is now compared to the theoretical ITA form of Fig. 3. This is a more complex exercise because the limit-cycle reached in Fig. 9 is strongly non-linear: the frequency shifts from 68 Hz in the linear zone to 50 Hz in the limit-cycle, as can be seen on a PSD in Fig. 11. Although the study of non-linearities is not the purpose of this work, this frequency drop is likely to be due to an effective increase of the mean flame length as its feet are periodically pushed upstream by high amplitude acoustic waves passing through the rim [START_REF] Mejia | [END_REF][START_REF] Kedia | [END_REF]. This would result in an effective increase of the flame delay τ (the mean flow rate being fixed) and a subsequent drop of the instability frequency f 1 ITA = 1/(2 τ ). Figure 11 shows that sub-harmonics of the first ITA mode are present in the frequency content of the limitcycle. They overlap higher order ITA modes (q > 1), which are most likely also triggered by non-linearities, although being linearly stable. More insight on the physical phenomena at play in intrinsic thermoacoustic instabilities are given in Section 7. The reduced frequency of the first ITA mode in the limit-cycle f 1 DNS, l-c = 50 Hz will be used in the remainder of this section for the study of mode structure, which will show that, despite this frequency shift, the mode structure reached in the DNS limit-cycle is very close to the theoretical ITA model of Fig. 3. Limit-cycle oscillations shown in Fig. 9 are reproduced for four periods in Fig. 12. As predicted by theory, and despite non-linear effects, the global heat release rate and pressure signals are perfectly in phase, whereas the velocity signal in the inlet duct is completely out of phase with pressure and heat release, proving that this specificity of ITA modes is well captured by the DNS. The amplitudes of the acoustic waves are also recorded at the inlet and outlet boundaries and plotted in Fig. 13, showing that the characteristic treatment of the anechoic boundaries is effective, even for very high amplitude oscillations: the upstream and downstream incoming waves A + 1 and A - 2 are zero at all times. The mode structure can also be recovered using pressure and velocity signals located along the symmetry plane of the burner (y = 0, see Fig. 1). The comparison between DNS and theoretical predictions 14 (Eq. ( 21) and Fig. 3) is shown in Fig. 14. The π phase shift of the acoustic velocity across the flame is perfectly recovered in the DNS (Fig. 14d), as well as the linear phase unwrapping of the acoustic pressure, which is characteristic of unidirectional planar waves (Fig. 14c). The amplitude of pressure fluctuations remains constant through the flame (Fig. 14a). The flame zone acts as an acoustic monopole radiating upstream and downstream with the same amplitude. The veloc-

ity jump S1/S2 (θ |F( f 1 ITA )| -1
) is correctly predicted despite slight discrepancies (Fig. 14b). The latter discrepancies may be attributed 14 For consistency reasons, the comparison is done by using the actual frequency of the limit-cycle oscillations f 1 DNS, l-c = 50 Hz, leading to the effective flame delay in limit-cycle τDNS, l-c = 10 ms used in Eq. ( 21). to non-linear effects intervening in limit-cycle which result in reducing the gain. A Flame Describing Function F (ω r , |u ′ 1 |) that takes into account the oscillation amplitude would be required for a more accurate estimation of the gain in limit-cycle [36]. As mentioned in Section 5, the validity of the adiabatic velocity jump θ adiab = 6.3 may also be questioned in this non-adiabatic setup. Finally, note that, in Fig. 14, theoretical predictions are not plotted in the flame zone since this region is not described in the model of Eqs. ( 3) and (4), which uses the compact flame assumption.

Description of physical phenomena responsible for ITA modes

The final part of this work is devoted to the analysis of the scenario of the ITA mode revealed by DNS. Indeed, at this point, even though characteristics of these modes, such as frequency, mode structure, or stability limits, match theoretical predictions, the phenomena responsible for this new kind of instability have not yet been described. Previous theoretical studies suggest that vortex shedding may be at the origin of ITA modes [3], while others mention a possible counter-reaction of the flame to its own production of acoustics [START_REF] Hoeijmakers | [END_REF].

This section shows that mode conversion [17,[37][38][39][40], i.e. the transformation of an acoustic wave into a convective wave is the mechanism that must be introduced into the analysis to understand the instability loop. In the present case the transformation takes place at the sharp corner where the flame is anchored, and where vorticity is generated by unsteady shear layers excited by acoustic waves. The second mechanism needed to elucidate the ITA instability is Kinematic Over-Restoration (KOR), which controls the second part of the instability loop. To gain more insight into these mechanisms, a complete ITA cycle is analyzed using the DNS data exposed in Section 6 for the case S 2 /S 1 = 2. Figure 15 shows four snapshots of the flame position during one cycle. The time evolution of the reference acoustic velocity and pressure in the chamber are plotted on the right hand-side of the figure.

The cycle can be described as follows: assume that a pair of vortices (visualized by Q-criterion isocontours in Fig. 15) is created on both sides of the flame at t = t 1 when the acoustic velocity is maximal at the dump plane (as observed in all pulsated jets [17]). These vortices push the flame anchoring points apart and a convective perturbation then propagates along the flame front (instant t 2 ). This perturbation takes approximately a time τ before the flame reaches its maximum elongation and the heat release rate is maximum at time t 3 . Kinematic Over-Restoration subsequently takes place. Standard kinematic restoration happens when a flame, pushed away from its stable position, naturally reverts to steady state [42,43]. Here, not only is the flame far from its equilibrium position at t = t 3 , but the inlet velocity is also minimal at this instant: these two combined factors lead to a very fast contraction of the flame between instants t 3 and t 4 . This is why the fluctuation of reaction rate becomes negative at time t 4 . Strong contraction of the flame results in negative dilatation and negative acoustic pressure, which in turn induces a velocity surge upstream from the flame. This positive inlet acoustic velocity again initiates mode conversion at the corner and creates a new pair of vortices, thereby closing the unstable loop at t = t 1 + 2τ . In this scenario, acoustics plays a role only between instants t 3 and t 4 to trigger a pair of vortices through mode conversion at the corner. No acoustic mode of the full system is involved. The only acoustic propagation involved in the ITA mode takes place between the flame zone and the chamber inlet plane, and is oriented upstream.

Conclusion

The understanding of thermoacoustic instabilities has recently evolved in a significant way through the theoretical derivation of Intrinsic ThermoAcoustic (ITA) instability modes [START_REF] Hoeijmakers | [END_REF]3] that do not depend on the acoustic eigenmodes of the combustor where the flame is located. These ITA modes are studied here using DNS of a laminar premixed flame stabilized in a dump combustor: most features predicted theoretically are captured by DNS with very good accuracy as soon as flames with sufficient confinement are considered (small cross-section ratio between injection duct and combustion chamber). These include frequency of the instability and spatial structure of the mode via pressure and velocity fields. It was shown that theoretical predictions of the stability and the frequency of ITA modes entirely rely on the determination of the Flame Transfer Function. A simple visual approach based on a Bode diagram of the FTF and a critical gain expression are given for the determination of these characteristics. This approach gives very good agreement with DNS results in terms of frequency, while showing some limitations regarding stability limits. However, results confirm that ITA modes are more unstable in confined combustion chambers. Finally, the mechanisms involved in the intrinsic feedback loop, i.e. mode conversion at the corners of the chamber inlet plane and Kinematic Over-Restoration, were described.

Intrinsic ThermoAcoustic instabilities may appear as a curiosity since their formal derivation requires anechoic boundaries. Thermoacoustic modes in real combustion chambers feature acoustic feedback through nozzles, compressors, turbines, etc., which are obviously non anechoic [START_REF] Duran | 19th AIAA/CEAS Aeroacoustics Conference[END_REF][START_REF] Marble | [END_REF]. However, the intrinsic feedback responsible for ITA modes can still play a role in these cases and interact with standard feedback mechanisms such as acoustic reflexion or entropyacoustic mixed modes [START_REF] Duran | 19th AIAA/CEAS Aeroacoustics Conference[END_REF]. Moreover, systems where ITA modes exist may react differently to changes of impedance at inlet or outlet compared to usual combustors. Increasing acoustic losses may actually make these systems more unstable, a property which is quite unexpected in the combustion instability community.

Fig. 1 .

 1 Fig. 1. Configuration of the 2D laminar premixed flame, colored by velocity magnitude. The steady flame position is indicated by a black iso-contour of heat release rate. Downstream and upstream propagating acoustic waves A + j and A - j are also represented in the inlet duct (i = 1) and the combustion chamber (i = 2).

Fig. 2 .

 2 Fig. 2. Stability map for ITA modes. Critical FTF gain threshold as a function of the cross-section area ratio Fc = f (S 2 /S 1 ) ( ) for a temperature jump θ = 6.3.

Fig. 3 .

 3 Fig. 3. Theoretical structure of an ITA mode: pressure () and velocity ( ). Left: modulus. Right: phase with different scales for pressure and velocity. The application case corresponds to the numerical setup described in Section 3.

Fig. 4 .

 4 Fig. 4. Flame transfer function of the flame of Fig. 1 for a confinement corresponding to S 2 /S 1 = 6 (case REF-6.0). DNS () and experimental results for the unconfined flame ( ). From left to right: modulus and argument of F against forcing frequency. The reference point P is located 1 cm upstream of the combustion chamber in the inlet duct.

  period via the flame delay τ (Eq. (11)).

Fig. 5 .

 5 Fig. 5. FTF gain and phase computed for various confinements S 2 /S 1 = 1.5 ( ), 2.0 ( ), 3.0 ( ) and 6 ( ).

Fig. 6 .

 6 Fig.6. Illustration of the matching instability criteria for the gain (left) and phase (right) using a typical FTF of premixed laminar flames.

Fig. 7 .

 7 Fig. 7. Stability map of ITA modes for various confinements. Critical gain |F| = Fc(S 2 /S 1 ) for θ = 6.3 (

S 2

 2 /S 1 f = ω r /(2π) [Hz]

Fig. 8 .

 8 Fig. 8. Frequency map of ITA modes for various confinements. Theoretical predictions based on the Flame Transfer Function (). Unstable ( ) and stable ( ) results from DNS. Ordinates are oscillation frequencies in the linear regime.

Fig. 9 .

 9 Fig. 9. Global heat release rate at the onset and limit-cycle of an unstable ITA mode, from DNS for case REF-2.0.

Fig. 11 .

 11 Fig. 11. Normalized power spectral density (PSD) of integrated heat release rate in the limit-cycle ( ) of the unstable ITA mode presented in Fig. 9 (case REF-2.0). Spectrum obtained with Hann windowing and zero-padding.

1 Fig. 12 . 2 Fig. 13 .

 112213 Fig. 13. Acoustic waves traveling at the inlet boundary x = -5 cm (left) and outlet boundary x = +10 cm (right) of the combustor, from DNS with S 2 /S 1 = 2 (case REF-2.0). Downstream propagating waves A + j ( ) and upstream propagating waves A - j

Fig. 14 .

 14 Fig. 14. Comparison between theory () and DNS ( ) for the spatial structure of the first ITA mode.

Fig. 15 .

 15 Fig. 15. One cycle of the ITA mode in limit-cycle on case REF-2.0. The black isoline corresponds to a high value of the Q criterion and marks intense vortices [41]. The heat release rate is indicated by red isolines. The background color corresponds to the axial velocity field. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

More generally, for any causal temporal signal g ′ (t), ĝ(ω) would denote the corresponding Laplace transform, defined for any complex frequency ω = ωr + iω i .

3 To take into account mean flow effects (low Mach number regime), a different formulation using mass flow rates should be used instead of volume flow rates in the velocity jump relation.

From an input-output point of view, the temporal notion of causality is respected in Crocco's model since the input is the reference velocity taken at a previous time u ′ 1 (tτ ). Only when the acoustic decomposition of the reference velocity is performed does the causality paradox appear.

The Nyquist criterion invokes the concepts of open and closed loop transfer functions in the context of intrinsic thermoacoustic feedback, which are not mentioned in this paper, but have been discussed by Emmert et al.[3]. It implies that F(ω) has no poles, which is true for the n -τ model (see Eq. (14)), but can be difficult to justify for arbitrary shaped FTFs. To overcome this problem, another interpretation of the Nyquist criterion in the context of thermo-acoustics was proposed by Kopitz & Polifke[11], in which the FTF is seen as a conformal mapping and its angle-conserving properties are used to find the sign of ω i .

7 Recalling that the frequency response F(ωr) is the restriction of the FTF at real frequencies, solving Eq. (16) for F(ωr) is only an approximation. This is relevant for most laminar premixed flames whose phase appears to be independent of ω i , as can be shown on a phase plot analysis at complex frequencies (see for instance Appendix B of[START_REF] Cuquel | Dynamics and nonlinear thermo-acoustic stability analysis of premixed conical flames[END_REF]).

Obviously, any given disturbance (even numerical) from the stable solution would trigger the instability, the stable solution obtained in the case with non-zero reflection coefficients becoming an unstable equilibrium state in the anechoic configuration. The gaussian impulse was used here to control the timing of onset of the instability for comparison with other cases.

Acknowledgements

This work was granted access to the high-performance computing resources of CINES under the allocation x20142b7036 made by Grand Equipement National de Calcul Intensif. The support of the European commission, through the ERC advanced grant INTECOCIS GA 319067, is acknowledged. The authors also thank the Center for Turbulence Research at Stanford University for hosting them to work on this project during the Summer Program of July 2014.

Appendix A. Derivation of the stability criterion of ITA modes using Rayleigh's criterion

Rayleigh's criterion states that the temporal variations of the period-averaged acoustic energy E results in the balance of the period-averaged source term due to combustion R and the periodaveraged acoustic fluxes leaving the domain through the boundaries F j . Consequently, unstable combustion-acoustics coupling will be found for [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]:

where the definition of the acoustic source term and fluxes are:

where ε j = ( -1) j because û j is an oriented quantity, and the star ( * ) symbol denotes the complex conjugate. Using the analogous of Eq. ( 20) with the interaction index n instead of the FTF modulus, the acoustic source term and fluxes can be expressed in terms of the geometry and flame parameters:

Finally, injecting Eq. (A.3) in Rayleigh's criterion (Eq. (A.1)), one can easily recover the criterion for instability of Eq. ( 12):


Bao Anh Nguyen

Quoc Dung
email: pqdung@hcmut.edu.vn

Phan Ho

Chi Minh

Minh Duy

Kien Long Nguyen

Olivier Nguyen

Pascal Durrieu

Maussion

Parameter Identification Method for a 3-phase Induction Heating System

Keywords: induction heating, multiphase, modeling, identification, parameter tuning

This paper describes a new method for the on-line parameter estimation of an induction heating system. Simulations and experiments are presented in order to measure its impedance matrix for more exact control in closed loop. In previous papers, various parameter identification methods including off-line methods were introduced and compared for current inverters. It has been demonstrated that parameter identification is necessary to achieve good control of the inductor currents. A "pseudo-energy" method for a simple and fast implementation is compared to a classical "V/I with phase shift" method. They are experienced on a reduced power 3-phase coupled resonant system supplied with voltage inverters with satisfying results.

I. INTRODUCTION

Nowadays, induction heating systems are widely used in industry because they provide safety, cleanness, better performance and higher efficiency when compared to the classical heating systems, convection heating systems and radiation heating system because the heat could be generated deep inside the material. They are usually used in metal industry for many applications such as heating, welding, melting, drying and merging [START_REF] Lucia | Induction Heating Technology and Its Applications: Past Developments, Current Technology, and Future Challenges[END_REF] [START_REF] Yilmaz | Medium-Frequency Induction Melting Furnace as a Load on the Power System[END_REF].

With the development of power electronics and microprocessors, higher temperature, more precise temperature profile and power could be achieved. The structure of a classical induction heating system consists of power converters, micro-controller and resonant circuit. The power converter configurations can be half-bridge, fullbridge, the single switch or multi-level converter [START_REF] Millán | Resonant inverter topologies for three concentric planar windings applied to domestic induction heating[END_REF]. Numerous applications can be found for domestic appliances in [START_REF] Lucia | Induction Heating Technology and Its Applications: Past Developments, Current Technology, and Future Challenges[END_REF][3] [START_REF] Lucía | Multiple-Output Resonant Matrix Converter for Multiple Induction Heaters[END_REF] for many years. The resonant circuit consists of resistor, inductor and capacitor. It is well known in induction heating, that the inductor is one of the most important parts of a resonant converter because it has to face high power, high frequency and high temperature. The inductance of inductor sets the inverter operating point and safe operating area and it is affected by the operating frequency and temperature. So, monitoring this parameter is very important to get the best performance and remains an issue.

Previous papers [START_REF] Souley | Optimization of the Settings of Multiphase Induction Heating System[END_REF] put in evidence that a precise control of the inductor currents in multi-phase systems was necessary to reach the reference temperature or power density profile in work piece. Multi-phase systems are good candidates to increase the heating power in industry applications. As a consequence, as the equivalent inductor depends on the load, system parameters are required for optimal controller tuning. As the load properties can vary (type of metal, thickness, form, temperature, width etc.) inverter loads will change a lot and on-line parameter identification becomes an issue. Some algorithms have been developed to monitor the impedance, they depend on the settings of estimation methods and may not always converge [START_REF] Zhe | Model Identification of the Continuous Casting Billet Induction Heating Process for Hot Rolling[END_REF]. A "pseudoenergy" method has been proposed [START_REF] Souley | Methodology to characterize the impedance matrix of multi-coil induction heating device[END_REF], but harmonics cause some errors in the calculation. All these algorithms can be qualified as "off-line", i.e. the parameters are calculated when the system does not operate. Another algorithm is "3D Numerical Modeling" which uses numerical method to identify the permeability, but it not used to calculate the impedance matrix [START_REF] Canova | Identification of Equivalent Material Properties for 3-D Numerical Modeling of Induction Heating of Ferromagnetic Workpieces[END_REF]. The performance will be increased if the impedance is monitored on-line while the system is still working. Different on-line methods have already been developed but they are limited to one-phase systems and the impedance matrix is rather complex [START_REF] Jimenez | Analysis and Implementation of FPGA-Based On-line Parametric Identification Algorithms for Resonant Power Converters[END_REF] [START_REF] Lucia | FPGA Implementation of a Gain-Scheduled Controller for Transient Optimization of Resonant Converters Applied to Induction Heating[END_REF]. Particle swarm optimization is used in [START_REF] Dominguez | Load identification of domestic induction heating based on Particle Swarm Optimization[END_REF] to achieve parameter identification of a system load in domestic induction heating, which equivalent single series resistanceinductance highly depends on temperature, pot material, etc. But PSO are known to be time consuming. Parameter identification could be based on an analytical state-space model of an electrical equivalent circuit and multi-model reset observers for domestic induction cookers in [START_REF] Paesa | Adaptive Simmering Control for Domestic Induction Cookers[END_REF]. A large amount of unknown pots (150 indeed) has been tested. This observer needs parameter settings, the hybrid observer behaves as a proportional-integral observer but the integral term is reset according to a specific reset map. Parameters extraction of the electrical equivalent model of a coupled double concentric coil for induction heating purposes is provided with Fourier transform in [START_REF] Carretero | First harmonic equivalent impedance of coupled inductive loads for induction heating applications[END_REF] but with high sampling rate. Integration requires several periods to improve evaluation accuracy of the parameters of the 2x2 matrix. Besides, as the "pseudo-energy" method [START_REF] Souley | Methodology to characterize the impedance matrix of multi-coil induction heating device[END_REF] and another one based on the calculation of some [I inv ,V]/I terms have also already been proposed for multiphase current inverters [START_REF] Nguyen | Modeling and parameter identification of a multi-phase induction heating systems[END_REF], in an off-line implementation. This paper deals with their on-line adaptation in classical single board computer dedicated to control, with low sampling rate, without specific measurement equipments and at a reduced computational cost. The final objective is to apply the pseudo-energy method to a multiphase induction heating system, with voltage inverters in order to identify its impedance matrix and to properly set the current controllers. The whole system and concept is firstly checked with simulation results and verified with experimental results on a reduced power system. The test bench is limited to the electrical part (no thermal part in here), composed of a reduced power 3-phase coupled resonant system supplied with voltage inverters.

II. THREE PHASE INDUCTION HEATING SYSTEM

A previous paper [START_REF] Souley | Optimization of the Settings of Multiphase Induction Heating System[END_REF] focused on a 3-phase induction heating system composed by three concentric inductor coils which are arranged face to face in a transverse flux and a disk plate i.e. the load to be heated. It has been shown on this device that parameter estimation was necessary and possible. In the present work, a reduced element system is proposed to emulate the electrical part of a three phase induction system with voltage control and with a resonant circuit in a serial association. The test bench is composed of six voltage inverters with a common dc source but only the first three inverters are used here for simplicity reasons and validation of the concept. The three inverters supply a three phase coupled load as showed in Fig. 1. Each phase includes a full-bridge converter, a resistor, a capacitor and a serial inductor situated on a common core, as it can be seen in Fig. 2. A mathematic model of this structure of coils is built in Matlab in order to study its theoretic impedance matrix, which is an also basement to compare with other methods. This model also allows to calculate the impedance matrix with different numbers of turns of each inductor or air gaps between the coils. The process includes two stages: calculating reluctances by equivalent magnetic circuit presented in Fig. 3. and calculating impedance matrix by the V/I method with phase shift. According to (1) to ( 6), when phase 1 is supplied by a sinusoidal voltage V S , inductor current I 1 is calculated from V S and Z 11 . Equivalent magnetic circuit will help to calculate fluxes ࢥ 1 , ࢥ 21 and ࢥ 31 which flow through inductor coil 1, coil 2 and coil 3 respectively. This leads to the induced voltages V 21 and V 31 across inductor 2 and inductor 3, which allow calculating Z 21 , Z 31 (R 21 , L 21 , R 31 , L 31 ). The other terms, Z 22 , Z 12 , Z 32 and Z 33 , Z 13 , Z 23 are identified by the same process on phase 2 and phase 3.

2 e R i ii i N L (1) ii ii ii Z R j L Z (2) S i ii V I Z (3) e R i i i i I N I (4) ij i ij V N i ZI (5) ij ij j V Z I (6)
L ii : self inductance of inductor i, N i : number of turns of inductor i, R ei : equivalent reluctance of case i, ࢥ i : flux through coil I,

As showed in Fig. 4, the power electronic part is composed of six independent voltage invertors with 6 separate control signals. This inverter consists of six independent full-bridge converters (special integrated circuit) with dead-time and thermal shutdown implemented inside. The outputs of the inverter legs can be connected in parallel or series. Six protection systems are used against overcurrent or short-circuit protection. Technical parameters of this reduced power test bench are expressed as max DC voltage: 48 V, wide range of switching frequency: up to 100 kHz, and max current: 5A Capacitor bank has been calculated in order to compensate the reactive power on the three phases. This leads to C = [3.74 µF, 3.37µF, 3.75µF]. The switch control signals can be generated either from a pulse generator or from Matlab/Simulink with a digital control board. Moreover, these control signals are isolated from the power stage for safety reasons. For example, inverter voltage n°1 and current are presented in Fig. 5, for example while the other currents can be seen in Fig. 6. In Fig. 3, the load consists of three inductors organized in series configuration and a ferrite core as work piece coupled. A matrix model of the system is given in [START_REF] Souley | Methodology to characterize the impedance matrix of multi-coil induction heating device[END_REF] where sinusoidal currents I 1 , I 2 and I 3 feed the three coils through the three inverters. The so-called "impedance matrix" Z in (8)( 9) and ( 10) carries all the information on the state of the load via its parameters.

ܸ ଵ ഥ ܸ ଶ ഥ ܸ ଷ ഥ = ܴ ଵ + ܮ݆ ଵ ߱ ܯ݆ ଵଶ ߱ ܯ݆ ଶଵ ߱ ܴ ଶ + ܮ݆ ଶ ߱ ܯ݆ ଵଷ ߱ ܯ݆ ଶଷ ߱ ܯ݆ ଷଵ ߱ ܯ݆ ଷଶ ߱ ܴ ଷ + ܮ݆ ଷ ߱ ൩ ܫ ଵ ഥ ܫ ଶ ഥ ܫ ଷ ഥ (7)
ܴ , ܮ : self-resistance and inductance for inductor I, ܯ : mutual inductance between phase i and j.

ܼ ҧ = ܼ ଵଵ തതതത ܼ ଵଶ തതതത ܼ ଵଷ തതതത ܼ ଶଵ തതതത ܼ ଶଶ തതതത ܼ ଶଷ തതതത ܼ ଷଵ തതതത ܼ ଷଶ തതതത ܼ ଷଷ തതതത (8)
ܼ ప,ప തതതത = ܴ + ܮ݆߱ (9)

ܼ ప, തതതതത = ܯ݆߱ , (10)

III. IDENTIFICATION BY V/I METHOD WITH PHASE SHIFT

In this method, all the three capacitors have been removed for the load to be supplied directly by a sine wave from a frequency generator. Consequently, only resistors, self-inductors and coupling terms are identified. In the first experiment, phase 1 is supplied with a sinusoidal voltage while phase 2 and phase 3 are opened. The necessary values to be measured are V 1rms , V 2rms , V 3rms , I 1rms and phase shifts between V 1 /I 1 , V 2 /I 1 , V 3 /I 1 in Fig. 7. Then, the impedance parameters Z 11 , Z 21 , Z 31 are calculated by [START_REF] Dominguez | Load identification of domestic induction heating based on Particle Swarm Optimization[END_REF].

i ij ij ij ij ij j V Z R jX R j L I Z (11)
Similarly, measurements on phase 2 and phase 3 will give Z 12 , Z 22 , Z 32 and Z 13 , Z 23 , Z 33 . It is easy to understand that this method does not lead to accurate. Indeed, I measurements that could be noisy, will lead to division by a false term and false impedance values, while phase measurements via zero crossing are often difficult to achieve, particularly in an on line automatic mode. In the pseudo-energy method, the capacitors are not removed in order to stay in the resonant mode. The identification process is composed of three successive steps during which the 3 phases are supplied independently and successively. This method was first described in [START_REF] Souley | Methodology to characterize the impedance matrix of multi-coil induction heating device[END_REF] in an off-line mode thanks to a specific and complex apparatus (Rogowski coils and Le Croy oscilloscope). It will be first implemented in the off-line mode in real time single board controller, with a reduced sample frequency and 12 bits AD converters. Then, the on-line mode will be presented in section V.

In the first step, phase 1 is supplied by the square voltage generated by the corresponding inverter while phase 2 and phase 3 are opened (their inverters control signals are not activated). Identification of the terms of the first column of matrix Z is achieved. Then, the second step consists in supplying phase j while I 2 and I 3 are null, in order to measure the parameters of column 2 in the Z matrix. Finally, step 3 is the same for phase 3. At resonant frequency, a sinusoidal current appears in phase i. The measured values to calculate R ii , R ij , R ik , X ii , X ij and X ik , i.e. terms in the first column of [START_REF] Canova | Identification of Equivalent Material Properties for 3-D Numerical Modeling of Induction Heating of Ferromagnetic Workpieces[END_REF], are the voltages across the RL parts of phase i, j and k which are V i , V j , V k , respectively, and current in phase i, named I i . Current in phase j and k must be equal to zero in order to cancel the coupling terms in phase i and in [START_REF] Souley | Methodology to characterize the impedance matrix of multi-coil induction heating device[END_REF]. Assuming that the voltages and currents are sinusoidal, the following equations can be written. Starting from the same principle as in the classical method, the RMS values of the currents [START_REF] Paesa | Adaptive Simmering Control for Domestic Induction Cookers[END_REF] and the voltages are first determined and cos(ij) is calculated by [START_REF] Carretero | First harmonic equivalent impedance of coupled inductive loads for induction heating applications[END_REF] with n the number of integration periods.

dt i T . n 1 I nT 0 2 rms ³ (12) j i j i I . V i . v ) cos( ² ¢ M (13) ij nT 0 j i j i P dt i . v T . n 1 i v ² ¢ ³ (14) 2 j ij ij I P R (15) 2 j 2 ij 2 ij j i I P S X ( 16
) j i ij I . V S (17) 2 2 i i i P S Q (18)
Of course, accuracy of cos(ij) and the other terms will be increased with the number of integration periods. The pseudo-power average is given by ( 14). The coupling terms involving i and j quantities with ij are called "pseudopower" because they are the product of currents and voltages that do not necessarily coexist in the same circuit. They represent power in terms of unity but may have no physical meaning. The mutual resistance is obtained by dividing the "pseudo active power" by the square of the corresponding RMS inductor current as described in [START_REF] Nguyen | Robustness of a Resonant Controller for a Multiphase Induction Heating System[END_REF]. The calculation of the capacitance or reactance (16) requires the calculation of "pseudo apparent power" (17) which is calculated as the product of the supplied RMS current and the voltage induced. In fact, as it can be easily seen in Fig 8, voltages are not sinusoidal. Consequently, all the measured values have to be low-pass filtered through first order digital filters at the fundamental frequency (1500Hz). As depicted in Fig. 8 where phase 2 and 3 are open circuited, current in phase 1 is perfectly sinusoidal while the other currents in phase 2 and 3 are equal to zero. The parameter identification runs properly. Fig. 10 gives the bloc diagram which is implemented on the DSP board to measure R 11 , R 21 , R 31 , X 11 , X 21 and X 31 parameters during the first identification step, when phase 1 is supplied and currents should be null in the other phases. The importance of this last point will be discussed later. It will be shown that the system must operate at reduced power just to ensure the currents in phase 2 or 3 will not exist due to the back emf which is generated by the coupling effects between phase 1 and 2 or phase 1 and 3. In that case, currents will circulate through the reverse diodes of the bridges. When the control signals of the inverter 2 is OFF, even if the phase 1 voltage is maximum as in Fig. 9, the induced voltage on phase 2 is not higher than the supply voltage. Thus there is no current flowing through the inverter 2 diodes and this leads to good identification results, such as listed in Table IV. Finally, it is worth noticing that these rather good identification results are obtained with a reduced sample period, i.e. only 12 times the fundamental frequency. This consideration will help to reduce the computational burden and is positive for real time implementation of this identification method.

Moreover, capacitors may vary with temperature or may differ from the supplier values given in the data sheet. This method will also provide estimation of C, the capacitor bank as listed in Table II. Voltages across them could be obtained through (19) assuming as previously stated, that voltages across the RL parts are measured with dedicated sensors and filters. Once each inverter voltage can be deduced from the dc bus voltage and the duty cycles, and not directly measured, capacitors derive from (20). Of course, these voltages across C i are far from sinus and also need low-pass filtering. Fast identification of the system parameters becomes an objective of dramatic importance in the very particular case of rolling plates. Their characteristics such as type, dimensions, temperature, etc. when placed under the induction generator, can suddenly change. It is particularly true when the industrial process needs to heat two different types of work pieces one after another. Assuming that a warning signal is sent during (or before) the change, the three identification steps can be launched successively, just during the duration of the change, i.e. time for the inductor crossing. Under the assumption of a 1m/s speed of the rolling piece and because each identification step is only few periods long, 100 periods i.e. 66ms for example, the work piece displacement would be no more than 6.6cm, just the dimension of the width of the soldering strip between the two work pieces! As a consequence, the resulting temperature disturbance will certainly be negligible. As seen in Fig. 11, when the gate signal at the beginning of the identification step, switches phase 2 and phase 3 from ON to OFF, the corresponding currents decrease down to zero. In order to reduce the estimation error, the number of identification periods T off must be high compared to the transients. Moreover, identification must start after the transients. A survey on this value is given in Table III. A convenient value is 100 times the system fundamental period which will lead to 7.2% error.

RLi invi Ci V V V (19) > @ > @ ci i i V I C Z (20)
Because the system is supposed to work at different operating points and with some resonant frequency variations, additional experimental results are provided in Table V with satisfying results. It has been shown in [START_REF] Nguyen | Robustness of a Resonant Controller for a Multiphase Induction Heating System[END_REF] that resonant controllers which have been successfully used in this type of application are robust enough to cope with misalignment and parameter estimation error. Then, control control of the inductor currents will be performed correctly.

VI. CONCLUSION

This paper describes a new method for parameter estimation of a multiphase induction heating system. Simulations and experiments were presented in order to measure its impedance matrix for more exact control in closed loop. The "pseudo-energy" method is applied in comparison with the classical V/I method. It achieves good results on a reduced power 3-phase resonant system with voltage inverters, even with a limited sample frequency.

Future work will deal with the on-line application of this method on a 6 phase system and on a real full power induction heating system with transverse flux and concentric coils. Moreover, heating magnetic steel with a non-linear behavior and the corresponding parameter identification will also be an important issue in the future.

Fig. 1 .

 1 Fig. 1. System PSIM model with voltage inverters, serial capacitors, coupled-inductors and sensors

Fig. 2 .

 2 Fig. 2. Three coupled coils

Fig. 3 .

 3 Fig. 3. Three coupled coils and their equivalent magnetic circuits

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. The 6-phase voltage inverter

Fig. 7 .

 7 Fig. 7. Voltage V1 (CH1), induced voltages V2 (CH2), V3 (CH3) and current I1 (CH4) IV. OFF-LINE IDENTIFICATION BY PSEUDO ENERGY METHOD

Fig. 8 .

 8 Fig. 8. Phase 1 and 2 voltages and currents during the identification step on phase 1 with modulated angle on inverter 1 and phases 2 and 3 open circuited. CH1:V1; CH2:I1; CH3:V2 and CH4:I2

Fig. 9 .

 9 Fig. 9. Phase 1 and 2 voltages and currents during the identification step on phase 1 with reduced current in inverter 1. CH1: V1; CH2: I1; CH3 :V2 and CH4 :I2

Fig. 10 .

 10 Fig. 10. Bloc-diagram for the identification process on phase 1

Fig. 11 .

 11 Fig. 11. Command signal of identification process (CH1), phase 2 voltage V2 (CH4) and phase 2 current I2 (CH2)

TABLE I PARAMETERS

 I

			OF THE IMPEDANCE MATRIX
	i	j	1	2	3
	1	R () L (mH)	2.27 2.93	0 0.086	0 0.027
	2	R () L(mH)	0 0.086	2.44 3.15	0 0.171
	3	R () L (mH)	0 0.027	0 0.171	2.74 2.86

TABLE III INFLUENCE

 III OF TOFF OF THE IDENTIFICATION RESULTS

	toff	Time (ms) diference
	100.T0	66.7	-7.2%
	20.T0	13.3	-25.2%
	5.T0	3.3	-31.8%