3.1. Impedance measurements and model derivation

Measurements of $\mathcal {Z}$ and $\mathcal {A}$ are made using the experimental set-ups respectively shown in figures 5(a) and 5(b). Acoustic forcing is applied with loudspeakers and the multi-microphone method (Schuermans et al. Reference Schuermans, Bellucci, Guethe, Meili, Flohr and Paschereit2004) is used to reconstruct the amplitude and phase of the forward and backward acoustic Riemann invariants $f$ and $g$, with which the reflection coefficient $\mathcal {R}$ and the specific impedance or admittance at a reference plane are deduced. The specific impedance of the cavity opening $\mathcal {Z}$ and the specific admittance of the deep cavity $\mathcal {A}$, which are measured with this method, are respectively presented in figures 6(a) and 6(b) for several bulk flow velocities $U$, and in figures 6(c) and 6(d) for several cavity lengths $L$. The white circles correspond to the specific impedance $\mathcal {Z}$ without flow. In that case, the specific resistance $\textrm {Re}(\mathcal {Z})$ is rather constant and it is equal to approximately 0.25 for the considered frequency range. This can be explained by considering the analogy between the travelling acoustic waves $f$ and $g$ in the side-branch cavity (see figure 5a) and the ones travelling toward and from an idealized compact area expansion in a duct. In the latter configuration, the amplitude reflection and transmission coefficients are respectively given by $\mathcal {R}=(\varepsilon -1)/(\varepsilon +1)$ and $\mathcal {T}=2\varepsilon /(\varepsilon +1)$, and the specific impedance is $\mathcal {Z}=\varepsilon$, where $\varepsilon$ is the area ratio at the sudden area expansion. In the present configuration, one can approximate the effective area expansion as the ratio between the cavity cross-section $WH$ and twice the wind channel cross-section $2H^{2}$, because acoustic energy originating from the cavity is transmitted in both the upstream and downstream directions of the channel. This leads to $\varepsilon \approx W/2H = 0.24$, which is very close to the measured specific acoustic resistance. In contrast with the real valued $\mathcal {Z}=\varepsilon$, the measured specific reactance is not zero. This is due to the fact that the inertia of the air at the cavity opening is not accounted for. The inertia of this attached air mass can be modelled as a fictive length of the cavity opening, based on the momentum balance $\rho \delta HW s \hat {u} = HW \hat {p}$, where $s=\textrm {i}\omega$ is the Laplace variable ($\omega =2{\rm \pi} f$ is the angular frequency), and the so-called end correction $\delta$ represents the fictive length of the cavity opening. From this balance, one has $\textrm {Im}(\mathcal {Z})=\omega \delta /c$, and, as explained in § 6.7 of Rienstra & Hirschberg (Reference Rienstra and Hirschberg2012), this end correction is of the order of the hydraulic radius of the cavity opening, which is approximately 25 mm in the present geometry. Based on the linear trend of the measured specific reactance of the opening without flow (white circles in figure 6b), one can deduce that $\delta \simeq 12$ mm.

Figure 5. (a) Experimental set-up for measuring the reflection coefficient $\mathcal {R}$, which is the ratio of the forward and backward acoustic Riemann invariants $f$ and $g$, and the specific impedance $\mathcal {Z}$ of the interface between the cavity and the channel along which the shear layer develops. (b) Experimental set-up for measuring the specific admittance of the deep cavity $\mathcal {A}$.

Figure 6. (a,b) Specific resistance $\textrm {Re}(\mathcal {Z})$ and reactance $\textrm {Im}(\mathcal {Z})$ of the cavity opening along which the shear layer develops for several bulk flow velocities $U$. The black lines correspond to the fits based on (3.1). The red dotted line in (a) highlights the presence of frequency ranges for which the resistance is negative, which means that reflected waves $g$ exhibit higher amplitude than incident waves $f$ (necessary condition for an aeroacoustic instability in the configuration of figure 1). (c,d) Modulus and phase of the specific admittance of the deep cavity for different lengths $L$. The black lines correspond to the fits based on (3.2).

In presence of the shear layer, the specific impedance of the cavity opening is significantly affected. One can see in figures 6(a) and 6(b) that the specific resistance and reactance both exhibit a frequency dependence oscillating around the values obtained without channel flow, which is typical of the impedance of side branches subject to grazing flow, e.g. Karlsson & Åbom (Reference Karlsson and Åbom2010). This is due to the response of the first longitudinal aerodynamic mode to the incident acoustic wave of complex amplitude $f$. Further insights into the nonlinear aerodynamic response of this shear layer when it is subject to incident acoustic waves are provided in the study of Boujo et al. (Reference Boujo, Bauerheim and Noiray2018). One can see in figure 6(a) that, for each velocity $U$, there exists a resistance minimum, which corresponds to the eigenfrequency $f_1$ of this aerodynamic eigenmode. The frequency of the minimum of $\textrm {Re}(\mathcal {Z})$ increases linearly with the flow velocity, scaling with the Strouhal number $St_1= f_1W/U\approx 0.4$, which has been already observed in several works on the topic, e.g. Dequand et al. (Reference Dequand, Hulshoff and Hirschberg2003). This corresponds to flow perturbations originating from the upstream corner that travel at approximately $0.4U$ across the opening during one acoustic oscillation cycle. The advection speed of the perturbations is roughly equal to the mean value of the velocity across the shear layer, which ranges from very low velocities in the cavity to the bulk velocity $U$ in the mean channel (see figure 4). It is important to note that the specific resistance minimum becomes negative for velocities exceeding $65\ \textrm {m}\ \textrm {s}^{-1}$. In fact, when the reflection coefficient satisfies $|\mathcal {R}|=|g/f|>1$, then $\textrm {Re}(\mathcal {Z})<0$, which is a necessary condition for self-sustained aeroacoustic oscillations in the configuration presented in figure 1. This occurs (i) if the time and spatially averaged projection of the unsteady component of the Lamb vector $(\boldsymbol {\omega }\times \boldsymbol {v})'$ onto the acoustic field is positive ($\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {v}$ is the vorticity), which corresponds to acoustic energy production according to Howe's energy corollary (Howe Reference Howe1980), and (ii) if this acoustic energy production exceeds the radiation losses in the wind channel. Readers interested by the space–time evolution of the unsteady component of the Lamb vector in a similar configuration (self-sustained aeroacoustic oscillations of a bottle whose neck is subject to a grazing flow) can refer to the recent work of Boujo et al. (Reference Boujo, Bourquard, Xiong and Noiray2020). In the present study, the measured specific impedance is fitted using the following second-order transfer function

(3.1)\begin{equation} \mathcal{Z}(s) = \frac{\hat{p}}{\rho c \hat{u}} = n\frac{s^{2} + 2m s + \omega_l^{2}}{s^{2} +2 d s + \omega_r^{2}}, \end{equation}

with $n$ the gain, $m$ and $d$ the damping coefficients and $\omega _l$ and $\omega _r$ the left and right angular frequencies associated with the acoustic feedback of the shear layer. Approximating the impedance with this second-order transfer function has two advantages: as explained below, (i) it can be nicely fitted to the experiments, and (ii) it leads to a system of coupled oscillators for describing the aeroacoustic instability. However, it has the drawback of having parameters that cannot be directly linked to the physics of the shear layer dynamics. Besides, one can note that Howe (Reference Howe1997) proposed a four-pole approximation for the conductivity of a rectangular aperture subject to low Strouhal flow, as an alternative to his physics-based model, and whose coefficients were calibrated to best reproduce the latter model.

