Data assimilation in a nonlinear time-delayed dynamical system

TL;DR

This paper introduces a data assimilation framework for nonlinear, time-delayed thermoacoustic models to improve their quantitative predictive capability, demonstrated through a horizontal Rijke tube experiment.

Contribution

It is the first application of data assimilation to nonlinear thermoacoustic models, enhancing real-time calibration with experimental data.

Findings

Decomposing pressure signals into a finite number of modes improves state estimation.

Mode decomposition is especially effective during nonlinear modal interactions.

The approach enables better quantitative predictions of thermoacoustic oscillations.

Abstract

When the heat released by a flame is sufficiently in phase with the acoustic pressure, a self-excited thermoacoustic oscillation can arise. These nonlinear oscillations are one of the biggest challenges faced in the design of safe and reliable gas turbines and rocket motors. In the worst-case scenario, uncontrolled thermoacoustic oscillations can shake an engine apart. Reduced-order thermoacoustic models, which are nonlinear and time-delayed, can only qualitatively predict thermoacoustic oscillations. To make reduced-order models quantitatively predictive, we develop a data assimilation framework for state estimation. We numerically estimate the most likely nonlinear state of a Galerkin-discretized time delayed model of a horizontal Rijke tube, which is a prototypical combustor. Data assimilation is an optimal blending of observations with previous state estimates (background) to produce optimal initial conditions. A cost functional is defined to measure the statistical distance between the model output and the measurements from experiments; and the distance between the initial conditions and the background knowledge. Its minimum corresponds to the optimal state, which is computed by Lagrangian optimization with the aid of adjoint equations. We study the influence of the number of Galerkin modes, which are the natural acoustic modes of the duct, with which the model is discretized. We show that decomposing the measured pressure signal in a finite number of modes is an effective way to enhance state estimation, especially when nonlinear modal interactions occur during the assimilation window. This work represents the first application of data assimilation to nonlinear thermoacoustics, which opens up new possibilities for real-time calibration of reduced-order models with experimental measurements.

Figures

Equations

$$\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x}$$

$$\frac{\partial p}{\partial t}+\frac{\partial u}{\partial x}+\zeta p-\dot{\mathcal{Q}}{\delta }_D(x-x_f)$$

$$\dot{\mathcal{Q}}\equiv \beta \text{Poly}(u_f(t-\tau )),$$

$$u(x,t)=\sum _{j=1}^{N_m}{\eta }_j(t)\cos (j\pi x),$$

$$p(x,t)=\sum _{j=1}^{N_m}\left(\frac{{\dot{\eta }}_j(t)}{j\pi }\right)\sin (j\pi x),$$

$${\text{F}}_{1j}$$

$${\text{F}}_{2j}$$

$$J$$

$$J_{bg}^a$$

$$J_{bg}^b$$

$$J_{bg}^c$$

$$J_{obs}^a$$

$$J_{obs}^b$$

$${\nabla }_{{\text{x}}_0}(J)={\nabla }_{{\text{x}}_0}(J_{bg})+\sum _{i=1}^{N_{obs}}{\nabla }_{{\text{x}}_0}(J_{obs,i})=0,$$

$${\text{F}}_3$$

$${\text{G}}_{1j}$$

$${\text{G}}_{2j}$$

$$\left[a,b\right]=\frac{1}{T}\underset{0}{\overset{T}{\int }}abdt$$

$$\mathcal{L}\equiv J_{bg}+J_{obs,i}+\sum _{j=1}^{N_m}{\mathcal{L}}_j-\left[\bar{\xi }(t),{\text{F}}_3\right],$$

$${\mathcal{L}}_j\equiv -\left[\frac{{\xi }_j}{j\pi },{\text{F}}_{1j}\right]-\left[{\nu }_j,{\text{F}}_{2j}\right]-b_{1j}{\text{G}}_{1j}-b_{2j}{\text{G}}_{2j},$$

$$\left[\frac{\partial \mathcal{L}}{\partial \mathbf{x}},\delta \mathbf{x}\right]\equiv \underset{ϵ\to 0}{\lim }\frac{\mathcal{L}(\mathbf{x}+ϵ\delta \mathbf{x})-\mathcal{L}(\mathbf{x})}{ϵ}.$$

