Bayesian inference of physics-based models of acoustically-forced laminar premixed conical flames

Abstract

We perform twenty experiments on an acoustically-forced laminar premixed Bunsen flame and assimilate high-speed footage of the natural emission into a physics-based model containing seven parameters. The experimental rig is a ducted Bunsen flame supplied by a mixture of methane and ethylene. A high-speed camera captures the natural emission of the flame, from which we extract the position of the flame front. We use Bayesian inference to combine this experimental data with our prior knowledge of this flame’s behaviour. This prior knowledge is expressed through (i) a model of the kinematics of a flame front moving through a model of the perturbed velocity field, and (ii) a priori estimates of the parameters of the above model with quantified uncertainties. We find the most probable a posteriori model parameters using Bayesian parameter inference, and quantify their uncertainties using Laplace’s method combined with first-order adjoint methods. This is substantially cheaper than other common Bayesian inference frameworks, such as Markov Chain Monte Carlo. This process results in a quantitatively-accurate physics-based reduced-order model of the acoustically forced Bunsen flame for injection velocities ranging from $1.75\,\mbox{m/s}$ to $2.99\,\mbox{m/s}$ and equivalence ratio values ranging from $1.26$ to $1.47$, using seven parameters. We use this model to evaluate the heat release rate between experimental snapshots, to extrapolate to different experimental conditions, and to calculate the flame transfer function and its uncertainty for all the flames. Since the proposed model relies on only seven parameters, it can be trained with little data and successfully extrapolates beyond the training dataset. Matlab code is provided so that the reader can apply it to assimilate further flame images into the model.

[a]organization=Department of Mechanics, Mathematics and Management, Politecnico di Bari, addressline=via Re David, 200, city=Bari, postcode=70125, country=Italy

[b]organization=Department of Engineering, University of Cambridge, addressline=Trumpington Street, city=Cambridge, postcode=CB2 1PZ, country=UK

Novelty and Significance Statement

Thermoacoustic systems tend to be extremely sensitive to small parameter changes, which makes them difficult to model a priori from existing models in the literature. This means, however, that thermoacoustic models tend to be easy to train using data-driven methods because, with well-chosen experiments, their parameters can be easily observed from experimental data. This paper presents a novel use of Bayesian inference to combine experimental measurements, numerical simulations, and prior knowledge about flame behaviour. We outline our methodology and demonstrate its effectiveness using a laminar premixed Bunsen flame. Our approach yields a quantitatively-accurate physics-based model that predicts the expected value and uncertainty bounds of the flame transfer function between velocity and heat release rate perturbations. The proposed model contains only seven physical parameters, which is fewer parameters than non-physics-based models, and can therefore be trained on relatively little data. We also illustrate how the trained model effectively extrapolates beyond the training dataset. Our numerical code and experimental data are open access.

1 Introduction

The efficiency of the mechanism driving thermoacoustic oscillations depends strongly on the phase difference between the heat release rate (h.r.r) and pressure oscillations [1]. In turn, this phase difference usually depends strongly on the flame’s dynamics [2] and on the acoustics of a combustion chamber, particularly if time delays in the flame are similar to or greater than the period of acoustic chamber modes [3]. Many factors affect the flame dynamics and combustion chamber acoustics, meaning that a combustor’s thermoacoustic behaviour tends to be extremely sensitive to small changes [4]. For the same reason, the outputs of faithful thermoacoustic models are sensitive to small changes in models and their parameters. On the negative side, this makes thermoacoustic systems difficult to model a priori from existing models in the literature. On the positive side, this makes thermoacoustic models easy to train from experimental data because their parameters tend to be easily observable from experimental data. With well-chosen experiments, we can therefore (i) tune the parameters of candidate models and (ii) compare candidate models against each other and select the one with most evidence, given the experimental data [5, 6]. Also on the positive side, this extreme sensitivity means that thermoacoustic systems can often be stabilized by making small changes, which is attractive in industrial settings. The challenge is to model systems accurately and to determine stabilizing modifications as quickly and accurately as possible.

The most influential component of a thermoacoustic system is the mechanism linking acoustic velocity and pressure fluctuations at the base of the flame to subsequent h.r.r. fluctuations in the body of the flame. (In most settings, only the acoustic velocity fluctuations are influential.) This mechanism is difficult to model or simulate [7], so researchers often rely on experimental measurements of the h.r.r. as a function of the acoustic velocity. The fluctuating natural emission from the flame is not, however, a reliable measurement of the fluctuating h.r.r.[8]. Alternatives, such as PLIF to identify reaction zones [9] are more accurate but are technically difficult and, in large systems, impractical. This paper seeks to circumvent the above problems by combining experimental measurements with numerical simulations of a flame model. The flame model is physics-based and qualitatively-accurate. It gives, amongst other things, the flame front dynamics under harmonic forcing. The physics-based model has several parameters, whose expected values and uncertainties are inferred from the flame position, as found from experimental natural emission measurements. We then use the model to calculate the h.r.r., and its uncertainty, as a function of the acoustic velocity perturbations.

In previous work, Juniper & Yoko [5] applied this process to a hot wire Rijke tube. From several candidate models in the literature, they selected those with the most evidence, given the data, and created a quantitatively accurate physics-based model of this thermoacoustic system. In subsequent work on the same system, they used Bayesian optimal experimental design [10] to evaluate the most informative experiments. Following this, they extended the method to infer the flame response of laminar flames from pressure measurements [11]. This study revealed that, in some cases, the uncertainty in the inferred flame responses could be large when the source of information was pressure measurements alone. In this paper we revisit the flames studied in [11] in order to assess an alternative source of information about the flame response: the dynamics of the flame front position. We use a high-speed camera to record the natural emission of the flames in steady conditions and then during a forced limit cycle. We assimilate this data into a physics-based model of a conical Bunsen flame from a previous study [12]. This model contains sufficient parameters to describe the forced experimental conical flame and is differentiable with respect to its parameters. We can therefore obtain the first derivatives of the model outputs with respect to model parameters and use these gradients to (i) find the most probable parameter values, given some data, and (ii) quantify our uncertainty in those parameters and the resulting uncertainty in the model predictions. This process creates a quantitatively-accurate physics-based model of the h.r.r. fluctuations as a function of the velocity perturbation. Unlike the output of many other data-driven approaches, the resulting model is interpretable, trustworthy and extrapolatable.