Here, for each bulk flow velocity, optimization of the parameters of our second-order model is performed in order to find the best fit of the measured specific impedance over the frequency range of interest. The results are presented in figures 6(a) and 6(b), where the solid lines show the best fits. As a further step, these optimized model parameters can be linked to the system parameters ($U$, $W$ and $c$) in the form of non-dimensional numbers. First, very good estimates of the gains $n$ can be obtained using the relationship $n_0=n /M^{2}=11.8$, where $M=U/c$ is the Mach number. Second, the parameters $m$ and $d$ providing the best impedance predictions can be very well approximated using $d_1={d}W/U=0.273$ and $m_1=m W M^{3}/c =2.75 \times 10^{-4}$. Third, the left and right angular frequencies $\omega _l=2{\rm \pi} \,f_l$ and $\omega _r=2{\rm \pi} \,f_r$ can be deduced from the following Strouhal numbers $St_l=\,f_l W/U=0.375$ and $St_r=\,f_r W/U=0.461$. It is important to mention that the values obtained here for $n_0$, $d_1$, $m_1$, $St_l$ and $St_r$ cannot be generalized because they depend on the detail of the side-branch geometry and on the turbulent boundary layer thickness in the channel upstream of the cavity. It is, for instance, expected that they would differ for other shapes of the cavity opening, e.g. with round corners, or for a rough wall surface. Also, this model of the specific impedance of the cavity opening captures the acoustic feedback of the shear layer eigenmodes from 600 to 1300 Hz only. In that regard, one can refer to the work from Boujo et al. (Reference Boujo, Bauerheim and Noiray2018) for a computation of these eigenmodes in a side-branch configuration of the same opening width $W$ and the same duct height $H$ ($D$ in their paper), and subject to a turbulent flow of bulk velocity of $56\ \textrm {m}\ \textrm {s}^{-1}$, which corresponds to the yellow dots in figure 6(a). The minimum of the specific resistance in the present work is at approximately 750 Hz, which may be attributed to a 3-D shear layer mode that somehow corresponds to mode 1 (also 750 Hz) in § 3.3 of Boujo et al. (Reference Boujo, Bauerheim and Noiray2018). In fact, one cannot fully compare the latter numerical analysis based on the incompressible linearized Navier–Stokes equations (LNSE) with the present experiments. Indeed, while the incompressible assumption in the work of Boujo et al. (Reference Boujo, Bauerheim and Noiray2018) is valid because the shear layer is compact in this frequency range, there are three differences to keep in mind: (i) the LNSE analysis is two-dimensional (2-D) only; (ii) the 2-D mean flow for the LNSE was obtained from 3-D LES of a slice of 10 mm thickness with periodic conditions at the side boundaries, which differs from the finite spanwise extension of 62 mm with no-slip side boundaries of our experiment; (iii) the boundary layer thickness upstream of the opening in the LES may also differ from the one in our experiment.

Now, having found a suitable model for the specific impedance of the cavity opening, one focuses on the modelling of the cavity's specific admittance at $y=0$, whose measured modulus and phase are shown in figures 6(c) and 6(d). The specific admittance $\mathcal {A}$ is governed by the quarter wave resonances of the deep cavity, which occur at frequencies $f_n=(2n-1)c/4L$, i.e. $\omega _n L/c={\rm \pi} /2 \mod {\rm \pi}$. The specific admittance of a non-dissipative closed duct, below its cut-on frequency (pure one-dimensional acoustic propagation), is $\mathcal {A}=i\tan \omega L/c$. It was shown in § 2.1.2 of the paper from Bourquard & Noiray (Reference Bourquard and Noiray2019), that $\mathcal {A}$ asymptotically behaves, around $\omega _n$, as $-\gamma s/(s^{2}+\omega _n^{2})$, with $\gamma =2c/L$. This approximation provides an explicit formulation of quarter wave type resonances as a second-order harmonic oscillator with equivalent mass $\rho LWH/2$, i.e. half of the mass of air in the deep cavity, and with equivalent stiffness $[(2n-1)^{2}{\rm \pi} ^{2}/8] K_0$ where $K_0=\rho c^{2} (WH/L)$ is the stiffness associated with the bulk compression of an air column of length $L$ and cross-sectional area $WH$. Therefore, it is natural to consider the simple transfer function

(3.2)\begin{equation} \mathcal{A}(s) = \frac{\rho c \hat{u}}{\hat{p}} = -\frac{\gamma s}{s^{2} + 2\alpha s + \omega_a^{2}}, \end{equation}

where $\alpha$ is the acoustic damping in the cavity and $\omega _a$ is the angular frequency of the three-quarter wave resonance of the deep cavity. Using the measured specific admittance (coloured dots in figures 6c and 6d) and setting $\gamma =2c/L$ and $\omega _a= 3 {\rm \pi}c /2 L$, it is found that a damping $\alpha \simeq 40\ \textrm {rad}\ \textrm {s}^{-1}$ provides an excellent match between the measured $\mathcal {A}$ and the above transfer function model (solid lines in the figure) for the range of deep cavity lengths considered in this work.

Combining (3.2) and (3.1), and expressing these transfer functions in the time domain, one obtains the following system of differential equations for the acoustic pressure $p$ and the acoustic velocity $u$ at the cavity opening:

(3.3)\begin{equation} \left.\begin{gathered} \ddot{u} +2\alpha \dot{u} + \omega_a^{2} u = -\gamma \dot{p}\\ \ddot{p} + 2d \dot{p} + \omega_r^{2} p = n (\ddot{u} +2m \dot{u} + \omega_l^{2} u). \end{gathered}\right\}\end{equation}

Using the first equation to express the second time derivative of the acoustic velocity $\ddot {u}$, the system can be rewritten as

(3.4)\begin{equation} \left.\begin{gathered} \ddot{u} + 2\alpha \,\dot{u} + \omega_a^{2}u = -\gamma \dot{p},\\ \ddot{p} + 2\beta \,\dot{p} + \omega_r^{2} p = \mu \dot{u} + \sigma u, \end{gathered} \right\} \end{equation}

with $\beta =(2d+n\gamma )/2$, $\mu =2n(m-\alpha )$ and $\sigma =n(\omega _l^{2}-\omega _a^{2})$. This system of oscillators with resistive and reactive coupling depends on a set of a parameters, whose values are directly linked, as described above, to the physical parameters $U$, $L$ and $W$. It will now be used to predict the aeroacoustic stability of the deep cavity subject to grazing turbulent flow.

3.2. Linear stability analysis

Using $\boldsymbol {x}=(u, p)^\textrm {T}$, the system is expressed in the following matrix form:

(3.5)\begin{equation} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\ddot{\boldsymbol{x}} + \begin{bmatrix} 2\alpha & \gamma \\ -\mu & 2\beta \end{bmatrix} \dot{\boldsymbol{x}} +\begin{bmatrix} \omega_a^{2} & 0 \\ -\sigma & \omega_r^{2} \end{bmatrix} \boldsymbol{x} = 0. \end{equation}

The linear stability depends on the sign of the real part of the system's eigenvalues $\lambda$, which are the roots of the characteristic polynomial

(3.6)\begin{equation} \lambda^{4} + \lambda^{3} 2(\alpha+\beta) + \lambda^{2} (\omega_a^{2}+\omega_r^{2} + 4\alpha \beta + \gamma \mu) + \lambda (2\alpha\omega_r^{2} + 2\beta\omega_a^{2} +\gamma \sigma)+ \omega_r^{2}\omega_a^{2}=0. \end{equation}

Based on the previously identified scaling laws for the parameters of the linear model, the polynomial roots are computed for a range of bulk flow velocity $U$ and deep cavity length $L$. These roots are two pairs of complex conjugate eigenvalues, of which only the ones with positive imaginary part, i.e. positive angular frequency, are presented in figure 7. Several comments are now made about this figure.

Figure 7. Eigenvalues $\lambda$ of the aeroacoustic system for a range of $U$ and $L$. These eigenvalues are the roots of (3.6). (a) Map of the real part (growth rate) of the largest eigenvalue $\lambda _m$ as a function of $L$ and $U$. The linear stability limit is drawn in red, and the line where $\omega _a=\omega _r$ in white. (b) Value of $\textrm {Re}(\lambda )$ for $L=250$ mm as a function of $U$ (black lines). The stability limit is drawn as a dashed red line, and the coloured circles correspond to the eigenvalues with the largest real part, at the values of $U$ considered in figure 6. (c) Map of the imaginary part (frequency) of the most unstable eigenvalue $\lambda _m$. (d) Value of $\textrm {Im}(\lambda )$ for $L=250$ mm as a function of $U$ (black lines), superimposed with the acoustic mapping of figure 2. The black dashed lines correspond to the frequencies of each oscillator of the coupled system: $\omega _a$ (horizontal line) and $\omega _r$ (linearly increasing with $U$). The coloured circles correspond to the most unstable eigenvalues $\lambda _m$ at the values of $U$ considered in figure 6. (e) Prediction of the most unstable eigenvalue for the velocities considered in figure 6. (f) Corresponding experimental spectra for $L=250$ mm and for these velocities.

Firstly, the predicted linear stability of the aeroacoustic system for varying $L$ and $U$ is presented in figures 7(a) and 7(c). The former and the latter respectively show the real and imaginary parts (linear growth rate and oscillation frequency) of the most unstable of the two eigenvalues, which is denoted $\lambda _m$. The red line indicates the stability limit $\textrm {Re}(\lambda _m)=0$. It shows that, for $L<270$ mm and $U>65\ \textrm {m}\ \textrm {s}^{-1}$, the system is linearly unstable around the white line, indicating coincidence of the resonance frequency of the deep cavity $\omega _a$, which governs the oscillator equation for the acoustic velocity, and the resonance frequency of the shear layer $\omega _r$, which governs the oscillator equation for the acoustic pressure.