Full text

1 Nonlinear time-delayed thermoacoustic model

We investigate the acoustics of a resonator excited by a heat source, which is a monopole source of sound. The main assumptions of the reduced-order model are [7]: (i) the acoustics evolve on top of a uniform mean flow; (ii) the mean-flow Mach number is negligible, therefore the acoustics are linear and no flow inhomogeneities are convected; (iii) the flow is isentropic except at the heat-source location; (iv) the length of the duct is sufficiently larger than the diameter, such that the cut-on frequency is high, i.e., only longitudinal acoustics are considered; (v) the heat source is compact, i.e., it excites the acoustics at a specific location, $x_{f}$; (vi) the boundary conditions are ideal and open-ended, i.e., the acoustic pressure at the ends is zero. Under these assumptions, the non-dimensional momentum and energy equations read, respectively [5]

[TABLE]

where $u$ is the acoustic velocity; $p$ is the acoustic pressure; $t$ is the time; $x$ is the axial coordinate of the duct; $\delta_{D}(x-x_{f})$ is the Dirac delta distribution at the heat source location, $x_{f}$; $\zeta$ is the damping factor, which models the acoustic energy radiation from the boundaries and thermo-viscous losses; and $\dot{\mathcal{Q}}$ is the heat release rate (or, simply, heat release). The heat release, $\dot{\mathcal{Q}}$, is modelled by a nonlinear time delayed law [10]

[TABLE]

where $\tau$ is the time delay; $\beta$ is the strength of the heat source; and $\textrm{Poly}(\tilde{u}(\tilde{t}))=a_{1}\tilde{u}^{5}(\tilde{t})+\cdots+a_{5}\tilde{u}(\tilde{t})$. The time delay and strength of the heat source are the two key parameters of a reduced-order model for the flame [3]. Physically, $\tau$ is the time that the heat release takes to respond to a velocity perturbation at the flame’s base; while $\beta$ provides the strength of the coupling between the heat source and the acoustics. Velocity, pressure, length and time are nondimensionalized as in [5]. The set of nonlinear time-delayed partial differential equations (1)-(2) provides a physics-based reduced-order model for the nonlinear thermoacoustic dynamics. Owing to the assumptions we made, the model can only qualitatively replicate the nonlinear thermoacoustic behaviour. In this paper, we propose a Lagrangian method to make a qualitative model quantitatively more accurate any time that reference data can be assimilated. Such data can come, for example, from sensors in experiments or time series from high-fidelity simulations.

1.1 Numerical discretization with acoustic modes

We use separation of variables to decouple the time and spatial dependencies of the solution. The spatial basis on to which the solution is projected consists of the natural acoustic modes. When decomposed on the natural acoustic eigenfunctions, the acoustic velocity and pressure read, respectively

[TABLE]

where the relationship between $\eta_{j}$ and $\dot{\eta}_{j}$ has yet to be found and $N_{m}$ is the number of acoustic modes considered. This discretization is sometimes known as the Galerkin discretization [12]. The state of the system is provided by the amplitudes of the Galerkin modes that represent velocity, $\eta_{j}$, and those that represent pressure, $\dot{\eta}_{j}/(j\pi)$. The damping has modal components, $\zeta_{j}=C_{1}j^{2}+C_{2}\sqrt{j}$, where $C_{1}$ and $C_{2}$ are damping coefficients [6, 9, 1, 5, 8]. In vector notation, $\bm{\eta}\equiv(\eta_{1},\cdots,\eta_{N})^{\textrm{T}}$ and $\bm{\dot{\eta}}\equiv(\dot{\eta}_{1}/\pi,\cdots,\dot{\eta}_{N}/(N\pi))^{\textrm{T}}$. The state vector of the discretized system is the column vector $\textbf{x}\equiv(\bm{\eta};\bm{\dot{\eta}})$. The governing equations of the Galerkin modes are found by substituting (4)-(5) into (1)-(3). Equation (2) is then multiplied by $\sin(k\pi x)$ and integrated over the domain $x=[0,1]$ (projection). In so doing, the spatial dependency is removed and the Galerkin amplitudes are governed by a set of nonlinear time-delayed differential equations

