A reduced-order mean-field synchronization model for thermoacoustic systems

The combustor comprises a settling chamber, a mixing duct, a combustion chamber, and a circular bluff body (see Figure 1a). Air passes through the settling chamber, where the flow irregularities are removed, and enters the mixing tube. The fuel (liquified petroleum gas, 40% propane and 60% butane) enters the mixing tube via radial injection ports on the central shaft. Both air and fuel mix in the mixing tube and flow into the combustion chamber.

The length of the combustion chamber is 1100 mm with a square cross-section of dimensions 90 $\times$ 90 mm². The combustion chamber near the bluff body has a backward-facing step. The reactants enter the combustion chamber from this side and are ignited with the help of a spark plug. The other end of the combustion chamber embodies a decoupler of dimensions 1000 $\times$ 500 $\times$ 500 mm³, which isolates the combustion chamber from the ambient fluctuations outside. The flame is stabilized by the circular bluff body having a thickness of 10 mm and a diameter of 47 mm. The bluff body is located 32 mm downstream of the backward-facing step. A fixed fuel flow rate of 48 SLPM is maintained, with the airflow rate varying from 752 SLPM to 1446 SLPM. This leads to a variation of the Reynolds number $Re$ from 3.9 $\times$ $10^{4}$ to 6.4 $\times$ $10^{4}$ and the equivalence ratio $\phi$ from 1 to 0.52. Here, the Reynolds number is given as $Re=Vd/\nu$, where $V$ is the mean velocity of the flow, $d$ is the bluff-body diameter and $\nu$ is the kinematic viscosity. The equivalence ratio is defined as the ratio of the actual fuel-air ratio to the stoichiometric fuel-air ratio.

The heat release rate fluctuations are measured using the fluctuations in chemiluminescence intensity. [31, 32] These intensity fluctuations are acquired using a high-speed camera equipped with a CH^∗ filter and a 100 mm Carl-Zeiss lens. The images were acquired at a sampling rate of 2 kHz at a resolution of 520 $\times$ 520 pixels. The noise effects were removed by coarse-graining the images by combining 10 $\times$ 10 pixels. A detailed description of the dump combustor setup and the measurements can be found in Sudarsanan et al. [33].

Figure 1(b) shows a schematic of the annular combustor. The annular combustor comprises 16 burner tubes, each 30 mm in diameter and 150 mm long. A technically premixed fuel-air mixture enters the settling chamber through twelve equally spaced inlet ports. The fuel used was liquified petroleum gas. A honeycomb mesh in the settling chamber removes the non-uniformities in the flow. The flow is then guided uniformly in the 16 burner tubes through a hemispherical divider. Each burner tube has a swirler on the top for flame stabilization and a flame arrestor at the bottom for preventing flashback. Swirlers mounted on each burner comprise a central shaft (15 mm diameter) and six guide vanes, each making a 60^∘ angle with the shaft axis (see Fig. 1c).

The combustion chamber consists of two concentric ducts with inner and outer diameters of 300 mm and 400 mm, respectively. The length of the outer duct is 400 mm, and the inner duct is 200 mm. A non-premixed pilot flame (see Fig. 1d) ignites the main combustor. The pilot flame is later extinguished after flame stabilization. The airflow rate is fixed at 1400 SLPM. The fuel flow rate is increased from 40 to 48 SLPM ($\phi$ = 0.82 to 0.98) in the forward direction and decreased from 48 to 36 SLPM ($\phi$ = 0.98 to 0.73) in the reverse direction. Considering the fixed airflow rate, which is much larger than the fuel flow rate, the Reynolds number is calculated to be $Re_{d}\approx 8600$ using the exit diameter of the burner tube ($d=15$ mm).

The CH^∗ chemiluminescence intensity fluctuations of swirling flames were captured using a high-speed camera. These intensity fluctuations are used to obtain the heat release rate oscillations. Images were acquired at a resolution of 1280 $\times$ 800 pixels at a sampling rate of 2 kHz. Only one-half of the combustor backplane (the shaded portion in Fig. 1d) was imaged with the help of an air-cooled mirror 1 m above the combustor. The camera lens used was Nikkon AF Nikkor with 70 - 210 mm $f$/4 to $f$/5.6. Additional information about the annular combustor setup can be found in Roy et al. [12], and Singh et al. [34].