The paper is structured as follows: in section 2 we explain the experiments and the image processing; in section 3 we report the key feature of the model and how we generate the gradients of the model output with respect to its parameters; in section 4 we explain the Bayesian inference methodology and Laplace’s approximation; in section 5.1 and 5.2 we show how we assimilate the experimental data to create the quantitatively-accurate physics-based model of the flame; in section 5.3 we show how we employ the model to predict the flame transfer function and its uncertainty; in section 5.4 we show how the model is capable of extrapolating outside the training dataset; in section 6 we discuss the conclusions and outline some future work.

2 Experiments

The experimental observations consist of high-speed footage of a conical flame under steady and acoustically-forced conditions. The footage is processed to isolate the position of the flame front.

2.1 Data Collection

The experimental configuration is a laminar premixed conical flame inserted into a vertical duct, as illustrated in figure 1. The lower end of the duct is fixed to a plenum chamber, through which co-flow air is supplied. The upper end is open to the atmosphere. The duct is a 0.8 $\mathrm{m}$ long section of quartz tube with an internal diameter of 75 $\mathrm{mm}$. Quartz is used to allow for optical access to the flame.

Figure 1: Diagram of the experimental rig.

The burner is a 0.85 $\mathrm{m}$ long section of brass tubing with an internal diameter of 14 $\mathrm{mm}$. The outlet of the burner is fitted with a nozzle that is chosen such that the system could become thermoacoustically unstable. At the injection plane, the nozzle diameter is 9.35 $\mathrm{mm}$. The burner is fuelled by a mixture of methane and ethylene over a wide range of equivalence ratios and fuel mass flow rates, which allow us to study flames with a wide range of thermoacoustic responses. The reactants are metered using a set of mass flow controllers, from which they flow through a choke plate, decoupling the supply lines from the acoustic fluctuations in the rig.

We use a high-speed camera to record the flame under both steady and perturbed conditions. If the system is linearly stable, we achieve the perturbed condition by using a loudspeaker to harmonically force the flame near the first resonant frequency of the system, $f=230\mbox{Hz}$. If the system is self-excited, we achieve the steady condition using active control with a phase-shift amplifier.

We study 20 flames, which we select to produce a wide range of flame responses. We parameterize the flames based on (i) the convective time delay, $\tau_{c}=L_{f}/\bar{u}$, which is the time taken for a perturbation travelling at the bulk velocity in the burner tube, $\bar{u}$, to traverse the length of the flame, $L_{f}$, and (ii) the mean heat release rate of the inner cone, $\bar{q}$. We group the 20 flames into five sets of four flames each, and select the composition of each flame such that the time delay is constant within a group and the mean heat release rate varies. The convective time delays range from 9.5 ms to 16 ms, and the mean heat release rates range from 375 W to 600 W. These flames produce a thermoacoustic behaviour ranging from strongly damped, to neutral, to strongly driven. The flame properties are summarized in table 1 in A, and illustrated in figure 2.

2.2 Data Processing

The raw footage is processed in order to (i) isolate the flame front position and (ii) produce a Euclidian distance field with the flame front as the zeroth level set. The images are first averaged to improve the signal-to-noise ratio. For the steady flames, this involves averaging over 200 frames. For the unsteady flames, it involves phase-averaging over 20 frames captured at 10 phase angles. The averaged images are straightened and centred so that the symmetry axis of the flame is vertical and in the centre of the frame. This is necessary to undo any misalignment in the camera setup. To do this, we exploit the symmetry of the flame and find the rigid transformation that maps the image to its mirror. From this, we can determine the transformation that moves the symmetry axis of the flame to be vertical and centred in the frame. Once the flame is centred, we perform an Abel deconvolution to undo the line-of-sight integration of the natural luminosity of the flame. We use an implementation of the ‘onion-peeling’ deconvolution [13]. This process produces a clearly defined flame front but accumulates noise towards the centre line of the image. We perform image segmentation on the Abel deconvolved image using Chan-Vese segmentation [14]. This separates the image signal from the background and allows us to interrogate the connectivity map of each cluster of pixels. We discard weakly connected clusters to denoise the image, leaving just the flame front as the dominant image structure. Finally, we perform a Euclidian distance transform on the binary image of the isolated flame front. This produces a matrix of the same dimensions as the original image, where each element contains the normal distance between the corresponding pixel and the flame front. We use this in our data assimilation process to quantify the discrepancy between the predicted and observed flame front positions.

Figure 2: Processed steady flame images from the 20 flames. Images are artificially coloured according to the approximate convective time delay. Flames with equal mean heat release rates are overlaid.

3 Reduced-Order Model of a laminar premixed Bunsen flame

3.1 Model description

The flame model is detailed in Ref. [12] and can be expressed compactly through:

$$\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),t;\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}),\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}(0)=\mathbf{x}_{0}}},$$ (1)

where $\mathbf{x}(t)=\left(\mathbf{r}(t),\mathbf{z}(t)\right)$ is the state vector defining the flame front position in terms of radial coordinates, $\mathbf{r}$, and longitudinal coordinates, $\mathbf{z}$, at a given time, $t$, and $\mathbf{x_{0}}$ is the initial state. The vector $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ contains the physics-based parameters, while $\mathbf{f}$ is the physics-based nonlinear operator encapsulating the flame front dynamics. The main feature of the model is that it has been designed to be differentiable with respect to the state and the parameters. The heat release rate $q(t)$ can also be expressed in compact form through:

$$q(t)=g(\mathbf{x}(t);\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}),$$ (2)

where $g$ is a nonlinear operator that links the heat release rate $q(t)$ to the flame front position and shape defined by $\mathbf{x}(t)$. The model parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ are: (i) the flame aspect ratio $\beta=\beta(\bar{u},S_{L})$, which depends on the bulk velocity in the burner tube $\bar{u}$ and the unstretched laminar flame speed $S_{L}$; (ii) the nondimensional Markstein length $\mathcal{M}$, which depends on the thermal expansion parameter, the effective Lewis number of the mixture, the Zel’dovich number and the thermal conductivity of the mixture [15]; (iii) the shape parameter $\mu$, which linearly combines a uniform and a parabolic mean velocity profile [16]; (iv) the wavelength of the harmonic perturbation velocity field $K$; (v) the amplitude of the acoustic forcing $\varepsilon_{V}$; (vi) the amplitude of the flame base oscillations $\lambda$; (vii) the initial phase of the flame base oscillations $\phi_{0}$. These seven parameters are sufficient to describe the flame front dynamics qualitatively. All lengths, including the Markstein length, are normalised by the nozzle radius $R=4.675$ $\mathrm{mm}$.

3.2 Obtaining periodic solutions and their sensitivities using adjoint methods

All flames are perturbed by a harmonic acoustic forcing with period $T$. We consider the following vector-valued function:

$$\mathbf{F}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}};\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}})=\mathbf{x}(T)-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}},\quad\text{where }\mathbf{x}(T)=\int_{0}^{T}\mathbf{f}(t,\mathbf{x};\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}})\,\mbox{d}t+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}.$$ (3)

For a given parameter set $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$, a solution $\mathbf{x}(t)$ of equation (1) with an initial state ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}$ is periodic with period $T$ if:

$$\mathbf{F}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}};\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}})=\mathbf{0}.$$ (4)

We solve equation (4) using the trust-region dogleg algorithm and by providing the partial derivative of $\mathbf{F}$ with respect to the initial state ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}$

$$\partialderivative{\mathbf{F}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}=\partialderivative{\mathbf{x}(T)}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}-\mathbf{I},$$ (5)

where $\mathbf{I}$ is the identity matrix, and $\partial{\mathbf{x}(T)}/\partial{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}$ is the partial derivative of the state at time $T$ with respect to changes in the initial state ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}$, which we compute using adjoint methods. We can also obtain the gradient of the periodic solution ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}(t)}}$ with respect to the parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$, by differentiating equation (3) with respect to $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$:

$$\begin{split}\derivative{\mathbf{F}}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}=&\partialderivative{\mathbf{x}(T)}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}+\partialderivative{\mathbf{x}(T)}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}\derivative{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}-\derivative{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}=\mathbf{0}\\ \implies&\derivative{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}=-\left(\partialderivative{\mathbf{x}(T)}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}-\mathbf{I}\right)^{-1}\partialderivative{\mathbf{x}(T)}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}},\end{split}$$ (6)

and then by differentiating:

$$\derivative{\mathbf{x}(t)}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}=\partialderivative{\mathbf{x}(t)}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}+\partialderivative{\mathbf{x}(t)}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}\derivative{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbf{x}_{0}}}}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}.$$ (7)

Therefore we compute the gradients of the h.r.r. with respect to the parameters by differentiating (2):

$$\derivative{q(t)}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}=\partialderivative{g(\mathbf{x}(t);\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}})}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}+\partialderivative{g(\mathbf{x}(t);\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}})}{\mathbf{x}(t)}\derivative{\mathbf{x}(t)}{\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}}.$$ (8)

This quantity is essential for quick calculation of the uncertainties of the heat release rate and therefore the flame transfer function.

4 Bayesian data assimilation

Bayesian data assimilation provides us with a set of tools that allow us to (i) infer the most probable parameters of a physics-based model, given some data, (ii) quantify the uncertainty in the model parameters, and estimate the systematic uncertainty in the model and data, (iii) rank several candidate models to select the best model, given the data, and (iv) identify optimal experiment designs and sensor placements. In this paper, we use points (i) and (ii), so these are described in detail. A demonstration of model ranking can be found in [11], and a demonstration of optimal experiment design can be found in [10].

In this section, we introduce the following notation to outline the Bayesian inference methodology: the data $\mathbf{D}$ contains experimental observations of the flame position in steady and/or unsteady configurations; the model $\mathcal{H}$ encodes a physics-based reduced-order model, such that for a given set of parameters, $\mathbf{a}$, the model $\mathcal{H}(\mathbf{a})$ gives a prediction of the data $\mathbf{D}$.

4.1 Parameter inference

When performing parameter inference, we assume that the model is structurally correct and we infer its most probable parameters, $\mathbf{a}_{\mathrm{MP}}$. The model is rarely free of structural error, however, and we will revisit this assumption later. We encode our level of uncertainty in the parameter values through a probability distribution, which we denote $P(\bullet)$. Using any prior knowledge we have about the unknown parameters (which may be none at all), we propose a prior probability distribution over the parameter values, $P(\mathbf{a}|\mathcal{H})$. We then assimilate the data, $\mathbf{D}$, by performing a Bayesian update on the parameter values:

$$P(\mathbf{a}|\mathbf{D},\mathcal{H})=\frac{P(\mathbf{D}|\mathbf{a},\mathcal{H})P(\mathbf{a}|\mathcal{H})}{P(\mathbf{D}|\mathcal{H})}$$ (9)

The quantity on the left-hand side of equation (9) is the posterior probability of the parameters, given the data. It is generally computationally intractable to calculate the full posterior, because this requires integration over a large parameter space. The integral typically cannot be evaluated analytically, and requires thousands of model evaluations to compute numerically. At the parameter inference stage, however, we are only interested in finding the most probable parameters, which are those that maximize the posterior. We therefore use an optimization algorithm to find the peak of the posterior without evaluating the full distribution. This process is made computationally efficient by (i) assuming that the experimental uncertainty is Gaussian distributed, and (ii) choosing the prior parameter distribution to be Gaussian. Assumption (i) is reasonable for well-designed experiments in which the uncertainty is dominated by random error, which is typically Gaussian distributed. For assumption (ii) we note that the choice of prior is often the prerogative of the researcher, and we are free to exploit the mathematical convenience offered by the Gaussian distribution. The correlations between model parameters are rarely known a-priori, so an independent Gaussian distribution is often used for the prior, as in this paper. When finding the most probable parameters, we neglect the denominator of the right-hand side of equation (9), because it does not depend on the parameters. It is then convenient to define a cost function, $\mathcal{J}$, to be the negative log of the numerator of equation (9), which we minimize. In this study, the data, $\mathbf{D}$, consists of a set of experimental observations of the flame position at various time frames and experimental configurations. In its most general form, the cost function for a single image of a flame is:

$$\begin{split}\mathcal{J}(\mathbf{a})=&-\log{\{P(\mathbf{D}|\mathbf{a},\mathcal{H})P(\mathbf{a}|\mathcal{H})}\}\\ =&\frac{1}{2}(\mathbf{s}(\mathbf{a})-\mathbf{z})^{T}\mathbf{C}_{ee}^{-1}(\mathbf{s}(\mathbf{a})-\mathbf{z})\\ &+\frac{1}{2}(\mathbf{\mathbf{a}-\mathbf{a_{p}}})^{T}\mathbf{C}_{aa}^{-1}(\mathbf{\mathbf{a}-\mathbf{a_{p}}})+\text{const.}\end{split}$$ (10)

