Parameter identification of two coupled oscillator model for pure intrinsic thermoacoustic instability

Roeland Wildemans 1

Viktor Kornilov Department of Mechanical Engineering, Eindhoven University of Technology Eindhoven, The Netherlands

Ines Lopez Arteaga Department of Mechanical Engineering, Eindhoven University of Technology Eindhoven, The Netherlands KTH Royal Institute of Technology, Department of Engineering Mechanics, Marcus Wallenberg Laboratory Stockholm, Sweden

ABSTRACT In this paper, a nonlinear phenomenological model of two coupled oscillators is presented, which is able to describe the rich nonlinear behaviour stemming from an unstable pure Intrinsic Thermo- Acoustic (ITA) mode of a simple combustion system. In a bifurcation analysis of a pure ITA mode, it was observed that for increasing mean upstream velocity the flames loose stability through a supercritical Hopf bifurcation and subsequently exhibit limit cycle, quasi-periodic and period-2 limit cycle oscillations. The quasi-periodic oscillations were characterized by low frequent amplitude and frequency modulation. It is shown that a phenomenological model consisting of two coupled oscillators is able to reproduce qualitatively all the di ff erent experimentally observed regimes. This model consists of a nonlinear Van der Pol oscillator and a linear mass-spring-damper oscillator, which are nonlinearly coupled to each other. Furthermore, a parameter identification of the model parameters is executed, which shows a good quantitative match between the model and the experimental data. Hence, the presented phenomenological dynamical model accurately describes the nonlinear self-excited acoustic behaviour of premixed flames.

1. INTRODUCTION

Thermo-acoustic instabilities pose a severe challenge in the development of clean burning techniques using lean fully premixed laminar flames. These instabilities can typically be divided in two types, namely (i) modes with a duct acoustic origin and (ii) modes with a burner / flame Intrinsic Thermo- Acoustic (ITA) origin [1]. The latter type of mode is caused by an internal feedback loop independent of the up- and downstream acoustics as shown by Bomberg et al. [2] and Emmert et al. [3]. The

1 r.wildemans@tue.nl

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

existence of an unstable pure ITA mode, which is not a ff ected by the acoustics of the setup, was demonstrated experimentally by Wildemans et al. [4]. This setup provides close to anechoic boundary conditions that enables the decoupling of the ITA mode from the acoustics. It was shown that the self-excited ITA mode reveals a rich variety of nonlinear oscillations. An experimental bifurcation analysis of the pure ITA mode was conducted by Wildemans et al. [5]. For increasing values of the bifurcation parameter, chosen as the mean velocity in the burner deck holes, subsequently a stable fixed point, limit cycle oscillations, quasi-periodic oscillations and period-2 limit cycle oscillations were observed. Besides the bifurcation analysis, a phenomenological model was proposed that was able to reproduce all observed attractors. Accurate models of the flame dynamics are an important tool in understanding and mitigating thermo-acoustic instabilities. Since the flames exhibit nonlinear behaviour, nonlinear models are required to accurately describe the flame dynamics. In most of these models, the heat-release rate is described by a static, nonlinear function of the (delayed) velocity perturbations, e.g, [6–9]. These models are capable of predicting the frequency and amplitude of self-excited thermo-acoustic oscillations originated from the acoustic modes of the system. Most of these models require the acoustics to produce self-excited oscillations. However, in case of an anechoic acoustic embedding of the burner / flame, the feedback through the acoustic modes is not present. Therefore, most of these models are not able to describe the nonlinear behaviour of the ITA modes. For pure ITA modes to be present, flame-flow feedback should be incorporated in the model to establish the ITA feedback loop. In network models, this feedback loop is incorporated due to its structure, see Bomberg et al. [2] and Emmert et al. [3]. Nonlinear flame models that are identified by externally forcing the flames, e.g., [9, 10], could easily be incorporated into network models. However, since these models are identified based on the acoustic forced dynamics of the flames, it remains unclear if they are capable of producing the rich variety of nonlinear behaviour originated from the self-excited ITA mode. The main goal of this paper is to investigate if the phenomenological model, as introduced in Wildemans et al. [5], is able to quantitatively describe the nonlinear dynamics of a pure ITA. Thereto, the nonlinear dynamical model of two coupled oscillators is investigated in detail to understand the exhibited oscillations. Furthermore, a parameter identification is conducted to show that, besides the qualitative match, there is a quantitative match between the experimental time-series and the model output. The structure of the paper is as follows. In the first section, the experimental setup is introduced. Next, the main results of the experimental bifurcation analysis are presented. In Section 4, the model of coupled oscillators is introduced. Thereafter, the parameter identification procedure is discussed. In Section 6, the results of the parameter identification are shown and discussed. Finally, conclusions are drawn in the last section.

