4.1. Step 1: Determination of conditional probability vortex shedding time period

The evolution of the joint probability density function, ${\mathscr{f}}$ , for $\varGamma$ and $\xi$ in the elapsed time $t^{\prime}$ is given by the following Fokker–Planck equation (corresponding to (3.1)–(3.5)):

(4.1) \begin{align} \frac {\partial {\mathscr{f}}}{\partial t^{\prime}}+\left (\frac {(1+u)^2}{2}+\xi \right ) \frac {\partial {\mathscr{f}}} {\partial \varGamma }&=\frac {1}{\epsilon _n}\frac {\partial \left (\xi {\mathscr{f}}\right )} {\partial \xi }+\frac {1}{\epsilon _n} \frac {\partial ^2 \left (A_n^2 {\mathscr{f}}\right )} {\partial \xi ^2}, \end{align}

(4.2) \begin{align} {\mathscr{f}}(\varGamma ,\xi ,t^{\prime}=0)&=\delta (\varGamma )\frac {e^{-\left ({\xi }^2/2{A_{n0}^2}\right )}}{\sqrt {2\pi } A_{n0}}\quad \textrm {initial condition}, \end{align}

(4.3) \begin{align} {\mathscr{f}}(\varGamma =\varGamma _{\textit{sep}},\xi ,t^{\prime})&=0\quad \textrm {boundary condition}, \end{align}

(4.4) \begin{align} {\mathscr{f}}(\varGamma \rightarrow -\infty ,\xi ,t^{\prime})&=0\quad \textrm {boundary condition}, \end{align}

(4.5) \begin{align} {\mathscr{f}}(\varGamma ,\xi \rightarrow \pm \infty ,t^{\prime})&=0\quad \textrm {boundary condition}, \end{align}

(4.6) \begin{align} \varGamma _{\textit{sep}}=\frac {1+u(t^{\prime})}{2}&=\frac {1+A\sin (\omega t^{\prime}+\phi )}{2}. \\[9pt] \nonumber \end{align}

The joint probability ${\mathscr{f}}(\varGamma ,\xi ,t^{\prime})$ is read as the probability that the vortex circulation and noise takes a value between $\varGamma \,\&\, {\textrm{d}}\varGamma$ and $\xi \,\&\, \xi +{\textrm{d}}\xi$ , respectively at a time $t^{\prime}$ , given its initial distribution at time $t^{\prime}=0$ . The above-mentioned Fokker–Planck equation is first order in $t^{\prime},\,\varGamma$ and second order in $\xi$ , requiring one initial, one boundary and two boundary conditions in $t^{\prime},\,\varGamma$ and $\xi$ , respectively. Here, $\varGamma$ takes a deterministic value of zero at the end of an $(m-1)$ th vortex shedding event. Therefore, the dependence of ${\mathscr{f}}$ on $\varGamma$ at $t^{\prime}=0$ is a Dirac-delta function placed at the origin, $\delta (\varGamma )$ .

At this point, one can choose various functional dependencies of ${\mathscr{f}}$ on $\xi$ at $t^{\prime}=0$ . Since the present investigation requires the stationary statistics of the shedding events, we choose the stationary probability distribution of $\xi$ as the initial condition for ${\mathscr{f}}$ . The stationary p.d.f. is Gaussian: ${e^{-({\xi }/{2A_{n0}})^2}}/{\sqrt {2\pi } A_{n0}}$ (zero mean and standard deviation $A_{n0}$ ). A similar choice is made by Palleschi & De Rosa (Reference Palleschi and De Rosa1992).

For $t^{\prime}\gt 0$ , the circulation of the $m$ th vortex increases according to (3.1)–(3.2). It increases until $\varGamma =\varGamma _{\textit{sep}}$ , leading to the $m$ th shedding. Since $\varGamma$ is reset to zero immediately after the $m$ th shedding, the corresponding probability, ${\mathscr{f}}(\varGamma =\varGamma _{\textit{sep}},\xi ,t^{\prime})=0$ (4.3). Boundary conditions (4.4)–(4.5) emphasise that ${\mathscr{f}}$ vanishes at $\pm \infty$ in the $\varGamma ,\,\xi$ domain, as appropriate.

The boundary condition (4.3) is time-varying due to $\varGamma _{\textit{sep}}$ , which poses difficulty in the numerical solution procedure. Since the initial dependence of ${\mathscr{f}}$ is Gaussian (4.2), it is possible to obtain a series solution in terms of Hermite polynomials (similar to Frisch & Nowakowski Reference Frisch and Nowakowski1993). The above-mentioned two points prompt us to perform a change of independent variables as follows: $\hat t=t,\ \hat \varGamma =1-{\varGamma }/{\varGamma _{\textit{sep}}},\ \hat \xi = {\xi }/A_{n0}$ . The transformation renders (4.1)–(4.5) to become the following:

(4.7) \begin{align} \frac {\partial {\mathscr{f}}}{\partial \hat t}+\left (\frac {(1-\hat \varGamma )\,{\textrm{d}}u/{\textrm{d}}\hat t-(1+u)^2}{1+u}-\frac {2A_{n0}\hat \xi }{(1+u)} \right ) \frac {\partial {\mathscr{f}}} {\partial \hat \varGamma }&=\frac {1}{\epsilon _n}\frac {\partial \big(\hat \xi {\mathscr{f}}\big)} {\partial \hat \xi }+\frac {1}{\epsilon _n}\left (\frac {A_n}{A_{n0}}\right )^2\frac {\partial ^2 {\mathscr{f}} } {\partial \hat \xi ^2}, \end{align}

(4.8) \begin{align} {\mathscr{f}}(\hat \varGamma ,\hat \xi ,\hat t=0)=\delta (\hat \varGamma -1)\frac {e^{-\big (\hat \xi ^2/2\big )}}{|1+A\sin \phi |\sqrt {\pi /2} A_{n0}}&\quad \textrm {initial condition}, \end{align}