where: $\mathbf{s}(\mathbf{a})$ and $\mathbf{z}$ are column vectors of the predicted and measured flame front positions respectively; $\mathbf{C}_{ee}$ is the covariance matrix describing the experimental uncertainty; $\mathbf{a}$ and $\mathbf{a}_{p}$ are column vectors of the current and prior parameter values respectively; $\mathbf{C}_{aa}$ is the covariance matrix describing the uncertainty in the prior, and const. is a constant that arises from the Gaussian pre-exponential factors, which does not impact the maximum posterior parameter estimation $\mathbf{a}_{\mathrm{MP}}$. In this study, we evaluate the discrepancy between the model predictions and the experimental observations by projecting the predicted flame front position onto the Euclidian distance field produced from the corresponding frame of the high speed footage. This produces a column vector in which each entry is the normal distance between the predicted and observed flame front position at each point in the model.

We see from equation (10) that assuming Gaussian distributions for the data and prior reduces the task of parameter inference to a quadratic optimization problem. We solve this optimization problem using the BFGS algorithm, where the gradient information is provided using adjoint methods.

The parameter inference process is illustrated in figure 3 for a simple system with a single parameter, $a$, and a single observable variable, $z$. In (a) we show the marginal probability distributions of the prior, $p(a)$, and the data, $p(z)$. The prior and data are independent, so we construct the joint distribution, $p(a,z)$ by multiplying the two marginals. In (b), we overlay the model predictions, $s$, for various values of $a$. Marginalizing along the model predictions yields the true posterior, $p(a|z)$. This is possible for a cheap model with a single parameter, but exact marginalization quickly becomes intractable as the number of parameters increases. In (c) we plot the cost function, $\mathcal{J}$, which is the negative log of the unnormalized posterior. We show the three steps of gradient-based optimization that are required to find the local minimum, which corresponds to the most probable parameters, $a_{\mathrm{MP}}$.

Figure 3: Illustration of parameter inference on a simple univariate system. (a) the marginal probability distributions of the prior and data, $p(a)$ and $p(z)$, as well as their joint distribution, $p(a,z)$ are plotted on axes of parameter value, $a$, vs observation outcome, $z$. (b) the model, $\mathcal{H}$, imposes a functional relationship between the parameters, $a$, and the predictions, $s$. Marginalizing along the model predictions yields the true posterior, $p(a|z)$. This is computationally intractable for even moderately large parameter spaces. (c) instead of evaluating the full posterior, we use gradient-based optimization to find its peak. This yields the most probable parameters, $a_{\mathrm{MP}}$

4.2 Uncertainty quantification

Uncertainty quantification can be split into two steps: (i) quantifying the parametric uncertainty and propagating it to the model prediction uncertainty, and (ii) estimating the systematic and structural uncertainty in the experiments and model predictions. We will deal with these separately.

4.2.1 Parameteric uncertainty

Once we have found the most probable parameters, $\mathbf{a}_{\mathrm{MP}}$, we estimate the uncertainty in these parameter values using Laplace’s method [17, 18, 5]. This method approximates the posterior as a Gaussian distribution with a mean of $\mathbf{a_{\mathrm{MP}}}$, and a covariance given by the Hessian of the cost function:

$$\begin{split}\mathbf{C}_{aa}^{\mathrm{MP}-1}&\approx\frac{\partial^{2}\mathcal{J}}{\partial a_{i}\partial a_{j}}\end{split}$$ (11)

The Hessian of the cost function is estimated by the BFGS algorithm during the minimization process. This estimate is based on the change in the value and gradient of the cost function over the previous iteration.

The accuracy of the Laplace approximation depends on the functional dependence between the model predictions and the parameters. This is shown graphically in figure 4 for three univariate systems. In (a), the model is linear in the parameters. Marginalizing a Gaussian joint distribution along any intersecting line produces a Gaussian posterior, so Laplace’s method is exact. In (b), the model is weakly nonlinear in the parameters. The true posterior is skewed, but the Gaussian approximation is still reasonable. This panel also shows a geometric interpretation of Laplace’s method: the approximate posterior is given by linearizing the model around $\mathbf{a}_{\mathrm{MP}}$, and marginalizing the joint distribution along the linearized model. In (c), the model is strongly nonlinear in the parameters, so the true posterior is multi-modal and the main peak is highly skewed. In this case, the gradient-based optimization algorithm will only find a single local minimum, which will depend on the choice of initial condition for the optimization. Furthermore, we see that the posterior estimated using Laplace’s method is a poor approximation of the true posterior. This problem can be avoided by reducing the extent of the nonlinearity captured by the joint distribution by (i) shrinking the joint distribution by providing more precise prior information or more precise experimental data, or (ii) re-parameterizing the model to reduce the strength of the nonlinearity [18, Chapter 27].

Figure 4: Illustration of uncertainty quantification for three univariate systems, comparing the true posterior, $p(a|z)$ to the approximate posterior from Laplace’s method, $p(a|z)_{L}$. (a) the model is linear in the parameters, so the true posterior is Gaussian and Laplace’s method is exact. (b) the model is weakly nonlinear in the parameters, the true posterior is slightly skewed, but Laplace’s method yields a reasonable approximation. (c) the model is strongly nonlinear in the parameters, the posterior is multi-modal and Laplace’s method underestimates the uncertainty

In many cases, the use of this approximate inference framework approximation will be justifiable, given the substantial reduction in computational cost compared to sampling methods, which are the only viable alternative for constructing the posterior. In previous work, we compared the computational cost of our framework to two sampling approaches [11]. The comparison was done on a computationally cheap thermoacoustic network model. Applying our framework to this model, we can compute the posterior probability of five unknown parameters in under 5 seconds on a single core on a laptop. The same inference problem takes 35 CPU hours running on a workstation when solved with Markov Chain Monte Carlo, and 22 CPU hours when solved with Hamiltonian Monte Carlo.

4.2.2 Uncertainty propagation

To quantify the uncertainty in the model predictions due to the parametric uncertainty, we linearise the model around $\mathbf{a}_{\mathrm{MP}}$ and propagate the parameter uncertainties through the linearised model [19]. The uncertainty in the model prediction is given by:

$$\mathbf{C}_{xx}=\mathbf{J}^{T}\mathbf{C}_{aa}\mathbf{J}$$ (12)

where $\mathbf{C}_{xx}$ represents the covariance in the model predictions, and $\mathbf{J}$ denotes the Jacobian matrix of the model predictions with respect to the parameters $\mathbf{a}$. The marginal uncertainty in each model prediction, $(\sigma_{x_{i}})^{2}$, is given by the diagonal entries of $\mathbf{C}_{xx}$.

4.2.3 Systematic uncertainty