2. EXPERIMENTAL SETUP

Pure ITA modes are observed when the burner / flame configuration is decoupled from the up- and downstream acoustics of the setup. This is typically realized by using anechoic up- and downstream mu ffl ers. In Wildemans et al. [4] a setup is introduced that provide close to anechoic conditions in a wide frequency range. This setup is depicted in Figure 1 and consists of three main parts: the burner and an up- and downstream duct with their mu ffl ers (denoted by G and H, respectively). A 1 mm thick perforated brass burner deck is used that consists of 109 holes with a diameter of 2 mm and a pitch of 4 . 5 mm as shown in Figure 1. This burner deck is installed on a water cooled burner deck holder. Using mass flow controllers the equivalence ratio and the mean upstream velocity of the mixture can easily be controlled. The equivalence ratio in this paper is kept constant at φ = 0 . 70. The mean upstream velocity in the holes of the burner deck is varied and denoted by V . The setup is equipped with three sensors. A Constant Temperature Anemometer probe (B) is used to measure the velocity oscillations just upstream of the burner deck. At the downstream side, a Photomultiplier equipped

H

C

120 380

Burner deck

B A

D

4 . 5

E F

375

290

∅ 2

50

∅ 50

G

260

Figure 1: Schematic overview of experimental setup: A - burner deck, B - CTA sensor, C - PMT sensor, D - pressure transducer, E - loudspeaker inlet, F - mixture inlet, G - upstream mu ffl er and H - downstream mu ffl er. All dimensions are in mm.

with a OH ∗ light filter (310 ± 10 nm) is used to measure the OH chemiluminescence as an indicator for the heat-release rate. Furthermore, a pressure transducer (D) is used to measure the pressure oscillations in the upstream duct. It is possible to force the flames with a loudspeaker (E). However, in this paper only self-excited oscillations of the flames are considered. Hence, all experimental result discussed in this paper are obtained with a closed loudspeaker inlet. The experimental setup is discussed in more detail in [4,5].

3. EXPERIMENTAL BIFURCATION ANALYSIS

The setup reveals self-excited oscillations when the mean velocity in the burner deck holes is at least 0 . 50 m / s. It is shown by Wildemans et al. [5] that these self-excited oscillations originate from an unstable pure ITA mode. Furthermore, it was observed that the characteristics of the self-excited oscillations depend on the mean velocity in the burner deck holes. When this velocity is varied, di ff erent nonlinear oscillations are exhibited by the flames. In Wildemans et al. [5] the nonlinear dynamics of the pure ITA mode are experimentally investigated by means of an bifurcation analysis. In this bifurcation analysis the mean velocity in the burner deck holes, V , is used as a bifurcation parameter. In Figure 2 the bifurcation diagram is depicted. The dots represents the local maxima of the time- series of the acoustic velocity just upstream of the burner deck. The bifurcation diagram reveals four distinct regimes. The first regime, denoted by I, represents a stable flame, since the velocity perturbation u ′ = 0. When the bifurcation parameter increases to V = 0 . 50 m / s, the system bifurcates to limit cycle oscillations through a super critical Hopf bifurcation. Limit cycles are characterized by a single dot in the bifurcation diagram. However, the observed local maxima reveal a small spread due to noise. Regime III is observed in the range V = 1 . 50 − 2 . 40 m / s and is characterized by a spread of local maxima. This signature indicates a time-series with an amplitude modulation and the oscillations are classified as quasi-periodic oscillations. This regime is reached via a Neimark Sacker bifurcation. The time-series of the last regime, IV, reveal two distinct local maxima. This indicates a period-2 limit cycle which is reached via a period doubling bifurcation. Based on the bifurcation analysis it is observed that the pure ITA mode reveals a rich variety of di ff erent nonlinear oscillations. Hence, a proper model should be able to describe all observed regimes. In the next section, a model of two coupled oscillators is introduced that produces all four regimes.

