Self-excited vortex-acoustic lock-in in a bluff body combustor

Abstract

We investigate self-excited vortex-acoustic lock-in in a bluff body combustor. The bluff body not only anchors the flame but also sheds vortices. Mutual interaction between vortex shedding and the combustor acoustic field leads to lock-in. A lower-order model for vortex shedding is coupled to the acoustic equations to obtain a set of discrete dynamical maps, which are then solved numerically. Lock-in is identified when a definite phase relationship between vortex shedding and acoustic field remains unaltered with time. The common frequency during (phase) lock-in is either close to the natural acoustic or vortex shedding frequencies, accordingly termed A- or V-lock-in, respectively. Instability and amplitude suppression occur with most parts of the A- and V-lock-in regions, respectively. Furthermore, we relate the lock-in and pre-lock-in regimes observed in our previous experiments (Singh & Mariappan, Combust. Sci. Technol., 2019, pp. 1–29) to A- and V-lock-in phenomena, respectively. Among the two lock-in phenomena, combustion instability is favoured by $1:1$ A-lock-in, whose boundaries in the parameter space can be obtained from a forced response of the vortex shedding process. Finally, we discuss the bifurcations leading to lock-in.

1. Introduction

Recent combustor designs focus on low-emission and fuel-flexible operations. Large amplitude self-sustained detrimental oscillations, termed combustion or thermoacoustic instability, occur during such operations (Lieuwen & Yang 2005; Khalil & Gupta 2013). Conventional combustors and afterburners employ swirlers (Huang & Yang 2009) and bluff bodies (Hertzberg, Shepherd & Talbot 1991) to anchor the flame. Both the geometries support vortex shedding (Schadow & Gutmark 1992). Several experimental (Poinsot et al. 1987; Emerson & Lieuwen 2015) and theoretical (Matveev & Culick 2003) investigations highlight the role of vortex shedding in thermoacoustic instability. Vortical structures in combustors are observed due to two mechanisms: (i) (absolute) wake instability in an unexcited flow (Monkewitz 1988) and (ii) longitudinal (Poinsot et al. 1987; Candel 1992) and transverse (Rogers & Marble 1956; Kushwaha et al. 2017) acoustic excitations. In-turn vortices perturb the heat release rate by affecting the flame surface area (Emerson & Lieuwen 2015), equivalence ratio (Lieuwen & Zinn 1998) etc. Heat release rate fluctuations generate the acoustic field, thus closing the thermoacoustic feedback.

In vortex shedding combustors two self-excited systems (oscillators) are involved. (i) Vortex shedding behind flame stabilizers, causing heat release rate fluctuations of the flame at the frequency of shedding (Poinsot et al. 1987). (ii) Acoustic field of the combustor can be excited in two ways. The first one is due to the thermoacoustic feedback with an unsteady flame, irrespective of whether the flow is laminar (for example, Kabiraj & Sujith 2012) or turbulent (for example, Nair, Thampi & Sujith 2014). In the second way occurring only in turbulent combustors, the acoustic field can be excited by the turbulent flow fluctuations, even without an active unsteady flame (p.684 of Pawar et al. 2017). In general, the two systems oscillate at their own frequency. However, at certain operating conditions, they interact, leading to lock-in (Poinsot et al. 1987; Singh & Mariappan 2019). During (phase) lock-in, both the systems oscillate at a common frequency and evolve with a constant phase difference (Pikovsky, Rosenblum & Kurths 2003; Balanov et al. 2008). Often, the common frequency is close to the natural duct acoustic frequency (Poinsot et al. 1987). Therefore, many investigations replaced the acoustic feedback with a unidirectional velocity excitation to study the hydrodynamic behaviour, both in non-reacting (Li & Juniper 2013a) and reacting (Emerson & Lieuwen 2015) flow configurations.

In practical systems involving a set of can-annular combustors, the mutual interaction between each of the self-excited thermoacoustic systems (STS) leads to lock-in or desynchronization depending on the strength and time delay in coupling (Moon et al. 2020). Furthermore, the interaction leads to amplification or amplitude death (Hyodo, Iwasaki & Biwa 2020) of pressure oscillations even when the individual uncoupled combustors are stable or unstable, respectively (Jegal et al. 2019). Specifically, if one or more can-combustors have a weaker thermal power, the interaction leads to new regions of pressure amplifications (Doranehgard, Gupta & Li 2022) in the weaker combustor. The above dynamics depend on the individual system parameters, strength and type of coupling between the combustors. Overall, lock-in is frequently encountered in practical combustors, reiterating its importance.

A systematic investigation on the effects of external excitation in relation to lock-in on a low-density (non-reacting) jet experiencing absolute instability at a frequency $(f_n)$ was investigated by Li & Juniper (2013a). Excitation frequency $(f_f)$ and amplitude ($A_f$) were varied. On forcing in the vicinity of frequency ratio $f=f_f/f_{n}=1$, a V-shaped lock-in region, stemming at $f=1$ occurs in the excitation amplitude–frequency plane. Inside the lock-in region, the jet oscillates at the forcing frequency. Two bifurcations to lock-in are observed as the excitation amplitude is increased. (i) Saddle-node bifurcation: when $f$ is close to 1, $f_n$ is pulled towards $f_f$, keeping the response amplitude almost constant. A theoretical analysis on a physically relevant lower-order model for vortex shedding does indicate the above bifurcation to lock-in (Britto & Mariappan 2019). (ii) Inverse Nemark–Sacker bifurcation: when $f$ is farther from 1, the response amplitude at $f_n$ decreases (no movement in $f_n$) and vanishes completely during lock-in.

Apart from the two bifurcations, two kinds of asymmetries are observed in the lock-in boundary about $f=1$. (i) Asymmetry in the geometry of the boundary. Li & Juniper (2013a) (non-reacting flow) and Li & Juniper (2013b) (reacting flow) reported the occurrence of lock-in at lower excitation amplitudes in the $f<1$ region than in the $f>1$ region and vice versa, respectively. (ii) Asymmetry in the response amplitude of the system. In the non-reacting case (Li & Juniper 2013a), where the response is measured as velocity fluctuations, the response is stronger and weaker than the unexcited case at $f<1$ and $f>1$ regions, respectively. Heat release rate fluctuations $\dot q'$, measured as the response in the reacting case (Li & Juniper 2013b), indicated an amplification/reduction (compared with the unexcited case) in $f<1/f>1$ regions. Therefore, lock-in can be associated with both amplification and reduction of the response of the system. On a similar line, investigations by Emerson & Lieuwen (2015) reported a reduction in $\dot {q}'$ of a bluff body stabilized flame when excited at the natural (reacting) hydrodynamic frequency. Furthermore, when a STS is forced far away from its natural frequency, asynchronous quenching occurs, where suppression in the amplitude of pressure oscillations is observed (Balusamy et al. 2015; Guan et al. 2019a; Mondal, Pawar & Sujith 2019).

The above investigations focus on the forced synchronization of a thermoacoustic system whose natural state is periodic self-excited oscillation. However, a thermoacoustic system exhibits a wide range of dynamics from periodic, quasiperiodic to chaotic oscillations (Kabiraj, Sujith & Wahi 2012; Kashinath, Waugh & Juniper 2014). In this interest, forced synchronization of an aperiodic (quasiperiodic and chaotic systems) thermoacoustic system was first investigated by Kashinath, Li & Juniper (2018). The introduction of harmonic external forcing on the aperiodic system destroys the aperiodic attractor leading to the emergence of order in the system (period-1 or higher-periodic oscillations). A quasiperiodic thermoacoustic system that naturally contains two incommensurate frequencies first undergoes partial synchronization (Guan et al. 2019b) wherein the weaker natural mode is locked in with the forcing system before complete synchronization. The synchronization results in a period-1 limit cycle, where the system oscillates at the forcing frequency. Furthermore, strong forcing near the weaker natural mode is sufficient for complete synchronization resulting in the most significant reduction of thermoacoustic oscillations through asynchronous quenching.

In summary, the investigations in the literature can be categorized into two baskets. (i) Study of unidirectional forcing on a self-excited system (oscillator). The system can be a hydrodynamic (Li & Juniper 2013a; Emerson & Lieuwen 2015) or a thermoacoustic (Balusamy et al. 2015; Kashinath et al. 2018) system. Extensive studies on the characterization of lock-in boundary (Li & Juniper 2013a,b,c), exploration of an open-loop control strategy through asynchronous quenching were performed (Guan et al. 2019a; Mondal et al. 2019; Kushwaha et al. 2022). (ii) Study of lock-in due to the interaction between two self-excited systems. Self-excited thermoacoustic interaction in a vortex shedding combustor falls into this category. Combustor acoustic field and vortex shedding process form the systems. Existing literature (for example, Poinsot et al. 1987; Chakravarthy et al. 2007; Singh & Mariappan 2019) dealing with the mutual interaction of two self-excited systems shows that vortex-acoustic lock-in is accompanied by large amplitude oscillations.

The following questions arise concerning vortex shedding combustors. (1) It is unclear whether suppression of instability (similar to unidirectionally forced oscillators) exists in thermoacoustic systems involving two coupled self-excited oscillators. Recently, mutual synchronization between two STS has been studied using experiments, and a comparison is presented with two coupled Van der Pol oscillators (Guan et al. 2021). They have reported that the coupled model exhibits several behaviours observed in the experiments: quasiperiodicity, lock-in, amplitude death, etc. Although they have considered swirling flow that can exhibit vortex shedding, they did not explicitly address the existence and influence of vortex shedding in their systematic study. (2) Is the common frequency during lock-in always close to the combustor acoustic mode?. We make an attempt to address the above questions through theoretical analysis and performing comparisons with our experimental results (Singh & Mariappan 2019). Furthermore, dynamical bifurcations leading to lock-in are identified.

We consider a Rijke tube combustor with the flame stabilized by a bluff body: a geometry experimentally investigated in Singh & Mariappan (2019). The bluff body supports vortex shedding, modulating the heat release rate of the flame. Pressure time series show frequency peaks associated with acoustic field and vortex shedding modes. As the reactant flow rate varies, the modes interact and vortex-acoustic lock-in is observed. In theory, the acoustic field is modelled by a wave equation, with the heat source appearing from the vortex shedding process. The lower-order model for vortex shedding from Matveev & Culick (2003) is used. The continuous-time domain equation is converted to discrete maps, following the procedure discussed in our earlier papers (Britto & Mariappan 2019, 2021). Lock-in is identified as a phase synchronization between vortex shedding events and acoustic velocity fluctuations. For the parameter range associated with the experiments, we predict the occurrence of two lock-in orders: $1:1$ and $2:1$. We show that the former and latter lead to combustion instability and amplitude suppression, respectively. Furthermore, the common frequency during lock-in is close to natural acoustic and vortex shedding modes for $1:1$ and $2:1$ lock-in, respectively.

The rest of this paper is divided into three sections. Governing equations for vortex shedding and the acoustic field, followed by the construction of the dynamical map, are discussed in § 2. Section 3 describes the numerical results, identifying the occurrence of phase lock-in and the associated bifurcations. An interpretation of instability and suppression of vortex-driven thermoacoustic oscillations, in relation to the experimental results, are presented § 4. Section 5 ends the paper with the salient conclusions.

2. Theoretical formulation