(4.9) \begin{align} {\mathscr{f}}(\hat \varGamma =0,\hat \xi ,\hat t)=0&\quad \textrm {boundary condition}, \end{align}

(4.10) \begin{align} {\mathscr{f}}(\hat \varGamma \rightarrow \infty ,\hat \xi ,\hat t)=0&\quad \textrm {boundary condition}, \end{align}

(4.11) \begin{align} {\mathscr{f}}(\hat \varGamma ,\hat \xi \rightarrow \pm \infty ,\hat t)=0&\quad \textrm {boundary condition}. \\[9pt] \nonumber \end{align}

We use Galerkin method to solve  (4.7). Orthonormal Hermite functions $(\varphi _p,\,\psi _p$ of $p$ th degree) in both $\hat \varGamma$ and $\hat \xi$ are used as basis functions. They satisfy the boundary conditions, (4.9)–(4.11). The solution is expanded in the following form:

(4.12) \begin{align} {\mathscr{f}}\big (\hat \varGamma , \hat \xi , \hat t \big )=\sum _{\substack {r=0,1,2,\ldots R\\p=1,3,5,\ldots P}}a_{rp}\left (t\right )\varphi _r\big (\hat \xi \big ) \psi _p\big (\hat \varGamma \big ), \quad t\geqslant 0,-\infty \lt \hat \xi \lt \infty ,\,0\leqslant \hat \varGamma \lt \infty, \end{align}

where $\varphi _p (\hat \xi )=H_p (\hat \xi )e^{-\hat \xi ^2/2}/\sqrt { (\sqrt \pi 2^pp!) )},\,\psi _p (\hat \varGamma )=H_p (\hat \varGamma )e^{-\hat \varGamma ^2/2}/\sqrt { (\sqrt \pi 2^{p-1}p!) )}$ , with $\int _{-\infty }^{\infty }\varphi _r(x) \varphi _s(x)\,{\textrm{d}}x=\delta ^k_{rs},\,\int _{0}^{\infty }\psi _p(x) \psi _q(x)\,{\textrm{d}}x=\delta ^k_{pq}$ . Here, $\delta ^k_{rs}$ is the Kronecker-delta function (superscript $k$ is to differentiate the symbol from the Dirac-delta function used earlier). Additionally, $H_r$ is the $r$ th degree physicist’s Hermite polynomials. The Galerkin expansion satisfies the boundary conditions (4.9)–(4.11). Here, $R,\,P$ indicate the total number of Galerkin modes. Applying Galerkin projection on (4.7)–(4.8), we obtain a first-order ordinary differential equation for the time evolution of $a_{sq}$ , along with its initial condition,

(4.13) \begin{align} \frac {{\textrm{d}} a_{sq}}{{\textrm{d}}\hat t} \!=\! \sum _{\substack {r,p\\}}a_{rp} \! \left [\! I^{\alpha }\frac {\left (1+u\right )^2-{\textrm{d}}u/{\textrm{d}}\hat t}{1+u}+I^{\beta }\frac {{\textrm{d}}u/{\textrm{d}}\hat t}{1+u}+I^{\gamma }\frac {2A_{n0}}{(1+u)}+\frac {I^{\delta }+\zeta (1+u)^2I^{\varUpsilon }/4}{\epsilon _n}\! \right ]\!, \\[-28pt] \nonumber \end{align}

(4.14) \begin{align} a_{sq}(\hat t=0)=\frac {\sqrt {2}H_q(1)e^{-1/2}\delta ^k_{s0}}{|1+A\sin \phi |A_{n0}\pi 2^{s+q-1}s!q!} ,\end{align}

where $I^{\alpha },\,I^{\beta },\ldots, I^{\varUpsilon }$ are the Galerkin integrals, which are functions of the Hermite functions of degrees $p,q,r,s$ . Their functional forms are given in Appendix A. They are evaluated a priori to the simulations.

Equation (4.13) is time-marched to obtain the joint probability ${\mathscr{f}}$ . The conditional probability ${\mathscr{p}}$ for the circulation is then evaluated as ${\mathscr{p}}(\varGamma ,t)=\int _{\xi =-\infty }^{\infty }{\mathscr{f}}\,{\textrm{d}}\xi$ . Determination of the conditional probability of vortex shedding time period $\rho (\tau /\phi )$ is posed as a first passage time problem, when $\varGamma$ reaches $\varGamma _{\textit{sep}}$ for the first time. An expression is as follows: $\rho (\tau /\phi )=- (\textrm{d}/{{\textrm{d}}t})(\int _{-\infty }^{\varGamma _{\textit{sep}}(t)}{\mathscr{p}}\, {\textrm{d}}\varGamma )_{t=\tau }$ (details from Mariappan Reference Mariappan2022). After applying the coordinate transformation and Galerkin ansatz, the following expression for $\rho$ is obtained in terms of the Galerkin coefficients:

(4.15) \begin{equation} \rho (\tau /\phi )=-\frac {A_{n0}}{2}\sum _{\substack {r,p}}\left [(1+u)\frac {{\textrm{d}}a_{rp}}{{\textrm{d}}t}+\frac {{\textrm{d}}u}{{\textrm{d}}t}a_{rp}\right ]I^{\nu }_rI^{\varsigma }_p .\end{equation}

The dependence on $\phi$ occurs through the initial condition (4.8). In the next step, we assemble the shedding events to obtain a spike train and its corresponding long-time statistics.

4.2. Step 2: Assembly of vortex shedding events and identification of stochastic lock-in

