## Abstract

This paper presents a data-driven identification framework with the objective to retrieve a flame model from the nonlinear limit cycle. The motivation is to identify a flame model for configurations, which do not allow the determination of the flame dynamics: that is commonly for industrial applications where (i) optical access for nonintrusive measurements of velocity and heat release fluctuations are not feasible and (ii) unstable combustion is monitored via multiple pressure recordings. To demonstrate the usefulness of the method, we chose three test cases: (i) a classical Rijke tube; (ii) an experiment of a laminar flame (EM2C case), (iii) and a high-pressure, turbulent premixed flame (German Aerospace Center (DLR) case). The procedure is as follows: First, acoustic network models for the three cases are generated for which the in-house software taX is employed. Next, the acoustic network models are embedded in an optimization framework with the objective to identify flame parameters that match the experimental limit cycle data: based on the instability frequency and pressure amplitudes, we formulate physical constraints and an objective function in order to identify the flame model parameters gain *n*_{opt} and time delay *τ*_{opt} in the nonlinear regime. We demonstrate for the three cases that the identified flame parameters reproduce the unstable combustion processes and highlight the usefulness of the method for control purposes.

## 1 Introduction

Design requirements for a combustion chamber are prescribed by [1]:

thermal efficiency given by the overall pressure ratio;

adherence of emission rates such as NO

_{x}and CO;stability limits regarding flashback, flame blow-off, and combustion instabilities (CIs).

These counteracting features require a tradeoff for the best design solution, e.g., lean premixed flames reduce emission rates but are, at the same time, more prone to trigger CIs as lean flames are more sensitive to flow perturbations [2]. Typically, CIs appear at a stage, where the design process regarding emissions and efficiency is terminated. This means that major revisions of the design for a stable combustor will cause high cost and this is one reason why CIs represent a bottleneck in industrial design of low-emission combustion systems.

Combustion instabilities arise due to a constructive feedback loop between fluctuations of heat release, acoustics, and flow field [3–5]. In this situation, the energy released by the flame is transferred to the acoustic field [6,7] so that pressure fluctuations grow in time and reach high levels that enhance fatigue and, in extreme cases, risk of engine failure.

*n*and a time delay

*τ*at an angular frequency $\omega =2\pi f$, known as the

*n*–

*τ*model [10,11]

Unfortunately, the flame response is not known a priori and has to be obtained for each flame placed in a combustion chamber. This can be done computationally [12–15] or experimentally [16–21]. Typically, in experiments, the acoustic velocity $u\u2032$ is measured in the upstream flame vicinity (by, e.g., hot-wire or laser Doppler anemometry). The acoustic velocity modulation is induced by loudspeakers, which have to be mounted within the experimental test rig. Heat release fluctuations $Q\u02d9\u2032$ are measured, in the simplest case, by a photomultiplier, for which optical access has to be ensured to integrate the light intensity spatially inside the combustion chamber.

Determination of the FTF is thus a necessary but costly procedure. It is mostly done for well-equipped and controllable lab-scale configurations, which operate at atmospheric pressures and simplified inflow conditions that allow installation of an upstream forcing device. Only very few industrial applications are designed for FTF determination. For practical reasons, industrial configurations feature some acoustic pressure sensors mounted in the plenum or the combustion chamber but have by no means a flow forcing device or optical access to the combustion chamber.

The motivation of the paper is the following: Suppose that pressure time series of an industrial device exhibit CIs. Now, the objective is to eliminate the unstable combustion process, but due to high cost of a trial-and-error redesign (which is not even guaranteeing to find a stable configuration), the optimal approach would be to employ a numerical model of the prototype, which is computationally cheap and allows fast stability analysis. This can be done, e.g., by an acoustic network model [22–25], which computes, if an appropriate flame model is included, the stability of the system. However, for most industrial systems, the FTF is not known, so that no flame model is provided for the acoustic network model, and, hence, stability analysis cannot be performed.