Secondly, a subset (for $L=250$ mm) of the predicted system eigenvalues $\lambda$ is presented in figures 7(b) and 7(d). In the former, one of the eigenvalues has a significantly smaller real part for the considered range of $U$, which implies that it is much more stable than the other. The eigenvalue with the largest real part crosses the complex plane imaginary axis when $U\simeq 65\ \textrm {m}\ \textrm {s}^{-1}$, i.e. the system becomes linearly unstable beyond this bulk flow velocity. In figure 7(d), the excellent match between the peak frequency of the overlayed power spectral density and the frequency of the least-stable eigenvalue indicates that the present coupled oscillators model performs very well. Moreover, in contrast with the phenomenon of frequency lock-in in flow-induced vibration problems, which are also modelled as coupled oscillators (De Langre Reference De Langre2006; Mohany et al. Reference Mohany, Arthurs, Bolduc, Hassan and Ziada2014; Shoshani Reference Shoshani2018; Dolci & Carmo Reference Dolci and Carmo2019), the frequencies of the eigenvalues $\textrm {Im}(\lambda )$ do not merge in the range of $L$ and $U$ for which the present aeroacoustic system is linearly unstable. Indeed, the black lines corresponding to the root loci are repelled from the intersection point of the natural frequencies of the two oscillators (dashed black lines). The fundamental topological difference of the coupled-oscillator root loci between flow-induced vibration problems and aeroacoustic instabilities of deep cavities subject to grazing flow originates from the differences in stability and coupling nature of the coupled-oscillator model. In the case of the flow-induced vibration problems, the hydrodynamic oscillator is linearly unstable and it is typically modelled as a van der Pol oscillator (De Langre Reference De Langre2006), the mechanical oscillator is linearly stable, and the coupling between these two oscillators is usually purely reactive. In the present case of aeroacoustic instabilities of deep cavities, which we model in this section with the system (3.4), both oscillators are linearly stable and the system can become linearly unstable because of the presence of both resistive and reactive coupling terms.

Finally, in figure 7(e), the most unstable eigenvalue $\lambda _m$ is plotted in the complex plane for the velocities considered in figure 6) and for $L=250$ mm. These eigenvalues can be compared with the corresponding experimental power spectral densities for the same velocities $U$ that are presented in figure 7(f). The good agreement in terms of frequency and the sharpening of the peak for $U>65\ \textrm {m}\ \textrm {s}^{-1}$ again contribute to the linear model validation.

4.1. Describing function analysis

Impedance measurements of the cavity opening are performed with the set-up shown in figure 5(b) for a range of forcing amplitudes. This forcing amplitude is deduced from the multi-microphone method and corresponds to the amplitude at a pressure antinode. The results of these measurements for $U=74\ \textrm {m}\ \textrm {s}^{-1}$ are presented in figures 8(a) and 8(b), respectively showing $\textrm {Re}(\mathcal {Z})$ and $\textrm {Im}(\mathcal {Z})$. For increasing amplitude, there is a monotonic decrease of the deviation of the specific impedance from the one without flow, which shows that the shear layer responds with less strength to the acoustic forcing. The underlying mechanism has been presented by Boujo et al. (Reference Boujo, Bauerheim and Noiray2018) for $U=56\ \textrm {m}\ \textrm {s}^{-1}$ and is in line with previous work on the subject: as the forcing amplitude grows, coherent Reynolds stresses thicken the mean shear layer, which reduces the potential for perturbation amplification at the forcing frequency. As a consequence, the range across which the real part of the impedance is negative reduces progressively as the amplitude increases. It is seen that, beyond a forcing amplitude of 200 Pa, $\textrm {Re}(\mathcal {Z})$ is positive for the whole frequency range. Noticeably, in the range 800 to about 1150 Hz, for which the shear layer produces acoustic energy, this contribution is less energetic than the radiation losses to the wind channel, and consequently, $\textrm {Re}(\mathcal {Z})>0$ and the modulus of the reflection coefficient $|\mathcal {R}|$ is lower than 1.

Figure 8. (a) Real and (b) imaginary part of the specific impedance of the cavity opening for different acoustic forcing amplitudes and for $U=74\ \textrm {m}\ \textrm {s}^{-1}$. The white circles correspond to the specific impedance without flow.

As in § 3.1, for each forcing amplitude, the parameters of the transfer function given in (3.1) are optimized to best fit the measured specific impedance $\mathcal {Z}$ over the frequency range presented in figures 8(a) and 8(b). The solid lines in this figure show these best fits. This optimization has also been performed with measurements of the specific impedance for the same set of bulk flow velocities as in figure 6. For each forcing amplitude, the optimized model parameters ($d$, $m$, $n$, $\omega _l$ and $\omega _r$) were linked to the system parameters $U$, $W$ and $c$ with the scaling laws presented in § 3.1. The non-dimensional numbers, with which these optimized model parameters can be deduced, are presented as coloured dots in figure 9 for several acoustic forcing amplitudes, with the same colour code for the bulk flow velocity $U$ as in figure 6.

Figure 9. Non-dimensional numbers, which link the system parameters $U$, $W$ and $c$ to the optimized model parameters $\omega _l$, $\omega _r$, $d$, $m$ and $n$ as functions of the acoustic forcing amplitude. The symbols are coloured with the same colour code for the bulk flow velocity $U$ as in figure 6. The dashed lines show the scaling laws used in the model of § 4.2. (a,b) Respectively show $St_l = \omega _l W / 2 {\rm \pi}U$ and $St_r = \omega _r W / 2 {\rm \pi}U$ that are assumed independent of the forcing amplitude in § 4.2. (c) Shows $dW/U$ and the scaling law used in § 4.2: $dW/U= d_1 + d_2 |p|$. (d) Shows $mWM^{3}/c$ and the scaling law used in § 4.2: $mWM^{3}/c= m_1 + m_3 p^{2}$. (e) Shows $n_0=n/M^{2}$ that is assumed constant in § 4.2.

Based on these data, it is possible to perform a describing function analysis for predicting the amplitude of the aeroacoustic limit cycle, as was done by Noiray et al. (Reference Noiray, Durox, Schuller and Candel2008) in the case of a thermoacoustic instability. In the present situation, it is possible to predict the limit cycle amplitude and frequency as functions of the cavity length $L$, as in the work of Noiray et al. (Reference Noiray, Durox, Schuller and Candel2008), and thanks to the model and the identified scaling laws for its parameters, one can also construct bifurcation diagrams as functions of the bulk flow velocity in the channel $U$. This is now performed for $L=250$ mm and $U=74\ \textrm {m}\ \textrm {s}^{-1}$: the eigenvalues of the coupled system (3.4) are computed using the optimized model parameters. For each of the forcing amplitudes, the real and imaginary parts of the most unstable eigenvalue $\lambda _m$ are displayed in figures 10(a) and 10(b). When the amplitude increases, the linear growth rate of the system monotonically decreases while the oscillation frequency does not significantly vary. The linear growth rate is positive at very low amplitude (approximately $20\ \textrm {rad}\ \textrm {s}^{-1}$) and, according to the describing function framework, the aeroacoustic system becomes marginally stable when the real part of $\lambda _m$ vanishes, which in the present case occurs at approximately 150 Pa. This predicted limit cycle amplitude is now compared to the actual self-sustained oscillation amplitude, which is measured for $L=250$ mm and $U=74\ \textrm {m}\ \textrm {s}^{-1}$ using the experimental set-up presented in figure 1. This comparison can be done with figures 10(c) and 10(d), which respectively show the band-pass filtered acoustic signal measured at a pressure antinode (only the positive acoustic pressure fluctuations are shown on this figure), and the probability density function (p.d.f.) of the signal's envelope. One can draw the following conclusions from figure 10: First, this describing function analysis provides a realistic estimate of the limit cycle amplitude (approximately 150 Pa), but it cannot be used for a quantitative prediction of the most probable amplitude (approximately 250 Pa). Second, the deterministic and frequency-domain describing function framework cannot capture the significant random fluctuations of the aeroacoustic limit cycle amplitude that are induced by the forcing from the intense turbulence in the channel. The next section aims at filling this gap by considering a nonlinear time-domain analysis of the problem and the turbulent stochastic forcing will be accounted for in § 5.

Figure 10. (a) Real and (b) imaginary parts of the most unstable eigenvalue from system (3.5) as functions of the acoustic amplitude for $L=250$ mm and $U=74\ \textrm {m}\ \textrm {s}^{-1}$. The eigenvalues are the roots of (3.6) in which the optimized model parameters $\omega _l$, $\omega _r$, $d$, $m$ and $n$ (dark blue circles in figure 9) were set. (c) Acoustic pressure filtered around the limit cycle frequency as done in figure 2(d) (grey line, only $p>0$ is shown), and corresponding envelope $B$ (thick black line), for $L=250$ mm and $U=74 \ \textrm {m}\ \textrm {s}^{-1}$. (d) Probability density function of the acoustic signal envelope $\mathcal {P}(B)$.