The conditional probability $(\rho )$ of the $m$ th vortex shedding time is determined solely by the $m-1$ th shedding phase $\phi$ . Therefore, given a probability distribution $(\chi^{\boldsymbol{0}})$ of the first vortex shedding phase, one can evolve the probability distribution of the subsequent shedding phases $( \chi^{\boldsymbol{1}},\,\chi ^{\boldsymbol{2}},\,\chi^{\boldsymbol{3}},\ldots )$ . The successive probability distributions are related by $\chi^{\boldsymbol{m+1}} = \boldsymbol{T}\chi ^{\boldsymbol{m}}$ , where $\boldsymbol{T}$ is the transition probability matrix. The matrix acts as a dynamical operator, which transitions the shedding phase probability distribution from one ( $m$ th) to the next ( $m+1$ th), shedding instance. An expression of the matrix requires the introduction of additional variables. Therefore, readers are directed to refer to § 2.2 of Mariappan (Reference Mariappan2022). In long time $t\rightarrow \infty$ , the shedding phase probability, $\chi ^{\boldsymbol{s}}=\chi ^{\boldsymbol{m\rightarrow \infty}}$ , reaches a stationary distribution, $\chi^{\boldsymbol{s}}$ . Here, $\chi^{\boldsymbol{s}}$ is the first eigenvector of $\boldsymbol{T}$ , having an eigenvalue of 1 (Plesser & Geisel Reference Plesser and Geisel1999). The corresponding probability of the stationary vortex shedding time period, $\rho _s(\tau )=\int _0^{2\pi }\rho (\tau /\phi )\,{\textrm{d}}\phi$ , is calculated. The second and subsequent eigenvalues of $\boldsymbol{T}$ decide the route taken by the system on its way to $\chi^{\boldsymbol{s}}$ . Their qualitative change from complex to real values determines the onset of stochastic lock-in (s-lock-in). We restrict ourselves to 1 : 1 s-lock-in in this paper as our previous investigation (Mariappan Reference Mariappan2022) indicated that it favours combustion instability. Therefore, the region of the phase of the second eigenvalue equal to zero indicates the s-lock-in (Tateno & Jimbo Reference Tateno and Jimbo2000).

Having discussed the solution procedure, the next step would be to identify the parameters and their numerical ranges, which are of interest to academic and industrial combustor configurations. There are five parameters in the model (3.1)–(3.2): (i) amplitude $A$ ; (ii) frequency $f$ of the harmonic excitation; (iii) amplitude $A_{n0}$ ; (iv) correlation time $\epsilon _n$ and (v) multiplicative noise factor $\zeta$ of the OU noise. Their ranges are determined in the following section.

7.1. Effects of additive OU noise

Figure 2 shows the p.d.f. ( $\rho _s$ ) of the shedding period $\tau$ for pure additive noise. Results for increasing noise amplitude $A_{n0}$ are depicted along the columns. We make three observations. The first one is broadening the $\rho _s$ p.d.f. with the increase in $A_{n0}$ for a given $\epsilon _n$ : an expected behaviour.

The second observation is the spreading of the p.d.f. (standard deviation) for increasing $\epsilon _n$ , for a given noise amplitude $A_{n0}$ (individual panels). It is nicely visible in the simulation for large noise amplitude (panel $c$ ). The reason can be explained with the time domain simulation, shown in figure 3. Time evolution of the noise $\xi$ (panels a–c) and $\varGamma$ (panels d–f) for three different $\epsilon _n$ at $A_{n0}=0.1$ are plotted. The sudden drop in $\varGamma$ marks the shedding instances. At the lowest $\epsilon _n=0.01$ (panels a,d), there is almost no relaxation time for $\xi$ . Therefore, $\xi$ switches between positive and negative values quickly in comparison to the natural vortex shedding period. Therefore, the peak in $\rho _s$ is close to 1.0. At an intermediate $\epsilon _n=1.5$ (panels b,e), the larger relaxation time allows $\xi$ to stay on the positive or negative value side for longer. It allows quicker or slower, respectively, shedding sequences, leading to the spread in $\rho _s$ . The sequence appears to have intermittent fluctuations in the shedding frequency. It occurs when the natural shedding frequency is close to $\epsilon _n$ : $\epsilon _n\approx 1$ . The $\xi$ relaxes even slower for the largest $\epsilon _n=10$ (panels c,f), leading to a rise in the spread of $\rho _s$ (figure 3).

The third observation (related to the second one) is that the spread is skewed towards a lower shedding period for increasing $\epsilon _n$ and $A_{n0}$ . The skewed spread makes the p.d.f. non-Gaussian. Such p.d.f.s indicate the presence of intermittent fluctuations (Reuther & Kähler Reference Reuther and Kähler2020). The tail regions on either side of the p.d.f. $\rho _s$ captures the occasional events with either small or large vortex shedding period. Such intermittent fluctuations in the shedding period are observed in the time domain simulation of the circulation $\varGamma$ (figure 3 e). Therefore, the information regarding the transient dynamics (such as intermittent fluctuations) indeed appears in the p.d.f.

The mode (maximum probability of the shedding time period) of $\rho _s$ (white dashed curve in figure 2) indicates the gradual drop. It is observed even at $\zeta =0$ , indicating its cause is $\epsilon _n$ . An explanation in this regard is provided for the coloured OU noise. A formal solution in the absence of multiplicative noise ( $\zeta =0$ ) for (3.2) is as follows:

(7.1) \begin{equation} \xi (\tau )=\xi _0e^{-\tau /\epsilon _n}+\frac {A_{n0}}{\sqrt {2\epsilon _n}}e^{-\tau /\epsilon _n}\int _0^{\tau } e^{s/\epsilon _n}\eta _s\, {\textrm{d}}s, \end{equation}

where $\eta _s$ is the Gaussian white noise in the independent variable $s$ . Here, $\xi _0$ is the initial condition sampled from a Gaussian distribution of zero mean and $A_{n0}$ (noise amplitude of $\xi$ ) standard deviation for each realisation. Inserting (7.1) in (3.1) and performing a formal integration leads to the time evolution of $\varGamma$ . The expression is valid between the shedding instances,

(7.2) \begin{align} \varGamma (\tau ) = &\int _0^\tau \frac {(1+u(\tau ^*))^2}{2}\,{\textrm{d}} \tau ^* +\xi _0 \epsilon _n\left (1-e^{-\tau /\epsilon _n}\right ) \nonumber\\&+\, \frac {A_{n0}}{\sqrt {\epsilon _n/2}}\int _0^\tau e^{-\tau ^*/\epsilon _n}\left (\int _0^{\tau ^*}e^{s/\epsilon _n}\eta _s \,{\textrm{d}}s \right )\,{\textrm{d}}\tau ^* .\end{align}