The idea of this study is to bridge this gap by estimating an appropriate flame model of the unstable system based on experimental measurements. We propose a method that is based on an optimization procedure, which bounds the solution by minimization to nonlinear inequalities. To do so, it uses the pressure data measured inside the plenum and combustion system together with an acoustic network model to give an estimate of the FTF in its most simple form by a gain and a time delay of the flame at the instability frequency (*n*–*τ* model). Although this flame model comes along with simplifications, it is sufficient to estimate the flame response during limit cycle [26]. We present an identification methodology, which allows combining a network model with experimental data to retrieve a data-driven *n*-*τ* model for the limit cycle. To demonstrate the robustness of the identification strategy, three cases are presented that cover a wide range of applications: (i) a classical Rijke tube; (ii) an experiment of a laminar flame (EM2C case); and (iii) a high-pressure, turbulent premixed flame (German Aerospace Center (DLR) case).

The paper is structured as follows: Sec. 2 describes briefly the employed acoustic network model followed by Sec. 3, where the identification procedure is presented, in which the acoustic network model is embedded. Next, this methodology is applied to a Rijke tube (Sec. 4), a laminar flame (Sec. 5) and finally to a turbulent flame operated at high-pressure (Sec. 6).

## 2 Acoustic Network Models as Fast Design Tools

*f*and upstream

*g*within the subsystems and are defined as

Acoustic velocity $u\u2032$ and pressure fluctuations $p\u2032$ are linked via characteristic acoustic wave amplitudes *f* and *g*. From Eq. (2), primitive variables can be computed via $u\u2032=f\u2212g$ and $p\u2032/\rho c=f+g$. At inlet and outlet boundaries, the characteristic waves express the reflection coefficient as $Ri=g/f$ and $Ro=f/g$, respectively. The subsystems are written as a set of homogeneous equations, where the characteristics wave amplitudes are unknowns. The determinant of the interconnected system delivers the stability properties of the system. Since taX, the in-house acoustic network software, utilizes a generic state-space structure, the matrices coefficients are frequency independent and constant, so that the numerical solution of the eigenvalue problem becomes linear. As a result, a highly flexible tool for thermoacoustic stability analysis is available, in which all system parameter such as geometrical lengths, acoustic boundary conditions, mean flow effects, acoustic losses, and the flame model parameter can be tested within seconds.

The workflow for regular stability analysis is shown in Fig. 1. It determines the growth rate of each eigenfrequency of the system and gives access to the spatial distribution of acoustic pressure and velocity fields. However, in this case, the flame element (FTF or *n*-*τ* model) has to be provided by an external work process using experimental or computational methods. The identification procedure for the FTF estimate presented in this work turns the workflow around: based on the measured pressure recordings, a data-driven flame model given as *n*_{opt} and *τ*_{opt} is identified. To do so, the computationally efficient taX software is integrated into an identification framework using a nonlinear optimization procedure based on the interior-point method for the inverse problem presented in Sec. 3.

## 3 Identification strategy

The objective function $fobj(x)$ describes the output that has to be minimized (or maximized). This is done by variation of the control variables, here summarized as *x*. We have also introduced equality constraints denoted as *h*(*x*), which bound the solution space during optimization. The problem discussed previously can be treaded in this framework of a constrained nonlinear optimization and is described next.

*n*and time delay

*τ*, that reproduce together with the network model the experimental behavior of the unstable combustor within a small error margin. Therefore, the experimental pressure data are used: first, the target instability frequency

*f*

_{exp}is introduced into the objective function

*f*

_{taX}is the instability frequency of the network model. Note that in this work, only the first unstable mode is considered. This very simple objective function minimizes the error in frequency. Besides the objective function (Eq. (5)), physical constraints are introduced, that make use of the amplitude of the pressure data

where *i* refers to the number of pressure measurement. The physical constraints are not limited to pressure recordings as acoustic velocity fluctuations can also be used (as will be shown in Sec. 5). The optimization strategy is subdivided into two optimization loops (Fig. 2): a local loop evaluates the network model for each parameter set of the flame model while the global loop [30] ensures a restart of the randomly distributed flame parameters: a necessary step to avoid local minima of the objective function [31]. The description of the local optimization loop is as follows:

- The initial, randomly chosen flame model parameters
*n*_{ini}and_{t}*τ*_{ini}, denoted as the design parameters, are read by the network model. Typically, the values of both parameters have different orders of magnitude ($O(1)$ for the gain and $O(10\u22122)$ for the time delay) and units. It is then convenient to introduce normalized design parameters, which are normalized by their minimal and maximal parameter range_{t}$nnorm=ninit\u2212nminnmax\u2212nmin$(7)Minimal and maximal values are given in Table 1. Similarly, the time delay is normalized. With the normalized parameters, it is possible to introduce a minimum step width