Theoretical analysis is performed on the experimental configuration of Singh & Mariappan (2019). A premixed flame (liquefied petroleum gas and air) anchored behind a circular disc bluff body shedding vortices is confined in a Rijke tube (refer to figure 1a). The confinement sustains the acoustic field. At low air flow rates, frequency peaks associated with vortex shedding and acoustic modes are observed. The peaks are wide apart, with the frequency of the vortex shedding mode lower than that of the acoustic mode. As the flow rate increases, the frequency of the vortex shedding mode approaches the acoustic mode, eventually culminating at vortex-acoustic lock-in (figure 8a). A sequence of well-defined regimes was identified during the above dynamic transition (discussed in § 4).

Figure 1. Schematic of vortex shedding behind a circular disc bluff body, placed in a Rijke type combustor of length $\tilde {l}_a$. (a) Vortex shedding/convection process and steady-state CH$^*$ chemiluminescence (representative of flame heat release rate) are illustrated in the left and right half of the bluff body, respectively. Chemiluminescence contours are from our earlier experimental investigation (Singh & Mariappan 2019). Here $\tilde {x}_{sep}$ represents the axial location of vortex shedding. (b) The axial variation of transversely integrated heat release rate shows a peak at $\tilde {x}_c$. The vortex convection distance, $\tilde {l}_c=\tilde {x}_c-\tilde {x}_{sep}$. (c) Schematic representation of shedding (blue dots) and burning (red dots) time instances of vortices. Subscripts in $t$ indicate the vortex number, while the superscripts $s$ and $c$ specify shedding and burning events, respectively.

A schematic of the experimental geometry and our theoretical representation are illustrated in figure 1. As said earlier, vortex shedding combustors have two individual oscillators: vortex shedding and the acoustic field. The following subsections discuss their modelling. The presented modelling framework is similar to Matveev & Culick (2003) and the bifurcation analysis follows our previous investigation (Singaravelu & Mariappan 2016).

2.1. Modelling vortex shedding process

The boundary layer flowing over a sharp corner (separation edge) creates a linearly unstable shear layer, which transitions (Bénard von-Kármán instability) to large coherent vortical structures. Such structures are observed behind the wake of the bluff body. The separation edge of the bluff body is located at $x_{sep}$ (figure 1a). The evolution of the shedding process is described by an integrate-and-fire type model for the vortex circulation $\tilde {\varGamma }_m$ as follows (Matveev & Culick 2003):

(2.1a)\begin{gather} \frac{\text{d}\tilde{\varGamma}_m}{\text{d} \tilde{t}}=\frac{1}{2}\tilde{u}_{sep}^2, \end{gather}

(2.1b)\begin{gather}\tilde{\varGamma}_{sep}=\frac{\tilde{d}}{2St} \tilde{u}_{sep}. \end{gather}

Here $({\cdot })_{sep}=({\cdot })(x_{sep},t)$ and $\tilde {u}_{sep}=\tilde {u}(\tilde {x}_{sep},\tilde {t})=\bar {\tilde {u}}+\tilde {u}_{sep}'$ is the total flow velocity just upstream of the separation edge. Here $\tilde {d}$ represents the characteristic length of the bluff body: in this case its diameter. Strouhal number of the flow configuration is $St$. Furthermore, $\ \tilde{}\ $, $\ \bar{}\ $ and $'$ indicate dimensional, steady and fluctuating quantities, respectively. Subscript $m$ is associated with the $m$th vortex. At the beginning ($t=0$), $\tilde {\varGamma }$ is set to zero. It increases according to (2.1a), till its numerical value equals $\tilde {\varGamma }_{sep}$ (2.1b). At this moment, the $m$th vortex is shed with the circulation $\tilde {\varGamma }_{sep}$ and $\tilde {\varGamma }_{m+1}$ is reset to zero for the evolution of the next $(m+1)$th vortex. In the absence of acoustic perturbations $(\tilde {u}_{sep}'=0)$, (2.1a) and (2.1b) are solved to obtain the steady-state vortex shedding frequency, $\tilde {f}_{s0}= \bar {\tilde {u}} St /\tilde {d}$. The presence of acoustic fluctuations $(\tilde {u}_{sep}')$ alters the instantaneous shedding frequency. The unsteady component, $\tilde {u}_{sep}'(t)$, appears from the solution of the acoustic wave equation (discussed in § 2.2).

The bluff body anchoring the flame is enclosed in a Rijke-type combustor of length $(\tilde {l}_a)$ that facilitates acoustic feedback. The length $\tilde {l}_a$ is kept a constant value (1 m) in Singh & Mariappan (2019). We therefore use the acoustic length $(\tilde {l}_a)$ and time $(\tilde {l}_a/\bar {\tilde {c}})$ scales to non-dimensionalize (2.1a) and (2.1b). Steady-state flow velocity $(\bar {\tilde {u}})$, which is assumed to be spatially constant, is taken as the reference scale for velocity fluctuations. Therefore, we have the following:

(2.2a–e)\begin{equation} \varGamma_m=\frac{\tilde{\varGamma}_m}{\bar{\tilde{u}} \tilde{l}_a};\quad \varGamma_{sep}=\frac{\tilde{\varGamma}_{sep}}{\bar{\tilde{u}}\tilde{l}_a}; \quad u'=\frac{\tilde{u}'}{\bar{\tilde{u}}};\quad t=\frac{\tilde{t} \bar{\tilde{c}}}{\tilde{l}_a};\quad d=\frac{\tilde{d}}{\tilde{l}_a}. \end{equation}

Consequently, the non-dimensional governing equations for vortex shedding become