In most cases, experimental data will contain some systematic uncertainty, and models will contain some structural uncertainty. These uncertainty sources cannot be quantified a-priori, and are often referred to as “unknown unknowns”. We can, however, construct a total covariance matrix, $\mathbf{C}_{tt}$, which encodes the total uncertainty due to (i) the known experimental uncertainty, (ii) the unknown systematic experimental uncertainty, and (iii) the unknown structural model uncertainty. We can then estimate this total covariance from the posterior discrepancy between the model and the data. This must be done simultaneously with parameter inference, because the posterior parameter distribution depends on the total uncertainty in the model and data. We therefore replace $\mathbf{C}_{ee}$ with $\mathbf{C}_{tt}$ in equation (10), and estimate the total uncertainty by simultaneously minimizing $\mathcal{J}$ with respect to $\mathbf{a}$ and $\mathbf{C}_{tt}^{-1}$.

We begin by calculating the derivative of $\mathcal{J}$ with respect to $\mathbf{C}_{tt}^{-1}$, assuming that the observed variables are uncorrelated, and keeping in mind that the normalizing constant, $K$, depends on $\mathbf{C}_{tt}$:

$$\mathcal{J}=\frac{1}{2}(\mathbf{s}(\mathbf{a})-\mathbf{z})^{T}\mathbf{C}_{tt}^{-1}(\mathbf{s}(\mathbf{a})-\mathbf{z}) +\log\left(\sqrt{(2\pi)^{k}|\mathbf{C}_{tt}|}\right)$$ (13)

$$+\frac{1}{2}(\mathbf{a}-\mathbf{a_{p}})^{T}\mathbf{C}_{aa}^{-1}(\mathbf{a}-\mathbf{a_{p}}) +\log\left(\sqrt{(2\pi)^{k}|\mathbf{C}_{aa}|}\right)$$

$${\frac{\partial\mathcal{J}}{\partial\mathbf{C}_{tt}^{-1}}=\frac{1}{2}(\mathbf{s}(\mathbf{a})-\mathbf{z})(\mathbf{s}(\mathbf{a})-\mathbf{z})^{T}\circ\mathbf{I}-\frac{1}{2}\mathbf{C}_{tt}}$$ (14)

where $\mathbf{I}$ is the identity matrix, and $\circ$ denotes the Hadamard product. For a given set of parameters, the most probable $\mathbf{C}_{tt}$ sets equation (14) to zero. This gives the estimate:

$${\mathbf{C}_{tt}=(\mathbf{s}(\mathbf{a})-\mathbf{z})(\mathbf{s}(\mathbf{a})-\mathbf{z})^{T}\circ\mathbf{I}}$$ (15)

which is the expected result that the total variance in the model and data is the square of the discrepancy between the model predictions and the data. Although we cannot directly identify the source of the unknown uncertainty because the experimental and model uncertainties cannot be disentangled, the inferred total uncertainty can assist the researcher with identifying potential error sources. For example, if the unknown error in a single sensor is unexpectedly large, this could indicate a faulty sensor or bad installation. If the unknown error at a certain experimental operating condition is large, this could prompt the researcher to repeat that experiment. If the unknown error grows with one of the input variables, the researcher might investigate the model to see if any important physical phenomena has been neglected.

5 Creating the quantitatively-accurate physics-based model

In this study, we perform Bayesian data assimilation in two steps. In the first step, we assimilate the flame front position data of each flame individually. Therefore, for each experimental configuration (i.e. for each combination of mean flow velocity $\bar{u}$ and equivalence ratio $\phi$), we determine the most probable set of model parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=(\mathcal{M},\mu,\beta,K,\varepsilon_{V},\lambda,\phi_{0})$, which have been introduced in section 3. This process does not create a general model that can predict the behaviour of an arbitrary flame, but is useful for assessing the influence of the model parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ on the model predictions. After completion of this step, we need to introduce a set of hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$, such that $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}(\phi,\bar{u};\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})$. We then simultaneously assimilate the flame front position data of all flames. This allows us to determine the most probable set of hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$ so that, for each combination of mean flow velocity $\bar{u}$ and equivalence ratio $\phi$, the model $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}(\phi,\bar{u};\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})$ provides the best prediction of the physics-based model parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ and, consequently, of the flame front dynamics. In figure 5, we show a diagram describing the methodology followed in this study. The end goal is to create a general model that provides the flame transfer function, along with quantified uncertainty bounds, for any given pair of equivalence ratio and bulk velocity in the burner tube.

Figure 5: Illustration showing the methodology for constructing a quantitatively-accurate reduced-order model. The general model depends on the hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$ and comprises two components: the physics-based reduced order model, which depends on the physical parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$, and a linking model $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}(\phi,\bar{u};\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})$ that establishes the relationship between these physical parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ and the hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$. The model is fully differentiable and provides the flame front dynamics and the flame transfer function. By comparing the model’s predictions of the flame front with experimental images, we refine our hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$ through Bayesian updates. The resulting model predicts the flame transfer function, with confidence bounds, for any combination of equivalence ratio $\phi$ and the mean flow velocity $\bar{u}$. This model is both interpretable and capable of extrapolation.

5.1 Step 1: assimilation of each flame individually

First, we assimilate the steady-state snapshots in order to infer $\mathcal{M}$, $\beta$ and $\mu$. Second, we use these inferred values as priors for the second step, in which we assimilate the unsteady flame images and infer all parameters simultaneously: $\mathcal{M}$, $\beta$, $\mu$, $K$, $\varepsilon_{V}$, $\lambda$, $\phi_{0}$. Here we show the process for flame 17 in Table 1, which has an equivalence ratio $\phi=1.34$, a flow rate of $\bar{u}=2.86\,\mbox{m/s}$ and an equal volume fraction of methane $(\mathrm{CH_{4}})$ and ethylene $(\mathrm{C_{2}H_{4}})$.

5.1.1 Steady-state data assimilation