Figure 2: Bifurcation diagram of pure ITA mode with bifurcation parameter V .

4. MODEL OF COUPLED OSCILLATORS

A first step in system identification is the selection of a model structure that is capable of describing the nonlinear dynamics of the flames of all four regimes. Many dynamical models exhibit limit cycle and period-2 limit cycle oscillations. However, the regime with quasi-periodic oscillations is special due to its distinct peak at the low modulation frequency as shown by Wildemans et al. [5]. Mostly, the frequency spectrum of systems that produce quasi-periodic oscillations does not show a distinct peak at the low modulating frequency, e.g., the model in Kuznetsov et al. [11]. However, the dynamical system as presented in Anishchenko et al. [12], is able to reproduce this distinct peak at the low modulation frequency. Furthermore, in Wildemans et al. [5] it was shown that this system is able to reproduce all observed regimes of the bifurcation diagram. Therefore, this model provides a promising model structure to describe the flame dynamics. In this section we discuss the important features of this model. The model represents two nonlinearly coupled oscillators and is given by:

˙ x =  m − dx 2  x + y − x ϕ, (1)

˙ y = − x , (2) ˙ z = ϕ, (3) ˙ ϕ = − γϕ − gz + γ Φ ( x ) , (4)

with

 1 , x > 0 0 , x ≤ 0 . (5)

Φ ( x ) = I ( x ) x 2 , I ( x ) =

The first two equations represent the Rayleigh equation, which can easily be trans formed to the famous Van der Pol equation. This oscillator is unstable at small values of x ( x < √ m / d ) and produces self-excited oscillations. At lar ge values of x , the oscillations will be damped and the oscillation amplitude saturates ( x > √ m / d ). Hence, the amplitude of the oscillations is determined by the ratio m / d . The second oscillator represents a damped oscillator, which is equivalent to a linear mass- spring-damper system. Herein, γ is the damping constant and g the spring constant. Both oscillators are nonlinearly coupled to each other. The coupling in the first oscillator, − x ϕ , a ff ects the damping through the first time derivative of the state of the second oscillator. Hence, when x and ϕ are in phase, negative damping is added and the amplitude of the first oscillator will increase. The coupling in the second oscillator, Φ ( x ), is only non-zero when x > 0. This coupling will destabilize the

second oscillator and enables it to produce oscillations, which will in turn a ff ect the first oscillator. These two coupled oscillators are able to generate complex dynamical oscillations due to their mutual interaction. The Rayleigh equation, as presented in Equations 1 and 2, is rewritten into a Van der Pol equation to obtain better fitting results. Furthermore, all terms in the dynamical system are parametrized to increase the flexibility in describing the experimental data. This results in the following set of di ff erential equations:

˙ x 1 = x 2 , (6)

˙ x 2 =  α − δ x 2 1  x 2 − β x 1 − µ 1 x 2 x 4 , (7)

˙ x 3 = x 4 , (8) ˙ x 4 = − ε x 4 − ζ x 3 + µ 2 Φ ( x 2 ) . (9)

In these equations, β and ζ are the sti ff ness parameters of both oscillators and mainly determine the oscillation frequency. The function Φ ( x ) is given by Equation 5. This system has an equilibrium at x 0 = [0 0 0 0] T , which stability is given by the eigenvalues of the Jacobian matrix of the system evaluated at the equilibrium point. The eigenvalues of the Jacobian matrix are given by:

λ 1 , 2 = α

2 ± 1

α 2 − 4 β, (10)

p

2

λ 3 , 4 = − ε

ε 2 − 4 ζ. (11)

2 ± p

Based on these eigenvalues, it is concluded that the model converges to a stable equilibrium point x 0 when α < 0 and ε > 0. Hence, for these parameter settings the first regime is obtained. In order to find the parameters for which the model reveals the oscillations of the other regimes, a parameter identification is conducted.

5. PARAMETER IDENTIFICATION