Using (7.2), an analytical expression for the ensemble average (angular brackets) of $\varGamma ^2$ can be written as

(7.3) \begin{equation} \big \langle \varGamma ^2(\tau )\big \rangle =\left [\int _0^\tau \frac {(1+u(\tau ^*))^2}{2}\,{\textrm{d}} \tau ^*\right ]^2+ 2A_{n0}^2\epsilon _n\left (\epsilon _n e^{-\tau /\epsilon _n}-\epsilon _n+\tau \right)\!,\quad \tau \gt 0. \end{equation}

The following inputs are used in the derivation of the above-mentioned expression: $\langle \xi _0\rangle =0,\,\langle \eta _s\eta _r\rangle =\delta (s-r)$ . Deterministic and stochastic ( $S_c$ ) contributions form the first and second terms in the $\langle \varGamma ^2\rangle$ expression. Additive noise appears as an additive correction $(S_c)$ to $\langle \varGamma ^2\rangle$ . A contour plot of $S_c$ with the shedding period $\tau$ for the range of $\epsilon _n$ is shown in figure 4 for a given noise amplitude $A_{n0}$ . Here, $S_c$ increases with $\epsilon _n$ . Therefore, $S_c$ quickens the circulation $\varGamma$ to reach the critical value $\varGamma _{\textit{sep}}$ . It results in a reduction in the maximum probability of the shedding time period (in an ensemble average sense). Therefore, the mode of $\rho _s$ shows a decreasing trend with $\epsilon _n$ . The trend becomes stronger, as $S_c \propto A_{n0}^2$ . The decrease in the mode of $\rho _s$ plays a crucial factor in the shape of the lock-in boundary, when excited with the harmonic signal (§ 7.3).

Figure 4.

Contours of the stochastic term $S_c$ in the $\tau {-}\epsilon _n$ plane. The parameters are $A_{n0}=0.4,\,\zeta =0$ ).

Figure 5.

Variation of statistical measures of $\rho _s$ : mode ( $\tau _{mode}$ ) and standard deviation ( $\tau _{\textit{sd}}$ ) for various $\epsilon _n$ and $\zeta$ . Three noise amplitudes $A_{n0}=0.02,\,0.2$ and 0.4 are considered.

Note that the mode $\tau _{mode}$ in figures 2 and 5 is above 1 (natural shedding frequency) by ${\sim} 6\,\%$ for low $A_{n0}$ and $\epsilon _n$ . We expect it to be close to 1. We believe that the discrepancy is due to the 5 % variation cutoff in $\rho _s$ as the convergence criteria to decide the number of Galerkin modes $R,\,P$ (Appendix B).

7.2. Effects of multiplicative OU noise

The influence of the multiplicative noise in the shedding period statistics is presented in figure 5. Three representative noise amplitudes in the small ( $A_{n0}$ = 0.02, panels a,d), medium ( $A_{n0}$ = 0.2, panels b,e), and large ( $A_{n0}$ = 0.4, panels c,f) values are chosen for illustration. Mode ( $\tau _{mode}$ ) and standard deviation of the probability density $\rho _s$ are plotted in panels (c) and (f), respectively. The varied parameters are the correlation time ( $\epsilon _n$ ) and the multiplicative noise factor ( $\zeta$ ).

Generally, the mode, $\tau _{mode}$ (panels $a{-}c$ ) hovers around the natural vortex shedding frequency (1.0). A qualitative similarity in the features of $\tau _{mode}$ and $\tau _{\textit{sd}}$ distribution is observed for various $A_{n0}$ . For instance: (i) the presence of a hump in $\tau _{mode}$ (panels a–c) and (ii) the increase of low $\tau _{\textit{sd}}$ regions with $\epsilon _n$ and $\zeta$ . It indicates the underlying role of noise is the same for all $A_{n0}$ . Nevertheless, we notice quantitative differences. (i) In general, the mode $\tau _{mode}$ decreases with $A_{n0}$ , indicating the ensemble shedding period is quickened. A plausible reason is the factor $A_{n0}^2$ hastening $\langle \varGamma ^2\rangle$  (7.3) to reach the shedding circulation. (ii) Similarly, the standard deviation $\tau _{\textit{sd}}$ increases with higher noise levels ( $A_{n0}$ values).

In the low- $\zeta$ case, the noise amplitude $A_n$ is almost constant ( $\approx A_{n0}$ ), leaving the effect of $\epsilon _n$ dominant. As briefed earlier, the monotonic decrease in the $\tau _{mode}$ results from the stochastic contribution $S_c$ . As one increases $\zeta$ , the noise intensity $A_n$ fluctuates according to the circulation $\varGamma$ . Since the dependence of $A_n^2$ with $\varGamma$ is quadratic (3.2), there is always an increase in the noise amplitude. As per the expression for $S_c$ , one expects a reduction in the shedding period with $\zeta$ for a given $\epsilon _n$ . Although we observe the above-mentioned trend for low noise level $A_{n0}=0.02$ (panel $a$ ) at small $\epsilon _n$ , we observe a hump (non-monotonic behaviour, marked in panels a–c). We hypothesise that the hump is created due to the multiplicative nature of noise.