The steady flame front position depends on $\beta$, $\mathcal{M}$ and $\mu$. We can include prior knowledge about the parameter $\beta=\beta(\bar{u},S_{L})$ because it depends on the unstretched flame speed $S_{L}$, which can be estimated using the open source chemical kinetics code Cantera [20], and the mean flow velocity $\bar{u}$, which can be measured during the experiments. Furthermore, we can estimate the value of the Markstein length by using Matalon’s formulation [15]. Following this process for flame $17$, we assign a prior estimate to the unstretched laminar flame speed of $S_{L}=0.4036\,\mbox{ms}^{-1}$ and a dimensional Markstein length of $\tilde{\mathcal{M}}=0.2579\,\mbox{mm}$, corresponding to $\mathcal{M}=\tilde{\mathcal{M}}/R=0.0552$. To construct the parameter covariance matrix, we must assign confidence bounds to these values, which is naturally a subjective process. In this case, we assign an uncertainty of $6\sigma_{S_{L}}=10\%S_{L}$ and $6\sigma_{\mathcal{M}}=100\%\mathcal{M}$ to indicate that we are relatively more confident in our ability to predict laminar flame speed than Markstein length. The only prior knowledge available for the velocity shape parameter $\mu$ is that it can vary between $0$ and $1$, so we set a prior with large uncertainty for this parameter. From the steady-state flame snapshot, we can infer the maximum a posteriori probability estimate of the aspect ratio, $\beta$, the nondimensional Markstein length, $\mathcal{M}$, and the velocity shape parameter, $\mu$. Figure 6(a-b) shows the prior and posterior probability distributions of the model parameters. The contours show $1,2$ and $3$ standard deviations from the prior $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{p}$ (blue) and posterior $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{\mathrm{MP}}$ (red). We see that the experimental observations have reduced the uncertainties in the model parameters. Figure 6(c) shows the model prediction $\pm 3$ standard deviations overlaid on the steady flame image. The largest uncertainties are found at the flame tip, which is the region that is most sensitive to the parameters. The model prediction uncertainty is found by solving equation (12).

Figure 6: (a-b) Probability distributions of parameters in the steady case: prior (blue) and posterior (red). Contours denote 1, 2, and 3 standard deviations from the prior $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{p}$ and posterior $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{\mathrm{MP}}$ estimations. Parameters include the flame aspect ratio $\beta$, which is a function of the unstretched laminar flame speed $S_{L}$ and the bulk velocity in the burner tube $\bar{u}$, the shape parameter for the velocity field $\mu$, and the nondimensional Markstein length $\mathcal{M}$. Both panels represent the same quantities but panel (a) is scaled to the prior distributions, while panel (b) is scaled to the posterior distributions. (c) Comparison between model predictions before (red line) and after (blue line) assimilating the experimental steady-state flame image. Confidence intervals, spanning 3 standard deviations, are illustrated by the red and blue shadings overlaying the steady flame image.

5.1.2 Unsteady-state (acoustically forced)

We repeat the process with the images of the acoustically forced flame. In this case, we set the values of the three steady parameters using the information gained from the previous step, and we set large uncertainties in the prior values of $K$, $\varepsilon_{V}$, $\lambda$ and $\phi_{0}$. Figure 7(a-b) shows the prior and posterior distribution of the parameters before and after assimilating the unsteady flame images. The results show that the uncertainty in the values of $K$, $\varepsilon_{V}$, $\lambda$ and $\phi_{0}$ has been reduced significantly through the assimilation of this data, and the values of $\mathcal{M}$ and $\mu$ have shifted slightly from the values found during the assimilation of the steady data, without exceeding the 3-standard deviation bounds. Figure 7(c) shows the model prediction $\pm 3$ standard deviations, plotted against the unsteady flame images before (blue) and after (red) assimilating the experimental video footage. As in the steady-state case, the largest uncertainties are found at the flame tip. We see that the model matches the data well at each frame. With this, we have produced a digital twin of the flame, which we can use to estimate additional quantities of interest such as the unstretched laminar flame speed $S_{L}$, the Markstein length $\mathcal{M}$, the convective speed of velocity perturbations $K$ and the heat release rate, which were not directly measured in the experiments. Furthermore, because the model is probabilistic, we can also quantify our uncertainty in these estimates.

Figure 7: (a-b) Probability distributions of parameters for the acoustically forced case: prior (blue) and posterior (red). Contours depict 1, 2, and 3 standard deviations from the prior $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{p}$ and posterior $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{\mathrm{MP}}$ estimations. Both panels represent the same quantities; however, panel (a) is scaled to the prior distributions, while panel (b) is scaled to the posterior distributions. Parameters include the flame aspect ratio $\beta$, which is a function of the unstretched laminar flame speed $S_{L}$ and the bulk velocity in the burner tube $\bar{u}$, the shape parameter for the velocity field $\mu$, the nondimensional Markstein length $\mathcal{M}$, the nondimensional velocity perturbation convective speed $K$, and the amplitude of velocity perturbations $\varepsilon_{V}$. (c) Comparison of model predictions before (red line) and after (blue line) assimilating the experimental video footage. The confidence interval, spanning 3 standard deviations, is indicated by red and blue shadings overlaid on the experimental flame images.

For the example of flame 17, the unstretched laminar flame speed and the Markstein length are estimated to be $S_{L}=0.4191\,\mbox{m/s}$ and $\mathcal{M}=0.2430\,\mbox{mm}$ compared to the value of $S^{c}_{L}=0.4036\,\mbox{m/s}$ and $\mathcal{M}^{c}=0.2579\,\mbox{mm}$ obtained from Cantera and the value of the convective speed is estimated to be $K=0.8874$ compared to $K=1$ in Schuller et al. [21], $K=0.9$ in Kashinath et al. [22] and $K=0.83$ in Orchini & Juniper [23]. In this case, we were able to provide reasonable prior information for $\mathcal{M}$ and $S_{L}$. Similar values are, however, inferred even if the prior is inaccurate and uncertain. This methodology could therefore be used to cheaply and easily estimate combustion properties of fuels for which data are not available.

The proposed model contains few parameters, allowing it to be trained on a relatively small amount of data. Furthermore, the data required to train the model is easy to collect from a simple experiment because only snapshots of the natural emission of the flame are required.

Once the most probable parameter values are inferred and the flame front position has been predicted, we can estimate the heat release rate using equation (2). The resulting heat release rate estimation is shown in figure 8. Panel (a) shows the model prediction (red lines) with a confidence of 3 standard deviations (red shading) overlaid on the experimental observations (grayscale). We also plot the predictions between two observations (blue lines and shading). Panel (b) shows the h.r.r. perturbation $q^{\prime}$ normalised by its mean value $\bar{q}$, computed using the model. The black dots correspond to the time frames observed in the experiments. Figure 8 highlights the advantages of the Bayesian inference approach: using experimental flame images we know the flame front position at ten different time points during a period of oscillation; we use this data to infer the seven parameters of a physics-based reduced order model, which gives the flame front position and the heat release rate at any point in time with their uncertainties.

Figure 8: Diagram showing the concept of using Bayesian inference to create a quantitatively accurate model capable of providing physical information not directly available from experiments. The experimental data consists of ten frames during one period of oscillation. Panel (a) shows the model prediction (red line) with a confidence interval of 3 standard deviations (red shading) plotted on top of the experimental flame images, as in 7. We also use the model to infer the flame position (blue lines) and its uncertainty (blue shading) at moments between observations. Panel (b) shows the h.r.r. perturbation normalised by its mean value, computed using the model. The black dots correspond to the timeframes observed in the experiments