4.2. Time-domain model of coupled oscillators

The starting point of this section is the time-domain model (3.4) of the two coupled linear oscillators developed in § 3. One now aims at incorporating into this model relevant nonlinear terms on the basis of the specific impedance measurements presented in the previous section. One can deduce from the model parameters, which were optimized to reproduce the specific impedance for different forcing amplitudes and flow velocities, scaling laws that not only depend on the system parameters, but also on the acoustic amplitude. These scaling laws are presented as dashed lines in figure 9. Considering that $St_l$, $St_r$ and $n_0$ do not vary significantly with the acoustic forcing, and that they are respectively equal to 0.375, 0.461 and 11.8, the same simple laws as in § 3.1 are used in the next sections: $\omega _l=2{\rm \pi} St_l U/W$, $\omega _r=2{\rm \pi} St_r U/W$ and $n=n_0M^{2}$.

On the other hand, one can see in figures 9(c) and 9(d) that the values of $dW/U$ and $m W M^{3} /c$ significantly depend on the acoustic amplitude, and they can be respectively approximated by a linear and quadratic regressions (dashed lines). Therefore, the model parameters $m$ and $d$ are respectively deduced from $dW/U={d}_1+{d}_2|p|$ with $d_1=0.273$ and $d_2=8.03 \times 10^{-4}\ \textrm {Pa}^{-1}$, and from $m W M^{3} /c = m_1+m_3 p^{2}$ with $m_1=2.75 \times 10^{-4}$ and $m_3=3.90 \times 10^{-10}\ \textrm {Pa}^{-2}$. These scaling laws are now incorporated into the time-domain model which becomes:

(4.1)\begin{equation} \left.\begin{gathered} \ddot{u} + 2\alpha \dot{u} + \omega_a^{2}~u = -\gamma \dot{p},\\ \ddot{p} + 2(\beta_1 + \beta_2|p|) \dot{p} + \omega_r^{2} p = (\mu_1 + \mu_3 u^{2}) \dot{u} + \sigma u, \end{gathered} \right\}\end{equation}

with

(4.2)\begin{align} \left.\begin{gathered} \beta_1=\dfrac{U}{W}d_1+\dfrac{\gamma M^{2}}{2}n_0, \quad \beta_2=\dfrac{U}{W}d_2,\quad \mu_1=2 n_0 M^{2} \left(\dfrac{c}{WM^{3}}m_1-\alpha\right),\quad\gamma=\dfrac{2c}{L}\\ \mu_3=(\rho c)^{2}\dfrac{2c}{WM}n_0m_3, \quad \sigma=n_0 M^{2}\left( \left[2{\rm \pi}\dfrac{U}{W}St_l\right]^{2} - \omega_a^{2} \right), \quad \omega_r=2{\rm \pi}\dfrac{U}{W}St_r. \end{gathered}\right\} \end{align}

It can be noted that (i) the increase of the acoustic pressure amplitude $p$ leads to an increase of the effective damping coefficient ($\beta _1 + \beta _2|p|$) of the second oscillator, which is associated with the lower receptivity of the thickened shear layer at higher amplitude, and (ii) the increase of the acoustic velocity amplitude $u$ yields an increase of the effective resistive coupling coefficient $(\mu _1 + \mu _3 u^{2})$.

This nonlinear deterministic model of coupled oscillators for describing the aeroacoustic dynamics of the deep cavity depends on a set of constants ($n_0$, $d_1$, $d_2$, $m_1$, $m_3$, $St_l$, $St_r$ and $\alpha$) that were quantified from measurements, and on the system parameters $\rho$, $U$, $c$, $L$, $W$ and the Mach number $M=U/c$. One also recalls that the angular frequency of the three-quarter wave resonance is $\omega _a= 3 {\rm \pi}c /2 L$. Before augmenting this model in § 5 with stochastic forcing from turbulence in order to explain the random fluctuations of the limit cycle amplitude, the amplitude and phase equations of this deterministic model are derived and analysed in the next section.

4.3. Amplitude and phase equations

The averaging procedure of Krylov & Bogoliubov (Reference Krylov and Bogoliubov1936) is now applied to the system of coupled oscillators (4.1) in order to obtain first-order differential equations for the oscillation amplitude and phases. One assumes that oscillations occur at $\omega$ and that this angular frequency satisfies $\omega \simeq \omega _a \simeq \omega _r$ and $\omega _a + \omega _r \simeq 2 \omega$, where $\omega _a$ and $\omega _r$ correspond to the natural angular frequency of the two coupled oscillators. The following ansatze for the acoustic velocity and pressure are used:

(4.3)\begin{equation} \left.\begin{gathered} u=A \cos (\omega t + \varphi_A) = \tfrac{1}{2}( a \ \textrm{e}^{\textrm{i}\omega t} + a^{*}\, \textrm{e}^{-\textrm{i}\omega t}) \quad \text{with} \ a=A\ \textrm{e}^{\textrm{i}\varphi_A}, \\ p=B \cos (\omega t + \varphi_B) = \tfrac{1}{2}( b \ \textrm{e}^{\textrm{i}\omega t} + b^{*}\, \textrm{e}^{-\textrm{i}\omega t}) \quad \text{with}\ b=B\ \textrm{e}^{\textrm{i}\varphi_B}, \end{gathered} \right\} \end{equation}

where $A$ represents the velocity amplitude and $B$ the pressure amplitude, and $\varphi _A$ and $\varphi _B$ their phases. This averaging approach can only be applied if the oscillation amplitudes and phases vary slowly with respect to the acoustic period, which is satisfied in the present problem, and which corresponds to damping and coupling terms that are small compared to the inertial and stiffness terms of the two oscillator equations. Assuming that the first derivative of the acoustic velocity can be written as

(4.4)\begin{equation} \dot{u}=\dfrac{\textrm{i} \omega}{2}( a \ \textrm{e}^{\textrm{i}\omega t} - a^{*}\ \textrm{e}^{-\textrm{i}\omega t}), \end{equation}

which implies, as explained by Balanov et al. (Reference Balanov, Janson, Postnov and Sosnovtseva2009), that $\dot {a} \ \textrm {e}^{\textrm {i}\omega t} + \dot {a}^{*}\ \textrm {e}^{-\textrm {i}\omega t} =0$, one can express the second time derivative of the acoustic velocity as

(4.5)\begin{equation} \ddot{u}=\textrm{i}\omega \dot{a}\ \textrm{e}^{\textrm{i}\omega t} - \dfrac{\omega^{2}}{2}( a\ \textrm{e}^{\textrm{i}\omega t} + a^{*}\, \textrm{e}^{-\textrm{i}\omega t}). \end{equation}

Following the same procedure for $p$, substituting these expressions into the first oscillator equation of (4.1), multiplying by $\textrm {e}^{-\textrm {i}\omega t}/\textrm {i}\omega$, integrating over one cycle, dividing by $\textrm {e}^{\textrm {i}\varphi _A}$ and taking the real and imaginary parts of the equation yields

(4.6)\begin{equation} \left. \begin{gathered} \dot{A}= - \alpha A - \dfrac{\gamma}{2} B \cos(\varphi_A - \varphi_B), \\ \dot{\varphi}_A=\dfrac{\omega_a^{2} - \omega^{2}}{2\omega} + \dfrac{\gamma}{2} \dfrac{B}{A} \sin(\varphi_A - \varphi_B). \end{gathered} \right\} \end{equation}

Similar treatment of the second oscillator equation in (4.1) yields

(4.7)\begin{align} \left.\begin{gathered} \dot{B}= - \beta_1 B - \beta_2 \dfrac{4}{3{\rm \pi}} B^{2} + \left(\dfrac{\mu_1}{2} + \dfrac{\mu_3}{8} A^{2}\right) A \cos(\varphi_A - \varphi_B) + \dfrac{\sigma}{2\omega}A\sin(\varphi_A - \varphi_B),\\ \dot{\varphi_B}=\dfrac{\omega_r^{2} - \omega^{2}}{2\omega} + \left(\dfrac{\mu_1}{2} + \dfrac{\mu_3}{8} A^{2}\right) \dfrac{A}{B} \sin(\varphi_A - \varphi_B) - \dfrac{\sigma}{2\omega}\dfrac{A}{B}\cos(\varphi_A - \varphi_B). \end{gathered}\right\} \end{align}

Defining the phase difference $\phi =\varphi _A - \varphi _B$ yields the following system of three coupled equations for the slowly varying amplitudes and the phase of the coupled oscillators:

(4.8)\begin{equation} \left.\begin{gathered} \dot{A} = - \alpha A - \dfrac{\gamma}{2} B\cos\phi, \\ \dot{B} = - \beta_1 B - \beta_2 \dfrac{4}{3{\rm \pi}} B^{2} + \left(\dfrac{\mu_1}{2} + \dfrac{\mu_3}{8} A^{2}\right) A \cos\phi + \dfrac{\sigma}{2\omega}A\sin\phi, \\ \dot{ \phi} = \varDelta - \left( \dfrac{\gamma}{2} \dfrac{B}{A} - \left(\dfrac{\mu_1}{2} + \dfrac{\mu_3}{8} A^{2}\right) \dfrac{A}{B} \right) \sin\phi + \dfrac{\sigma}{2\omega}\dfrac{A}{B}\cos\phi, \end{gathered}\right\} \end{equation}

with $\varDelta = (\omega _a^{2} - \omega _r^{2})/2\omega \simeq \omega _a - \omega _r$ the detuning constant. The fixed points of this system of first-order differential equations are found by searching for amplitudes and phase difference that cancel the right-hand side expressions of (4.8) and that thus lead to $\dot {A}=\dot {B}=\dot {\phi }=0$. For a few combinations of $L$ and $U$, this search is initialized from the amplitudes $A$ and $B$ and phase difference $\phi$ of limit cycles that were computed by integrating the second-order system (4.1). The fixed points of the slow-flow system (4.8) are then obtained for the full range of $L$ and $U$ by continuation of the local minimum search from neighbouring solutions. The results are presented in figure 11, which shows the acoustic velocity amplitude $A$, the acoustic pressure amplitude $B$ that is afterward compared with experimental data and their phase difference $\phi$ at conditions featuring aeroacoustic limit cycles. These results show that, according to this deterministic model, the deep cavity subject to turbulent grazing flow undergoes supercritical Hopf bifurcations only. This is further illustrated in figures 11(a) and 11(b) for three values of the bulk velocity $U$. One can see that over the entire range of combinations of $L$ and $U$ leading to limit cycles, the phase difference between the acoustic pressure $p$ and the acoustic velocity $u$ is nearly constant and slightly below $-{\rm \pi} /2$, which is typical of standing mode oscillations in quarter wave resonators. It is important to stress that the maps predicting the stability border and the limit cycle amplitudes in figure 11 for a wide range of bulk velocity $U$ and cavity length $L$, were obtained by using the scaling laws presented in figure 9 as dashed lines. These scaling laws allow reliable prediction of the specific impedance $\mathcal {Z}$ of the cavity opening geometry investigated in this work, for broad ranges of bulk velocity $U$ and acoustic amplitude.

Figure 11. (b,d) Limit cycle amplitudes $A$ and $B$ and (f) phase difference $\phi$ as functions of $L$ and $U$. They are fixed points of the slow-flow equations (4.8) and they were obtained by searching the zeros of the right-hand side of these equations. The origin is the only fixed point in the white region delimited by the stability border presented in figure 7; (a), (c) and (e) cuts showing the limit cycle amplitude for $U=67$, 70 and $74\ \textrm {m}\ \textrm {s}^{-1}$ as functions of $L$.

In figure 12(a), the limit cycle amplitude $B$, which is predicted from the slow-flow system (4.8) for $U=74\ \textrm {m}\ \textrm {s}^{-1}$ and for several cavity lengths $L$, is compared to the p.d.f. of the band-pass filtered acoustic pressure envelope measured with the experimental set-up shown in figure 1. The overall agreement between model predictions and most probable amplitude from the experiments is very good, albeit the former displays (i) a smoother dependence on the cavity length $L$, (ii) with bifurcation points shifted by a relatively small offset of approximately 10 mm (approximately 5 % of the range of cavity length displayed in this figure). The former difference (i) is due to the fact that the actual specific impedance is not as smooth as the second-order transfer function adopted in this work to model it (see figure 12(b) showing $\textrm {Re}(\mathcal {Z})$ for $U=74\ \textrm {m}\ \textrm {s}^{-1}$). The latter small difference (ii) may be reduced by using the values of the parameters obtained from the best fits of the specific impedance at $U=74\ \textrm {m}\ \textrm {s}^{-1}$ (dark blue dots in figure 9), which yield, at low acoustic amplitude, the solid black line in figure 12(b), instead of using the more general, but slightly less accurate at a given $U$, scaling laws defining the model parameters for broad ranges of velocity and acoustic amplitude (dashed lines in figure 9).

Figure 12. (a) The p.d.f.s of the envelope of the band-pass filtered acoustic pressure for several cavity lengths $L$, compared to the predicted bifurcation diagram of the slow-flow system (4.8), which shows the limit cycle amplitude $B$ for $U=74 \ \textrm {m}\ \textrm {s}^{-1}$ as a function of $L$. For each $L$, the p.d.f. is normalized by its maximum $\mathcal {P}_{m}$. The one in the red rectangle is also shown in figure 10(d). (b) Comparison between the measured and modelled specific resistance of the opening for $U=74\ \textrm {m}\ \textrm {s}^{-1}$. This is a close-up view of data presented in figure 8 in order to highlight the smoothness of the modelled $\textrm {Re}(\mathcal {Z})$.

So far, we have seen that good overall predictions of the stability borders and limit cycle amplitudes can be obtained from the deterministic model (4.1) and its corresponding slow-flow dynamics (4.8) over the full range of cavity lengths $L$ and duct velocities $U$, and this with a few scaling laws for defining all the model parameters as function of $U$, $W$ and $c$. However, the stochastic dynamics of the acoustic pressure and the resulting p.d.f.s cannot be predicted with these deterministic approaches, which motivates the developments presented in the next section.

5.1. Intermittency modelled with randomly forced coupled oscillators

The Gaussian-like p.d.f.s of the presumably aeroacoustically stable cases ($L=200$ and 270 mm in figure 13) can justify incorporating additive stochastic forcing to the deterministic model (4.1). Indeed, linearly stable oscillators that are subject to additive white noise forcing display such Gaussian-like p.d.f.s. Furthermore, the p.d.f.s of the acoustic pressure for $L=220$ and 230 mm are typical of marginally stable oscillators and weakly nonlinear self-sustained oscillators subject to additive white noise forcing (Boujo & Noiray Reference Boujo and Noiray2017). This random additive forcing at the side-branch cavity opening can be attributed to the broadband noise generated by the highly turbulent flow in the wind channel and its air supply line. Therefore, we add a Gaussian additive white noise $\xi (t)$ of intensity $\varGamma _\xi$ to the right-hand side of the equation for $u$ in the system of coupled oscillators (4.1).

Now, when $L$ is decreased from 270 mm to 260 and then to 250 mm, the base of the p.d.f. thickens. This indicates that the mean acoustic pressure amplitude increases, which can be explained by the crossing of the supercritical Hopf bifurcation that was discussed in the previous sections. However, these two p.d.f.s present distinct features that clearly differ from the ones of weakly nonlinear self-oscillators subject to additive random forcing only. Indeed, there is a central peak remaining, which is the marker of high probability of low acoustic amplitude periods. In fact this is the signature of intermittency between low- and high-amplitude acoustic oscillations.

As a first remark, the intermittency at play in the present system differs from the one found in systems that are governed by combined subcritical-Hopf and saddle-node bifurcations, and that are subject to purely additive stochastic forcing. A thermoacoustic example of the latter systems was investigated by Bonciolini et al. (Reference Bonciolini, Ebi, Boujo and Noiray2018) for quasi-steady and finite-rate ramping of the bifurcation parameter. These systems exhibit, for a range of bifurcation parameter values, two basins of attraction between which intermittent jumps can be triggered by the additive random forcing, with resulting p.d.f.s having also a prominent central peak. In contrast, the deterministic model derived in the previous sections on the basis of specific impedance measurements clearly show that the present aeroacoustic system features a supercritical-Hopf bifurcation. Consequently, the only possible explanation for the observed intermittency is the presence of parametric noise in the system, which can be linked to the theoretical investigation of Mohamad & Sapsis (Reference Mohamad and Sapsis2015).

Therefore, considering the intense turbulence of the channel flow as well as the low-frequency 3-D dynamics of the recirculation region, we propose to include a stochastic component to our bifurcation parameter $U$, which leads to several terms with stochastic multiplicative forcing. One can note that the low-frequency random fluctuations of the flow at the opening of the deep cavity probably share several features with the ones that were investigated by Basley et al. (Reference Basley, Pastur, Lusseyran, Soria and Delprat2014) in the case of shallow cavities. Their detailed analysis shows that the broadband slow fluctuations of the recirculating flow in shallow cavities, which result from centrifugal instabilities, interfere with the (fast) Rossiter modes. In the present deep cavity configuration, similar 3-D structures may alter the strength of the aeroacoustic coupling at time scales that are long compared to the period of the three-quarter wave acoustic eigenmode. We therefore express the bulk velocity in our low-order model (4.1) as