In Wildemans et al. [5] it is shown that the model as presented in Equations 6-9 is able to produce all observed regimes in the bifurcation analysis. In order to find a quantitative match between the model and the experimental data, the parameters of the model should be identified for the di ff erent regimes. Thereto, a cost function is required to evaluate a certain parameter set. A natural choice for a cost function is to compare the time-series of a certain model output with an experimental time-series. A similarity is observed between the evolution of the state x 2 and the experimental time-series of the acoustic velocity in the limit cycle regime. Both time-series reveal a relatively fast rise and a relatively slow decay. This feature of the signal shape is also observed in CFD simulations of a pure ITA mode by Courtine et al. [13] and Silva et al. [14]. Therefore, the state x 2 is chosen as the output of the model that should describe the experimental acoustic velocity. More specific, the time evolution of the state x 2 should match the experimental time-series data as close as possible. Hence, the following cost function should be minimized

N t X

 ˜ y j − h ( x ( t j , θ ) , t j , θ )  2 . (12)

J ( θ ) =

j = 1

Herein, ˜ y j is the j th sample of the experimental reference signal, N t is total number of samples, h is the output of the model given by h ( x , t ) = x 2 ( t ) and θ is the vector with the unknown parameters. The unknown parameters are all model parameters including the initial conditions of the states x 1 , x 3 and x 4 . Hence, θ is given by

θ = h α β δ µ 1 ε ζ µ 2 x 1 (0) x 3 (0) x 4 (0) i T . (13)

Note that the initial condition of x 2 is equal to the first value of the experimental time-series due to the definition of the cost function. Since the experimental time-series originate from self- excited oscillations and the experimental system is not externally forced, this parameter identification is executed on output data only. The cost function is minimized by using a limited memory quasi-Newton algorithm, namely the constrained limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS-B) [15]. For the minimization of the cost function, as given by Equation 12, the algorithm requires the derivative of the cost function with respect to the unknown parameters θ . This derivative is obtained by using the forward sensitivity of the states x i and the output y with respect to the unknown parameters θ . To compute these sensitivities, it is assumed that a nonlinear dynamical system in general can be written as ˙ x = f ( x ( t , θ ) , t , θ ) and y ( t ) = h ( x ( t , θ ) , t , θ ) with the initial conditions x ( t 0 , θ ) = x 0 ( θ ). The gradient of the cost function with respect to the parameter θ k is given by

N t X

∂ J ∂θ k = − 2

 ˜ y j − h ( x ( t j , θ ) , t j , θ )  s y j k , (14)

j = 1

with the output sensitivity with respect to the parameter θ k given by

∂θ k  ∇ x h j  ∂ x j

∂θ k + ∂ h j

s y j k = ∂ h

∂θ k . (15)

For a given model structure, ∇ x h and ∂ h /∂θ k are known functions which can depend on the states. Hence, these functions can depend on the evolution of the states over time. In the above equation, ∂ x j /∂θ k is the sensitivity of the states with respect to θ k . This state sensitivity is obtained by computing the partial derivative of the system equation with respect to the parameter θ k :

∂ ˙ x j ∂θ k  ∇ x f j  ∂ x j

∂θ k + ∂ f j

∂θ k = ˙ s x j k =  ∇ x f j  s x j k + ∂ f j

∂θ k . (16)

In this equation, ∇ x f and ∂ f /∂θ k are known functions and can depend on the states of the system. Equation 16 resembles an ODE and can be solved by forward time integration in parallel with the simulation of the model, Equations 6-9. The cost function is minimized as follows. First, the algorithm is initialized by initializing the parameter vector θ . Then the system equations and the state sensitivities are integrated with this parameter set. Based on this integration, the cost function value and its gradients are calculated. Based on these values, the L-BFGS-B algorithm computes new values for the parameter vector θ and the previous steps are repeated until the cost function reaches a minimum. Since the state sensitivities depend on the derivative of the the system equations, Equations 6-9, to the states ( ∇ x f ), the system equations should be continuous in the states. In the original model equations, the function Φ ( x ) is not continuous due to the Heaviside function. Therefore, this function is finally substituted by the following approximation:

˜ Φ ( x ) = 1 1 + e − 40 x x 2 . (17)

5.1. Experimental data and initialisation of parameter identification The parameter identification in this paper is based on the normalized acoustic velocity just upstream of the burner deck, u ′ / ¯ u , that is measured with a CTA sensor. This data is chosen for two main reasons. First, the shape of this signal reveals a remarkable similarity with the shape of the time-series of one of the model states, namely x 2 . Second, the CTA signals are less corrupted by noise compared to the pressure transducer and PMT sensor signals, which makes them suitable for the parameter identification. Before the experimental data is provided to the optimization algorithm, several pre- processing steps are taken to improve the results. The experimental data contains noise, which