[TABLE]

where $u_{f}(t-\tau)=\sum_{j=1}^{N_{m}}\eta_{j}(t-\tau)c_{j}$; $s_{j}\equiv\sin(j\pi x_{f})$ and $c_{j}\equiv\cos(j\pi x_{f})$. The labels $\textrm{F}_{\bullet}$ are introduced for the definition of the Lagrangian (Sec. 3.2). Because the Galerkin expansions (4)-(5) are truncated at the $N_{m}$-th mode, we obtain a reduced-order model of the original thermoacoustic system (1)-(2) with $2N_{m}$ degrees of freedom (6)-(8). The reduced-order model is physically a set of $2N_{m}$ time-delayed oscillators, which are nonlinearly coupled through the heat release term. In the following sections, we employ 4D-Var data assimilation to improve the accuracy of such a reduced-order model111Although different from our study, it is worth mentioning that a study that combined 4D-Var data assimilation with reduced-order models of the Navier-Stokes equations based on proper orthogonal decomposition can be found in [4]..

2 Data assimilation as a constrained optimization problem

The ingredients of data assimilation are (i) a reduced-order model to predict the amplitude of thermoacoustic oscillations, which provides the so-called background state vector $\mathbf{x}^{bg}$ at any time, $t$ (red thick line in Fig. 1); (ii) data from external observations, $\mathbf{y}^{i}$, which is available from high-fidelity simulations or experiments at times $t^{i}$ (grey diamonds in Fig. 1); and (iii) an educated guess on the parameters of the system, which originates from past experience. The uncertainties on the background solution and observations are here assumed normal and unbiased. $\mathbf{B}$ and $\mathbf{R}$ are the background and observation covariance matrices, respectively, which need to be prescribed. For simplicity, we assume that $\mathbf{B}$ and $\mathbf{R}$ are diagonal with variances $B$ and $R$ (i.e., errors are statistically independent). A cost functional is defined to measure the statistical distance between the background predictions and the evidence from observations. First, we want the state of the system to be as close as possible to the observations. Second, we do not want the improved solution to be “too far” away from the background solution. This is because we trust that our reduced-order model provides a reasonable solution. Mathematically, these two requirements can be met, respectively, by minimizing the following cost functional

[TABLE]

where $N_{obs}$ is the number of observations over the assimilation window $[0,T]$. $\mathbf{M}$ is a linear measurement operator, which maps the state space to the observable space (see Eqs. (4)-(5)). Moreover, $||\textbf{x}_{0}-\textbf{x}_{0}^{bg}||_{\textrm{{B}}}^{2}\equiv\left(\textbf{x}_{0}-\textbf{x}_{0}^{bg}\right)^{\textrm{T}}\textrm{{B}}^{-1}\left(\textbf{x}_{0}-\textbf{x}_{0}^{bg}\right)$. Likewise, $||\mathbf{M}\mathbf{x}^{i}-\textbf{y}^{i}||_{\textrm{{R}}}^{2}\equiv\left(\mathbf{M}\mathbf{x}^{i}-\textbf{y}^{i}\right)^{\textrm{T}}\textrm{{R}}^{-1}\left(\mathbf{M}\mathbf{x}^{i}-\textbf{y}^{i}\right)$.

The objective of state estimation is to improve the prediction of the state, x, over the interval $[0,T]$, by reinitializing the background initial conditions, $\textbf{x}_{0}^{bg}$, with optimal initial conditions. These optimal initial conditions are called analysis initial conditions, $\textbf{x}_{0}^{analysis}$, which are the minimizers of the cost functional (9). We start from a background knowledge of the model’s initial conditions, $\textbf{x}_{0}^{bg}$, which is provided by the reduced-order model when data is not assimilated. By integrating the system from $\textbf{x}_{0}^{bg}$, we obtain the red trajectory in Fig. 1, $\textbf{x}^{bg}(t)$. The set of observations is assumed to be distributed over an assimilation window at some time instants. Pictorially, the analysis trajectory corresponds to the green thin line in Fig. 1, which is the minimal statistical distance between the background initial condition (magenta thick arrow) and observations (blue thin arrows). This algorithm is known as 4D-Var in weather forecasting [2]. State estimation enables the adaptive update of reduced-order thermoacoustic models whenever data is available.