*δ*in the parameter space to accelerate the optimization process (Table 1). Run the network model with initial design parameters. The frequency, spatial pressure amplitude distribution, and growth rate of the first unstable eigenmode

*f*_{taX}, $p\u2032taX$ and*α*are computed, respectively.The computed frequency

*f*_{taX}is used to evaluate the objective function (Eq. (5)). For errors of ±20 Hz, the optimization loop is stopped and a new set of parameters is used. This is an a arbitrary chosen error margin that ensures that only the instability frequency of interest is considered.The physical constraints defined in Eq. (6) are evaluated. For errors

*ϵ*larger than the defined error margin (Table 1) between the experimental and network model pressure amplitude $p\u2032\u2009exp\u2009$ and $p\u2032taX$, the solution is abandoned. We use at least two constrains for which at least two experimental data sets have to be provided. This standard optimization problem with nonlinear equality constraints is solved using the interior-point method proposed by Ref. [32], which is provided by matlab.Depending on the error of the objective function and the constraints, the local optimization can be repeated. The optimization tolerance on the objective function and the constraints is set to $O(10\u22121)$. Once a sufficient parameter set has been obtained that fulfills the objective function and the physical constraints, the local optimization is stopped.

Termination of the local optimization initiates the global optimization loop, which uses a randomly chosen set of new design parameters

*n*_{ini}and_{t}*τ*_{ini}. The global optimization process is repeated 40 times in order to cover a wide range in the predefined min. and max. flame parameter for gain and phase. This iterative process is repeated, until the global minimum is found. Note that we have tested the number of global optimization runs and found no significant change in the optimization results. All relevant identification parameter are given in Table 1._{t}

Parameter | Value |
---|---|

Range for gain n | $0.1\u22123.0$ |

Range for time delay τ (ms) | $0.5\u22125.0$ |

Min. step width δ | $1\xd7102$ |

Optimization tolerance ϵ | $1\xd710\u22121$ |

Global optimizations G_{n} | 40 |

Parameter | Value |
---|---|

Range for gain n | $0.1\u22123.0$ |

Range for time delay τ (ms) | $0.5\u22125.0$ |

Min. step width δ | $1\xd7102$ |

Optimization tolerance ϵ | $1\xd710\u22121$ |

Global optimizations G_{n} | 40 |

## 4 The Rijke Tube—A Pedagogical Example

### 4.1 Setup and Reference Data.

The Rijke tube is an open duct in which a heated gauze (or a thin flame) is placed (Fig. 3). The flame model is given by *n *=* *2 and $\tau =0.353\u2009ms$. The configuration dimensions and the flame parameter are arbitrarily chosen and do not correspond to a particular experiment. Hence, the reference data are synthetic and are used to apply the identification method by way of example.

This Rijke tube is constructed within the network model and includes the elements shown in Fig. 4. Now, the acoustic modes of the system can be computed with taX. The first mode of the system is at 378 + j27 Hz and is unstable. Pressure fluctuations are extracted at the probe locations 1 and 2 (*P*_{1} = 25.7 Pa and *P*_{2} = 47.96 Pa in Fig. 3). The instability frequency and the pressure amplitudes serve as reference data for the identification procedure performed in Sec. 4.2.

### 4.2 Identification of the Model Parameters.

In the foregoing section, we have generated the reference data of the Rijke tube by employing the regular work flow (Fig. 1). Now, the objective is to retrieve the flame parameters (*n *=* *2 and $\tau =0.353\u2009ms$) with the identification framework (Sec. 3). To do so, the target frequency is set to $f=378\u2009Hz$ and the physical constraints (Eq. 6) are given by *P*_{1}= 25.7 Pa and *P*_{2} = 47.96 Pa. The optimization settings are given in Table 1.