Table 1: Relation between parameters of scaled and non-scaled model.

Parameter ˜ α ˜ β ˜ δ ˜ ε ˜ ζ

Expression T re f α T 2 re f β T re f δ T re f ε T re f ζ

might deteriorate the obtained results, therefore the data is filtered. Filtering the highly nonlinear experimental time-series with a low pass filter, will filter out the higher harmonics, that alters the nonlinear oscillations. Therefore, the experimental data is filtered with a nonlinear filter as described in Schreiber [16]. The experimental data is sampled with a frequency of 3000 Hz. It was found out that an increased sampling frequency resulted in a better parameter identification. With a higher sampling frequency, the typical oscillations over time are more precisely covered in the reference data. The sampling frequency of the filtered experimental data is increased with a factor 4 by using the Modified Akima cubic Hermite interpolation (makima) function in MATLAB. The dominant frequencies in the experimental time-series are in the range of 100 − 200 Hz, indicating that the response of the model should be in the order of milliseconds. Since the oscillation frequency of the model is mainly determined by the masses and spring sti ff ness values, this indicates that they should di ff er by roughly three orders of magnitude. Large di ff erences between parameters will make the minimization of the cost function with a gradient based algorithm a challenging task. To avoid these problems, the time variable of the model is scaled to a dimensionless time τ such that the model parameters will be of the same order. The dimensionless time is given by τ = t / T ref , in which T re f describes a typical timescale of the oscillating flames. A Van der Pol oscillator with mass and spring values of O (1) and a strength of the nonlinear damping that provides limit cycle oscillations with a comparable shape as in the experimental data, has a typical period time of τ = 5. Several characteristic times can be used to obtain a suitable reference time. For example, in Doehner et al. [10] the restoration time is used, which is the typical time for a flame to recover to its initial position after it is perturbed. In this paper, the period time of a typical limit cycle is used. At a velocity of V = 1 . 00 m / s, the limit cycle frequency is 150 . 3 Hz with a characteristic time period of t = 6 . 7 ms. Hence, the reference time to slow down the flame dynamics is given by

τ = 6 . 7 × 10 − 3

T ref = t

5 1 751 . 5 s . (18)

Applying this time scaling to the model in Equations 6-9 will result in a model of the same structure with the parameters as given in Table 1. Furthermore, the time of the experimental data is scaled similarly to match the new dimensionless time of the model equations. After all preprocessing steps, the unknown parameter vector θ should be initialized as a starting point for the optimization algorithm. Initially, the unknown parameters are all set to unity. However, with this initialization the algorithm is not always capable to find a minimum that provides a satisfactory match between the model output and the experimental data. In these cases, multiple initializations should be tried to obtain a satisfactory parameter set.

6. RESULTS

In this section, the results of the parameter identification for the di ff erent regimes are discussed. Based on the time-scaled data, the scaled parameters are determined by the L-BFGS-B algorithm. These parameters are transformed back according to the relations in Table 1 and the identified systems are simulated and compared with the experimental time-series.

6.1. Regime II: limit cycle oscillations After the supercritical Hopf bifurcation, limit cycles are observed in the range of V ∈ [0 . 50 , 1 . 50] m / s. A typical response at V = 1 . 00 m / s is selected for the parameter identification of the limit cycle

10 0

1

Experimental data Model response

Experimental data Model response

0.8

10 -2

0.6

0.4

10 -4

0.2

PSD

0

10 -6

-0.2

-0.4

10 -8

-0.6

0 0.005 0.01 0.015 0.02 0.025 Time [s]

0 100 200 300 400 500 600 700 Frequency [Hz]

(a)

(b)

Figure 3: Limit cycle oscillations at V = 1 . 00 m / s. In (a) the time-series are shown and in (b) the corresponding PSDs.

oscillations. In Figure 3a the time-series from the experiment and model ( x 2 ) are depicted. This result shows that the model is very capable of describing the limit cycle oscillations. The frequency and the shape of the limit cycle are predicted accurately by the model. The limit cycle is characterized by a bend close before the minima, which is also predicted by the model. The Power Spectral Density (PSD) of both time-series is depicted in Figure 3b. The oscillation frequency and its higher harmonics are well captured by the model. Especially the dominant frequency and the first harmonic reveal a similar strength. The second and third harmonic are slightly less strongly present in the model response. The corresponding model parameters are given in Table 2. As expected, the parameters of the scaled model are O (1).