3 Data assimilation for nonlinear thermoacoustic dynamics

We propose a set of cost functionals to perform data assimilation with the thermoacoustic model introduced in Sec. 1. We also introduce the formalism to perform Lagrangian optimization, in which adjoint equations enable the efficient calculation of the gradients of the thermoacoustic cost functionals with respect to the initial state.

3.1 Thermoacoustic cost functionals

Crucial to data assimilation is the definition of the cost functionals $J_{bg}$ and $J_{obs}$. Three physical thermoacoustic cost functionals are proposed and compared to reproduce different scenarios. For the background error

[TABLE]

For the observation error

[TABLE]

where $x_{m}$ is the location where the measurement is taken and $t_{obs}^{i}$ is the instant at which the $i-$th observation is assimilated. On the one hand, by using $J_{bg}^{a}$ and $J_{obs}^{a}$ the analysis solution is optimized against the background pressure at $t=0$ and the measured pressure at $t=t_{obs}^{i}$, $i=1,\dots,N_{obs}$. Physically, this means that the acoustic pressure is the model output, $p(0)_{bg}$, and the observable from the sensors, $p_{obs}^{i}$. On the other hand, $J_{bg}^{b}$ and $J_{obs}^{b}$ constrain the pressure modes. Physically, this means that every pressure mode provided by the background solution is a model output, $\left(\dot{\eta}_{j0}/(j\pi)\right)_{bg}$, and it is assumed that the modes of the acoustic pressure,$\left(\dot{\eta}_{j}/(j\pi)\right)^{(i)}_{obs}$, can be calculated from the sensors on the fly. For the background cost functional, we also defined $J_{bg}^{c}$ as a norm of the modes, which does not have a corresponding observation cost functional because the spatial dependency is not explicit.

To attain a minimum of $J$, a necessary condition is that the gradient vanishes, i.e.,

[TABLE]

where $\nabla_{{\textbf{x}_{0}}}$ is the gradient with respect to the initial conditions. There exists $\textbf{x}_{0}$ such that $\nabla_{{\textbf{x}_{0}}}(J)=0$ because of the convexity of the cost functionals in the neighbourhood of a local extremum. To compute $\nabla_{{\textbf{x}_{0}}}(J_{bg})$ and $\nabla_{{\textbf{x}_{0}}}(J_{obs,i})$, we use calculus of variation (Sec. 3.2). The Lagrange multipliers, also known as adjoint, or dual, or co-state variables (Sec. 3.3), provide the gradient information with respect to the initial state.

3.2 Lagrangian of the thermoacoustic system

The governing equations and their initial conditions are rewritten in the form of constraints, F, which hold over time intervals, while G are the constraints that hold for a specific time only, i.e., $t=t_{0}$. Together with equations (6)-(8) and by defining the auxiliary variable $\bar{\eta}(t)\equiv u_{f}(t-\tau)$, they read

[TABLE]

The constraints for the initial conditions read

[TABLE]

By defining an inner product

[TABLE]

where $a$ and $b$ are arbitrary functions, the Lagrangian of the nonlinear system can be written as

[TABLE]

where each $\mathcal{L}_{j}$ is

[TABLE]

where $\bar{\xi}$, $\xi_{j}/j\pi$, $\nu_{j}$ and $b_{\bullet j}$ are the Lagrange multipliers, or adjoint variables, of the corresponding constraints. Because we wish to derive the adjoint equations for the cost functional $J_{obs,i}$, we consider the time window to be $T=t_{obs}^{i}$.

3.3 Adjoint equations

We briefly report the steps to derive the evolution equations of the Lagrange multipliers (adjoint equations) [7]. First, the Lagrangian (21) is integrated by parts to make the dependence on the direct variables explicit. Secondly, the first variation is calculated by a Fréchet derivative

[TABLE]

Thirdly, the derivatives of (21) are taken with respect to the initial condition of each variable of the system, $\partial\mathcal{L}/\partial\textbf{x}_{0}$. These expressions will be used later to compute the gradient. Finally, to find the set of Lagrange multipliers that characterizes an extremum of the Lagrangian, $\mathcal{L}$, variations with respect to $\delta\textbf{x}$ are set to zero. The adjoint equations and their initial conditions are derived by setting variations of $\delta\eta_{j}$, $\delta\left(\dot{\eta}_{j}/(j\pi)\right)$ and $\delta\bar{\eta}$ to zero over $[0,T]$.