(2.3a)\begin{gather} \frac{ \text{d}{\varGamma}_m(t)}{\text{d} {t}} =M\frac{1}{2}\left[1+u'_{sep}(t)\right]^2, \end{gather}

(2.3b)\begin{gather}\varGamma_{sep}(t)=\frac{d}{2 St} \left[1+u'_{sep}(t)\right], \end{gather}

where $M$ is the steady-state flow Mach number. In the absence of acoustic perturbations ($u'=0$), vortices are shed at the natural frequency, $f_{s0}=M St /d$. Hence, the vortex shedding model represents a self-sustained oscillator. The shed vortex is assumed to travel at a constant velocity: a fraction $(\alpha )$ of the steady-state velocity ($\alpha =0.4-0.6$, Dotson, Koshigoe & Pace 1997). The assumption greatly simplifies the analytical expressions (refer (2.11)). The next step is to model the acoustic fluctuations and their source.

2.2. Modelling acoustic field of the combustor

The non-dimensional acoustic momentum and energy equations (at low Mach number and constant steady-state density) in the presence of unsteady heat release fluctuations are (Juniper & Sujith 2018)

(2.4a)\begin{gather} \frac{\partial u'}{\partial t}+\frac{\partial p'}{\partial x}=0, \end{gather}

(2.4b)\begin{gather}\frac{\partial p'}{\partial t}+\frac{\partial u'}{\partial x}= \frac{\gamma-1}{\gamma M} \dot{q}'(x,t), \end{gather}

where $p'(x,t)=\tilde {p}'/\gamma M \bar {\tilde {p}}$, $u'(x,t)=\tilde {u}'/\bar {\tilde {u}}$ and $\dot {q}'(\tilde {x},t)=\dot {\tilde {q}}'/\bar {\tilde {p}} \bar {\tilde {c}}$ are the pressure, velocity and heat release rate fluctuations per unit length, respectively. Furthermore, $\gamma$ is the specific heat capacity ratio and $\bar {\tilde {p}}$ is the spatially uniform steady-state pressure. In the present paper, $\dot {q}'(x,t)$ occurs due to perturbations in the flame caused by the movement of vortices.

In diffusion flames, the mixing rate of fuel and oxidizer linearly depends on the circulation of the vortex (Cetegen & Sirignano 1990). The interaction of vortices with flames enhances reaction rate (Cetegen & Basu 2006). The enhancement is proportional to the circulation of the vortex (Peters & Williams 1989). Therefore, we assume that the heat release rate fluctuation $\dot {q}'(x,t)$ is directly proportional to the circulation of the shed vortex $\varGamma _{sep}$. Furthermore, as the vortex travels over the flame, $\dot {q}'(x,t)$ occurs non-locally both in space and time (figure 1b). Modelling it significantly complicates the analysis. Therefore, we assume that $\dot {q}'$ occurs instantly in space and time when the vortex reaches the location $(x_c)$ of the maximum steady-state heat release rate. Such a spatially compact assumption is accepted as the ratio of the flame zone to acoustic length is ${\sim }10^{-3}$ as reported in Singh & Mariappan (2019).

Given the above arguments, the heat release rate fluctuation $\dot {q}'(x,t)$ is modelled by

(2.5)\begin{equation} \dot{q}'(x,t)= \beta \sum_m \varGamma_m \delta (t-t_m^c) \delta (x-x_c),\end{equation}

where $\beta$ is a proportionality factor, termed the heat release rate coefficient; $\delta$ is the Dirac-delta function, representing the compactness of heat release rate fluctuations both in space and time; $t_m^c$ is the time instant, when the $m$th vortex reaches $x_c$.

With the assumption of constant convection velocity, all the shed vortices take the same travel time ($\tau _t$) to move from $x_{sep}$ to $x_c$, given by $\tau _t=l_c/M\alpha$, where $l_c=x_c-x_{sep}$. The location of $x_c$ is obtained from the experiments (Singh & Mariappan 2019) as the axial location of maximum transversely integrated steady-state heat release (figure 1a,b). Later from § 2.3, we show that $\tau _t$ determines the order of the dynamical map formed by coupling the evolution equation for vortex shedding and the acoustic field.

We use the Galerkin technique to transform the set of partial differential equations, (2.4a) and (2.4b), to a second-order ordinary differential equation, by using the following mode shapes for $u'$ and $p'$:

(2.6a,b)\begin{equation} u'(x,t)=\frac{U_n(t)}{{n {\rm \pi}}} \frac{\text{d}\sigma_n(x)}{\text{d}\kern0.7pt x},\quad p'(x,t)= P_n(t) \sigma_n(x). \end{equation}

Here, $\sigma _n(x)$ is the $n$th mode shape, satisfying the boundary conditions of the duct. Since the Rijke tube employed in Singh & Mariappan (2019) is open to the atmosphere, we choose an acoustically open-open boundary condition: $\sigma _n(x)=\sin ({n {\rm \pi}} x)$, where $n$ is a positive integer representing the acoustic mode number. Motivated by the experimental observation, where the first acoustic mode is dominant, we choose a one-mode approximation by setting $n=1$ in the current study. The one-mode approximation also promotes analytical tractability. Although we set $n=1$, for the simulations in this paper, we derive the dynamical equations for a general $n$. On substitution of (2.6a,b) in (2.4a) and (2.4b), a second-order differential equation for $U_n$ is obtained,

(2.7)\begin{equation} \ddot{U}_n+2\zeta_n \omega_n \dot{U}_n+({n {\rm \pi}})^2 U_n={-}2 {n {\rm \pi}} \frac{(\gamma-1)}{\gamma M} \beta \sin ({n {\rm \pi}} x_c) \sum_m \varGamma_m \delta(t-t_m^c). \end{equation}

The above non-dimensional equation represents a kicked harmonic oscillator (Matveev & Culick 2003; Tuwankotta & Ihsan 2019). The Galerkin coefficients for pressure are obtained as $P_n=-\dot {U}_n/{n {\rm \pi}}$. Furthermore, a linear damping term ($\zeta _n \omega _n$) has been introduced in (2.7). The damping coefficient $\zeta _n$ is dependent on the frequency: $\zeta _n=(c_1 n+c_2 \sqrt {1/n})/2{\rm \pi}$, where $c_1$ and $c_2$ correspond to the end and boundary layer losses, respectively (Sterling & Zukoski 1991; Matveev & Culick 2003).

Due to the damping term ($\zeta _n$), the acoustic velocity and pressure fluctuations vanish when the heat release rate due to vortex shedding is absent in the combustor. Hence, the model for the acoustic field does not represent a self-sustained oscillator. As mentioned in the introduction (§ 1), acoustic modes are excited in turbulent combustors, which can be modelled as an additional noise term in the above (2.7) (Sujith & Unni 2020). However, the addition leads to the loss of analytical tractability (§ 2.3) and interpretation (§ 2.5) of instability and amplitude during lock-in. Therefore, we avoid including a noise term in (2.7). Given the above, the theory of mutual synchronization to study lock-in in the present thermoacoustic system cannot be applied in a strict sense. However, we still borrow the concept of phase synchronization to define lock-in in the current case. The reason is as follows. Acoustic fluctuations follow the equation of a (damped) harmonic oscillator, indicating a definite centre of rotation in the phase space ($u'\unicode{x2013}p'$ plane, refer to figure 3a). Therefore, a generalized phase $(\psi _A)$ can be defined using one of the dynamical variables: $U_n$. In turn, a definite relationship between $\psi _A$ and the phase of vortex shedding events allows us to identify lock-in (details in § 3). In the present formulation, acoustic fluctuation continues to exist (although damped) at its natural frequency all through the simulation, as it is generated by the kicks of the (self-sustained) vortex shedding process.

In the following section we convert the continuous-time dynamical system to a discrete dynamical map connecting $U_n,~P_n$ between two successive burning and earlier shedding events. The map allows us to identify and study the bifurcation characteristics of lock-in.

2.3. Construction of a dynamical map

The kick event (right-hand side of (2.7)) occurs at the vortex burning instant ($t_m^c$). It adds acoustic energy instantly, leading to a sudden increase (kick) in pressure fluctuations (velocity fluctuations are continuous), i.e. (Matveev & Culick 2003)

(2.8a)\begin{gather} U_n(t_{m+}^c)=U_n(t_{m-}^c), \end{gather}

(2.8b)\begin{gather}P_n(t_{m+}^c)=P_n(t_{m-}^c)+P_{kick,m}, \end{gather}

(2.8c)\begin{gather}P_{kick,m}=\left(\frac{\gamma-1}{\gamma}\right) \frac{\beta}{f_{s0}} \sin({n {\rm \pi}} x_c)\left[ 1+ \cos({n {\rm \pi}} x_{sep}) U_n(t_{m}^s) \right], \end{gather}

where subscripts $m-$ and $m+$ indicate the time instants just before and after the $m$th kick, respectively. The amplitude of the $m$th kick ($P_{kick}$) depends on the circulation strength of the $m$th vortex through the term $U_n(t_m^s)$, the heat release coefficient $\beta$ and the natural shedding frequency $f_{s0}$.

Between the two successive kicking instances ($t_{m}^c$ and $t_{m+1}^c$), the right-hand side of (2.7) is zero, forming a damped harmonic oscillator. Therefore, in between the kicks, the solution of (2.7) is given by

(2.9a)\begin{gather} U_n(t)=\exp({- {n {\rm \pi}} \zeta_n \tau_{m}^c})\left[ r \sin\left(\omega_d \tau_{m}^c + \alpha_1 \right)\right], \quad t_{m+}^c< t< t_{(m+1)-}^c, \end{gather}

(2.9b)\begin{gather}P_n(t)=\exp({- {n {\rm \pi}} \zeta_n \tau_{m}^c})\left[ r \sin\left(\omega_d \tau_{m}^c + \alpha_2 \right)\right], \quad t_{m+}^c< t< t_{(m+1)-}^c , \end{gather}

where $\omega _d={n {\rm \pi}} \sqrt {1-\zeta _n^2}$. $\tau _{m}^c=t-t_{m+}^c$ indicates the elapsed time, since the burning of the $m$th vortex. The above expression for $U_n$ and $P_n$ is valid till the burning of the $(m+1)$th vortex. An $m$th kick event characterized by the jump in $P_n(t_{m}^c)$, governs the subsequent evolution of $U_n$ and $P_n$ till the $(m+1)$th kick event. Here $\alpha _1, \alpha _2$ and $r$ represent the initial phases and amplitude of $U_n$ and $P_n$, respectively. They are determined by the Galerkin coefficients at the end of the $m$th kick. They serve as an initial condition for the evolution of $U_n$ and $P_n$ after the $m$th kick. We have

(2.10a)\begin{gather} r= \sqrt{A_{1}^2+B_{1}^2}, \quad \alpha_{1,2}=\arctan{\frac{A_{1,2}}{B_{1,2}}}, \end{gather}

(2.10b)\begin{gather}A_1= U_n(t_{m+}^c), \quad B_1=\frac{\zeta_n U_n(t_{m+}^c)-P_n(t_{m+}^c)}{\sqrt{1-\zeta_n^2}}, \end{gather}

(2.10c)\begin{gather}A_2=P_n(t_{m+}^c), \quad B_2= \frac{U_n(t_{m+}^c)-\zeta_n P_n(t_{m+}^c)}{\sqrt{1-\zeta_n^2}}. \end{gather}

Combining the expressions for the evolution of $U_n,~P_n$ in between ((2.9a) and (2.9b)) and across ((2.8a)–(2.8c)) the kicks, a map relating the Galerkin co-coefficients just after the burning of the $(m-1)$th and $m$th vortices are obtained as

(2.11a)\begin{align} &U_n(t_{m+}^c)=\exp({-{n {\rm \pi}} \zeta_n \Delta t_m^c}) \nonumber\\ &\quad\times \left[U_n(t_{(m-1)+}^c) \cos(\omega_d \Delta t_m^c)+\frac{\zeta_n U_n(t_{(m-1)+}^c)-P_n(t_{(m-1)+}^c)}{\sqrt{1-\zeta_n^2}}\sin(\omega_d \Delta t_m^c)\right], \end{align}

(2.11b)\begin{align} &P_n(t_{m+}^c)=\exp({- {n {\rm \pi}} \zeta_n \Delta t_m^c})\nonumber\\ &\quad\times \left[P_n(t_{(m-1)+}^c) \cos(\omega_d \Delta t_m^c)+\frac{ U_n(t_{(m-1)+}^c)-\zeta_n P_n(t_{(m-1)+}^c)}{\sqrt{1-\zeta_n^2}}\sin(\omega_d \Delta t_m^c)\right] \nonumber\\ &\quad+P_{kick,m}, \end{align}

where $\Delta t_m^c=t_m^c-t_{m-1}^c$ is the time period between the $m$th and $(m-1)$th vortex burning events. Since the vortex convection velocity is assumed to be constant (§ 2.1), $\Delta t_m^c=\Delta t_m^s$, indicating that the burning and shedding time periods are the same for an $m$th vortex. Here $\Delta t_m^s$ is calculated by solving (2.3), using the expression for $u'_{sep}$ from (2.9a) (details are in next subsection).

It should be noted that the amplitude of the $m$th kick ($P_{kick,m}$) depends on $U_n(t_m^s)$ (refer to (2.8c)). A (constant) time lag between the burning $(t_m^c$) and shedding $(t_m^s$) instances of the $m$th vortex occurs due to the vortex convection time ($\tau _t$). Considering $k_m$ vortices are shed in the time interval $\tau _t$, the initial condition (2.10) for $U_n(t_m^s)$ (during the evolution of the $m$th vortex) depends on the burning of the $(m-k_m)$th vortex. Here $k_m$ is the greatest integer that satisfies the condition $\sum _{j=0}^{k_m} (\Delta t_{m-j}) \le \tau _t$. Therefore, $P_{kick,m}$ depends on the velocity fluctuations determined by the burning of the $(m-k_m)$th vortex. Therefore, (2.11a) and (2.11b) form a $k_m$th order coupled dynamic map, where $k_m$ depends on the $m$th vortex. An expression for $U_n(t_m^s)$ is given as

(2.12)\begin{align} &U_n(t_{m+}^s)= \exp({- {n {\rm \pi}} \zeta_n \Delta t_{m-k_m}^s}) \nonumber\\ &\quad \left[U_n(t_{(m-k_m)+}^c) \cos(\omega_d \Delta t_{m-k_m}^s)+\frac{\zeta_n U_n(t_{(m-k_m)+}^c)-P_n(t_{(m-k_m)+}^c)}{\sqrt{1-\zeta_n^2}}\sin(\omega_d \Delta t_{m-k_m}^s)\right], \end{align}

where $\Delta t_{m-k_m}=t_m^s-t_{m-k_m}^c$ is the time interval between the shedding of the $m$th vortex and burning of the $(m-k_m)$th vortex. To summarize, (2.8c), (2.11) and (2.12) form a coupled $k_m$th order map. Here $k_m$ depends on the shedding time periods of vortices between the $(m-k_m)$th and $m$th vortices. The unknown shedding time periods $\Delta t_m$s (and, therefore, $k_m$) are found by solving another map, discussed in the next subsection.

2.4. Dynamical map to determine shedding instances

Acoustic velocity fluctuations at the separation edge, $u_{sep}'= U_n(t) \cos ({n {\rm \pi}} x_{sep})$, determines the vortex shedding time instance $t_m^s$ (therefore, the instantaneous shedding period, $\Delta t_m^s$) and its circulation $\varGamma _{sep}$. We substitute the expression for $u'_{sep}$ in (2.3), evaluating the integral to obtain a map relating successive vortex shedding instances $t_{m-1}^s$ and $t_{m}^s$,

(2.13)\begin{equation} \int_{t_{m-1}^s}^{t_m^s} \left[1+U_n(t) \cos({n {\rm \pi}} x_{sep})\right]^2\, \text{d} t - \frac{1}{f_{s0}} \left[1+U_n(t_m^s) \cos({n {\rm \pi}} x_{sep})\right]=0. \end{equation}

The above integral must be evaluated in chunks, as $r$ and $\alpha _{1}$ in the expression for $U_n$ (refer to (2.9a) and (2.10a)) changes after a burning event. Consider $b$ number of burning events occur between $t_{m-1}^s$ and $t_m^s$. The integral is divided into $b+2$ chunks as

(2.14)\begin{align} &\int_{t_{m-1}^s}^{t_m^s} \left[1+U_n(t) \cos({n {\rm \pi}} x_{sep})\right]^2 \,\text{d} t=\int_{t_{m-1}^s}^{t_{m-k_m}^c} \left[1+U_n(t) \cos({n {\rm \pi}} x_{sep})\right]^2 \,\text{d} t \nonumber\\ &\qquad + \sum_{j=m-k_m}^{m-k_m+b}\left( \int_{t_{j}^c}^{t_{j+1}^c}\left[1+U_n(t) \cos({n {\rm \pi}} x_{sep})\right]^2\,\text{d}t\right) \nonumber\\ &\qquad + \int_{t_{m-k_m+b}^c}^{t_{m}^s}\left[1+U_n(t) \cos({n {\rm \pi}} x_{sep})\right]^2 \,\text{d}t, \end{align}

where $j=m-k_m,m-k_m+1,\ldots (m-k_m+b)$ is an integer, with $0\le b\leq (k_m-1)$. The value of $b$ is the greatest integer that satisfies the condition $\sum _{b=1}^{k_m-1} ( \Delta t_{m-k_m+b}^c) \le \Delta t_m^s$, where $\Delta t_{m-k_m+b}^c=t_{m-k_m+b}^c-t_{m-k_m+b-1}^c$ and $\Delta t_m^s=t_m^s-t_{m-1}^s$ are the burning and shedding time periods of the $(m-k_m+b)$th and $m$th vortices, respectively. Now in each individual chunk, $\alpha _{1}$ (in $U_n$) remains the same, allowing us to perform the integral analytically. The expression for the implicit map (2.14) is given in (B1).

Figure 2 shows the shedding and burning instances, indicating the extraction of the instantaneous vortex shedding time period $\Delta t_m^s$, when there is no burning event or the presence of a single burning event between two successive shedding instances. The evolution/reset of $\varGamma _m$ (black curve) and shedding circulation $\varGamma _{sep}$ (blue curve) are illustrated. Vortex shedding and burning events are marked using solid blue circle and red cross symbols, respectively. We set at $t=0$, $\varGamma _0=U_n=P_n=0$. The circulation $\varGamma _0$ increases till it reaches $\varGamma _{sep}$ (at $t^s_1$), leading to the shedding of the first vortex. The shed vortex travels at a constant velocity and burns at time $t_1^c=\tau _t+ t_1^s$ (red cross symbol), where $\tau _t=l_c/(M\alpha )$ is the vortex travel time. Subsequently, at the separation edge, the circulation of the second vortex begins to increase. The initial condition (after the shedding of the first vortex) for $(U_n,P_n)$ still remains $(0,0)$ as there is no burning event till the shedding time $t_2^s< t_1^c$. In this case, the burning event does not occur till the shedding of the third vortex. Vortex shedding occurs at its natural time period, $\Delta t_{s0}=d/M St$.

Figure 2. Evolution of circulation strength $\varGamma _m$ and circulation threshold $\varGamma _{sep}$. Blue solid circle and red cross symbols represent vortex shedding and burning events, respectively. Arrows indicate the evolution and reset of $\varGamma _m$ as the vortices are formed and shed. Burning (shown as a dashed red vertical line) of the first vortex $(t_1^c)$ occurs during the evolution of the fourth vortex ($k_4=3$). Therefore, $\varGamma _m$ shows a slope discontinuity at $t_1^c$.

During the formation of the fourth vortex, the burning event of the first vortex occurs. The burning event (marked as a red vertical dashed line at $t_1^c$) creates the kick given by the map equations, (2.11) and (2.12). It therefore alters the subsequent (after $t_1^c$) evolution of $U_n,P_n$ through the initial condition (2.11). To identify the shedding instance of the fourth vortex, the integral in (2.14) is separated into two sub-integrals with limits, $t_3^s-t_{1-}^c$ and $t_{1+}^c-t_4^s$. In the first integral, $U_n(t),P_n(t)$ are evaluated using the initial condition $U_n(0)=P_n(0)=0$. For the second integral, the initial condition becomes $U_n(t_{1+}^c), P_n(t_{1+}^c)$, which in turn are evaluated from (2.11) and (2.12). The discontinuity in $P_n$ at $t_{1}^c$ translates as a slope discontinuity in the evolution of $\varGamma _4$ (figure 2). After the first kick, the evolution of $\varGamma _4$ and, therefore, the instantaneous vortex shedding time period is altered from the natural period. The alteration occurs for subsequent shedding events due to the presence of successive kicks.

Geometric and flow parameters used in the study are presented in table 1. The parameters $x_{sep},x_c,d,\gamma$ and $St$ are chosen from the experimental configuration of Singh & Mariappan (2019). The vortex convection coefficient ($\alpha$) is assumed to be 0.6 (Dotson et al. 1997). We use the values of damping coefficients $c_1,c_2$ from Matveev & Culick (2003). In Singh & Mariappan (2019) the flow rate is increased to cause a rise in the natural vortex shedding frequency while keeping the length of the Rijke tube (acoustic length) constant. Since the occurrence of vortex-acoustic lock-in strongly depends on the relative values of the natural acoustic and vortex shedding frequency, we vary the frequency ratio $f=f_{a0}/f_{s0}$ ($\,f_{a0}=n/2$ is the natural frequency of the $n$th acoustic mode) in our simulations. At this point, the heat release coefficient $\beta$ cannot be evaluated from the experiments or existing literature. We, therefore, choose $\beta$ as another control parameter. Simulations are performed by varying $(f,\beta )$ in the range $(0.4\unicode{x2013}2.5,\,0\unicode{x2013}0.75)$. The results are discussed in relation to the experiments in §§ 3 and 4.

Table 1. Simulation parameters chosen for the study.

Before we end, the geometry of the evolution of the obtained dynamical map in phase space is studied. It illustrates the possibility that kicks during the burning of vortices can lead to both an increase and decrease in the amplitude of the oscillation, which in turn is related to Rayleigh criteria.

2.5. Phase space representation

A phase space representation of the coupled map equations presented in § 2.4 is shown in figure 3. Acoustic velocity and pressure fluctuations govern the successive vortex shedding and burning events. Therefore, the dynamics of the thermoacoustic system are represented in $U_n,P_n$. Following (2.9) in between the kicks, the phase difference between $U_n,P_n$ equals $\alpha _2-\alpha _1=\alpha _s-{\rm \pi} /2=\arctan (\zeta _n/\sqrt {1-\zeta _n^2})-{\rm \pi} /2$. Eliminating $t$ from (2.9), leads to $P_n^2+U_n^2-2 P_n U_n \sin \alpha _s =r^2 \exp ({-2 {n {\rm \pi}} \zeta _n \tau _m^c}) \cos ^2 \alpha _s$, representing the phase space trajectory. In the absence of damping, $\alpha _s=0$, the above equation represents a circle, with $\alpha _2-\alpha _1=-{\rm \pi} /2$, indicating that the trajectories evolve in the clockwise direction. For $\alpha _s={\rm \pi} /2 \text { and } 3{\rm \pi} /2$ $(\zeta _n\neq 0)$, the above equation represents a straight line ($P_n=\pm U_n$). The phase trajectory is an ellipse for other values of $\alpha _s$ (Cundy & Rollett 1989). The damping term $\zeta _n$ causes an exponential decrease in the size of the ellipse with time. A typical phase trajectory is shown in panel (a).

Figure 3. (a) Schematic illustration of a typical phase space trajectory in $U_n\unicode{x2013}P_n$ plane. The trajectory is a damped rotated ellipse. Vortex shedding and burning events are represented using blue solid circle and red cross symbols, respectively. Panels (b,c) illustrate the effect of kick events on the amplitude of pressure and velocity fluctuations. Red and green elliptical trajectories indicate an increase (instability) and decrease (amplitude suppression) of the oscillation amplitude resulting from a kick event. Arrows indicate the direction of the evolution.

A discontinuous change in the trajectory occurs during the burning events, which can lead to the growth (instability, panel b) or decay (amplitude suppression, panel c) of the ellipse, depending on the phase of the burning event. Shedding and burning events are represented using solid blue circle and red cross symbols, respectively. During a kick, $U_n$ is continuous, while $P_n$ experiences a discontinuous jump, marked by a horizontal red line in figure 3. The Matveev & Culick (2003) model for vortex shedding is valid only in the absence of a reverse flow at $x_{sep}$. It restricts $U_n<1$, which follows in all our simulations. Due to the forms of (2.3b) and (2.8c), $\varGamma _m$ (therefore, $\dot q'$, refer to (2.5)) and the kick $P_{kick,m}$ in $P_n$ are always positive.

Initially, at $t=0$, $U_n=P_n=0$, implying the trajectory begins at the origin (panel a). At the end of the first kick at $t_1^c$, $P_n$ suddenly jumps to a positive value. Till the next (second) kick (at $t_2^c$), the trajectory moves on the damped ellipse and vortices are shed (solid blue circles in panels b,c) before the second kick event. Depending on the phase of the trajectory at $t_2^c$, the second kick can increase (orange trajectory in panel b) or decrease (green trajectory in panel c) the oscillation amplitude. The above argument is the representation of Rayleigh criteria. As said already, $\dot q'$ is always positive. The Rayleigh index is positive when the pressure fluctuation just before the kick at $x_c$ is positive, $P_n(t_m^{c-})>0$. Therefore, instability and amplitude suppression occurs when the trajectory is in the first/fourth and second/third quadrants, respectively, during a kick. Using the above idea, we show that the current model exhibits both instability and amplitude suppression, followed by a qualitative comparison with the experiments in § 4.2.

3. Identification and characterization of lock-in

For a given $f$ and $\beta$ value, we use the numerical algorithm presented in Britto & Mariappan (2019) to calculate the vortex shedding time period. Numerical simulations are performed for 2000 shedding events. However, only the last 500 sheddings are used for the analyses (to discard the transients). Figure 4 shows the iterates of the instantaneous shedding frequency $1/\Delta t_m^s$ (top panels) and contours of power spectral density (PSD) of the corresponding $p'(x=0.36,t)$ (bottom panels) for (a,e) $f=0.6$, (b,f) $0.8$, (c,g) $0.9$ and (d,h) $1.01$ across $\beta =0-0.75$. Pressure fluctuations are computed at the same axial location, $x=0.36$, as that of the experiments (Singh & Mariappan 2019).

Figure 4. (a–d) Iterates of the instantaneous vortex shedding frequency, $f_{m}=1/\Delta t_{m}^s$ for various $\beta =0-0.75$ at $f=0.6,0.8,0.9$ and $1.1$. Plots (e–h) represent the contour that illustrates the amplitude of PSD of pressure fluctuations $p'$ evaluated at $x=0.36$. The range of kick strength $\beta$, where $f_{m}$ is close (less than 5 % deviation) to the natural vortex shedding frequency $f_{s0}$ and first acoustic mode $f_{a0}$ of the duct are illustrated using green and red double-headed arrows, respectively.

In the absence of thermoacoustic feedback, $\beta =0$, shedding occurs at $f_{s0}$. At the lowest considered frequency ratio $f=0.6$ (panel a), the period-1 solution for the instantaneous shedding frequency, $1/\Delta t_m^s$, occurs for $\beta <0.65$ (marked as a red vertical dashed line). In the next section we show that the period-1 solution corresponds to 1:1 lock-in. The numerical value of $1/\Delta t_m^s$ is close to the natural shedding frequency, $f_{s0}=0.83$. Hence, we identify the region as V-lock-in (‘V’ indicates vortex shedding mode). A higher-periodic solution occurs for $\beta >0.65$. The corresponding PSD contours of $p'(t)$ (panel e) reveal a single peak and multiple peaks for period-1 ($\beta <0.65$) and higher-periodic ($\beta >0.65$) solutions, respectively. Since pressure fluctuations (and not instantaneous vortex shedding frequency) are measured in the experiments (Singh & Mariappan 2019), PSD contours serve as a bridge in the comparative study (refer to § 4).

As we move to the next higher frequency ratio $f=0.8$ (which is closer to 1), at low kick strengths $\beta =0-0.33$ (panel b), iterates of $\Delta t_m^s$ are period-1 (V-lock-in), with the shedding occurring near to its natural vortex frequency $(f_{s0}=0.625)$, similar to that of panel (a). As the kick strength increases, the iterates exhibit a transition at $\beta =0.33$, where the shedding frequency jumps close to the frequency of the first acoustic mode $f_{a0}=0.5$, still retaining the period-1 feature (called A-lock-in, ‘A’ indicating the acoustic mode). At large kick strengths $\beta >0.63$, the shedding transitions out of a period-1 orbit and exhibits higher-periodic behaviour. In the PSD (panel f), the two transitions at $\beta =0.33$ and $\beta =0.63$ are reflected as a shift in the frequency of the dominant peak and the formation of multiple peaks, respectively.

As we move $f$ further closer to 1, ($\,f=0.9$, panels c,g), the range of $\beta$ where A-lock-in occurs, widens. Comparing the last two cases, the transition from V-lock-in to A-lock-in is sudden and smooth for $f=0.8$ and $f=0.9$, respectively. Considering $f=1.1$ (present on the other side of $f=1$), one could not identify separate A- and V-lock-in regions. In all the panels (a–d), the instantaneous shedding begins from $f_{s0}$ (period-1 solution) at $\beta =0$. It approaches $f_{a0}$ (also a period-1 solution) as $\beta$ is increased. The approach is both smooth and sudden, depending on $f$. We show in subsection § 3.2, the reason for sudden/smooth transition. Although there is a shift in the dominant frequency in the PSDs (bottom panels), the peak values almost remain constant, indicating similarity with the frequency-pulling route to synchronization (Balusamy et al. 2015).

Increasing $\beta$ leads to a more effective thermoacoustic coupling, leading to more significant $u'$ fluctuations, affecting the shedding frequency. Therefore, there is a tendency for $1/\Delta t_m^s$ to approach the acoustic frequency as $\beta$ is increased. In the next subsection we show that the period-1 solution lies in the phase lock-in region. Unlike the case of external excitation discussed in Britto & Mariappan (2019, 2021), the shedding frequency associated with lock-in in the present self-excited case is near to either $f_{s0}$ (low $\beta$) or $f_{a0}$ (high $\beta$). Thus, two lock-in (A- and V-lock-in) regions are observed.

3.1. Phase lock-in

Acoustic field and vortex shedding are the two oscillators in the present vortex shedding combustor. The relative evolution of phase associated with the two oscillators allows one to define phase lock-in. Velocity fluctuation at the separation edge ($u_{sep}'$) is chosen as the flow variable for the former, while the shedding instances are chosen for the latter. The above two variables were also considered in our previous study (Britto & Mariappan 2021).

The instantaneous phase $(\psi _A(t))$ of $u_{sep}'(t)$ is obtained using the Hilbert transform (Balanov et al. 2008) $\psi _{A}(t)=\arctan ({H(u_{sep}'(t))}/{u_{sep}'(t)})$, where $H$ indicates the Hilbert transform. Furthermore, the phase instances for vortex shedding are obtained based on the following implication. The evolution of circulation strength $\varGamma _m$ is reset to zero at the end of every vortex shedding. Therefore, the shedding system can be considered to have completed one cycle ($2{\rm \pi}$), when a vortex is shed. Hence, the phase ($\psi _{V,m}$) at the $m$th vortex shedding instance $t_m^s$ is given by $\psi _{V,m}(t_m^s)=2m{\rm \pi}$ (the phase at $t=0$ is set arbitrarily to zero). As discussed in Britto & Mariappan (2021), a $p:q$ phase lock-in is said to occur when the phase difference, $\Delta \psi _m(t_m^s)=q\psi _{A,m}(t_m^s)-p\psi _{V,m}(t_m^s)$, evolves to a constant value. Therefore, a zero difference $\Delta (\Delta \psi _m)$ between the phase differences during successive shedding, $\Delta \psi _m - \Delta \psi _{m-1}=0$, determines the phase lock-in state.

Figure 5 shows the contour of $\Delta (\Delta \psi _m)$, corresponding to $1:1$ (panel a) and $2:1$ (panel b) phase relations between $\psi _A$ and $\psi _V$. Blue regions represent $\Delta (\Delta \psi _m)=0$, indicating phase lock-in. It is observed that $1:1$ and $2:1$ phase lock-in dominates the $\beta \unicode{x2013}f$ plane, in the range $f=0.5\unicode{x2013}1.5$ (region I) and $f=1.5\unicode{x2013}2.5$ (region II), respectively. Higher-periodic and quasiperiodic solutions (unlocked state) occur at large $\beta$, identified as region III. Near $f=1.5$ (where there is a junction of $1:1$ and $2:1$ lock-in boundaries), one observes multiple V-shaped boundaries forming higher-order lock-ins such as $3:2$, $5:3$ and $7:4$ (region IV). Furthermore, below $f=0.5$, we observe a $1:2$ phase lock-in state (marked in panel a). The region has large $\beta$ values in a small range of $f$. Both regions IV and V occur for non-zero $\beta$. Since $1:1$ and $2:1$ lock-in regions occupy a large part and occur in the frequency range explored in the experiments (Singh & Mariappan 2019), the rest of the paper is focused on studying regions I and II. In particular, we focus on understanding (i) the transition from V- to A-lock-in, as $f$ approaches one from either side, (ii) the bifurcations involved in the system and (iii) the interpretation of V- and A-lock-ins in relation to our experiments (Singh & Mariappan 2019).

Figure 5. Contours of the difference in the phase difference, $\Delta (\Delta \psi _m)$ across the $\beta \unicode{x2013}f$ plane. Panels (a,b) show regions of $1:1$ and $2:1$ phase relations, respectively, between vortex shedding and the acoustic field. The blue colour represents the phase lock-in regions, where $\Delta (\Delta \psi _m)=0$.

3.2. Transition from V- to A-lock-in and its relation to unidirectional forcing

The understanding of the transition from V- to A-lock-in is illustrated through figure 6. As discussed in § 3, V- and A-lock-ins are described by the closeness of the shedding frequency $1/\Delta t_m^s$ to $f_{s0}$ and $f_{a0}$, respectively. Panels (a,b) show the contour of $1/\Delta t_m^s$ in the $\beta \unicode{x2013}f$ plane for the frequency range $f=0.5\unicode{x2013}1.5$ ($1:1$ phase lock-in) and $f=1.5\unicode{x2013}2.5$ ($2:1$ phase lock-in region), respectively. The natural acoustic frequency $f_{a0}$ is kept at 0.5, while $f_{s0}$ is varied to achieve the desired range of $f$. Masked white regions exclude $1:1$ and $2:1$ phase lock-in regions.

Figure 6. Panels (a,b) show the contours of $1/\Delta t_m^s$ across $f=0.5\unicode{x2013}1.5$ and $f=1.5\unicode{x2013}2.5$, respectively. The black curve represents $1:1$ and $2:1$ A-lock-in boundaries. Green and pink boundaries indicate the occurrence of Neimark–Sacker (NS) and saddle-node (SN) bifurcations, respectively. Panels (c,d) shows the transition of $1/\Delta t_m^s$ (blue) and the shedding phase $\psi _{A,m}$ (pink) with the frequency ratio $f=0.5\unicode{x2013}1.5$ and $f=1.5\unicode{x2013}2.5$, respectively, at $\beta =0.2$. They illustrate the sudden and smooth transition between V- and A-lock-in regions. The horizontal orange dashed line indicates the frequency of the first acoustic duct mode $f_{a0}$ and its first subharmonic $f_{a0}/2$ in panels (c,d), respectively. Furthermore, the natural vortex shedding frequency $(f_{s0})$ is shown using green dashed lines. Insets in panel (c) show the behaviour of the first return map (between $\psi _{A,m} - \psi _{A,m-1}$, blue curve) across the $1:1$ A-lock-in boundary. The map intersects the $45^{\circ }$ dashed black line at the fixed points. Green and red solid circles represent stable and unstable fixed points, respectively. Panels (e,f) compare the $1:1$ and $2:1$ A-lock-in boundaries, extracted from the present coupled simulations (black curve) and the analytical expression (red dashed curve) described in Britto & Mariappan (2019).

To recollect, the transition from V- to A-lock-in is sudden ($\,f=0.8$, figure 4b) and smooth ($\,f=0.9,1.1$, figure 4c,d) depending on the closeness of $f$ to 1. Therefore, for classification, we identify A-lock-in when $1/\Delta t_m^s$ lies less than ${\pm }5\,\%$ (arbitrary threshold) of $f_{a0}$ and $f_{a0}/2$ for $1:1$ and $2:1$ phase lock-in, respectively. Otherwise, V-lock-in prevails. The black curve in panels (a,b) mark the A-lock-in boundary. From the contours, we observe that $1/\Delta t_m^s$ close to the A-lock-in boundary in the $f<1$ and $f>1$ sides show steep and shallow gradients, respectively.

To scrutinize for the presence of bifurcations, we investigate a representative case, $\beta =0.2$, where the instantaneous shedding time period $1/\Delta t_m^s$ and phase $\psi _{A,m}(t_m^s)$ of $U_n$ are plotted with $f$ (panels c,d). The dashed green line indicates the natural vortex shedding $f_{s0}$. Brown horizontal lines indicate $f_{a0}$ (panel c) and $f_{a0}/2$ (panel d) with which 1:1 and 2:1 A-(phase) lock-in relation occurs, respectively.

In the $f<1$ region we observe that the response ($1/\Delta t_m^s$) occurs at the natural vortex shedding frequency $f_{s0}$ (V-lock-in) when $f$ is outside the black curve. Across the black curve, two types of transition occur as the frequency ratio is gradually increased. (i) A sudden jump in the shedding frequency $1/\Delta t_m^s$ (blue curve) occurring from $f_{s0}$ (V-lock-in) to $f_{a0}$ (A-lock-in) is observed for $f<1$ (panel c). A similar jump in $\psi _{A,m} (t_m^s)$ occurs (pink curve in panel (c) shown using the right-hand side handle of the vertical axis). (ii) On the other hand, for $f>1$, $1/\Delta t_m^s$ and $\psi _{A,m}(t_m^s)$ exhibit a smooth transition (panel c) from $f_{a0}$ to $f_{s0}$ as $f$ is increased. Both the sudden (V- to A-lock-in) and gradual (A- to V-lock-in) occur in the case of $2:1$ lock-in (panels b,d) also.

The insets in panel (c) show the first return maps between $\psi _{A,m} - \psi _{A,m-1}$ (blue curve), which illustrate the variation of stable (green solid circle) and unstable (red solid circle) fixed points across the black boundary. The map intersects the $45^{\circ }$ line (dashed black line) at the fixed points. The number of fixed points remains the same (two), and their stability is unchanged. Therefore, no bifurcation occurs. However, on the $f<1$ side, the change in the shape of the first return map forces a strong movement in the location of the stable fixed point. The corresponding shedding phase, $\psi _{A,m}$, shows a sharp variation, leading to the sudden drop in the vortex shedding frequency $(1/\Delta t^s_m)$. On the other side $(f>1)$, the movement of the stable fixed point is smooth, leading to the smooth transition in $\psi _{A,m}$ and, therefore, the shedding frequency.

Later in § 4.2 we show that the A-lock-in phenomenon provides the most favourable condition for the occurrence of combustion instability. A comparison of the A-lock-in boundary from the present simulation is made with the boundary obtained from a unidirectionally forced vortex shedding system. The acoustic feedback is replaced by an external velocity forcing. An analytical expression for the A-lock-in boundary, described in our earlier work (Britto & Mariappan 2019) is available: $| ({2 u_{rms}^2 + 2 (1-{f}/{p} )})({2 \sqrt {2} u_{rms} f/p})|= 1$, where $p=1,2$ for $1:1$ and $2:1$ A-lock-in, respectively. Here $u_{rms}$ indicates the root-mean-square (r.m.s.) value of the acoustic velocity fluctuations measured at the vortex separation location $x_{sep}$. It is obtained for a given $\beta$ from the coupled simulation. Solving the analytical expression provides two values for $f$, within which A-lock-in occurs. Panels (e,f) show the comparison of the A-lock-in boundaries obtained from the coupled simulation (black curve, same as that in panels a,b) and unidirectionally forced analytical solution (red dashed curve) for $1:1$ and $2:1$ lock-in, respectively. A fair qualitative match in the shape of the boundary is observed. Therefore, the analytical solution can be used as a first design prescription for combustors to avoid A-lock-in and combustion instability.

3.3. Bifurcations leading to lock-in

White regions (termed as the unlocked state) shown in panels (a,b) of figure 6, exhibit vortex shedding time periods, either at higher-periodic or quasiperiodic orbits. Bifurcations occurring during the transition between the phase-locked state and the white region are investigated by sweeping along $f$ for a fixed $\beta$, chosen as ${\beta =0.65}$. For illustration, we choose two transitions occurring at $f=0.628$ and $f=0.831$ (figure 6a). The corresponding Feigenbaum diagrams for $\Delta t_m^s$, displaying the transitions from V-lock-in to the unlocked state $(f=0.6-0.65)$ and the unlocked state to the A-lock-in state $(f=0.76-0.86)$ are shown in figure 7(a,b), respectively. Grey regions spanning $f=0.600-0.627$ and $0.830-0.860$ in panels (a,b) indicate $1:1$ V- and A-lock-in, respectively. Two successive iterates of $\Delta t_m^s$ are plotted to form the phase portrait in the insets for a few specific values of $f$ (vertical dashed lines). The colours in the insets correspond to $f$ at the vertical dashed lines. The phase portrait illustrates the dynamics observed in the Feigenbaum diagrams.

Figure 7. Feigenbaum diagram for $\Delta t_m^s$ across (a) $f=0.6-0.65$ and (b) $f=0.76-0.87$. Niemark–Sacker (NS) bifurcation is observed in panel (a) wherein the shedding exhibits a gradual transition from a period-1 fixed point to a quasiperiodic orbit. The size of the orbit smoothly increases with $f$ as shown in the phase portraits at various values of $f$ in the inset. The chosen frequencies for illustration in the phase portrait are marked using vertical dashed lines having the same colours as in the phase portrait. The transition from a higher-periodic state (black colour in the inset) to $1:1$ A-lock-in state (period-1 orbit, red colour) occurs through a saddle-node (SN) bifurcation as shown in panel (b).

At low values of $f$ (panel a), shedding exhibits a period-1 V-lock-in state, where in the vortices are shed close to their natural shedding frequency $f_{s0}$. The phase portrait (say at $f=0.62$, red line) becomes a point. On increasing $f$ beyond 0.627, the iterates densely fill a finite space in the Feigenbaum diagram. The corresponding phase portraits (say at $f=0.637$, black line) is a closed curve, indicating a quasiperiodic orbit. Furthermore, the size of the orbit increases gradually with a further increase in $f$ (compare portraits of $f=0.637$ black, $f=0.64$ green, $f=0.65$ yellow). There is a birth of ergodic torus (quasiperiodic orbit) from the period-1 fixed point (V-lock-in), indicating the occurrence of Neimark–Sacker bifurcation (marked as NS in figure 6a).

The transition back to a lock-in state on the higher value side of $f$ is illustrated in panel (b). At $f=0.76$ (black), vortex shedding occurs at higher-periodic orbits, leading to discrete (finite) points in the phase portrait. The dynamics occur in a resonant torus. At $f=0.831$, the resonant torus suddenly collapses to become a period-1 fixed point (A-lock-in). The associated phase portrait (say at $f=0.84$, red) becomes a point. The sudden death of the resonant torus indicates the occurrence of saddle-node bifurcation. The corresponding bifurcation boundary is marked in pink with a mention of the saddle node in figure 6(a). In between, the frequency ratio mentioned in the panels of figure 7 ($0.645< f<0.757$), there is a transition from ergodic to the resonant torus, occurring via a saddle-node bifurcation (not shown).

A similar exercise to identify the bifurcation leading to 2:1 lock-in is performed. It should be noted that phase trapping, a behaviour commonly observed in synchronization regions near the onset of lock-in (Li & Juniper 2013c), is not observed in the current simulations. Phase trapping was also not spotted in the external velocity excitation investigation in our previous papers (Britto & Mariappan 2019, 2021). Since acoustic field equations used in the study are linear, phase trapping is not expected in the current paper. In summary, saddle-node and Neimark–Sacker bifurcations are observed during the transition to $1:1$ and $2:1$ phase lock-in. Their locations are marked in figure 6(a,b). The next section discusses our interpretation of the A- and V-lock-in regions in light of the experimental observations of Singh & Mariappan (2019).

4. Reinterpretation of the experimental observations with respect to A- and V-lock-in phenomena

The experiments used a premixed fuel–air mixture to establish a turbulent flame behind a circular disc bluff body, shedding vortices. The fuel flow rate is fixed at 1.75 slpm, while the air flow rate is chosen as the control parameter to take the system from an unlocked to a lock-in state. On the journey of increasing air flow rate, three distinguishable regimes of interest occur, namely: (i) no lock-in, (ii) pre-lock-in and (iii) lock-in. The waterfall plot of PSD corresponding to the three regimes are shown in figure 8(a) using blue (no lock-in), green (pre-lock-in) and red (lock-in) colours.

Figure 8. Comparison of waterfall plots of PSD between the experimental results (panel a) of Singh & Mariappan (2019) and the current simulations (panels (b,c) correspond to $\beta =0.2,0.65$, respectively). The PSD is obtained from the pressure time series at $x=0.36$. The letters V and A indicate the peaks associated with vortex shedding and acoustic modes, respectively. The number 2 before the letters V and A represent the first harmonics. The air flow rates and the corresponding frequency ratio $f$ are represented in each panel.

In panel (a) (experiments), for low air flow rates 18–26 slpm (blue colour), the PSD indicates the existence of two dominant frequencies: vortex shedding (identified with the letter ‘V’) and acoustic (identified with the letter ‘A’) mode frequencies. The frequencies are incommensurate. The phase space reconstruction from the measured acoustic pressure time series shows quasiperiodic orbits (discussed subsequently in figure 9a–j). Since the frequencies of the vortex shedding and acoustic modes are far apart, both modes do not influence each other. Vortex shedding occurs at its natural frequency $f_{s0}$, which increases linearly with the flow rate. The frequency of the first acoustic mode remains (almost) constant at $f_{a0}=0.5$. Since both vortex shedding and acoustic modes are present, we identify it as a no-lock-in regime or unlocked state.

Figure 9. Recurrence plots (a–c and g–i) and their corresponding Poincaré maps (d–f and j–l) generated from the experimental acoustic pressure time series for the flow rates 18–28 slpm. A transition from the no-lock-in to pre-lock-in regime is observed in the above flow rates. Saddle-node bifurcation causes the transition.

On further increasing the flow rate, the combustor advances to the pre-lock-in region (green colour), where the acoustic mode is suppressed: fundamental (V) and first harmonic (2V) of the vortex shedding frequency are present. Furthermore, as the flow rate increases, the vortex shedding frequency approaches the acoustic frequency $(f\rightarrow 0.5)$, leading to lock-in (red colour). During lock-in, fundamental (A) and first harmonic (2A) of the acoustic mode are only present (vortices shed at the acoustic frequency). In the lock-in regime large amplitude pressure oscillations (instability) are observed.

Bifurcations during the transition from one regime to another are examined using recurrence plots and Poincaré maps in figure 9. Flow rates involving the transition from no-lock-in (18 slpm) to pre-lock-in (28 slpm) regimes are shown. The panels in the first (a–c)/third (g–i) and second (d–f)/fourth ( j–l) columns show the recurrence plots and their corresponding Poincaré maps, respectively. Acoustic pressure data from the experiments are Fourier filtered, keeping only the first two dominant modes (based on the amplitude) to remove the noise from the data. Time delay and embedding dimension required to construct the recurrence plots are obtained using average mutual information (Fraser & Swinney 1986) and false nearest neighbours (Abarbanel et al. 1993), respectively. Furthermore, the threshold for the revisit of a phase trajectory to an earlier dynamical state is fixed at 20 % of the size of the attractor (Nair et al. 2014). The Poincaré map is constructed between successive local peaks ($p'_{max}(m)$, $m$ being the index of the peak) in the pressure time series.

Unequal spacing among the diagonals in the recurrence plots for flow rates 18–24 slpm indicates quasiperiodic behaviour. The corresponding Poincaré maps show closed curves. At 26 slpm, the spacing between the diagonal lines becomes constant, displaying periodic orbits. Two distinct points in the Poincaré map show the occurrence of period-2 orbits. At the next higher flow rate (28 slpm), the system transitions to a period-1 limit cycle. Overall, the phase space attractor transitions from an ergodic torus (18–24 slpm) to a resonant torus (26 slpm), culminating in a closed curve (28 slpm). Therefore, the transition from a no-lock-in to a pre-lock-in regime occurs through a saddle-node bifurcation. On the other hand, the period-1 limit cycle continues as the system transitions from a pre-lock-in to lock-in regime (not shown), indicating the absence of a bifurcation.

Experiments show the presence of pre-lock-in and lock-in regimes, while the model displays V-and A-lock-in regimes, all of them displaying period-1 orbits. Our focus is to (i) establish the connection between the regimes observed in the experiments and model, and (ii) demonstrate the occurrence of instability and amplitude suppression during lock-in using the Rayleigh index (§ 4.2).

4.1. Journey to lock-in and its relation to V-, A-lock-in

As mentioned earlier, simulation parameters shown in table 1 are chosen from Dotson et al. (1997), Matveev & Culick (2003) and Singh & Mariappan (2019). The range of air flow rate, 18–58 slpm, in the experiments translates to the frequency ratio in the range $f=2.7\unicode{x2013}0.88$ ($\,f_{s0}$ increases linearly with flow rate, while $f_{a0}$ is kept at 0.5). As $\beta$ cannot be determined, we perform simulations for two values (a low and a high value) of $\beta =0.2,0.65$. Their corresponding waterfall plots for PSD are shown in panels (b,c) of figure 8. Traversing the air flow rate from 18 to 64 slpm, corresponds to the movement of $f$ from right to left in figures 5 and 6(a,b).

The waterfall plots in figure 8(b,c) are read along with figures 5 and 6(a,b) to relate the dynamics to V- and A-lock-ins. For $\beta =0.2$ (figure 8b), at low flow rates (18–22 slpm, marked in green), dominant frequencies are associated only with the vortex shedding mode. The flow regime qualitatively resembles the pre-lock-in regime observed in the experiments (figure 8a). Furthermore, the V-lock-in region has a $2:1$ phase relationship (refer to figure 5). As the flow rate gradually increases (24–28 slpm), the frequency ratio enters the $2:1$ A-lock-in state (red colour in figure 8b). Here, the peak at $f_{s0}$ is absent. Pressure fluctuations occur at $f_{a0}$ and its sub- and super-harmonics. A further increase in the flow rate (30–46 slpm, green colour) transitions the system from a $2:1$ A-lock-in state to a $1:1$ V-lock-in state (absence of the acoustic mode: similar to pre-lock-in region). The $1:1$ A-lock-in occurs at 48 slpm (figure 8b), which resembles the lock-in regime observed in the experiments (figure 8a).

From the qualitative similarity of the presence of peaks in the PSD (in panel (a) with panels b,c), one can connect that pre-lock-in and lock-in regimes observed in the experiments occur due to V- and A-lock-in phenomena, respectively. The connection is further strengthened by the (i) presence of saddle-node bifurcation leading to the lock-in from the unlocked (no-lock-in) state and (ii) absence of bifurcation during the transition between V-and A-lock-in phenomena, observed in both experiments and simulations. However, the amplitudes of the higher harmonics are comparable to the corresponding fundamental peaks in the V- and A-lock-in regions (green and red plots in figure 8b,c). In contrast, the higher harmonics are highly damped in the experiments (figure 8a). The above behaviour is a limitation of the model, discussed further in § 4.4.

The model also does not capture the no-lock-in region of the experiments (blue colour curves of figure 8a). The turbulent combustion noise excites weakly the acoustic modes and, therefore, is recorded in the pressure fluctuations in the experiments. Since the current thermoacoustic model excludes noise, the simulations cannot obtain the no-lock-in region. Furthermore, far from $f=1$, say at $f=2.7$ (where the no-lock-in regime occurs in experiments, refer to figure 8), the thermoacoustic coupling is too weak in generating appreciable acoustic waves to appear in the PSD of figure 8(b,c).

Simulations at large $\beta =0.65$ (panel c) is qualitatively similar to that of $\beta =0.2$ (panel b). However, one observes unlocked states (24–34 slpm, black colour) in between V-lock-in regions. The oscillations are found to be either quasi or higher periodic in the unlocked state (not shown in the paper). In the subsequent section we study the occurrence of thermoacoustic instability and suppression, followed by a comparison with the experiments.

4.2. Amplitude suppression and instability

As indicated in § 2.5, the phase of the vortex burning event with respect the acoustic pressure determines whether the lock-in leads to instability (figure 3b) or amplitude suppression (figure 3c). The present section calculates the Rayleigh index (RI) and identifies the regions of instability and suppression in the $\beta \unicode{x2013}f$ plane. Later, the two regions are identified in the experiments. The Rayleigh index in the non-dimensionalized form becomes

(4.1)\begin{equation} \text{RI }=\frac{1}{T} \int_0^T \int_0^1 p'(x,t) \dot{q}'(x,t)\,\text{d}\kern0.7pt x\, \text{d} t. \end{equation}

The above integral is evaluated over the combustor length and a time $T$, comprising of many oscillation cycles. In the current model pressure and heat release rate oscillations are piecewise continuous functions of time ($t_m^c$) and space ($x_c$) (refer to (2.5) and (2.8)). At the kick events ($t_m^c$) the jump in pressure (given in (2.8c)) and heat release rate fluctuations are proportional to $\varGamma _{sep}=(1+u_{sep}')/2d St$. Since $|u_{sep}'|<1$ in the explored parameter space, $\varGamma _{sep}>0$. Therefore, the jumps in pressure $P_{kick,m}$ and heat release rate $\dot {q}'$ remain positive.

Figure 10(a) shows the contour of $\varGamma _{sep}$ in the $\beta \unicode{x2013}f$ plane. Black curves represent $1:1$ and $2:1$ A-lock-in boundaries. For the parameter values considered in table 1 and in the absence of thermoacoustic feedback ($u_{sep}'(t)=0$), $\varGamma _{sep}=\varGamma _{sep,wfb}$ equals $0.05$. The related burning events result in a constant heat release rate (across the kicks). The presence of acoustic feedback $(\beta \neq 0,|u_{sep}'(t)|<1)$ alters $\varGamma _{sep}$ away from $\varGamma _{sep,wfb}$. The value of $\varGamma _{sep}$ depends on the phase of the velocity fluctuations $u_{sep}'$ at the time of shedding. Since $q'(x,t) \propto \varGamma _{sep}$, the contour in figure 10(a) highlights the variation of the heat release rate $\dot q'$. Yellow/red and green regions represent heat release rates higher and lower than the heat release rate without feedback, respectively. Most higher and lower heat release rate regions overlap with the A- and V-lock-in regions, hinting at the possibility of instability and amplitude suppression, respectively. The observation is incomplete since the analysis should be performed in terms of the Rayleigh index, which determines the actual addition of acoustic energy into the system (discussed later). Note that a large heat release rate (red contour regions) occurs slightly to the left of $f=1$ (red vertical dashed line). The result reinforces the experimental observation of Emerson & Lieuwen (2015), who highlighted a reduction in the heat release rate when a bluff body stabilized flame is excited at its natural frequency $(f_f=f_{s0})$.

Figure 10. Contours of (a) $\varGamma _{sep}$ and (b) Rayleigh index (${\rm RI}_{fb}$) due to thermoacoustic feedback in the $\beta \unicode{x2013}f$ plane. Heat release rate is proportional to $\varGamma _{sep}(t)$ (panel a). Regions where $\varGamma _{sep}(t) < d/2St$ and $\varGamma _{sep}(t) > d/2St$, shown using green and yellow/red shades, indicate heat release lower and higher than the heat release without feedback, respectively. Red/yellow and blue shaded regions in panel (b) indicate instability (${\rm RI}_{fb} > 0$) and suppression (${\rm RI}_{fb} < 0$) of the thermoacoustic oscillations, respectively. A vertical pink line demarcates $1:1$ and $2:1$ V-lock-in boundaries. Panel (c) shows the phase portrait at two close-by parameters, which are on either side of the $1:1$ A-lock-in boundary.

From the above discussion, we observe the jump in heat release rate can be decomposed into two components: heat release rate in the absence ($q'_{wfb}$) of and due $(q'_{fb})$ to acoustic feedback. They are caused due to the steady-state ($\bar {u}=1$) and fluctuating velocity ($u_{sep}'(t)$) components, respectively. Their expressions are

(4.2)\begin{gather} \dot{q}'_{wfb}(x,t)= \frac{d \beta}{2 St} \sum_m \delta(t-t_m^c)~ \delta(x-x_c), \end{gather}

(4.3)\begin{gather}\dot{q}'_{fb}(x,t)=\frac{d \beta}{2 St} \sum_m{u_{sep}'(t_m^s)}\delta(t-t_m^c), \delta(x-x_c) \end{gather}

where $\dot q'=\dot {q}'_{wfb}+\dot {q}'_{fb}$. Note that $\dot {q}'_{wfb}$ is always positive, and therefore, acoustic waves are continuously generated, even in the absence of thermoacoustic feedback. We choose the state as the reference and explore the amplification or attenuation of acoustic waves over the reference due to the thermoacoustic feedback. We therefore calculate the Rayleigh index (${\rm RI}_{fb}$) associated with the heat release due to feedback $(\dot {q}'_{fb})$. Using (4.1), ${\rm RI}_{fb}$ is calculated by replacing $\dot {q}'$ with $\dot {q}'_{fb}$. After substituting the Galerkin expansion for $p'$ and performing the spatial integration, it results in the expression

(4.4)\begin{equation} {\rm RI}_{fb}=\frac{F_1}{T} \sum_{m=1}^{M_{tot}} \left[ u_{sep}'(t_m^s) \right] \int_0^T P_n(t)\delta(t-t_m^c)\,{\rm d}t, \end{equation}

where $M_{tot}$ is the total number of shed vortices in time $T$ (taken as 20 non-dimensional time) and $F_1=\beta d \sin ({n {\rm \pi}} x_c)/(2nSt)$. In the time integral $P_n(t)$ is discontinuous at $t_m^c$, where the Dirac delta $\delta (t-t_m^c)$ function also fires. As suggested by Coutinho, Nogami & Toyama (2009), the average of $P_n(t)$ $(=(P_n(t^c_{m+})+P_n(t^c_{m-}))/2)$ across the $m$th kick is used to evaluate the integral. Therefore, (4.4) results in the expression

(4.5)$$\begin{gather} {\rm RI}_{fb}=\frac{F_1 F_2 \cos({n {\rm \pi}} x_{sep})}{T} \left[\sum_m^{M_{tot}} U_n^2(t_m^s) \cos({n {\rm \pi}} x_{sep}) +\sum_m^{M_{tot}} U_n(t_m^s)\right] \nonumber\\ + \frac{F_1 \cos({n {\rm \pi}} x_{sep})}{T} \sum_m^{M_{tot}} U_n(t_m^s) P_n(t_{m-}^c), \end{gather}$$

where $F_2=\beta f \sin ({n {\rm \pi}} x_c) (\gamma -1)/n\gamma$.

Equation (4.5) implies that ${\rm RI}_{fb}$ (acoustic energy addition) depends on two factors, namely: (i) velocity fluctuations at the separation time instant $U_n(t_m^s)$, translating to the heat release rate; and (ii) pressure fluctuations just before the $m$th kick at the burning location $x_c$. Positive and negative values of ${\rm RI}_{fb}$ correspond to instability and amplitude suppression of the thermoacoustic oscillations, respectively.

Figure 10(b) represents the contour of Rayleigh index (${\rm RI}_{fb}$). Yellow/red and blue contours represent regions where ${\rm RI}_{fb}>0$ (instability) and ${\rm RI}_{fb}<0$ (amplitude suppression), respectively. Two regions of instability are observed: (i) inside the $1:1$ A-lock-in region and (ii) the upper portion of the $2:1$ A-lock-in region. The following general conclusions can be made regarding instability.

1. (a) At low kick strength ($\beta <0.15$), both $1:1$ and $2:1$ A-lock-in regions do not excite the system towards instability (because ${\rm }RI_{fb}\leq 0$, refer to panel b), although the corresponding $(\dot {q}'_{fb})$ (refer to panel a) is positive.

2. (b) Large values of positive ${\rm RI}_{fb}$ is observed for higher $\beta$ in the $1:1$ A-lock-in. Therefore, $1:1$ A-lock-in region is the most favourable region for the occurrence of instability (to be discussed further in figure 11).

3. (c) Near the A-lock-in boundary (for example, $f<1$ side), ${\rm RI}_{fb}$ changes its sign leading to a sudden switch between instability and suppression of oscillation (refer to panel c). Therefore, the evaluation of the A-lock-in boundary is crucial. As shown in § 3.2, the boundary approximately coincides with the lock-in boundary of the unidirectionally excited vortex shedding system (figure 6e,f). The analytical expression for the same is available in Britto & Mariappan (2019, 2021), which can serve as an essential design guideline.

4. (d) Throughout the V-lock-in region, ${\rm RI}_{fb} \le 0$ indicates amplitude suppression and, therefore, the most suitable lock-in state for the combustor operation to avoid thermoacoustic instability.

Figure 11. Comparison of Rayleigh index and r.m.s. value of pressure fluctuations ($P_{rms}$) between the experiments (Singh & Mariappan 2019) (panel a) and the current simulations (panels (b,c) for $\beta =0.2,0.65$, respectively). Left (blue colour) and right (red colour) handles of the vertical axes are used to present ${\rm RI}_{fb}$ ($\cos \theta$ in panel (a) only: $\theta$ is the phase angle between heat release rate and pressure fluctuations at the dominant frequency) and $P_{rms}$ (at $x=0.36$), respectively. Panel (d) shows a flowchart of the dynamical states encountered in experiments and numerical simulations. Arrows indicate increasing flow direction. The comparison is colour coded with respect to figure 8.

4.3. Comparison with Singh & Mariappan (2019)

Figure 11 compares the Rayleigh index and r.m.s. value of pressure fluctuations ($P_{rms}$) between the experiments (Singh & Mariappan 2019) (panel a) and the current simulations (panels (b,c) for $\beta =0.2,0.65$, respectively). Left (blue colour) and right (red colour) handles of the vertical axes are used to present ${\rm RI}_{fb}$ ($\cos \theta$ in panel (a) only: $\theta$ is the phase angle between heat release rate and pressure fluctuations at the dominant frequency) and $P_{rms}$ (at $x=0.36$), respectively.

Considering panel (a) (experiments), at the lowest air flow rate (18 slpm), $\cos \theta$ is negative, leading to low amplitude oscillations. The increase in the air flow rate gradually intensifies the positive thermoacoustic feedback between $p'$ and $\dot q'$ (positive slope of $\cos \theta$ close to the first vertical black dashed line, 26 slpm). Simultaneously, the amplitude of pressure oscillations also increases. The first transition occurs when the combustor enters the pre-lock-in region (above 26 slpm), where in the interaction between $p'$ and $\dot q'$ leads to negative feedback. It marks the event of suppression (shown using green downward arrows). The corresponding $P_{rms}$ decreases slightly. A further increase in the air flow rate (above 30 slpm) intensifies the interaction between $p'$ and $\dot q'$ (blue curve marches upwards and reaches positive values: black upward arrows indicate positive feedback), and the pressure oscillations begin to grow again (red curve). The oscillations are now locked-in with the acoustic mode of the combustor (above 38 slpm, second vertical dashed line). In the lock-in flow regime the combustor exhibits instability ($\cos \theta$ reaches its maximum), and large amplitude pressure oscillations are observed (red curve). Flame blow-off occurs beyond the lock-in regime (52 slpm), leading to a sudden drop in $p'$ and $\dot q'$. As discussed in § 4.1, pre-lock-in and lock-in regimes are associated with V- and A-lock-in phenomena, respectively, through the PSD, recurrence plots and Poincaré maps of acoustic pressure time series. It is imperative to examine the above association regarding the Rayleigh index.

In the current model (panels b,c), the qualitative variation of ${\rm RI}_{fb}$ and $P_{rms}$ with air flow rate is similar to that of the experimental observation (panel a). The A-lock-in regions (both $2:1$ and $1:1$) are accompanied by positive ${\rm RI}_{fb}$ leading to instability. The behaviour is associated with the lock-in regime of the experiments, where $\cos \theta$ is positive with the oscillations occurring close to the natural duct acoustic frequency ($1:1$ A-lock-in).

On the other hand, V-lock-in ($2:1$ and $1:1$) is accompanied by ${\rm RI}_{fb}<0$, indicating amplitude suppression (panels b,c). A similar drop to negative values in $\cos \theta$ (experiments, panel a) is observed in the pre-lock-in regime. It allows one to associate the regime with V-lock-in phenomena based on the Rayleigh index. In a nutshell, lock-in and pre-lock-in regimes observed in the experiments correspond to $1:1$ A-lock-in and $2:1$/$1:1$ V-lock-in regions of the simulations, showing instability and amplitude suppression, respectively. A comparison summary of the observed dynamical states between experiments and the simulations is represented as a flowchart in panel (d). The comparison is colour coded with respect to figure 8. Arrows indicate the direction of the increasing flow rate.

The vortical structures are shed in the wake region behind the bluff body/step, where the time-averaged velocity profile is lower than the upstream flow velocity. It slows the vortex convection speed compared with the upstream flow (Shanbhogue, Husain & Lieuwen 2009). Based on the previous experimental and numerical investigations (Rossiter 1964; Zaman & Hussain 1981; Zhou & Antonia 1992; Kook & Mongeau 2002), the vortex convection speed downstream of a bluff body or in a cavity is found to be between $0.4\unicode{x2013}0.95$ times the free-stream velocity. All the simulations in this paper are performed by assuming the coefficient of vortex convection speed to be $\alpha =0.6$. The convection of the vortex produces a time delay between the origin of its birth to its effect on the heat release rate. Thermoacoustic instability strongly depends on the time delay (Huber & Polifke 2009). We therefore re-ran the simulations with $\alpha =0.9$ (green solid curves) and presented the comparison with $\alpha =0.6$ (original results shown in red dashed curves) in figure 12. Panel (a) compares the $1:1$ A-lock-in boundary. As expected, the convection velocity does influence the lock-in boundary; however, the qualitative nature of the boundary remains the same. Furthermore, panels (b,c) show the calculated Rayleigh index and r.m.s. values of pressure. The frequency ratio and the corresponding A- and V-lock-in regions are illustrated in the figure. The convection velocity does not produce a qualitative difference in the curves. It indicates that the conclusions on amplitude suppression and instability hold for other values of vortex convection speed.

Figure 12. Comparison of vortex convection velocity $\alpha$, between $\alpha =0.6$ (red dashed curve) and $0.9$ (green solid curve). Panel (a) compares lock-in boundaries in the $\beta \unicode{x2013}f$ plane. The parameter $\alpha$ influences the A-lock-in boundary in the $f>1$ region. Panels (b,c) represent the comparison of Rayleigh index ($R_{fb}$) and r.m.s. values of pressure fluctuations ($P_{rms}$), respectively. The V- and A-lock-in regions are illustrated along with the corresponding frequency ratios.

4.4. Limitations of the model

The current model has its limitations, resulting in an inaccurate match with the experiments, which are explained through figures 8 and 13. The first one is the inability of the model to reproduce the no-lock-in region (blue curves, 18–26 slpm) of the experiments (figure 8). Acoustic fluctuations are generated due to flow turbulent fluctuations (Pawar et al. 2017). The exclusion of noise in the model leads to the cause. The second limitation is that the simulations (panels b,c) exhibit higher harmonics of comparable amplitude to the fundamental peak. The behaviour is absent in the experiments (panel a). For clarity, figure 13 directly compares the PSD between the experiments and the current model at two flow rates. The flow rates $36$ and $50$ slpm are chosen as they represent V- and A-lock-in, respectively. Altering the damping $c_1,c_2$ or heat release $\beta$ coefficients does not produce a qualitative change in the simulated PSD (not shown).

Figure 13. Comparison of the non-dimensional PSD between the experimental results of Singh & Mariappan (2019) (black colour) and the current simulations ($\beta =0.2$, red colour). The PSD is obtained from the pressure time series at $x=0.36$ for the flow rates 36 and 50 slpm where $1:1$ V- and A-lock-in are observed, respectively.

The following assumption has led to the above contrasting observations between the experiment and the model. The heat release is assumed to occur instantaneously when the vortex reaches the spatial location $x_c$ at the temporal instant $t_m^c$ (2.5). It leads to using the Dirac-delta ($\delta$) functions. The $\delta$ function in time excites all the modes with equal amplitudes. In reality, as the vortex convects along the flame, the heat release rate fluctuation $\dot q'$ is a continuous function of space and time. It should be noted that accurate modelling of $\dot q'$ distribution in space and time, such as that observed in the experiments (figure 1b), destroys the analytical tractability (§§ 2.3–2.5), which removes substantial findings of the paper.

Our main aim of the paper is to show the essential effects and the outstanding questions that exist in the literature of vortex-acoustic lock-in: (i) amplification or suppression of thermoacoustic instability (§ 4.2) and their analytical interpretation (§ 2.4), (ii) existence of lock-in regions where the common frequency can be either associated with acoustic (A-lock-in) or vortex shedding (V-lock-in) modes. Experimental comparisons are made for qualitative purposes only. Therefore, the presented results are at the expense of the model's assumptions (figure 8) in favour of the analytical tractability rather than the accuracy of the match with the experiments. The quantitative predictions of the model can certainly be enhanced by relaxing the above assumptions and using Bayesian inference (Juniper & Yoko 2022; Nóvoa & Magri 2022) to assimilate data from experiments or computational fluid dynamics, which will be taken as our future work.

5. Conclusion

The paper investigates the mutual interaction between vortex shedding and acoustic waves to study vortex-acoustic lock-in in a bluff body combustor. The former is modelled by an existing, widely used low-order model (Matveev & Culick 2003), while the latter is governed by the wave equation. Coupling between them is established through heat release rate fluctuations from vortex shedding, whose shedding characteristics are in turn, altered by the acoustic field. A set of $k_m$th order discrete dynamical maps relating pressure and velocity fluctuations of $k_m$ shedding events are obtained. The heat release coefficient $\beta$ and frequency ratio $f$ are the chosen control parameters, while the other parameters are selected from our previous experiments (Singh & Mariappan 2019).

Considering vortex shedding and the acoustic field as two coupled oscillators, the concept of phase lock-in is used to identify the region of vortex-acoustic lock-in. Two orders of lock-in, $1:1$ and $2:1$, are observed in the explored parameter space. During lock-in, the common oscillation frequency can lie close to the natural acoustic or vortex shedding frequency. Accordingly, the phenomena are termed A- and V-lock-in, respectively. By comparing our simulations with Singh & Mariappan (2019), we make the following salient conclusions. (i) The A- and V-lock-in phenomena are associated with the lock-in and pre-lock-in regimes reported in the experiments, respectively. (ii) Instability and amplitude suppression occur in most of the A- and V-lock-in regions, respectively. (iii) Saddle-node or Neimark–Sacker bifurcations occur at the onset of lock-in. (iv) Finally, we show that the boundary of the 1:1 A-lock-in regime, which is the most favourable condition for instability, can be determined from the lock-in boundary corresponding to the forced response of the vortex shedding process. An analytical expression for the forced lock-in boundary is provided in our earlier papers (Britto & Mariappan 2019, 2021), which can serve as a tool for the design of quiet combustors. All the conclusions of the study rely on a qualitative comparison between simulations and experiments. The analytical elegance/simplicity of the model is preferred over a quantitative match with the experiments. Future work will aim to use recent data assimilation methods developed for thermoacoustics to have a better quantitative match.