5.2 Step 2: assimilation of all flames into a general model

At this stage, we introduce a linking model $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}(\phi,\bar{u};\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})$, designed to determine the model parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ based on the equivalence ratio $\phi$ and bulk flow rate in the burner tube $\bar{u}$, given any specified set of hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$. Figure 9 shows the maximum a posteriori probability distribution of the model parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}_{\mathrm{MP}}$ assimilated in the first step, denoted by blue markers with error bars representing 3 standard deviations. Panel (a) shows the assimilated unstretched laminar flame speed ($S_{L}$) against the equivalence ratio $\phi$. Panel (b) shows the assimilated Markstein length ($\mathcal{M}$) against the equivalence ratio $\phi$. Panel (c) shows the assimilated shape parameter ($\mu$) against the bulk velocity ($\bar{u}$) in the burner tube. The parameter $\mu$ decreases linearly with increasing $\bar{u}$. Panel (d) shows the assimilated convective speed ratio ($K$) against the acoustic-forcing Strouhal number St. The plot shows that $K$ falls within the range of $0.8$ to $1$, with no discernible pattern with respect to St. The remaining parameters $\varepsilon_{V}$, $\lambda$, and $\phi_{0}$ (not shown) depend on experimental conditions and vary according to external acoustic forcing, both of which are independent of the flame properties.

Based on these observations, the model $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}(\phi,\bar{u};\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})$ is formulated as follows: (i) the laminar flame speed $S_{L}=S_{L}^{c}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{1}\phi+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{2}$ is equal to the laminar flame speed computed using Cantera $S_{L}^{c}$ plus a linear correction that depends on the equivalence ratio; (ii) the Markstein length $\mathcal{M}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{3}\mathcal{M}^{c}$ is equal to the Markstein length computed using Cantera multiplied by the parameter ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{3}$; (iii) the shape parameter $\mu={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{4}\bar{u}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{5}$ is assumed to vary linearly with the bulk flow velocity and therefore with the Reynolds number of the bulk flow in the burner tube; (iv) the value $K={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}_{6}$ is assumed to be constant across all flame cases. We keep $\varepsilon_{V},\lambda$ and $\phi_{0}$ equal to the values obtained in the previous assimilation step because the amplitudes of the acoustic forcing and of the flame base oscillations depend on the specific experimental conditions and are not general properties of the flame.

We then assimilate the most probable hyperparameter set, denoted as $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}_{\mathrm{MP}}$, by minimizing the cost function expressed by equation (10) for all flames at once. The results of this step are shown in red in figure 9. The red lines display the estimates of the physical parameters $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}$ as functions of the equivalence ratio $\phi$, the mean flow velocity $\bar{u}$, and the Strouhal number St, after the assimilation of the hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$. The three red contours show 1, 2, and 3 standard deviations from the posterior values.

Figure 10 displays the predicted flame front dynamics for all 20 flames in the library at a single timestep. The experimental observations are shown as grayscale images, and the model predictions are overlaid as red lines with confidence intervals of 1, 2, and 3 standard deviations (depicted as red shading). The numbers at the base of the flames indicate the corresponding flame numbers referenced in Table 1. In the supplementary material, we provide the animated version of this figure, showing the predictions for all timesteps of the experimental video footage.

Figure 9: Results of the second step of the data assimilation process. From the first step of the process, we obtain the trends of the physical parameters with respect to the equivalence ratio $\phi$ and the bulk velocity $\bar{u}$, which is proportional to the Reynolds number, Re, inside the burner tube. Then we propose a general model for the parameters and assimilate all the flame conditions together. The blue markers represent the posterior values from the first step of the data assimilation process and the error bars correspond to 3 standard deviations. The red dotted lines correspond to the posterior values from the second step and the shaded areas represent 1, 2 and 3 standard deviations. Panel (a) shows the assimilated values of the unstretched laminar flame speed $S_{L}$ plotted against the equivalence ratio $\phi$. Panel (b) shows the assimilated values of the dimensional Markstein length $\tilde{\mathcal{M}}$ plotted against the equivalence ratio $\phi$. Panel (c) shows the assimilated values of the velocity shape parameter $\mu$ plotted against the bulk velocity in the burner tube $\bar{u}$. Panel (d) shows the assimilated values of the convective speed ratio $K$ plotted against the acoustic perturbation Strouhal number St.

Figure 10: The general model predictions for all flames in the library are represented by the red line, with confidence intervals of 3 standard deviations (depicted as red shading), overlaid on the experimental flame images. The numbers at the base of the flames indicate the corresponding flame numbers referenced in Table 1. Additionally, in the supplementary material, we provide an animated version of this figure, showing predictions at all time steps of the experimental video footage.

5.3 Flame transfer function predictions

We introduce the flame transfer function $\mathcal{F}$ as:

$$\mathcal{F}(\mbox{St};\phi,\bar{u},\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})=\frac{\hat{q}/\bar{q}}{\hat{u}_{ref}/\bar{u}},$$ (16)

where $\hat{\star}$ denotes the Fourier component of $\star$ at the Strouhal number St; $q$ denotes the fluctuating heat release rate and $u_{ref}$ denotes the velocity fluctuations at a reference position. In this study, we define the reference velocity $u_{ref}$ as the component of the velocity field perpendicular to the flame front at its base, adhering to the definition provided by Orchini & Juniper [24]. The heat release rate $q$ is normalised by the mean heat release rate, while the velocity $u$ is normalised by the mean flow velocity in the burner tube $\bar{u}$. We use the calibrated general model obtained in the second step to predict the gain and phase of the flame transfer function at the excitation Strouhal number St for all flames in Table 1.

In figure 11(a), we present the model predictions along with their uncertainties in a polar plot. Each flame transfer function is coloured according to its corresponding flame group as in figure 2. It is worth noting that the flames within the same group, which are grouped based on convective time delay $\tau$, exhibit similar phase delays on the polar plot, as one would expect. We extend these findings by assuming that the hyperparameters $\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}}$ remain constant across different forcing frequencies. Figure 11(b) illustrates the predictions of flame transfer functions along with their uncertainties over a wide range of Strouhal numbers (St). In future research, we will evaluate the validity of this assumption by assimilating data collected over a range of Strouhal numbers. Should it prove untenable, we will enhance the complexity of the linking model $\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}=\mathbf{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{a}}}(\phi,\bar{u};\boldsymbol{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\alpha}}})$ to incorporate the impact of the forcing frequency, by following the same methodology presented in this paper.