6.2. Regime III: quasi-periodic oscillations The third regime is characterized by quasi-periodic oscillations and is experimentally observed for V ∈ [1 . 60 , 2 . 40] m / s. The response at V = 2 . 00 m / s is used as a characteristic time-series for the quasi- periodic regime. The experimental time-series and the model response are depicted in Figure 4a. The amplitude modulation as revealed by the experimental time-series is also visible in the model response. However, the amplitude modulation produced by the model is less strong, especially when the amplitude is minimal. The minimal amplitude of the oscillations is slightly overpredicted by the model. Furthermore, the dominant oscillation frequency is predicted accurately by the model as shown in the PSD in Figure 4b. Both time-series reveal a relatively similar frequency spectrum. However, the peak at the low modulation frequency is much higher in the experimental data. The many sidebands in the PSD at f ± nf L with n = 1 , 2 , ... indicate that frequency modulation is present in both time-series. To verify this frequency modulation, the instantaneous frequency of both signals is computed by using the Hilbert transform of both time-series. The Hilbert transform enables the calculation of instantaneous attributes of time-series, such as amplitude, phase and frequency. This transform extends a real valued signal x r ( t ) into an analytical signal x ( t ) = x r ( t ) + ix i ( t ) in which i represents the imaginary unit. The imaginary part of the signal x i is similar to the real part of the signal x r , but with a phase shift of 90 ◦ . The instantaneous phase and frequency are closely related, since the latter is the time rate of change of the first. The instantaneous phase is given by

! . (19)

φ ( t ) = arctan x i ( t )

x r ( t )

10 0

1.5 Experimental data Model response

Experimental data Model response

10 -2

1

0.5

10 -4

PSD

0

10 -6

-0.5

10 -8

-1

0 0.05 0.1 0.15 Time [s]

0 100 200 300 400 500 600 700 Frequency [Hz]

(a)

(b)

Figure 4: Quasi-periodic oscillations at V = 2 . 00 m / s. In (a) the time-series are shown and in (b) the corresponding PSDs.

220

Experimental data Model response

210

200

Frequency [Hz]

190

180

170

160

150

140

0 0.05 0.1 0.15 Time [s]

Figure 5: Instantaneous frequency in the quasi-periodic regime.

Hence, the instantaneous frequency is given by ω ( t ) = d φ/ d t . This derivative is numerically computed using total-variation regularization [17]. The instantaneous frequency of both time-series is depicted in Figure 5. These lines reveal a similar pattern for both time-series. However, the maxima of the model response are slightly shifted. The minimum and maximum frequencies are correctly predicted by the model. The model parameters are also given in Table 2. According to expectation, ˜ β is increased since this parameter determines the oscillation frequency of the Van der Pol oscillator and the oscillation frequency is higher for this regime compared to the limit cycle oscillations. Furthermore, it is observed that the damping in both oscillators is changed in terms of ˜ δ and ˜ ζ . The coupling parameters µ 1 and µ 2 are significantly increased for the quasi-periodic oscillations compared to the limit cycle oscillations. These changes in parameters indicate that many parameters influence the change in oscillations and there is not a clear relation between a single parameter and the experimental bifurcation parameter.

6.3. Regime IV: period-2 limit cycle oscillations The last regime is characterized by period-2 limit cycles, which repeats itself after two oscillations. This regime is experimentally observed for V ∈ [2 . 50 , 3 . 50] m / s. For the parameter identification, a

10 0

1.5

Experimental data Model response

Experimental data Model response

1

10 -2

0.5

10 -4

PSD

0

10 -6

-0.5

10 -8

-1

0 0.01 0.02 0.03 0.04 0.05 Time [s]

0 100 200 300 400 500 600 700 Frequency [Hz]

(a)

(b)

Figure 6: Period-2 limit cycle oscillations at V = 3 . 00 m / s. In (a) the time-series are shown and in (b) the corresponding PSDs.

Table 2: Model parameters of scaled model for three di ff erent regimes.

Limit cycle Quasi-periodicity Period-2 limit cycle

˜ α 0.1425 0.1418 0.1305 ˜ β 1.6274 2.8996 2.1279 ˜ δ 2.9151 1.2549 1.4922