The mass flow controllers (MFCs) used to regulate the flow rates of fuel and air in both the annular and dump combustor setups are Alicat Scientific MCR Series. These MFCs have an uncertainty of $\pm 0.8\%$ of the measured reading and $\pm 0.2\%$ of the full-scale reading. The maximum uncertainties in $Re$ and $\phi$ are $\pm 0.8\%$ and $\pm 1.6\%$, respectively. A piezoelectric transducer PCB103B02, having a sensitivity of 217.5 mV/kPa and an uncertainty of $\pm 0.15$ Pa was used to obtain the pressure fluctuations. The signal is recorded for a duration of 3 s at a sampling frequency of 10 kHz for both combustor configurations. The chemiluminescence images were obtained using the high-speed camera (CMOS Phantom V12.1) fitted with a CH^∗ filter. The fluctuations in global heat release rate were acquired using a photomultiplier tube (PMT) equipped with a CH^∗ filter. The CH^∗ filter used in both the camera and the PMT has a bandwidth of 435 $\pm$ 10 nm.

The linearized momentum and energy equations in a one-dimensional thermoacoustic system with the negligible effect of temperature gradient and mean flow, driven by an oscillatory heat release rate, are [35, 36]

where $\tilde{p}^{\prime}$ and $\tilde{u}^{\prime}$ are the mean subtracted fluctuations of the acoustic pressure and velocity, respectively. $\tilde{p}_{0}$ and $\tilde{\rho}_{0}$ are the pressure and density of the mean flow, $\tilde{z}$ denotes the axial distance along the duct, $\gamma$ is the ratio of specific heat capacities, and $\tilde{t}$ is time. $\dot{\tilde{q}}^{\prime}$ represents the unsteady heat release rate fluctuations per unit area. These fluctuations are assumed to be acoustically compact with concentration at $\tilde{z}=\tilde{z}_{f}$, represented by the Dirac delta function $\delta(\tilde{z}-\tilde{z}_{f})$. $\tilde{()}$ indicates that the terms are dimensional.

The above partial differential equations can be transformed into ordinary differential equations by employing the Galerkin technique [24]. In this method, the acoustic pressure and velocity fluctuations are expanded as a superposition of spatial basis functions (sine and cosine) with time-varying coefficients ($\eta(\tilde{t})$ and $\dot{\eta}(\tilde{t})$) that satisfy the given boundary conditions. We choose the basis functions satisfying the boundary conditions of an open-closed duct.

Here $\tilde{c}_{0}$ is the average sound velocity in the duct, $\tilde{k}_{j}=(2j-1)\pi/(2\tilde{L})$ is the wavenumber, $\tilde{\Omega}_{j}=\tilde{c}_{0}\tilde{k}_{j}$ represents the natural frequency for the $j^{th}$ mode of the system, and $\tilde{L}$ is the duct length. Substituting the expansions for $\tilde{p}^{\prime}$ and $\tilde{u}^{\prime}$ in Eq. (2), we get

The resulting equation is projected onto the basis functions by computing an inner product of Eq. (5) with $\cos(\tilde{k}_{j}\tilde{z})$. This involves multiplying both sides of Eq. (5) with $\cos(\tilde{k}_{j}\tilde{z})$ and integrating over the spatial domain from $0$ to $\tilde{L}$. As a result, the second-order ordinary differential equations for individual modes are given as

Here, we use $\tilde{c}_{0}=\sqrt{\gamma\tilde{p}_{0}/\tilde{\rho}_{0}}$ and the identity, $\int_{0}^{\tilde{L}}cos^{2}(\tilde{k}_{j}\tilde{z})d\tilde{z}=\tilde{L}/2$. We neglect the effects of higher modes and assume that the single-mode analysis captures the transitions reasonably well, [37] to obtain