3.4 Gradient-based optimization

The optimization loop consists of the following steps:

Integrate the system forward from $t=0$ to $t=T$ from an initial state $\textbf{x}_{0}$;

• 2)

Initialize the adjoint variables;

• 3)

Evolve the adjoint variables backward from $t=T$ to $t=0$;

• 4)

Evaluate the gradient using the adjoint variables at $t=0$.

Once the gradient is numerically computed, the cost functional can be minimized via a gradient based optimization loop. The conjugate gradient [11] is used to update the cost functional until the condition $\nabla_{\textbf{x}_{0}}(J)=0$ is attained to a relative numerical tolerance of $10^{-4}$. By using a gradient based approach, we find a local minimum of $J$. We verify that there is no other local minimum by computing $J=J(\textbf{x}_{0})$ in the vicinity of $\textbf{x}_{0}^{analysis}$.

4 Results

We validate the data-assimilation algorithm by twin experiments: The true state solution, $\textbf{x}^{true}(t)$, is produced by perturbing the unstable fixed point at the origin of the phase space222Strong constraint 4D-Var assumes that the model is perfect and the uncertainty is only in the initial conditions, therefore, the true trajectory can be a model output.; the background trajectory, $\textbf{x}^{bg}(t)$, is obtained by perturbing each true mode initial condition with Gaussian error with variance $B=0.005^{2}$; the $i-$th observation is produced by adding Gaussian error with variance $R=0.005^{2}$ to $\textbf{x}^{true}(t^{i}_{obs})$. The outcome of twin experiments is summarized by the error plots shown in Figs. 3 and 4. The cyan vertical line indicates the end of the assimilation window, the red thick line is the difference between the true pressure and the background pressure, the green thin line is the difference between the true pressure and the analysis pressure. First, it is shown how the number of computed acoustic modes affects the solution of the system. Secondly, we investigate the effects that the different cost functionals have on the analysis solution. Finally, we discuss the effects that the number of observations have on the analysis.

The parameters we use are $\beta=1$, $\tau=0.02$, $C_{1}=0.05$, $C_{2}=0.01$ and $(a_{1},a_{2},a_{3},a_{4},a_{5})=(-0.012,0.059,-0.044,-0.108,0.5)$ for the heat release term $\dot{\mathcal{Q}}$. The position of the heat source is $x_{f}=0.3$ and all the measurements are taken at $x_{m}=0.8$.

4.1 Remarks on the thermoacoustic nonlinear dynamics

The thermoacoustic model is a set of $2N_{m}$ nonlinearly coupled oscillators, which we initialize by imposing non-equilibrium initial conditions. We compare two solutions, using $N_{m}=3$ and $N_{m}=10$ in Figs. 2a and 2b, respectively. Higher modes are quickly damped out, thus, after a transient where strong nonlinear modal coupling occurs, the solution obtained with $N_{m}=10$ is qualitatively similar to the solution obtained with $N_{m}=3$. During the transient, if sufficient modes are computed, the dynamics are more unpredictable because of the intricate modal interaction. The twin experiments are performed with $N_{m}=10$, which provide a more accurate solution. It is shown that state estimation is markedly affected depending on whether we observe the system during the transient or at regime.

4.1.1 Effect of the observation error

We can simulate two main scenarios, depending on the choice of $J_{obs}$ (Sec. 3.1). Figure 3a is obtained using $J_{obs}^{a}$, i.e., by modelling observations about the pressure only. The analysis pressure error slightly deviates from zero in the assimilation window. When the forecast window starts, the analysis suddenly approaches the background again. Figure 3b is obtained using $J_{obs}^{b}$, therefore the observations contain information about every pressure modes. The forecast quality is considerably enhanced. The assimilation window is $T_{as}=0.4$, thus, the observations are obtained during the transient, which lasts up to $t\approx 2$, where the dynamics are more unpredictable due to nonlinear interactions between modes. As we will show in the next subsection, increasing $N_{obs}$ is not an effective strategy to improve the forecast during the transient when the pressure is observed.