(5.1)\begin{equation} U=\bar{U}(1+\chi), \end{equation}

where $\chi$ is an Ornstein–Uhlenbeck process obeying the Langevin equation

(5.2)\begin{equation} \dot{\chi}=-\dfrac{\chi}{\tau}+ \zeta_\chi, \end{equation}

with $\tau$ the correlation time of the bifurcation parameter fluctuations and $\zeta _\chi$ a Gaussian white noise of intensity $\varGamma _\chi / \tau ^{2}$. By incorporating this unsteady formulation of the bulk velocity $U$, with a mean value $\bar {U}$ and a standard deviation $\sigma _U$, into the system of coupled oscillators, the parameters $\beta _1$, $\beta _2$, $\mu _1$, $\mu _3$, $\sigma$ and $\omega _r$ become time dependent. Using the model coefficients given in § 4.2, considering $\bar {U}=74\ \textrm {m}\ \textrm {s}^{-1}$ and setting $\varGamma _\xi =5\times 10^{9} \ \textrm {m}^{2}\ \textrm {s}^{-6}$, $\varGamma _\chi =5\times 10^{-3}$ and $\tau =0.05$ s, which corresponds to $\sigma _U=\bar {U}\sqrt {\varGamma _\chi /2\tau }\approx 16\ \textrm {m}\ \textrm {s}^{-1}$, time-domain simulations of (4.1) are performed for $L=245$ and 250 mm. In figure 14, these simulations are compared to the experimental results obtained at $L=255$ and 260 mm. We consider this 10 mm cavity length offset for the comparison because, as can be seen in figure 12, the predicted bifurcation diagram is staggered by approximately 10 mm from the experimental data. As a reminder, the origin of this relatively small offset was discussed at the end of § 4.3. The parameters defining the additive and multiplicative stochastic forcing, i.e. $\varGamma _\xi$, $\varGamma _\chi$ and $\tau$, were empirically adjusted in order to best match the main features of the experimental time traces.

Figure 14. Comparison of time traces and p.d.f.s of the acoustic pressure $p$ and its envelope $B$ between experiments (a,b) and simulations using the coupled stochastic differential equations presented in § 5.1 (c,d). The considered mean channel velocity is $\bar {U}=74\ \textrm {m}\ \textrm {s}^{-1}$, and the cavity lengths are $L=255$, 260, 245 and 250 mm in (a), (b), (c) and (d) respectively. A 10 mm offset is chosen for this comparison because, as shown in figure 12, the predicted bifurcation diagram is staggered by approximately 10 mm from the experimental data.

One can see in figure 14 that the model reproduces very well the p.d.f.s of the acoustic pressure $p$ and its envelope $B$, and especially the central peak of $\mathcal {P}(p)$. The underlying intermittency between periods of low amplitude and bursts of high amplitude is thus well captured by the model.

It is important to mention that the intensity of the fluctuating component of the bifurcation parameter $U$ is high enough, such that the instantaneous bulk velocity alternately takes values in the range for which our model of the aeroacoustic system is linearly stable, and in the one for which it is linearly unstable. Another important point is that the correlation time $\tau$ of the bifurcation parameter fluctuations is long enough for the system to adapt to the variations of $U$. Otherwise, we cannot reproduce with this model the experimentally observed intermittency.

In order to shed light on the conditions leading to these intermittent high-amplitude bursts or low-amplitude periods, a simplified model is scrutinized in the following section.

5.2. Intermittency modelled with a randomly forced van der Pol oscillator

We have established in the previous sections (i) that the present aeroacoustic system is linearly unstable for a sufficiently large bulk flow velocity $U$ and for a range of cavity lengths $L$, (ii) that at the border of the region corresponding to the limit cycles (red line in figure 7a), only supercritical Hopf bifurcations occur and (iii) that intermittency is induced by parametric noise. To gain further insight into the conditions that lead to intermittency, rather than carrying on with the model of coupled oscillators, it is convenient to consider a simpler version of the supercritical Hopf bifurcation with parametric noise. Reckoning with the fact that one of the eigenvalues of (4.1) is very stable (see figure 7b), it can be shown that, in the vicinity of the Hopf bifurcation, our model of two coupled oscillators for describing the aeroacoustic system can be approximated by a single van der Pol oscillator subject to additive and multiplicative stochastic forcing. This van der Pol model is

(5.3)\begin{equation} \ddot{\eta}+ (-2[\bar{\nu}+\chi(t)]+\kappa \eta^{2})\dot{\eta} + \omega_0^{2} \eta = \xi(t), \end{equation}

where $\eta$ is a new state variable that represents the aeroacoustic oscillations, $\omega _0$ is the natural frequency of the oscillator, $\nu (t)=\bar {\nu }+\chi (t)$ is the instantaneous linear growth rate, $\kappa >0$ is the saturation constant and, as in the previous section, $\xi$ is a Gaussian additive white noise of intensity $\varGamma _\xi$, i.e. of autocorrelation $\left \langle \xi \xi _\theta \right \rangle =\varGamma _\xi \delta (\theta )$ with $\delta$ the Dirac delta function, and $\chi$ a coloured Gaussian multiplicative noise, with correlation time $\tau$ and equilibrium variance $\sigma _\nu ^{2}=\varGamma _\chi /2\tau$, and whose Langevin equation is (5.2).

Time marching of this equation is performed for various combinations of mean linear growth rate $\bar {\nu }$, multiplicative noise correlation time $\tau$ and standard deviation $\sigma _\nu$, in order to illustrate the conditions required for a high probability of intermittency. These conditions are set by two non-dimensional parameters: $\bar {\nu }/\sigma _\nu$, which is related to the probability of crossing the stability border, and $\bar {\nu }\tau$, which compares the characteristic growth or decay time of the oscillation amplitude with the correlation time of the instantaneous growth or decay rate.

The natural frequency of the oscillator is set to $f_0=\omega _0/2 {\rm \pi}= 1050$ Hz, which is typical of the present aeroacoustic instability. The mean linear growth rate $\bar {\nu }$ is successively set to $-30$, $-8$, $-4$, 4, 8 and $30\ \textrm {rad}\ \textrm {s}^{-1}$, which are also typical values of the present system in the vicinity of the supercritical Hopf bifurcation, as shown in § 3.2. The saturation constant is set to $\kappa =0.01\ \textrm {s}^{-1}$ for all the simulations such that the resulting oscillation amplitude is of the same order of magnitude as the acoustic pressure in the cavity. The intensity of the additive noise is set to $\varGamma _\xi =10^{12}$, such that the variance of $\eta$ for $\bar {\nu }=-30\ \textrm {rad}\ \textrm {s}^{-1}$ is similar to the variance of the band-pass filtered acoustic signal for $U\simeq 64\ \textrm {m}\ \textrm {s}^{-1}$ and $L=250$ mm. The correlation time $\tau$ and variance $\sigma _\nu ^{2}$ of the multiplicative random forcing $\chi$ are changed for each simulation in order to obtain converged p.d.f.s for $\bar {\nu }/\sigma _\nu =-10$, $-1$, $-0.2$, 0.2, 1 and 10, and for $|\tau \bar {\nu }|=0.1$, 1 and 10. The simulated time is 500 s for $|\tau \bar {\nu }|=10$ and $\bar {\nu }/\sigma _\nu =-0.2$ and 0.2, and 240 s for the other simulated cases, which is sufficiently long to obtain converged statistics.

The results of these temporal simulations are gathered in figure 15, which shows the converged p.d.f.s of the state variable $\eta$. In this system, the parametric noise adds upon the additive forcing, and its effect is scrutinized in the next paragraphs.

Figure 15. (a) Bifurcation diagram of the deterministic problem (right axis, red curve), and characteristic times of the stochastic problem (left axis). These characteristic times are (i) the inverse of the mean growth or decay rate of the oscillation amplitude $1/\bar {\nu }$ indicated as a solid blue line, (ii) the oscillation period $T=2{\rm \pi} /\omega _0$ indicated as a dashed blue line and (iii) the correlation time $\tau$ of the fluctuations of the instantaneous linear growth rate $\nu$. The coloured dots indicate the combinations of $\tau$ and $\bar {\nu }$ which were used for the temporal simulations of (5.3), whose p.d.f.s are presented in (b,c). (b) The p.d.f.s of the linear growth rate. These Gaussians have a mean $\bar {\nu }$ and a standard deviation $\sigma _\nu$. (c) The p.d.f.s of the oscillations for several combinations of $\bar {\nu }/\sigma _\nu$ and $|\tau \bar {\nu }|$. Rare events are associated to heavy tails indicating high-amplitude bursts (e.g. for $\bar {\nu }/\sigma _\nu =-0.2$ and $|\tau \bar {\nu }|=10$), or to a pronounced central peak indicating sporadic quiet periods (e.g. for $\bar {\nu }/\sigma _\nu =1$ and $|\tau \bar {\nu }|=10$). They reflect the intermittency of the system and they are very likely when $|\tau \bar {\nu }| \geqslant {O}(1)$ and $|\bar {\nu }/\sigma _\nu | \leqslant {O}(1)$.