where we introduce the term $\tilde{\alpha}\dot{\eta}(\tilde{t})$ for accounting the acoustic damping of the system [38] with $\tilde{\alpha}$ being the damping coefficient.

A model for the heat release rate $\dot{\tilde{q}}^{\prime}$ is required to obtain the acoustic field in the system. The mean field thermoacoustic model considers the flame response as a combined effect of a population of Kuramoto phase oscillators coupled with the acoustic field. [26, 27] The dynamical equation for the $j^{th}$ phase oscillator having a phase ($\theta_{j}$) at time $\tilde{t}$ and having a natural frequency ($\tilde{\omega}_{j}$) is

where $\hat{R}$ is the normalized amplitude of the acoustic variable ($\eta$) and $\Phi$ is its phase. In this model, the term $\tilde{K}\hat{R}$ represents the effective coupling strength of the oscillators with the acoustic field, which signifies that high-amplitude acoustic fluctuations enhance coupling. The individual contribution of each Kuramoto phase oscillator is summed to obtain $\dot{\tilde{q}}^{\prime}$ as

where $\tilde{q}_{0}$ represents the heat release rate amplitude of individual oscillators and N denotes the number of oscillators. Equation (9) indicates that the heat release rate fluctuations are high (low) when the oscillators are synchronized (unsynchronized). The level of synchrony is quantified by the order parameter $z$, defined as

where $\psi$ and $r$ are the phase and magnitude of the complex order parameter $z$. Substituting the heat release rate from Eq. (9) in Eq. (7), the equation governing the acoustic field can be written as

where $\tilde{\beta}=2\tilde{q}_{0}(\gamma-1)/\tilde{L}\tilde{p}_{0}$ is the measure of flame strength. Equations (8) and (11) are non-dimensionalized using

This leads to the following non-dimensionalized equations:

Assuming the fluctuations to be quasi-harmonic, we decompose the acoustic variable $\eta(t)$ as [39]

where $\Phi(t)$ and $R(t)$ are the phase and the envelope amplitude, which evolve slower as compared to the fast time scale $2\pi/\Omega_{0}$. Using this decomposition in Eq. (14) yields a second-order ordinary differential equation in $R$ and $\Phi$. Further, the sines and cosines can be expressed in complex form using Euler’s representation, and the fast time scales can be eliminated using the method of averaging [40]. Equating the imaginary and real parts of the resulting differential equation leads to (see Appendix A for detailed derivation)

By multiplying both sides of Eq. (10) by $e^{-i\Phi}$ and equating the imaginary and real parts of the resulting equation, we can express the summation of sines and cosines in terms of $r$ and $\psi$ as

Using Eqs. (18) and (19) in Eqs. (16) and (17) leads to

We now normalize the envelope amplitude ($R$) using the amplitude of limit cycle oscillations ($R_{LCO}$) at $K\to\infty$. To compute $R_{LCO}$, we transform Eqs. (20) and (21) to Cartesian coordinates using the following transformation:

The transformed equations become

The limit cycle solution is given by $\ddot{A}=\dot{A}=0$ and $\ddot{B}=\dot{B}=0$. For $K\to\infty$, the oscillators will be synchronized, and the magnitude of the order parameter will be $r=1$. Using these conditions in Eqs. (23) and (24) leads to the following:

By using Eq. (22) and expressing $A$ and $B$ from Eqs. (25) and (26) in terms of $R$, we obtain,

The acoustic amplitude is normalized as $\hat{R}={R}/{R_{LCO}}$ in Eqs. (20) and (21), resulting in the following equations,

The system of equations (13), (28) and (29) represent the mean-field thermoacoustic model and can be integrated in time to obtain the bifurcation diagrams.

The equations (13), (28) and (29) have a dimensionality of $N$, 2 and 2, respectively, resulting in an overall dimensionality of N+4, where N is the number of oscillators. For a large number of oscillators, the model becomes analytically intractable and the simulations become computationally expensive. In this section, we utilize the Ott-Antonsen method to derive a reduced-order model. To reduce the dimensionality of the system, a continuum limit of $N\rightarrow$ $\infty$ is considered, and a probability density function ($f$) is defined such that at any time $t$, $f(\theta,\omega,t)d\theta d\omega$ represents the fraction of oscillators having natural frequency between $\omega$ and $\omega+d\omega$ and phase between $\theta$ and $\theta+d\theta$. The probability density function satisfies the normalization,

