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 Reference Matveev and Culick2003):

(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 (Reference Singh and Mariappan2019). 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 Reference Dotson, Koshigoe and Pace1997). 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 Reference Juniper and Sujith2018)

(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 Reference Cetegen and Sirignano1990). The interaction of vortices with flames enhances reaction rate (Cetegen & Basu Reference Cetegen and Basu2006). The enhancement is proportional to the circulation of the vortex (Peters & Williams Reference Peters and Williams1989). 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 (Reference Singh and Mariappan2019).

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 Reference Singh and Mariappan2019) 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 (Reference Singh and Mariappan2019) 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 Reference Matveev and Culick2003; Tuwankotta & Ihsan Reference Tuwankotta and Ihsan2019). 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 Reference Sterling and Zukoski1991; Matveev & Culick Reference Matveev and Culick2003).

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 Reference Sujith and Unni2020). 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 Reference Matveev and Culick2003)

(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 (Reference Singh and Mariappan2019). The vortex convection coefficient ($\alpha$) is assumed to be 0.6 (Dotson et al. Reference Dotson, Koshigoe and Pace1997). We use the values of damping coefficients $c_1,c_2$ from Matveev & Culick (Reference Matveev and Culick2003). In Singh & Mariappan (Reference Singh and Mariappan2019) 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 Reference Cundy and Rollett1989). 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 (Reference Matveev and Culick2003) 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.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 Reference Britto and Mariappan2021).

The instantaneous phase $(\psi _A(t))$ of $u_{sep}'(t)$ is obtained using the Hilbert transform (Balanov et al. Reference Balanov, Janson, Postnov and Sosnovtseva2008) $\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 (Reference Britto and Mariappan2021), 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 Reference Singh and Mariappan2019), 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 Reference Singh and Mariappan2019).

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 (Reference Britto and Mariappan2019).

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 Reference Britto and Mariappan2019) 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 Reference Li and Juniper2013c), 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 Reference Britto and Mariappan2019, Reference Britto and Mariappan2021). 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 (Reference Singh and Mariappan2019).

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. (Reference Dotson, Koshigoe and Pace1997), Matveev & Culick (Reference Matveev and Culick2003) and Singh & Mariappan (Reference Singh and Mariappan2019). 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 (Reference Emerson and Lieuwen2015), 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 (Reference Coutinho, Nogami and Toyama2009), 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 (Reference Britto and Mariappan2019, Reference Britto and Mariappan2021), 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 Reference Singh and Mariappan2019) (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 (Reference Singh and Mariappan2019)

Figure 11 compares the Rayleigh index and r.m.s. value of pressure fluctuations ($P_{rms}$) between the experiments (Singh & Mariappan Reference Singh and Mariappan2019) (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 Reference Shanbhogue, Husain and Lieuwen2009). Based on the previous experimental and numerical investigations (Rossiter Reference Rossiter1964; Zaman & Hussain Reference Zaman and Hussain1981; Zhou & Antonia Reference Zhou and Antonia1992; Kook & Mongeau Reference Kook and Mongeau2002), 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 Reference Huber and Polifke2009). 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. Reference Pawar, Seshadri, Unni and Sujith2017). 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 (Reference Singh and Mariappan2019) (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 Reference Juniper and Yoko2022; Nóvoa & Magri Reference Nóvoa and Magri2022) to assimilate data from experiments or computational fluid dynamics, which will be taken as our future work.