The first necessary condition for intermittency is that the standard deviation $\sigma _\nu$ of the fluctuations of $\nu$ is sufficiently large to have $\nu$ penetrating in $\mathbb {R}^{+}$ (respectively $\mathbb {R}^{-}$) when $\bar {\nu }<0$ (respectively $\bar {\nu }>0$). Indeed, for negative (respectively positive) mean growth rate, if the instantaneous growth rate experiences deep excursions in $\mathbb {R}^{+}$ (respectively in $\mathbb {R}^{-}$), rare events in the form of amplitude bursts (respectively intermittent quiet period) can occur. This necessary condition is fulfilled for $|\bar {\nu }/\sigma _\nu | \leqslant {O}(1)$. It corresponds to the four middle columns of figure 15(c), where heavy tailed p.d.f.s. and pronounced central peaks can be observed.

However, the condition $|\bar {\nu }/\sigma _\nu | \leqslant {O}(1)$ is not sufficient for intermittency to occur. In fact, sporadic excursions of $\nu$ in $\mathbb {R}^{+}$ or $\mathbb {R}^{-}$ must last sufficiently long for the dynamic system to react to this change of stability, i.e. for having a system intermittently unstable (fat tailed p.d.f.) or intermittently stable (marked central peak). Therefore, there is a second necessary condition for observing intermittency: the characteristic relaxation time of the system $1/\bar {\nu }$ must be shorter than the correlation time of the parametric stochastic forcing $\tau$. This condition is satisfied when $|\bar {\nu }\tau | \geqslant {O}(1)$, which corresponds to the two upper rows in figure 15(c). When $|\bar {\nu }\tau |\gg {O}(1)$, the system is subject to slow random fluctuations of the instantaneous linear growth rate $\nu$, which corresponds to quasi-steady change of the bifurcation parameter. When $|\bar {\nu }\tau |={O}(1)$, rare events can still occur but the intermittence is less marked. Finally, if $|\bar {\nu }\tau |\ll {O}(1)$, the fluctuations of the instantaneous growth rate exhibit a correlation that is smaller than the characteristic relaxation time of the oscillation amplitude, i.e. they are too fast to allow the system to adapt to them and the instantaneous attractor cannot be reached.

In summary, the rare events happen (p.d.f.s with fat tails or marked central peak) when the correlation time of the fluctuations of the instantaneous linear growth rate is of the order of, or longer than, the inverse of the mean linear growth rate, and when the standard deviation of these fluctuations is larger than the mean linear growth rate. It is important to note that the model developed in this section can be used for a very wide range of phenomena that can be described by a van der Pol oscillator subject to additive and multiplicative coloured noise.

5.3. Amplitude dynamics: Langevin and Fokker–Planck equations

In order to complement the analysis of carried out in § 5.2, the amplitude equation associated with the stochastic differential equations (5.3) is derived by performing deterministic and stochastic averaging (Krylov & Bogoliubov Reference Krylov and Bogoliubov1936; Stratonovich Reference Stratonovich1967). The derivation is based on the fact that the aeroacoustic system is weakly nonlinear, which implies that the limit cycle is quasi-sinusoidal, and on the fact that the mean linear growth rate satisfies $\bar {\nu } \ll \omega _0$, which means that the amplitude and phase drift of the oscillation slowly vary with respect to the acoustic period $T=2{\rm \pi} /\omega _0$. It is therefore convenient to investigate the system dynamics using the coordinate system $(A, \varphi )$, with

(5.4a,b)\begin{equation} \eta=A \cos (\omega_0 t + \varphi)\quad\text{and}\quad\dot\eta=-A\omega_0\sin(\omega_0 t+\varphi), \end{equation}

and $A=\sqrt {\eta ^{2}+(\dot {\eta }/\omega _0)^{2}}$ and $\varphi =-\arctan (\dot {\eta }/\omega _0\eta )-\omega _0 t$. Furthermore, we make the assumption that the multiplicative Ornstein–Ulhenbeck noise $\chi$, which is governed by the Langevin equation (5.2), exhibits a correlation time $\tau$ that is significantly longer than the acoustic period. In other words, when $\omega _0\tau /2{\rm \pi} \gg 1$, $\chi$ can be considered as constant during one oscillation period and one can apply the averaging process described by e.g. Noiray (Reference Noiray2017). For most of the cases considered in § 5.2, $\tau$ is at least one order of magnitude longer than $T$ (see coloured circles in figure 15a) and this condition is well satisfied. The stochastic averaging procedure leads to a system of equations for the amplitude $A$, the phase drift $\varphi$ and the parametric noise $\chi$ of the form

(5.5)\begin{equation} \dot{\boldsymbol{Y}} = \boldsymbol{F}(\boldsymbol{Y}) + \boldsymbol{B}(\boldsymbol{Y}) \boldsymbol{N}. \end{equation}

In this equation, the dynamics of the system state vector $\boldsymbol {Y}=(A,\varphi ,\chi )^\textrm {T}$ is defined by a deterministic contribution $\boldsymbol {F}(\boldsymbol {Y})$ and by a stochastic forcing term $\boldsymbol {B}(\boldsymbol {Y}) \boldsymbol {N}$, where $\boldsymbol {N}=(\zeta _A, \zeta _\varphi , \zeta _\chi )^\textrm {T}$ is a vector of independent white Gaussian noises. The intensity of $\zeta _A$ and $\zeta _\varphi$ is $\varGamma _\xi / 2\omega _0^{2}$, and that of $\zeta _\chi$ is $\varGamma _\chi / \tau ^{2}$. The deterministic components of this coupled system of Langevin equations are

(5.6a–c)\begin{equation} F_A(\boldsymbol{Y})=\left(\bar{\nu} + \chi \right)A - \frac{\kappa}{8}A^{3} + \dfrac{\varGamma_\xi}{4\omega_0^{2} A},\quad F_\varphi(\boldsymbol{Y})=0,\quad F_\chi(\boldsymbol{Y})=-\dfrac{\chi}{\tau}, \end{equation}

and the stochastic components are given by

(5.7)\begin{equation} \boldsymbol{B}(\boldsymbol{Y}) \boldsymbol{N}= \left(\begin{array}{ccc}1 & 0 & 0 \\0 & A^{-1} & 0 \\0 & 0 & 1 \end{array}\right)\left(\begin{array}{c}\zeta_A\\ \zeta_\varphi\\ \zeta_\chi\end{array}\right). \end{equation}

One can notice that the equation for $A$ does not depend on $\varphi$ and one can therefore focus on the independent system of Langevin equations

(5.8)\begin{gather} \dot{A}=\left(\bar{\nu} + \chi \right)A - \frac{\kappa}{8}A^{3} + \dfrac{\varGamma_\xi}{4\omega_0^{2} A} + \zeta_A, \end{gather}

(5.9)\begin{gather}\dot{\chi}=-\dfrac{\chi}{\tau}+\zeta_\chi, \end{gather}

and write its corresponding two-dimensional Fokker–Planck equation

(5.10)\begin{equation} \dfrac{\partial \mathcal{P}}{\partial t} = -\dfrac{\partial }{\partial A} \left(F_A \mathcal{P}\right)-\dfrac{\partial }{\partial \chi} \left(F_\chi \mathcal{P}\right) + \dfrac{\partial^{2}}{\partial A^{2}} \left( \dfrac{\varGamma_\xi\mathcal{P}}{ 4\omega_0^{2}} \right) + \dfrac{\partial^{2}}{\partial \chi^{2}} \left( \dfrac{\varGamma_\chi \mathcal{P}}{ 2\tau^{2}} \right), \end{equation}

which describes the time evolution of the joint p.d.f. $\mathcal {P}(A,\chi; t)$. This Fokker–Planck equation is solved in the domain bounded by $A \in [0, 250]$ and $\chi \in [-200, 200]$ by using finite differences for the partial derivatives with respect to $A$ and $\chi$, and explicit Euler integration for the time derivative. A Dirichlet boundary condition $\mathcal {P}(A,\chi;t)=0$ is imposed for $A=0$, while the other boundaries are Neumann boundaries with no flux. The initial conditions $\mathcal {P}(A,\chi;0)=\mathcal {P}_{init}(A,\chi )$ are arbitrarily set for $\bar {\nu }\tau =0.1$ to the theoretical joint p.d.f. of the bivariate problem when the two random processes $A$ and $\chi$ are assumed independent, i.e. when $\chi$ is removed from (5.8)