which leads to,

Here, $g(\omega)$ is the time-independent natural frequency distribution of the oscillators. Since the number of oscillators is conserved, $f$ satisfies the continuity equation, [29]

where $v=\dot{\theta}$ is the angular velocity of the oscillators given by Eq. (13). The probability density function $f(\theta,\omega,t)$ can be expressed as a Fourier series expansion in $\theta$ as

where $a(\omega,t)$ is the Fourier series coefficient and $\overline{a}(\omega,t)$ is its complex conjugate. We assume $|a(\omega,t)|<1$ for ensuring convergence. This form of the probability density function, with the Fourier coefficients as a power series of $a(\omega,t)$, naturally connects the incoherent and partially synchronized states. [29] Substituting $f$ and $v$ from Eqs. (33) and (13) in Eq. (32) gives

In the continuum limit ($N\to\infty$), the order parameter defined in Eq. (10) is expressed as

Substituting the Fourier expansion of $f$ (Eq. 33) in the above equation leads to

The natural frequency distribution $g(\omega)$, which represents the density of oscillators with a natural frequency $\omega$, must be known to compute the above integral. Singh et al. [27] utilized the FFT of chemiluminescence data during the combustion noise state as the natural frequency distribution of the oscillators. This distribution can be effectively approximated using a weighted sum of Lorentzian distributions [41] as

by minimizing Kullback-Leibler (KL) divergence. Here, $\omega_{k}$ is the peak location, $\gamma_{k}$ is the half-width at half-maximum and $c_{k}$ is the weight for the $k^{th}$ Lorentzian distribution. $n$ denotes the total number of Lorentzian distributions used to approximate the experimental distribution. The natural frequency distribution satisfies the normalization, $\int_{-\infty}^{\infty}g(\omega)d\omega=1$ which leads to

Lorentzian distributions are chosen for the fitting because they contain simple poles, which makes the computation of integrals straightforward. Using the frequency distribution from Eq. (37), the integration in Eq. (36) is computed with the aid of the residue theorem. A closed contour (Fig. 2) is selected, consisting of the real $\omega$-axis and a semi-circle with a large radius ($\mathcal{R}\to\infty$) in the upper half of the complex $\omega$-plane ($\Im(\omega)>0$). The poles of the frequency distribution inside this contour are given by $\mathcal{Z}_{k}=\omega_{k}+i\gamma_{k}$. Using these poles in the residue theorem, the integration in Eq. (36) simplifies to

The complex Fourier coefficients $\overline{a}(\omega_{k}+i\gamma_{k},t)$ are expressed in polar forms as $r_{k}\exp(i\phi_{k})$ and substituted in Eq. (34). By comparing the imaginary and real parts on both sides of the resulting equation, we obtain the following set of equations

where,

The fixed-point solutions of the system are obtained by equating the time derivatives to zero. For the system of Eqs. (28), (29), and (40) - (42), the fixed point corresponds to the state where $\dot{r}_{k}$, $\dot{\phi}_{0k}$, $\dot{\hat{R}}$, $\ddot{\hat{R}}$, and $\ddot{\Phi}$ are zero. These equations are nonlinear and require numerical methods to solve them. We use the arclength continuation method, [42, 43] where the solution curve is traced iteratively. The initial guess is computed by taking a small step along the tangent of the curve for each iteration. This guess is then updated along a circle with a radius equal to the step length until convergence. Once the $r_{k}$ and $\phi_{k}$ at the fixed point solutions are computed, the order parameter can be obtained using Eq. (39). One key advantage of the arclength continuation method over the commonly used Newton-Raphson method is its ability to obtain all solutions, even when the function exhibits folds [44]. This enables us to capture the unstable limit cycle solutions.