The evolution of the optimization process (Fig. 5) shows the value of the objective function, which is similar to an error made on the target frequency. The size of the circle displays the error made on the physical constrains: large circles indicate unsatisfied physical constraints and are not considered. Solutions with large errors on the objective function feature also large errors on the constraints. Similar pattern are found for the laminar and the turbulent combustion systems presented later (not shown). A number of possible solutions fulfill the constraints with small errors for the objective function but only one minimum is found: the identified flame parameters are *n*_{opt} = 1.99 and $\tau opt=0.305\u2009ms$ within the defined optimization tolerances. Using the identified flame parameters with the regular work flow, the frequency of the first unstable mode is 377.8 + i22.6 Hz. The mode shape remains the same and no other mode is unstable in the stability map (not shown). Overall, the optimization framework reproduces in good agreement the reference flame parameter. The computational time is 9 min on a standard workstation using one Intel Core processor.

## 5 The EM2C Case—A Laminar Premixed Flame

The next test case is based on the experiments of Noiray et al. [34]. It is well-suited for the demonstration of our presented method since measurements for stable (Flame Transfer Functions) and unstable combustion (limit cycle frequencies and amplitudes for velocity and pressure fluctuations) are reported. This allows validating the flame model results obtained by the identification process against the measured FTF.

### 5.1 Setup and Operating Point.

The experimental rig consists of: (i) a cylindrical plenum tube, which allows adjusting its length by a movable piston and (ii) an exchangeable perforated plate anchoring a collection of Bunsen flames. A schematic view of the rig is presented in Fig. 6. The testing facility has a varying total length of $0.10\u2009m<L<0.75\u2009m$.

We focus on one of the three perforated plates, which were investigated by Noiray [35]: “Plate 1” has a plate thickness of $3\xd710\u22123\u2009m$ and features 420 perforations with each hole having a radius of $rp=1\xd710\u22123\u2009m$.

A hot wire (denoted as *V _{pl}*) is located 0.02 m upstream of the perforated plate and measures acoustic velocity fluctuations. A microphone (denoted as

*P*) records pressure fluctuations at a distance of 0.35 m downstream the combustion zone.

_{ff}The operating point is defined by a mean flow velocity in the plenum tube $u\xafpl=1.2\u2009m/s$ and an equivalence ratio of $\varphi =0.86$. The total thermal power of the methane/air-fueled combustor is 14 kW. Further relevant quantities are computed using the cantera software with a gri-mech 3.0 mechanism and summarized in Table 2.

Thermal power (kW) | 14 |

Mean flow velocity (m/s) | 1.2 |

Mean hole velocity (m/s) | 3.5 |

Equivalence ratio | 0.86 |

Pressure (Pa) | $1.03\xd7105$ |

Fresh gas temperature (K) | 300 |

Adiab. flame temperature (K) | 2082 |

Fresh gas density (kg/m^{3}) | 1.13 |

Burnt gas density (kg/m^{3}) | 0.16 |

Thermal power (kW) | 14 |

Mean flow velocity (m/s) | 1.2 |

Mean hole velocity (m/s) | 3.5 |

Equivalence ratio | 0.86 |

Pressure (Pa) | $1.03\xd7105$ |

Fresh gas temperature (K) | 300 |

Adiab. flame temperature (K) | 2082 |

Fresh gas density (kg/m^{3}) | 1.13 |

Burnt gas density (kg/m^{3}) | 0.16 |

### 5.2 Experimental Data.

In this section, we briefly summarize the experimental results of Noiray [35] and present the data used for the identification procedure.

Noiray [35] measured the FTF with linear and nonlinear excitation amplitude. In this work, we use the FTF for a nonlinear excitation amplitude of $u\u2032=0.34\u2009m/s$ (Fig. 7), which is close to the velocity fluctuation observed during limit cycle. It was already reported that the acoustic network model predictions vary significantly, if the linear FTF is employed [35].

For the already described operation condition (Table 2), time series for acoustic pressure and velocity fluctuations are shown in Fig. 8. During measurements, plate 1 was installed and the plenum length was adjusted to $L=0.46\u2009m$. We are interested in the mean limit cycle amplitudes, which will serve as input data for the identification procedure. The amplitude of velocity fluctuations during limit cycle upstream the perforated plate is about $u\u2032lc=0.35\u2009m/s$ and the pressure fluctuations measured in the far field are around 15 Pa. Spectral analysis displays a dominant peak at 522 Hz, which corresponds to a 3/4-wave mode of the plenum tube (not shown).

### 5.3 Acoustic Network Model.