(5.11)\begin{equation} \mathcal{P}_{init}(A,\chi)\propto \exp\left(-\frac{4\omega_0}{\varGamma_\xi}\left[\frac{\bar{\nu} A^{2}}{2}- \frac{\kappa A^{4}}{32}\right]\right)\exp\left(-\frac{\chi^{2}}{\varGamma_\chi/\tau}\right). \end{equation}

The converged steady solutions of this Fokker–Planck equation for $\bar {\nu }\tau =0.1$ and $\bar {\nu }\tau =1$ are presented in figures 16(i) and 16(j) respectively. Note that for $\bar {\nu }\tau =1$, the initial condition is not defined by (5.11), but it is set as the converged solution of the case $\bar {\nu }\tau =0.1$ in order to avoid numerical instabilities. Note also that these two cases correspond to the same set of parameters as the one used to generate the p.d.f.s presented in the second last column of figure 15(c). In this figure, for the sake of clarity, we show $\mathcal {P}(A,\nu )$, which is identical to $\mathcal {P}(A,\chi )$ in the coordinate system $(A,\nu =\bar {\nu }+\chi )$. From now on, we refer to stationary p.d.f.s when the time variable $t$ is not indicated as the argument of $\mathcal {P}$. The marginal p.d.f.s for $A$ and $\nu$ are also presented on the side and on the top of the joint distribution in figures 16(i) and 16(j). One can note that, instead of computing $\int _0^{\infty }\mathcal {P}(A,\nu )\ \textrm {d}A$ to find $\mathcal {P}(\nu )$, this marginal distribution could be directly found analytically: it is given by the Gaussian solution of the Fokker–Planck equation associated with the univariate Ornstein–Uhlenbeck process (5.9), which is proportional to the second exponential term in (5.11). One can clearly see that, when the correlation time $\tau$ of the parametric noise is increased from 12.5 ms to 125 ms while the mean and the standard deviation of the growth rate are kept constant ($\bar {\nu }=\sigma _\nu =8\ \textrm {rad}\ \textrm {s}^{-1}$), the peak of the joint p.d.f. contracts and adopts a more elongated shape that follows the bifurcation diagram of the deterministic system. This behaviour is due to an increased intermittency which is well captured by this Fokker–Planck formalism.

Figure 16. Time-domain simulations of the van der Pol oscillator (5.3) subject to additive and multiplicative stochastic forcing and stationary solutions of the Fokker–Planck equation (5.10) for $\tau \bar {\nu } = 0.1$ (a, b, c, d, i, k) and for $\tau \bar {\nu } = 1$ (e, f, g, h, j, l). (a,e) Time traces of $\nu$ which are governed by (5.2). (b,f) The p.d.f.s of $\nu$ deduced from the time-domain simulations of (5.2). (c,g) Time traces of $\eta$ and their envelope. (d,h) The p.d.f.s of $\eta$ and the envelope of $\eta$, which are deduced from the time-domain simulations of (5.3). (i,j) Stationary solutions $\mathcal {P}(A,\nu )$ of the Fokker–Planck equation (5.10) for $\tau \bar {\nu } = 0.1$ and for $\tau \bar {\nu } = 1$ respectively. The deterministic bifurcation diagram is superimposed as a dashed line, with $A=(8\nu / \kappa )^{1/2}$ for $\nu >0$. The marginal p.d.f.s $\mathcal {P}(A)=\int _{-\infty }^{\infty }\mathcal {P}(A,\nu )\ \textrm {d}\nu$ and $\mathcal {P}(\nu )=\int _0^{\infty }\mathcal {P}(A,\nu )\ \textrm {d}A$ are also shown as solid lines on the top and the side of the contour plot of $\mathcal {P}(A,\nu )$. (k,l) Comparison of $\mathcal {P}(\eta )$ obtained from the time-domain simulations of (5.2) and (5.3), and from the time-domain simulations of (5.10), respectively shown as dashed and solid lines.

Moreover, these p.d.f.s are compared to the ones deduced from time-domain simulations of the stochastically forced van der Pol oscillator governed by (5.3), and of its coloured multiplicative noise obeying equation (5.2). The simulated time traces of $\nu$ and $\eta$ are presented in figures 16(a), 16(c), 16(e) and 16(g) for $\bar {\nu }\tau =0.1$ and $\bar {\nu }\tau =1$. One can clearly see that the correlation time of the instantaneous growth rate is larger for the case of $\bar {\nu }\tau =1$. From these time traces, the histograms of $\nu$, $\eta$ and their envelope are computed to get the p.d.f.s shown in figures 16(b), 16(d), 16(f) and 16(h). One can see that for $\bar {\nu }\tau =1$, at $t\approx 8$ s, the instantaneous growth rate makes a sufficiently deep and sufficiently long excursion into $\mathbb {R}^{-}$ to induce a noticeable decrease of the oscillation amplitude in figure 16(g). This is an example of a sporadic quiet period, which is typical of intermittently stable dynamics, and which leads to an increased probability density function $\mathcal {P}(A)$ at low amplitude compared to the case $\bar {\nu }\tau =0.1$.

Finally, it is interesting to compare the p.d.f.s of $\eta$ obtained from the time-domain simulations of the van der Pol oscillator (5.3) and the ones deduced from the Fokker–Planck equations describing the slow-flow dynamics of this stochastically forced van der Pol oscillator. To that end, one has to deduce $\mathcal {P}(\eta )$ from the numerical solution of the Fokker–Planck equation (5.10), which requires a few steps. First, one assumes that the joint p.d.f. of the multivariate process $(A,\varphi ,\chi )$ is separable and can be written $\mathcal {P}(A,\varphi ,\chi )= \mathcal {P}(A,\chi )\mathcal {P}(\varphi )$ in order to deduce the univariate stationary p.d.f. for $\varphi$. By injecting this ansatz into the Fokker–Planck equation of the trivariate process, considering the stationary solution and the fact that $\mathcal {P}(A,\chi )$ satisfies (5.10), we can deduce that $\textrm {d}^{2} \mathcal {P}(\varphi )/\textrm {d}\varphi ^{2}=0$, which implies that $\mathcal {P}(\varphi )=c_1\varphi +c_2$ with $c_1$ and $c_2$ the integration constants. Considering that $\mathcal {P}(\varphi )$ must be periodic, one can deduce that $c_1=0$, which means that $\varphi$ is uniformly distributed between 0 and $2{\rm \pi}$. Second, we consider the projection of the oscillatory dynamics onto the slowly varying amplitude and phase that are used to derive the above analytical expressions. With the mapping $(A,\varphi )\rightarrow (\eta ,\dot \eta )$ given in (5.4a,b), one can write that $\mathcal {P}(\eta ,\dot \eta )\propto J^{-1}\mathcal {P}(A,\varphi )$, where $J=A\omega _0$ is the absolute value of the determinant of the Jacobian matrix associated with this mapping. Then, considering that $\mathcal {P}(\varphi )=1/2{\rm \pi}$, one can write that $\mathcal {P}(A,\varphi )\propto \int _{-\infty }^{\infty }\mathcal {P}(A,\chi )\ \textrm {d}\chi$, and therefore, that $\mathcal {P}(\eta ,\dot \eta )\propto A^{-1}\int _{-\infty }^{\infty }\mathcal {P}(A,\chi )\ \textrm {d}\chi$. Finally, the univariate p.d.f. for $\eta$ can be obtained from the marginal p.d.f. for $A$ that is deduced from the numerical solution of (5.10) by using $\mathcal {P}(\eta )=\int _{-\infty }^{\infty } \mathcal {P}(\eta ,\dot \eta ) \ \textrm {d}\dot \eta$. The p.d.f.s obtained in this way for $\bar {\nu }\tau =0.1$ and $\bar {\nu }\tau =1$ are shown in figures 16(k) and 16(l) as solid lines. As expected, there is an excellent overlap between these p.d.f.s and the ones from the histograms of the time-domain simulations of (5.3) shown as dashed lines. For the case shown in figure 16(l), the two necessary conditions for intermittency presented in § 5.2 are fulfilled. Indeed, one has $|\bar {\nu }/\sigma _\nu | \leqslant {O}(1)$ and $|\bar {\nu }\tau | \geqslant {O}(1)$, which leads to a significant change of the shape of $\mathcal {P}(\eta )$ compared to the case where $|\bar {\nu }\tau |=0.1$, and this change can be predicted using the Fokker–Planck description of the slow-flow dynamics.