A rigorous mathematical explanation is not possible for the following reasons. (a) Since $A_n$ is a function of the circulation $\varGamma$ , (3.1)–(3.2) become coupled. Therefore, a formal solution for the noise $\xi$ similar to (7.2) for pure additive noise is not possible. (b) Equation (3.2) becomes nonlinear in $\varGamma$ through the term $A_n$ . The above-mentioned reasons prevent us from obtaining an equation similar to (7.3), where the effect of noise on the mode of $\rho _s$ is discernible through $S_c$ . For multiplicative noise, a plausible heuristic argument is as follows. As mentioned in reason (a), both $\xi$ and $\varGamma$ evolve simultaneously, leading to a feedback loop (similar to heat release rate and acoustics in thermoacoustic interaction). The loop either accelerates or decelerates the vortex-shedding process, depending on $\epsilon _n$ . The noise correlation time $\epsilon _n$ (can be compared with the time delay in thermoacoustics) acts as a memory effect of the past noise value $\xi$ . Similar to thermoacoustic interaction, the acceleration (low- $\tau _{mode}$ regions) or deceleration (high- $\tau _{mode}$ regions) depends on $\epsilon _n$ non-monotonically, leading to the observed hump in $\tau _{mode}$ . The location of the hump occurs around a particular $\epsilon _n$ and it strongly depends on $A_{n0}$ . High noise levels dominate the effect of $\epsilon _n$ . Therefore, the region (figure 5 a−c) of accelerated vortex shedding (low $\tau _{mode}$ ) expands, pushing the hump towards the low $\epsilon _n$ values.

Panels (d–f) show the standard deviation ( $\tau _{\textit{sd}}$ ) of $\rho _s$ . As in the previous paragraph, one expects a monotonic rise in $\tau _{\textit{sd}}$ with $\zeta$ due to the quadratic nature of the multiplicative noise amplitude. However, the multiplicative noise affects the vortex dynamics non-trivially, leading to non-monotonic variation of $\tau _{\textit{sd}}$ too. In particular, we observe regions of low $\tau _{\textit{sd}}$ (blue region). The combination of the low $\tau _{\textit{sd}}$ region, along with the region of $\tau _{mode}$ (figure 2) close to the harmonic excitation time period, aids in determining the stochastic-lock-in region, described in the following subsection.

7.3. Effects of simultaneous OU noise and harmonic excitation

A deterministic version of the vortex-shedding model ( $A_{n0}=0$ ) represents a self-sustained oscillator with frequency 1, also termed a period-1 solution in the parlance of a dynamical system. When excited with a harmonic external forcing, the shedding frequency and its dynamical pattern change (from the period-1 solution). At appropriate excitation amplitudes and frequencies, the period-1 solution reappears (dynamic bifurcation occurs). This time, the shedding frequency is the exciting frequency. The phase difference between the shedding and the excitation signals is also maintained. The phenomenon is termed a 1 : 1 phase lock-in. In general, the shedding and excitation frequency ratio can be rational $q/p$ ( $p,q$ are integers). It is defined as $p:q$ lock-in. Detailed theoretical investigation in relation to thermoacoustics is given by Britto & Mariappan (Reference Britto and Mariappan2019, Reference Britto and Mariappan2021).

In the stochastic counterpart, lock-in is identified by the phase of the eigenvalues of the transition probability matrix $\boldsymbol{T}$ . The eigenvalues are ordered in descending order of their absolute values. The largest (first) eigenvalue of $\boldsymbol{T}$ always equals 1, leading to a stationary probability distribution $\rho _s$ of the vortex-shedding time period. The 1 : 1 stochastic lock-in occurs when the phase of the second eigenvalue becomes zero. The phase of the higher eigenvalues determines higher-order s-lock-in. Our earlier result (Mariappan Reference Mariappan2022) shows that 1 : 1 s-lock-in provides the most favourable condition for combustion instability in the presence of Gaussian white noise. It is also found to be true when re-interpreted with our previous experiments (Britto & Mariappan Reference Britto and Mariappan2023). Furthermore, since the present paper is focused on studying the effects of colour and multiplicative noise, we restrict our attention to 1 : 1 s-lock-in, for limiting the range of harmonic excitation parameters, $f,\,A$ . We discuss the shape of the s-lock-in boundary in the harmonic excitation parameters plane, $A{-}f$ , followed by the noise parameters space, $A_{n0},\,\epsilon _n,\,\zeta$ .

7.3.1. Lock-in boundary in the harmonic excitation parameters $(A,\,f)$ plane

Figure 6 shows the contours of the second eigenvalue phase $(\theta _2)$ in the harmonic excitation parameter space for three different noise amplitudes $A_{n0}$ (rows) and correlation times $\epsilon _n$ (columns). Additive noise is considered, $\zeta =0$ . The white curve indicates the s-lock-in boundary, inside which $\theta _2=0$ . At the noise amplitude $A_{n0}=0.02$ , the lock-in boundary remains almost the same for all the correlation times $\epsilon _n$ . The feature is due to the low noise level.

Figure 6.

Contour plots detailing the variation of the phase $\theta _2$ of the second eigenvalue of the transition probability matrix $\boldsymbol{T}$ for various $A,\,f$ parameter combinations at $\zeta =0$ . Continuous white and dashed green curves represent the 1 : 1 s-lock-in boundary for $\zeta =0$ and 0.2, respectively.

At the intermediate $A_{n0}=0.2$ (panels b,e,h) and high $A_{n0}=0.4$ (panels c,f,i) noise levels, the boundary shape strongly depends on $\epsilon _n$ . Notably, the boundary is shifted to the higher excitation frequency (mostly $f\gt 1$ ) for $A_{n0}=0.2$ . A reasonable connection between the lock-in boundary shape and the harmonically unexcited $\tau _{mode},\,\tau _{\textit{sd}}$ (figure 5) contours can be made as follows. The maximum probability of the shedding time period, $\tau _{mode}$ (frequency) decreases (increases) for higher values of $\epsilon _n$ . In this region, the excitation frequency $f$ coincides more often with $1/\tau _{mode}$ . Furthermore, the $\tau _{\textit{sd}}$ has lower values in the higher $\epsilon _n$ region (refer to $\zeta =0$ in figure 5). A combination of the above-mentioned two factors leads to a shift and change in the lock-in boundary shape for intermediate noise amplitudes, $A_{n0}=0.2$ .