˜ ε 0.0260 0.1762 0 ˜ ζ 1.4692 1.6437 0.5285

µ 1 0.4069 1.7803 0.7217

µ 2 0.9620 2.1183 0.0540

time-response at V = 3 . 00 m / s is selected. In Figure 6a it is shown that the model is able to describe this regime accurately. The frequency and the shape of the oscillations are corresponding with the experimental time-series. However, there is a small mismatch at the extrema. In Figure 6b the PSD of both time-series is depicted. The dominant oscillation frequency shows a good match. However, the sub- and higher harmonics are in general lower for the model response. This indicates that the shapes of the oscillations are not perfectly matching, which is also seen in Figure 6a. The identified model parameters are also given in Table 2. Compared to the quasi-periodic regime, the parameters ˜ α , ˜ β and ˜ δ does not reveal a significant change. However, the other parameters, which mainly a ff ects the second oscillator, have decreased significantly and are even lower compared to the values for the limit cycle regime. This might indicate that the change from quasi-periodic oscillations to period-2 limit cycles is mainly caused by the changed coupling and the changed influence of the second oscillator.

6.4. Discussion In this paper a phenomenological model of two coupled oscillators is proposed. The first oscillator resembles a Van der Pol oscillator, which is able to produce self-excited oscillations. The second oscillator is a standard mass-spring-damper oscillator. These two oscillators are nonlinearly coupled to each other. In Doehner et al. [10] also a model of two coupled oscillators is used to describe the nonlinear flame dynamics. That model consists of two identical Van der Pol oscillators which are linearly coupled. The model parameters are identified based on acoustically excited flames. This model describes the forced nonlinear behaviour of the flames accurately. However, it is unclear if

that model is also able to describe the nonlinear oscillations of the self-excited ITA modes. The nonlinear damping term x 2 ˙ x that characterizes the Van der Pol oscillator is also used in previous studies, e.g., [6–8]. However, in these studies a single oscillator is considered to describe both the flame and the acoustics. Overall, a model of two coupled oscillators provides a promising model structure for new nonlinear flame models due to its ability to describe both the forced and self-excited nonlinear behaviour of the flames. The presented model of two coupled oscillators is able to describe all observed regimes exhibited by the pure ITA mode. Especially the limit cycle oscillations are described very accurate. The model is also able to describe the quasi-periodic oscillations which are characterized by low frequent amplitude and frequency modulation. Both phenomena are captured by the model. In preliminary experiments it was observed that the discussed attractors are also observed for equivalence ratios slightly smaller and larger then φ = 0 . 70, indicating that the model is applicable for a broader range of operational conditions. Hence, this low-order phenomenological model provides insights in what type of nonlinear model is required to describe the intrinsic self-excited flame behaviour. This can guide the search for better nonlinear flame models. With these low-order models, insight in the (in) stability mechanisms can be obtained. These insights might be beneficial towards the development of stable combustors. Furthermore, low-order models are beneficial over high-order CFD models for designing control strategies to stabilize the nonlinear flames. The proposed model is a time-domain model, which can be directly used for time-domain simulations of thermo-acoustic systems. The proposed phenomenological model has also a few shortcomings. Firstly, it is yet unclear how this model is related to the heat-release rate produced by the flames. Therefore, this model is not yet a complete description of the flame dynamics. Furthermore, this model is a phenomenological model and thus mainly describes the observed phenomena. The downside is that there is no direct physical interpretation of the states and parameters of the model. An other remaining question is how the proposed model should be connected to an acoustic domain in order to simulate complete thermo-acoustic systems? A closely related question is if this model is able to describe more then the inherent ITA behaviour? Is the model also able to describe the flame response due to acoustic forcing of the flames?

7. CONCLUSION