Figure 11: Inferred flame transfer functions presented as (a) polar plots, which show gain on the radial axis and phase on the angular axis, and (b) Bode plots, which show (i) gain and (ii) phase against Strouhal number St. In (a) we plot flame transfer functions as shaded patches, which denote 1-3 standard deviations around the expected value. In (b) we plot the flame transfer functions as solid lines with a shaded patch representing the confidence bounds of 3 standard deviations, whereas the coloured dots represent the values corresponding to the experimental forcing frequency.

5.4 Extrapolation

The model in this paper contains a small number of parameters because it is physics-based. Consequently, the training process requires only a small amount of data. In this section, we show the model’s ability to predict the flame front dynamics of the 20 flames in our library (see Table 1) using data from just four flames. We train the model’s hyperparameters on flames 1, 6, 11, and 16 of Table 1 and predict the remaining flames with this model. These flames have the same convective time delay, but different powers. We therefore test the ability of the model to extrapolate outside the range of trained convective time delays. Figure 12 shows the predicted flame front dynamics for all flames at a single timestep. Red lines show the flames used for training. Blue lines show the flames used for testing. Red and blue shading denotes 1, 2, and 3 standard deviations from the predicted flame front position. All model predictions are overlaid onto the experimental flame images for comparison. The number at the base of each flames indicates the flame numbers in Table 1. In the supplementary material, we provide the animated version of this figure, showing the predictions for all timesteps of the experimental video footage. We now compare figure 12 with figure 10. When only one-fifth of the available data has been used, the predicted flame front positions are less accurate and have larger uncertainty. This is particularly clear for flame 7, which has a large acoustic amplitude and is close to pinch off in these figures. Indeed, the model shows significant sensitivity to the model parameters as the acoustic amplitude increases. Similarly, the prediction for flame 20 is less accurate because the flame length is much larger than the length of the four flames used for model training. Even with such a small amount of training data, however, the data still falls within the uncertainty bounds of the model predictions, even when the data is well outside the training range. We further note that figure 12 shows the timestep with the largest model misfit and uncertainties. In the supplementary material, we provide an animated version of this figure showing the predictions for all timesteps of the video footage.

Figure 12: Predicted flame front dynamics for all flames in the library at a single timestep. Red lines show flames used for training, while blue lines show flames used for testing. The shaded regions in red and blue indicate 1, 2, and 3 standard deviations from the optimal flame front estimation. Additionally, all model predictions are overlaid onto the experimental flame images for comparison. In the supplementary material, we provide an animated version of this figure, presenting predictions for all time steps of the experimental video footage.

6 Conclusions

In this study, we perform experiments on an acoustically forced, laminar premixed conical flame in a duct. We use a high-speed camera to record snapshots of the natural emission of the flame while steady and forced. We propose a physics-based reduced-order model that provides, amongst other things, the flame front dynamics under harmonic forcing. We combine the model output with the experimental images of the flame to infer the most probable model parameters. This process (i) turns a qualitatively-accurate model into a quantitatively-accurate model, and (ii) quantifies the uncertainty in the inferred model parameters and the model predictions.

The inference process produces a digital twin of the flame, which provides access to quantities that were not directly measured in the experiments. From observations of the steady flame, we can estimate combustion characteristics such as the laminar flame speed and the Markstein length. This can be used to characterise the combustion properties of a fuel for which data are not available. From observations of the perturbed flame, we can infer the fluctuating heat release rate as a response to the velocity perturbation, which is then used to calculate the flame’s thermoacoustic response through, for example, the flame transfer function. The proposed model contains few parameters and can therefore be trained on relatively little data. To demonstrate this we infer the model’s parameters using one-fifth of the available data and successfully predict the remaining flames in the dataset. The model is also differentiable with respect to the parameters meaning that it can be used as a design optimization tool.

In future work, the flame model will be included within a larger model of the thermoacoustic behaviour of the system, whose other parameters are tuned in the same way. The final result will be a quantitatively-accurate physics-based model of the flame and the thermoacoustic system that is interpretable, trustworthy, and extrapolatable.

The Matlab code developed in this study is made available to the reader, who is invited to modify and apply it to assimilate other flame image data into the flame model.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the research reported in this paper.

Data availability

The experimental data is openly available at http://doi.org/10.17863/CAM.106960. The Matlab code is available at https://github.com/aleGiannotta/Bayesian_Inference_Lam_Flame.git

Appendix A Properties of the flames studied in this paper

Table 1: Summary of the properties of the 20 flames studied. We show the average measured flow rates of air, methane ($\mathrm{CH_{4}}$) and ethylene ($\mathrm{C_{2}H_{4}}$), the equivalence ratio ($\phi$), the bulk velocity in the burner tube ($\bar{u}$), the measured flame length ($L_{f}$), and the mean heat release rate ($\bar{Q}$).

Flame No.	Air	$\mathrm{CH_{4}}$	$\mathrm{C_{2}H_{4}}$	$\phi$	$\bar{u}$	$L_{f}$ meas.	$\bar{Q}$
[-]	[ln/min]	[ln/min]	[ln/min]	[-]	[m/s]	[mm]	[W]
1	6.049	0.325	0.325	1.28	1.75	16.6	375
2	6.147	0.348	0.348	1.35	1.79	20.0	375
3	6.219	0.364	0.364	1.40	1.82	23.1	375
4	6.283	0.379	0.379	1.44	1.84	26.4	375
5	6.338	0.391	0.391	1.47	1.86	28.9	375
6	7.246	0.387	0.387	1.27	2.10	19.8	450
7	7.369	0.416	0.416	1.34	2.15	23.8	450
8	7.459	0.436	0.436	1.39	2.18	27.4	450
9	7.537	0.454	0.454	1.43	2.21	30.9	450
10	7.603	0.468	0.468	1.47	2.24	34.2	450
11	8.444	0.449	0.449	1.27	2.45	23.1	525
12	8.594	0.484	0.484	1.34	2.51	27.6	525
13	8.699	0.508	0.508	1.39	2.55	31.7	525
14	8.790	0.529	0.529	1.43	2.58	36.4	525
15	8.868	0.546	0.546	1.47	2.61	39.2	525
16	9.644	0.512	0.512	1.26	2.80	25.9	600
17	9.818	0.553	0.553	1.34	2.86	30.1	600
18	9.939	0.580	0.580	1.39	2.91	34.8	600
19	10.045	0.604	0.604	1.43	2.95	39.9	600
20	10.134	0.624	0.624	1.47	2.99	43.6	600