The lock-in region almost ceases to exist ( $\epsilon _n=0.77$ , panel $h$ in figure 6) and reappears ( $\epsilon _n=10$ , panel $i$ ) for the highest noise amplitude $A_{n0}=0.4$ . The cessation at $\epsilon _n=0.77$ can be attributed to the non-proximity between the values of $1/\tau _{mode}$ (figure 5 b) and $f$ . It is further fuelled by high values of $\tau _{\textit{sd}}$ at the regions where the above-mentioned proximity occurs. The reappearance of the lock-in region at $\epsilon _n=10$ (figure 6 i) is due to the low $\tau _{\textit{sd}}$ and proximity between $1/\tau _{mode}$ and $f$ .

Green dashed curves in figure 6 represent the lock-in boundary in the presence of multiplicative noise $\zeta =0.2$ . Since $\zeta$ has a substantial impact on $\tau _{mode}$ and $\tau _{\textit{sd}}$ (figure 5), we expect a significant difference in the lock-in boundary due to the multiplicative noise. Since the difference in $\tau _{mode}$ and $\tau _{\textit{sd}}$ between $\zeta =0$ and $0.2$ is significant for the lowest $\epsilon _n=0.01$ (figure 5), we observe a noteworthy difference in the lock-in boundaries in figure 6 (panels a,d,g). Furthermore, the same difference decreases with an increase in $\epsilon _n$ (figure 5), the lock-in boundary with multiplicative noise appears closer to the pure additive noise case for intermediate and high $\epsilon _n$ (last two columns of figure 6).

A practical observation from the above-mentioned exercise is that the s-lock-in location in the parameter space during harmonic excitation can be approximately identified from the (i) closeness of the mode of the shedding period to the harmonic excitation period and (ii) standard deviation of the shedding period, in the presence of only background turbulent fluctuations. The above-mentioned outcome is found to be common to both additive and multiplicative noise. Similar outcomes were used to detect an impending combustion instability (for example, Nair & Sujith Reference Nair and Sujith2014) and system identification (for example, Lee et al. Reference Lee, Kim, Gupta and Li2021) from the response of the system to turbulent fluctuations.

7.3.2. Lock-in boundary in the OU noise parameters $(A_{n0},\,\epsilon _n)$ plane

Figure 7 shows the occurrence of s-lock-in for various $\epsilon _n$ and $A_{n0}$ considering an additive noise $(\zeta =0)$ (panels a–d) and including a multiplicative noise (panels e–h). The exciting amplitude and frequency are the same: $A=0.25$ and $f=0.9$ . Lock-in is observed in the corresponding deterministic ( $A_{n0}=0$ ) case (refer to figure 6b of Britto & Mariappan Reference Britto and Mariappan2019).

Figure 7.

Contour plots detailing the variation of the phase $\theta _2$ (panels a,e) and the absolute value $|\kappa _2|$ (panels b, f) of the second eigenvalue of the transition probability matrix $\boldsymbol{T}$ . The next two columns represent the scaled mode $\tau _{mode} f$ of the p.d.f.s $\rho _s$ and the relative fluctuations $r_f$ . The harmonic excitation frequency and amplitude are $f=0.9$ and $A=0.25$ , respectively. The white dashed curve in all the panels represents the 1 : 1 s-lock-in boundary. The top and bottom row panels represent additive $\zeta =0$ and multiplicative noise $\zeta =0.2$ , respectively.

Considering the case of additive noise first (panels a–d), contours of $\theta _2$ are plotted. The dashed white curve indicates the s-lock-in boundary. At low noise amplitudes $A_{n0}\lt 0.09$ , s-lock-in continues to exist for all $\epsilon _n$ . It ceases to occur for large $\epsilon _n$ at higher noise amplitudes, due to the spread in the unforced shedding period distribution (figure 5). The cessation region widens with the rise in $A_{n0}$ . The white dashed curve (s-lock-in boundary) represents a stochastic bifurcation. It is defined as a transition point where the phase $\theta _2$ of the second eigenvalue of the transition probability matrix $\boldsymbol{T}$ becomes zero, as a parameter is varied. The transition of $\theta _2$ equalling zero is to identify 1 : 1 s-lock-in only, which is studied exclusively in this paper. Other orders of s-lock-in are identified when $\theta _2$ takes specific rational multiples of $\pi$ (Plesser & Tanaka Reference Plesser and Tanaka1997). Across a stochastic bifurcation, there is a qualitative change in the shedding period distribution $\rho _s$ (discussed in § 7.3.2 and shown in figure 8).

Figure 8.

P.d.f. of the shedding time period $\rho _s$ at a $\epsilon _n$ and $A_{n0}$ are shown in the first two columns: $\epsilon _n=0.77$ and $A_{n0}=0.15$ respectively. Vertical dashed lines locate the boundary of 1 : 1 s-lock-in. External excitation amplitude and frequency are the same as in figure 7.

Figure 7 $(b)$ shows the absolute value $(|\kappa _2|)$ contours of the second eigenvalue. Its value indicates the strength of the s-lock-in. Near the boundary (white dashed curve), $|\kappa _2|$ smoothly reduces between the s-lock-in and unlocked region. Panel ( $c$ ) provides a similar insight regarding the shedding period. The contour of the product of mode $(\tau _{mode})$ of $\rho _s$ and the excitation frequency $f$ equals 1 (green shade) in the s-lock-in region. Furthermore, the relative fluctuation $r_f$ (ratio of standard deviation $\tau _{\textit{sd}}$ and mean of $\rho _s$ ) is lower in the s-lock-in region and rises across the lock-in boundary. Large $r_f$ (panel d) leads to the unlocked state, although $\tau _{mode}f$ equals to 1 in a patch (marked in panel $c$ ) at high values of $A_{n0},\,\epsilon _n$ .

The presence of multiplicative noise (panels e–h) shrinks the s-lock-in region further (panel $a$ ). The phase $\theta _2$ takes values close to 0 (but not 0) at low $\epsilon _n$ . Hence, the system is close to s-lock-in: $\tau _{mode}f=1$ (panel g). Other than the above-mentioned differences, the qualitative features of $|\kappa _2|,\,\tau _{mode}f,\,r_f$ remain the same as those of the additive noise case.