The EM2C experiment is reconstructed as an acoustic network model (Fig. 9). The red circles correspond to the experimental positions of the hot wire position upstream the perforated plate and the pressure sensor in the far field (Fig. 6).

The boundary conditions are adopted from the experiment and are similar to Kaess et al. [36]. As Noiray [35] already discussed, the inlet can be considered as a rigid wall, so that *R*_{in} = 1. The outlet reflection coefficient *R _{o}* is frequency dependent and follows the theoretical considerations of [35]: the sound radiated by the collection of flames influences the acoustic impedance at the outlet. We use the model proposed by Noiray [35] which was successfully compared to measurements in his thesis.

### 5.4 Identification of the Flame Model.

Now, the network model is embedded into the identification framework (Sec. 3). The settings for the optimization procedure are identical to the Rijke tube and are provided in Table 1. The target frequency $f\u2009exp\u2009=522\u2009Hz$ together with the physical constraints on acoustic pressure and velocity (Fig. 8) are used for the identification procedure. As for the Rijke tube, one local minimum is found that fulfills the objective function and the physical constraints. The identified flame parameters are *n*_{opt} = 0.69 and $\tau opt=1.25\u2009ms$ and agree well with the experimental values *n _{t}* and

*τ*(Fig. 7). The computational time for the identification procedure is 38 min on a standard workstation using one Intel Core processor.

_{t}## 6 The DLR Case—A High-Pressure, Turbulent Premixed Flame

The DLR high-pressure, premixed, axisymmetric bluff-body stabilized burner is shown in Fig. 10. It consists of five main components; (1) a fuel-air mixing channel, (2) a tubular flow channel, (3) a cylindrical bluff body, equipped with active cooling and a conical tip (to accelerate the reactants and thereby reduce the risk of flashback), (4) an actively cooled adaptor mount and (5) an optically accessible combustion chamber (not shown). The dimensions of the bluff-body and reactant channel are similar to those in the bluff-body burner described in Balachandran et al. [37] with a central bluff-body face of 25 mm diameter, and reactant channel diameter of 35 mm, for a 50% blockage ratio. The central bluff-body contains internal cooling channels designed to protect it from the high thermal loads it experiences at elevated pressure conditions. In this study, compressed air was used as the bluff-body coolant.

This burner was mounted within an optically accessible, high-pressure test rig (DLR-HIPOT), which is capable of testing combustors up to 35 bars of pressure with a maximum thermal load of 300 kW. The test rig supplies up to 200 g/s of air, which can be preheated up to 800 K. It is mounted on a three-axis translation stage, to facilitate translation of the measurement position without adjustment or realignment of the optical diagnostics system.

### 6.1 Setup and Operating Point.

The rectangular combustion chamber has a length of 200 mm and is terminated by a choked nozzle. A schematic view of the rig is presented in Fig. 11. The testing facility has a total length of about 685 mm.

The operating point of the choked combustion chamber is defined by the bulk velocity and equivalence ratio (Table 3). The temperature of the fresh gases has been fixed to $T=300\u2009K$. The equivalence ratio is 0.67 and the mean pressure inside the chamber is 7 bar. Three pressure sensors are installed (red dots in Fig. 11), one in the plenum and two inside the combustion chamber.

Thermal power (kW) | 86.5 |

Mean flow velocity (m/s) | 10.1 |

Equivalence ratio | 0.67 |

Pressure (Pa) | $7.21\xd7105$ |

Fresh gas temperature (K) | 300 |

Adiab. flame temperature (K) | 2049 |

Fresh gas density (kg/m^{3}) | 7.91 |

Burnt gas density (kg/m^{3}) | 1.16 |

Thermal power (kW) | 86.5 |

Mean flow velocity (m/s) | 10.1 |

Equivalence ratio | 0.67 |

Pressure (Pa) | $7.21\xd7105$ |

Fresh gas temperature (K) | 300 |

Adiab. flame temperature (K) | 2049 |

Fresh gas density (kg/m^{3}) | 7.91 |

Burnt gas density (kg/m^{3}) | 1.16 |

### 6.2 Experimental Data of Unstable Combustion Process.

The operating point described above features large pressure oscillations (Fig. 12) in the entire test rig. The measured pressure fluctuations were highest in the plenum and lowest close to the nozzle. The frequency of the limit cycle is at 140 Hz with quasi-constant amplitude and shows no other dominant nonlinear frequency interaction (Fig. 13).