In this paper, a phenomenological model is presented that is able to describe the nonlinear self- excited nonlinear dynamics of pure ITA modes. This model consists of two nonlinearly coupled oscillators, where the first represents a Van der Pol oscillator and the second one a standard mass- damper-spring oscillator. In an experimental bifurcation analysis, it was observed that the fully premixed flames reveal subsequently four di ff erent attractors, namely a fixed point and limit cycle, quasi-periodic and period-2 limit cycle oscillations. For all four regimes, model parameters are identified such that the presented phenomenological model is able to qualitatively and quantitatively describe the experimentally observed nonlinear attractors. The regime with quasi-periodic oscillations is characterized by low frequency amplitude and frequency modulation, which phenomena are both captured by the model. Hence, this model provides a promising model structure to represent the internal nonlinear flame dynamics. The presented phenomenological model is able to describe the nonlinear flame dynamics accurately. However, it does not provide a full description of acoustic behaviour of the flames. Therefore, the model should be extended such that it can describe the heat-release rate oscillations as well. The extension of this model should be an interesting question for future research. Other interesting questions are how to couple this model to an acoustic domain and what is the physical meaning of the states and model parameters?

ACKNOWLEDGEMENTS

This publication is part of the project STAbLE (with project number 16315) of the research programme Open Technology Programme which is (partly) financed by the Dutch Research Council (NWO).

REFERENCES

[1] A. Orchini, C.F. Silva, G.A. Mensah, and J.P. Moeck. Thermoacoustic modes of intrinsic and acoustic origin and their interplay with exceptional points. Combustion and Flame , 211:83–95, jan 2020. [2] S. Bomberg, T. Emmert, and W. Polifke. Thermal versus acoustic response of velocity sensitive premixed flames. Proceedings of the Combustion Institute , 35(3):3185–3192, jan 2015. [3] T. Emmert, S. Bomberg, and W. Polifke. Intrinsic thermoacoustic instability of premixed flames. Combustion and Flame , 162(1):75–85, jan 2015. [4] R. Wildemans, V. Kornilov, P. de Goey, and I. Lopez Arteaga. Experimental proof of existence of pure intrinsic thermo-acoustic modes. In Proc. of the 10th European Combustion Meeting , pages 331–336, Naples, Italy, 2021. [5] R. Wildemans, V. Kornilov, P. de Goey, and I. Lopez-Arteaga. Nonlinear dynamics of pure intrinsic thermo-acoustic modes. Combustion and Flame , under review, 2022. [6] G. Bonciolini, A. Faure-Beaulieu, C. Bourquard, and N. Noiray. Low order modelling of thermoacoustic instabilities and intermittency: Flame response delay and nonlinearity. Combustion and Flame , 226:396–411, apr 2021. [7] F. Gant, G. Ghirardo, and M.R. Bothien. On the importance of time delay and noise in thermoacoustic modeling. Journal of Sound and Vibration , 501:116067, jun 2021. [8] P. Kasthuri, V.R. Unni, and R.I. Sujith. Bursting and mixed mode oscillations during the transition to limit cycle oscillations in a matrix burner. Chaos , 29(4):043117, apr 2019. [9] N. Tathawadekar, N.A.K. Doan, C.F. Silva, and N. Thuerey. Modeling of the nonlinear flame response of a Bunsen-type flame via multi-layer perceptron. Proceedings of the Combustion Institute , 38(4):6513–6520, 2021. [10] G. Doehner, M. Haeringer, and C.F. Silva. Nonlinear flame response modeling using coupled oscillators. In Proceedings of SoTiC: Industry meets Academia , Munich, Germany, 2021. [11] A.P. Kuznetsov, S.P. Kuznetsov, and N.V. Stankevich. A simple autonomous quasiperiodic self- oscillator. Communications in Nonlinear Science and Numerical Simulation , 15:1676–1681, july 2010. [12] V. Anishchenko, S. Nikolaev, and J. Kurths. Winding number locking on a two-dimensional torus: Synchronization of qausiperiodic motions. Physical Review E , 73:056202, may 2006. [13] E. Courtine, L. Selle, and T. Poinsot. DNS of intrinsic thermoacoustic modes in laminar premixed flames. Combustion and Flame , 162(11):4331–4341, 2015. [14] C.F. Silva, T. Emmert, S. Jaensch, and W. Polifke. Numerical study on intrinsic thermoacoustic instability of a laminar premixed flame. Combustion and Flame , 162(9):3370–3378, 2015. [15] R. Byrd, P. Lu, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained optimization. Journal of Scientific Computing , 16(5):1190–1208, 1995. [16] T. Schreiber. Extremely simple nonlinear noise-reduction method. Physical Review E , 47:2401– 2404, Apr 1993. [17] R. Chartrand. Numerical di ff erentiation of noisy, nonsmooth data. ISRN Appl. Math. , 2011, 01 2011.