The shrink in the s-lock-in region (at high $A_{n0}$ ) is observed at the location of a combination of low $\tau _{mode}$ with respect to the harmonic excitation time period $1/f=1.1$ and high standard deviation $\tau _{\textit{sd}}$ observed at $\zeta =0.2$ for zero harmonic excitation (figure 5). Therefore, we expect that multiplicative noise does not bring a new mechanism (dynamics), at least for s-lock-in. Intermittency or skewed $\rho _s$ distribution is observed for pure additive noise too, § 7.1.

Figure 8 shows the path the system takes as it crosses the s-lock-in boundary. It is illustrated through the contours of the p.d.f. $\rho _s$ for a given $\epsilon _n=0.77$ (panel a) and $A_{n0}=0.15$ (panel b). The other variables $A_{n0}$ and $\epsilon _n$ are varied respectively in the two panels. The product $\tau f$ is plotted for the p.d.f. random variable. For low noise amplitudes (panel a), where s-lock-in occurs, the $\rho _s$ distribution is narrow, peaking close to $\tau f=1$ (white horizontal line). The spread in $\rho _s$ increases as one approaches the s-lock-in boundary (vertical dashed line). Crossing the boundary leads to the peak shifting away from the excitation frequency. A similar situation occurs when $\epsilon _n$ (panel b) increases and for the case of multiplicative noise (not shown for brevity).

7.4. Effects of harmonic excitation frequency

The previous subsection dealt with considering a particular harmonic excitation frequency $f=0.9$ . The vortex shedding or acoustic frequency (essentially $f$ ) is often altered to keep the combustor in the unlocked state (Lieuwen Reference Lieuwen2012), thereby avoiding combustion instability. Therefore, we study the effects of noise parameters at various harmonic excitation frequencies $f$ .

Figure 9 shows the extent and strength of 1 : 1 s-lock-in in the $\zeta {-}\epsilon _n$ plane. Harmonic excitation frequency, $f$ , increases from the left to right panels. Panels (a)–(d) show the phase $\theta _2$ contours. The white curve in each panel indicates the boundary inside which $\theta _2=0$ and 1 : 1 s-lock-in occurs. In panels (a) and (b), where $f\lt 1$ , the size of the s-lock-in region increases as $f$ approaches 1. The genesis of the s-lock-in region occurs, where $\epsilon _n \sim 0.1$ (panel $a$ ). However, for $f\gt 1$ (panels $c,\,d$ ), it is pushed towards higher values of $\epsilon _n$ . The lock-in region also decreases as $f$ moves away from 1.

Figure 9.

Top and bottom row panels show the variation of the phase $\theta _2$ (normalised with $\pi$ ) and absolute $|\kappa _2|$ of the second eigenvalue of $\boldsymbol{T}$ in the $\zeta {-}\epsilon _n$ plane. The continuous white curve in each panel marks the boundary of the 1 : 1 s-lock-in region. Each panel represents the harmonic external excitation frequencies on either side of the natural vortex shedding frequency. Harmonic excitation and noise amplitudes are kept constant at $A=0.25$ and $A_{n0}=0.15$ , respectively. Four different harmonic excitation frequencies $f=0.8,\,0.9,\,1.1$ and 1.2 are considered.

The above-mentioned observation can be connected with the mode $(\tau _{mode})$ contours of the shedding frequency without harmonic forcing (figure 5 a), discussed previously in § 7.1. For $f\lt 1$ in figure 9(a, b), lock-in is observed in the region where $\tau _{mode}\gt 1$ in figure 5(a) (at $\epsilon _n \sim 0.1$ ). It indicates that lock-in occurs in the region where the maximum probability of the shedding period is close to the external harmonic excitation period. Using the same reasoning, one anticipates lock-in for $f\gt 1$ (figure 9 c, d) would occur for large and small (at least for $\zeta \gt 0.05$ ) values of $\epsilon _n$ . However, we observe lock-in only at large $\epsilon _n$ in figure 9(c, d). At small $\epsilon _n$ , the standard deviation $\tau _{\textit{sd}}$ (figure 5 b) is large leading to the disappearance of lock-in.

Figure 9(e–h) displays the $|\kappa _2|$ contours. The value of $|\kappa _2|$ denotes the strength of the s-lock-in. At $f=0.8$ (panel e), $|\kappa _2|$ takes relatively high values in the s-lock-in region. Its value further increases for $f=0.9$ and strengthens the s-lock-in region. However, in the $f\gt 1$ region, low values of $|\kappa _2|$ represent a weak s-lock-in. The analysis in the section concludes that the s-lock-in region is intense for $f\lt 1$ (excitation frequency is lower than the natural vortex shedding frequency) and it appears near $\epsilon _n \sim 0.1$ .

Figures 6–9 represent the s-lock-in region in various parameter space: $A,\,f,\,A_{n0},\,\epsilon _n$ and $\zeta$ . Lock-in is accompanied by combustion instability in many combustors (for example, Poinsot et al. Reference Poinsot, Trouve, Veynante, Candel and Esposito1987, Lieuwen Reference Lieuwen2012, Singh & Mariappan Reference Singh and Mariappan2021). Therefore, excluding the lock-in region from the above-mentioned figures serves as a safe operating map in the parameter space. Subsequently, it can also aid as a design guideline in the parameter selection of new combustors.

8.1. Contribution of stochastic resonance and s-lock-in

The results of stochastic resonance and distribution of $\varGamma _{\textit{sep}}$ in relation to combustion instability are shown in figure 10. Additive noise is considered in the first four panels. Panel (a) shows a patch of high normalised SNR in the $A_{n0}{-}\epsilon _n$ plane, indicating a strong stochastic resonance. Coloured noise promotes stochastic resonance, as reported by Nozaki, Collins & Yamamoto (Reference Nozaki, Collins and Yamamoto1999) for biological neurons. In the plot, stochastic resonance occurs at $A_{n0}=0.25,\,\epsilon _n=0.34$ . Many industrial gas turbine combustors operate in the above-mentioned turbulent noise parameter range (Hosseinalipour et al. Reference Hosseinalipour, Fattahi, Khalili, Tootoonchian and Karimi2020). Therefore, at first sight, stochastic resonance appears favourable for the occurrence of combustion instability in the practical combustors.