The system has a trivial solution $r_{n}=\hat{R}=0$, which cannot be numerically computed using the arclength continuation method due to ill-conditioning of the Jacobian matrix. To obtain the stability characteristics at these trivial fixed points, we transform the system to Cartesian coordinates using $\overline{a}(\omega_{k}+i\gamma_{k},t)=x_{k}+iy_{k}$ in Eq. (34). This transformation leads to

Equations (22), (23) and (24) can be normalized using the limit cycle amplitude ($R_{LCO}$) as

From Eq. (39), we have

Using Eqs. (46), (47) and (48), we obtain

Equations (43), (44), (49) and (50) have a trivial fixed-point solution at $x_{k}=y_{k}=\hat{A}=\hat{B}=0$, where the time-derivatives become zero. The stability of these fixed points can be evaluated by calculating the eigenvalues of the Jacobian at these points.

The stability of a fixed-point solution is determined by computing the eigenvalues of the Jacobian at that point. For a stable fixed point, all eigenvalues must have a negative real part. As we move along the solution branch, the eigenvalues vary smoothly, and some will cross the imaginary axis when the stability of the system changes. Therefore, the locations on the solution branch where the real part of any eigenvalue becomes zero correspond to the locations of bifurcation points [45].

Consider a dynamical system $\dot{X}=f(X)$, where $X=[x_{1},x_{2},...,x_{n}]^{T}$ is the state of the system at any time $t$ and $f=[f_{1},f_{2},...,f_{n}]^{T}$ is a differentiable function. The fixed-point solution $X^{*}$ is given as $f(X^{*})=0$. The Jacobian $J(X)$ of $f(X)$ and its eigenvalues $\lambda(X)$ at any state $X$ are given by

where $I$ is an identity matrix of order $n\times n$. The product of the real parts of the eigenvalues for any fixed-point solution $X^{*}$ is

If the real part of any eigenvalue corresponding to $X^{*}$ is zero, $F(X^{*})$ will be zero. Therefore, the fixed-point solutions satisfying $F(X^{*})=0$ correspond to the bifurcation locations.

We obtain the bifurcation plot by solving the system of equations (28), (29), and (40) - (42) for the fixed-point solutions. These equations require the parameters $c_{k}$, $\omega_{k}$ and $\gamma_{k}$ of the Lorentzian distributions in Eq. (37). For accurate results, we learn these parameters by minimizing Kullback-Leibler (KL) divergence between the experimentally obtained and model frequency distributions, using an optimization algorithm. The KL divergence, also known as relative entropy, is commonly used in information theory to quantify the similarity between two density functions. [46] Let $\mathbb{E}$ and $\mathbb{M}$ denote the experimental and model frequency distributions, then the KL divergence between them is expressed as

where $\mathbb{E}_{i}$ and $\mathbb{M}_{i}$ are probabilities corresponding to the frequency $f_{i}\in\{f_{1},f_{2},...,f_{\mathcal{N}}\}$. Here, $f_{1}=0$ and $f_{\mathcal{N}}$ is the maximum frequency. $\mathbb{E}_{i}$ is obtained for discrete frequencies $f_{i}$ from the FFT of the experimentally obtained heat release rate data. The model frequency distribution depends on the frequency $f_{i}$ and the vector of parameters $\mathbf{a}=[c_{1},c_{2},...,c_{n-1},\gamma_{1},\gamma_{2},...,\gamma_{n},\omega_{1},\omega_{2},...,\omega_{n}]$ as $\mathbb{M}_{i}=\mathbb{M}(f_{i},\mathbf{a})$, where $\mathbb{M}$ is given by Eq. (37). Here, $c_{n}$ is eliminated from $\mathbf{a}$ using Eq. (38). The parameter vector $\mathbf{a}$ is obtained by minimizing the KL divergence $\mathcal{D}$ using the gradient descent method [47]

where $\sigma_{L}$ is the learning rate. We use learning rates of $10^{-6}$ and $10^{-5}$ for frequency distributions of the annular and dump combustor, respectively. The gradient ($\nabla_{\mathbf{a}}\mathcal{D}$) of $\mathcal{D}$ with respect to $\mathbf{a}$ was evaluated using automatic differentiation. [48]

The authors have no conflict to disclose.