4.1.2 Effect of the background error

From a numerical standpoint, the background error, $J_{bg}$, acts as an observation at $t=0$. The analysis trajectory is an optimal blending of the information from the measurements and the previous educated model output. Therefore, we emphasize that the outcome of twin experiments is improved if the cost functionals of the background knowledge and observations are consistent. In the present framework, it means that $J_{bg}^{a}$ should be used with $J_{obs}^{a}$ and $J_{bg}^{b}$ should be used with $J_{obs}^{b}$. On the one hand, if the source of assimilated data favours the background knowledge (e.g. using $J_{bg}^{b}$ or $J_{bg}^{c}$ together with $J_{obs}^{a}$), the analysis trajectory will be closer to the background. On the other hand, if the source of assimilated data favours the observations (e.g. using $J_{obs}^{b}$ with $J_{bg}^{a}$), the analysis trajectory will be closer to the observations.

4.1.3 Effect of the number of observations

Generally speaking, the higher the number of observations the more the optimal solution will be similar to the true solution. This can be deduced by inspection of Figs. 4a and 4b. The value of $N_{obs}$ is increased from 50 to 250, over an assimilation window of 2.5 time units (the observations are not shown in the figures), resulting in a smaller error amplitude when more observations are available.

However, it is possible that the increasing of $N_{obs}$ will not result in a better state estimation. This happens if we measure the pressure during the transient, that is, using $J_{obs}^{a}$ and $T_{as}<2$. In the transient, the measured pressure is a combination of all modes. Therefore, the same pressure level is associated with different combinations of modes, hence, no useful information is assimilated to help determine the state of the system by knowing only the pressure (i.e., by using $J_{obs}^{a}$). Under these circumstances, the forecast quality will remain poor, as shown in Fig. 3a, regardless of the number of observations. When the assimilation window is extended up to 2.5 time units, as shown in Fig. 4, the pressure is observed also at regime, that is in $t\in[2,2.5]$, approximately. In this interval, the measured pressure signal is produced chiefly by the first three modes because higher modes are damped out. Therefore, the measured pressure becomes more effective information about the state to assimilate. At regime, as intuitively expected, increasing the observation frequency produces a better forecast.

We conclude that having poor information about the system’s state cannot be simply balanced by increasing the number of observations. Given a number of observations with their time distrbution, it is the synergy between an appropriate cost functional and the recognition of the type of dynamics (transient vs. regime) to be the key for a successful data assimilation.

5 Conclusions and future work

Preliminary thermoacoustic design is based on simplified and computationally cheap reduced-order models that capture the inception of thermoacoustic instabilities (linear analysis) and their saturation to finite amplitude oscillations (nonlinear analysis). We propose a Lagrangian method to make a qualitative reduced-order model quantitatively more accurate any time that reference data can be assimilated. Such data can come, for example, from sensors in experiments or time series from high-fidelity simulations. To test the method we perform a series of twin experiments with the thermoacoustic model of a resonator excited by a heat source (horizontal Rijke tube). When sufficient modes are computed, a clear distinction emerges between a transient state and the state at regime. The former is characterized by irregular dynamics due to the interaction between all modes, while at regime the dynamics are chiefly dominated by the first three modes. We find that, at regime, it is possible to enhance the forecast by assimilating data about the pressure. As intuitively expected, the higher the number of observations, the better the forecast accuracy. While testing the effectiveness of data assimilation during the transient, we find that it is not possible to improve the forecast by measuring the pressure only. Moreover, the quality of the forecast remains poor regardless of the number of observations. Therefore, we propose a more effective cost functional, which takes into consideration the spectral content of the measured signal to enable a successful state estimation also during the transient. In state estimation, we implicitly assume that the parameters are correct. However, this is rarely the case in thermoacoustics, where the parameters are uncertain and need to be optimally calibrated. Ongoing work includes simultaneous parameter and state estimation using Lagrangian optimization with state augmentation. This work opens up new possibilities for on-the-fly optimal calibration and state estimation of reduced-order models in thermoacoustics for applications in propulsion and power generation.