However, the 1 : 1 s-lock-in region (boundary shown by the continuous white curve) does not engulf the high stochastic resonance patch. A 1 : 1 s-lock-in is a period-1 dynamical solution for the shedding period $\tau$ , superimposed by stochastic fluctuations. The period-1 solution is centred around the excitation frequency (figure 7 c). An essential ingredient of a higher periodic solution begins to appear at the lock-in boundary (Balanov et al. Reference Balanov, Janson, Postnov and Sosnovtseva2008). Stochastic resonance does peak near the lock-in boundary.

8.2. Contribution of shedding circulation $\varGamma _{\textit{sep}}$

The probability $(\rho _{\varGamma _{\textit{sep}}})$ distribution of $\varGamma _{\textit{sep}}$ at three locations, low $A_{n0}=0.07$ , where s-lock-in occurs for all $\epsilon _n$ (panel b), medium $A_{n0}=0.25$ passing through the location of maximum stochastic resonance (panel c) and high $A_{n0}=0.38$ (panel d), where s-lock-occurs partly are plotted.

At low $A_{n0}$ , $\varGamma _{\textit{sep}}$ takes high probability (white regions) near $(1+A)/2$ for small $\epsilon _n\lessapprox 1$ , indicating a favourable condition for combustion instability. However, $\varGamma _{\textit{sep}}$ high probability switches to the lower side (near $(1-A)/2$ ) for $\epsilon _n\gtrapprox 1$ . Therefore, $\epsilon _n$ strongly influences the occurrence of combustion instability in the lock-in region.

In the medium $A_{n0}$ (panel c) case, a similar switch in $\varGamma _{\textit{sep}}$ probability occurs. The switch is closer to the location of maximum stochastic resonance. Therefore, stochastic resonance appears to play a favourable role for combustion instability. At the high $A_{n0}$ (panel d), high probability regions are both near $\varGamma _{\textit{sep}}=(1+A)/2$ and $(1-A)/2$ for $\epsilon \ll 1$ . Large $\epsilon _n$ values lead to high probability close to 0.5: the unperturbed shedding circulation.

Combining the above-mentioned observations on lock-in boundary and vortex circulation probability, we conclude that coloured noise promotes stochastic resonance significantly (approximately 10 times) compared with Gaussian white noise. The circulation ( $\varGamma _{v}$ ) at the time of vortex shedding is proportional to the unsteady heat release produced by the vortex (Matveev & Culick Reference Matveev and Culick2003). Since SNR is calculated based on $\varGamma _{v}$ , an approximate tenfold (figure 10 a) increase in the unsteady heat release is observed: a physical observable. Since the rise is significant compared with the Gaussian white noise (Mariappan Reference Mariappan2022), we expect coloured noise may play a crucial role in the occurrence of combustion instability through stochastic resonance. However, the above-mentioned argument is incomplete, as the occurrence of combustion instability also depends on the phase difference between the pressure and heat release rate oscillations (Rayleigh Reference Rayleigh1878). It requires two-way thermoacoustically coupled analysis (beyond the scope of the present paper). In any case, the tenfold rise in the unsteady heat release may still pose a threat to instability in practical combustors (Lieuwen & Yang Reference Lieuwen and Yang2005).

Adding the multiplicative noise (of highest $\zeta =0.2$ ), the normalised SNR dropped to ${\sim}1.5$ (panel $e$ ). Therefore, multiplicative noise does not promote stochastic resonance and, in turn, combustion instability. Other plots related to $\varGamma _{\textit{sep}}$ probability remain qualitatively similar to the additive noise case and therefore are omitted for brevity.

The highest multiplication factor $\zeta =0.2$ is chosen based on the simulations of Wilczek & Friedrich (Reference Wilczek and Friedrich2009) (refer to discussion in § 5). From figure 10, it is unclear whether the trend of normalised SNR monotonically decreases with $\zeta$ for intermediate values $(0\lt \zeta \lt 0.2)$ or completely suppresses the stochastic resonance at values $\zeta \gt 0.2$ . Performing simulations in the three-parameter space, $A_{n0},\,\epsilon _n$ and $\zeta$ , is computationally expensive. However, given the monotonic effects of $\zeta$ in the previous figures 5 and 9, we expect the same to hold true in stochastic resonance.

8.3. Possible passive control strategy

Based on the discussions made in §§ 7.3.1–8, it is possible to avoid s-lock-in by passively controlling the parameters $f,\,\epsilon _n$ and $\zeta$ . Other parameters are as follows. (i) Noise amplitude $A_{n0}$ and the correlation time $\epsilon _n$ are mostly determined by the incoming flow. They are known parameters. (ii) Acoustic oscillation amplitude $A$ is a result of thermoacoustic interaction. Therefore, the three parameters cannot be tuned for a passive control strategy. The author believes that there are two options available for passive control. (a) The easiest choice is to alter the frequency ratio $f$ , according to the s-lock-in maps shown in figures 6 and 9. It is achieved, either by tweaking the combustor (natural acoustic frequency, for example, Han et al. Reference Han, Han, Wang and Li2022) or burner (vortex shedding frequency, for example, Griffin Reference Griffin1985) geometries. (b) A second option is to alter the multiplicative noise factor $\zeta$ . Combustor liners can be built using acoustic metamaterials. They are inherently nonlinear in sound absorption (Gao et al. Reference Gao, Zhang, Deng, Guo, Cheng and Hou2022). Such nonlinear materials can be engineered to manipulate $\zeta$ . Although such materials are yet to be tested in the harsh environment of a combustor, future possibilities do exist.