Phase-averaged OH* images were recorded during the instability and were phase-locked with the pressure recordings. With this, the Rayleigh criterion is evaluated (Fig. 14). Pressure and heat release rate fluctuations are in phase (ca. 58 deg), meaning that the energy released by the combustion processes is transferred to the acoustic field of the combustor.

### 6.3 Acoustic Network Model.

The DLR test rig is reconstructed as a network model (Fig. 15). Since the inflow and outflow conditions in the experiments were choked, the acoustic boundary conditions at the inlet *R _{i}* and the outlet

*R*can be safely set to fully reflective. We take into account the two area changes at the flow straightener and bluff body passage. The flame is again represented via two parameter: the flame amplitude

_{o}*n*and time delay

*τ*.

### 6.4 Identification of the Flame Model.

The FTF of the DLR test rig cannot be measured as no flow forcing device is installed. Hence, the flame parameters are identified using the acoustic network model with the optimization strategy (Sec. 3). As for the foregoing test cases, we find a minimum objective function that satisfies the physical constraints. The constraints are given by two of the experimental pressure data (plenum and combustion chamber). For the flame parameters, we find *n*_{opt} = 1.85 and $\tau opt=2.6\u2009ms$, for which the predicted instability frequency is 142.4 + 68j Hz. A good match in instability frequency is found between experiments (Fig. 13) and the identified flame model results. The gain cannot be validated but its large gain may be reasonable in the low-frequency region. The time delay can be valid based on the Rayleigh criterion (Fig. 14). In previous work, we have shown that the Rayleigh criterion can give a good estimate of the time delay *τ* [38]. We use the same approach and retrieve a time delay based on a reconstructed velocity signal $u\u2032$ to be $\tau =2.8\u2009ms$. The time delay based on the reconstructed acoustic velocity is close to the optimized time delay $\tau opt=2.63\u2009ms$. The computational time for the identification procedure is 64 min on a standard workstation using one Intel Core processor.

## 7 Conclusion

This paper presents a data-driven identification framework with the objective to retrieve a flame model from experimental limit cycle measurements. To demonstrate the large applicability of the method, we chose three test cases: (i) a Rijke tube; (ii) an experiment of a laminar flame (EM2C case); and (iii) a high-pressure, turbulent flame (DLR case). The identification strategy retrieved the nonlinear flame model for all cases successfully, thereby satisfying the physical constraints for which experimental data were used. The computational time for all cases ranged from several minutes to maximum 1 h, which highlights the low computational cost of the identification strategy. However, as the complexity of the experiment increases, the number of discretizing acoustic elements increase and cause a longer computational time for the network model. This, in turn, increases the time for the identification strategy. One way to overcome long computational times is to parallelize the optimization framework, e.g., for the global optimization loops.

In a further step, this framework can be used for instability control purposes, e.g., the acoustic network model together with the data-driven flame model access the spatial distribution of the acoustic pressure so that passive control devices can be placed effectively.

This study is based on the most standard optimization framework and can be expanded to, e.g., stochastic frameworks, where the probability and, hence, the robustness of the identified flame parameters can be assessed.

## Funding Data

European Research Council (ERC) under the European Union's Horizon 2020 (Grant No. 682383; Funder ID: 10.13039/501100000781).

Air Force Office of Scientific Research (Grant No. FA9550-16-1-0044; Funder ID: 10.13039/100000181).

## Nomenclature

*c*=speed of sound

*f*,*g*=characteristic waves

*f*_{obj}=objective function

*G*=_{n}global optimization runs

*h*=_{i}constraints

*L*=length

*n*=gain

*P*=probe

- $p\u2032$ =
acoustic pressure fluctuation

- $Q\u2032\u02d9$ =
heat release rate fluctuation

*r*=radius

*R*=reflection coefficient

- $u\u2032$ =
acoustic velocity fluctuation

*V*=volume

*α*=growth rate

*δ*=step width

*ϵ*=optimization tolerance

*ρ*=density

*τ*=time delay

- $\varphi $ =
equivalence ratio

- $\phi $ =
phase

*ω*=angular frequency

## References

**160**(9), pp.