2.1. Flow domain

Figure 1.

Flow-domain sketch. The acoustically damping buffer zones of thickness $l_s$ are indicated by red colour. Inset: contour of mean velocity magnitude $\overline {u}$ with overlaid arrows indicating direction.

Figure 1. Long description

Panel A: A diagram showing the flow over a cavity with labeled sections including Gamma in, Gamma out, Gamma wall, and Gamma slip. The diagram includes dimensions such as H, W, L1, L2, D, and ls. Panel B: An inset showing a contour plot of mean velocity magnitude with arrows indicating the direction of flow. The velocity magnitude is color-coded from blue to red, with a scale ranging from 0 to 70 meters per second.

The two-dimensional flow domain consists of a straight duct that features a closed cavity side branch of rectangular shape, as shown in figure 1. The duct cross-section length is $D=6.2$ cm, while the cavity measures $W=3$ cm in width and its height $H$ is varied in this work between $18$ and $30$ cm, corresponding to $H/D \in [2.9,\,4.8]$ . Air flows inside the duct with a bulk velocity $U_b$ varied between 55 and $83\ \text{m s}^{-1}$ . This choice of parameters yields a set of flows characterised by small Mach numbers $M$ and large Reynolds numbers $Re$ , $U_b/c_0=M\in [0.16,\,0.24]$ , where $c_0=340\ \text{m s}^{-1}$ is the speed of sound and $U_b W/\nu =Re \simeq 1.5\times 10^{5}$ , where $\nu =1.5\times 10^{-5} \ \text{m}^2\,\text{s}^{-1}$ is the kinematic viscosity, while the Prandtl number $ \textit{Pr}=\rho _0 \nu c_p/ \kappa$ is around 0.7, where $\rho _0=1.2\,\text{kg m}^{-3}$ is the density, $c_p=10^3 \ \text{J kg}^{-1}\ \text{K}^{-1}$ the specific heat capacity and $\kappa =2.6\times 10^{-2} \ \text{J s}^{-1}\text{m}^{-1}\text{K}^{-1}$ the thermal conductivity.

2.2. Governing equations

The compressible NS equations are cast in non-dimensional form using the reference values presented above, as in the work of Meliga et al. (Reference Meliga, Sipp and Chomaz2010) and Fani et al. (Reference Fani, Citro, Giannetti and Auteri2018) on noise prediction from flow over a bluff body and a shallow cavity, respectively,

(2.1) \begin{align} \partial _t \varrho + \boldsymbol{\nabla }\boldsymbol{\cdot }(\varrho \boldsymbol{u})=0, \\[-28pt] \nonumber \end{align}

(2.2) \begin{align} \varrho \partial _t \boldsymbol{u} + \boldsymbol{N}(\varrho ,\boldsymbol{u},P)=0, \\[-28pt] \nonumber \end{align}

(2.3) \begin{align} \boldsymbol{N}(\varrho ,\boldsymbol{u},P)=\varrho (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u} + \frac {1}{\gamma M^2} \boldsymbol{\nabla }P - \frac {1}{Re}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }(\boldsymbol{u}), \\[-12pt] \nonumber \end{align}

where $\gamma$ is the ratio of specific heat capacities and $N$ denotes the non-dimensional momentum equation expressed in terms of the non-dimensional density $\rho$ , pressure $P$ , velocity $U$ , and shear stress $\tau$ . Furthermore, the non-dimensionalised energy equation reads

(2.4) \begin{align} \varrho \partial _tT + \boldsymbol{E}(\varrho ,\boldsymbol{u}, T)=0, \end{align}

with

(2.5) \begin{align} \boldsymbol{E}(\varrho ,\boldsymbol{u}, T) = \varrho (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })T + (\gamma -1)\varrho T\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} - \gamma (\gamma -1)\frac {M^2}{Re}\boldsymbol{\tau }(\boldsymbol{u}):\boldsymbol{d}(\boldsymbol{u}) -\frac {\gamma }{Pr Re}\boldsymbol{\nabla} ^2T . \end{align}

Closure is achieved by means of the equation of state $P = \varrho T$ , where the non-dimensional pressure $P$ is defined as $(p-p_0)/(\rho _{0}U_b^2)$ . The solution $\boldsymbol{q}=[\varrho ,\boldsymbol{u},T,P]^{\intercal }$ is decomposed into its time-averaged component $\overline {\boldsymbol{q}}(\boldsymbol{x})$ , coherent fluctuations $\widetilde {\boldsymbol{q}}(\boldsymbol{x},t)$ and turbulent fluctuations $\boldsymbol{q}'(\boldsymbol{x},t)$ following Hussain & Reynolds (Reference Hussain and Reynolds1970), where the different components can be isolated from the complete unsteady flow by phase averaging $(\langle \boldsymbol{\cdot }\rangle )$ and time averaging $(\overline {\boldsymbol{\cdot }})$

(2.6) \begin{align} \boldsymbol{q}=\overline {\boldsymbol{q}}+\widetilde {\boldsymbol{q}}+\boldsymbol{q}' , \\[-28pt] \nonumber \end{align}

(2.7) \begin{align} \overline {\widetilde {\boldsymbol{q}}+\boldsymbol{q}'}=0,\ \langle \boldsymbol{q}\rangle =\overline {\boldsymbol{q}} +\widetilde {\boldsymbol{q}},\ \langle \boldsymbol{q'}\rangle =0. \\[0pt] \nonumber \end{align}

Substituting the decomposition into the governing equations results in a coupled system of equations for the mean flow and the coherent fluctuations and omitting the terms $ \{ \mathcal{O}((q_i-\overline {q}_i)(q_j-\overline {q}_j)) | \ i=1,2,3,4, \ \ j=1,3,4 \}$ (i.e. omitting all terms involving products of two or more fluctuating quantities except for terms involving products of up to two velocity fluctuation components), which is a reasonable approximation for low-Mach-number flows, considering that fluctuations of density and temperature vanish in the incompressible limit. The approximation is also motivated by the study of Boujo et al. (Reference Boujo, Bauerheim and Noiray2018), where frequency responses computed with incompressible LNSE around a LES mean flow closely match the full LES solution of the flow field in response to external forcing applied to the present flow configuration. Notable discrepancies arise at relatively high forcing amplitudes, i.e. when the cavity forcing velocity approaches 3 %–5 % of the mean flow velocity. The conservation of mass, momentum and energy, together with the ideal gas law as closure, are given by

(2.8a) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }(\overline {\rho }\overline {\boldsymbol{u}})= 0, \\[-28pt] \nonumber \end{align}

(2.8b) \begin{align} \boldsymbol{N}(\overline {\varrho },\overline {\boldsymbol{u}})=-\boldsymbol{\nabla }\boldsymbol{\cdot } \big(\overline {\rho } \overline {\widetilde {\boldsymbol{u}} \widetilde {\boldsymbol{u}}} + \overline {\rho }\overline {\boldsymbol{u}' \boldsymbol{u}'} \big), \\[-28pt] \nonumber \end{align}

(2.8c) \begin{align} \boldsymbol{E} \big(\overline {\rho },\overline {\boldsymbol{u}},\overline {T} \big) =- \gamma (\gamma -1) \frac {M^2}{Re} \big[\overline {\boldsymbol{\tau }(\boldsymbol{u}'):\boldsymbol{d}(\boldsymbol{u}')} + \overline {\boldsymbol{\tau }(\widetilde {\boldsymbol{u}}):\boldsymbol{d}(\widetilde {\boldsymbol{u}})}\big ], \\[-28pt] \nonumber \end{align}

(2.8d) \begin{align} \overline {P} - \overline {\varrho } \overline {T} = 0. \\[0pt] \nonumber \end{align}

Accordingly, the LNSE in terms of the coherent fluctuations $\widetilde {\boldsymbol{q}}$ take the form

(2.9a) \begin{align} \partial _t\widetilde {\rho } + L_C(\overline {\rho },\overline {\boldsymbol{u}})\widetilde {\rho }=0, \\[-28pt] \nonumber \end{align}

(2.9b) \begin{align} \overline {\varrho }\partial _t\widetilde {\boldsymbol{u}}+\boldsymbol{L}_N(\overline {\rho },\overline {\boldsymbol{u}})\widetilde {\boldsymbol{u}}= -\boldsymbol{\nabla }\boldsymbol{\cdot } \big(\overline {\rho } \widetilde {\widetilde {\boldsymbol{u}} \widetilde {\boldsymbol{u}}} + \overline {\rho }\widetilde {\ \boldsymbol{u}' \boldsymbol{u}'} \big), \\[-28pt] \nonumber \end{align}

(2.9c) \begin{align} \hspace {3.5cm} \hspace {-2cm} \overline {\varrho }\partial _t \widetilde {T} +L_E(\overline {\rho },\overline {\boldsymbol{u}},\overline {T})\widetilde {T}= - \gamma (\gamma -1) \frac {M^2}{Re} \big[\widetilde {\boldsymbol{\tau }(\widetilde {\boldsymbol{u}}):\boldsymbol{d}(\widetilde {\boldsymbol{u}})} + \widetilde {\boldsymbol{\tau }(\boldsymbol{u}'):\boldsymbol{d}(\boldsymbol{u}')} \big], \\[-28pt] \nonumber \end{align}

(2.9d) \begin{align} \widetilde {P}-\widetilde {\rho }\overline {T}-\overline {\rho }\widetilde {T}=0. \\[0pt] \nonumber \end{align}

The expressions for the linear operators $L_C$ , $\boldsymbol{L}_N$ and $L_E$ can be found in Appendix A. The steady mean flow equations and the equations governing the evolution of the coherent fluctuations can be rewritten compactly using solution-vector notation. Grouping the right-hand sides into the nonlinear forcing terms $[\overline {\boldsymbol{F}}_{\textit{turb}}(\boldsymbol{q}'),\ \widetilde {\boldsymbol{F}}_{\textit{turb}}(\boldsymbol{q}')]^{\intercal }$ and $[\overline {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}}), \ \widetilde {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}})]^{\intercal }$ which cover the effects of turbulent and coherent fluctuations, respectively, and introducing the triple decomposition into (2.8)–(2.9) yields

(2.10a) \begin{align} \boldsymbol{N}_{\textit{tot}}(\overline {\boldsymbol{q}})=\overline {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}})+\overline {\boldsymbol{F}}_{\textit{turb}}(\boldsymbol{q}') , \\[-28pt] \nonumber \end{align}

(2.10b) \begin{align} \partial _t\boldsymbol{B}(\overline {\boldsymbol{q}}) \boldsymbol{\widetilde {q}} +\boldsymbol{L}_{\textit{tot}} (\overline {\boldsymbol{q}},\widetilde {\boldsymbol{q}})\widetilde {\boldsymbol{q}}=\widetilde {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}})+\widetilde {\boldsymbol{F}}_{\textit{turb}}(\boldsymbol{q}') , \\[-12pt] \nonumber \end{align}

where $\boldsymbol{B}(\overline {\boldsymbol{q}})=\text{diag}(1,\overline {\rho },\overline {\rho },0)$ is the mass matrix of the system and $\boldsymbol{L}_{\textit{tot}}=[L_c,\boldsymbol{L}_N,L_E]^{\intercal }$ .The effects of turbulence on the mean flow ( $\overline {\boldsymbol{F}}_{\textit{turb}}(\boldsymbol{q}')$ ) are captured by a $k$ - $\omega$ shear stress transport (SST) turbulence model which introduces the turbulent viscosity (TV) $\mu _t=\rho k/\omega$ , where $k$ is the turbulent kinetic energy and $\omega$ the specific turbulent dissipation rate. The mean Reynolds stresses $\overline {\boldsymbol{u}'\boldsymbol{u}'}$ are linked to the mean traceless strain rate (Boussinesq approximation) through

(2.11) \begin{align} -\rho \overline {\boldsymbol{u}'\boldsymbol{u}'}=\mu _t(\boldsymbol{x})\left [ \boldsymbol{\nabla }\overline {\boldsymbol{u}} + (\boldsymbol{\nabla }\overline {\boldsymbol{u}})^{\intercal } -\frac {2}{3}(\boldsymbol{\nabla }\boldsymbol{\cdot }\overline {\boldsymbol{u}})\boldsymbol{I}\right ] -\frac {2}{3}\rho \overline {k} \boldsymbol{I} , \end{align}

with $\boldsymbol{I}$ being the identity matrix. To link the coherent component of the turbulent Reynolds stresses $\widetilde {\boldsymbol{u}'\boldsymbol{u}'}$ to the traceless strain rate from the coherent fluctuations, the Boussinesq approximation can be used as in the work of Boujo et al. (Reference Boujo, Bauerheim and Noiray2018) and Yim et al. (Reference Yim, Meliga and Gallaire2019)

(2.12) \begin{align} -\rho \widetilde {\boldsymbol{u}'\boldsymbol{u}'}=\mu _t(\boldsymbol{x})\left [\boldsymbol{\nabla }\widetilde {\boldsymbol{u}} + (\boldsymbol{\nabla }\widetilde {\boldsymbol{u}})^{\intercal } -\frac {2}{3}(\boldsymbol{\nabla }\boldsymbol{\cdot }\widetilde {\boldsymbol{u}})\boldsymbol{I}\right ] -\frac {2}{3}\rho \widetilde {k} \boldsymbol{I}. \end{align}

Note that another approximation is introduced, as the turbulence model itself is not linearised, i.e. the coherent fluctuations $\widetilde {k}$ and $\widetilde {\omega }$ are not included in the self-consistent approach. Other works, such as the one of Meliga et al. (Reference Meliga, Pujals and Serre2012), have included the interaction between a coherent TV fluctuation, $\widetilde {\mu _t}$ , and $\widetilde {\boldsymbol{q}}$ . However, previous studies of the present flow configuration by Boujo et al. (Reference Boujo, Bauerheim and Noiray2018) showed that the LNSE method remains accurate despite this approximation, based on comparison with LES.

The governing (2.10) can then be recast as

(2.13a) \begin{align} \boldsymbol{N}_{\textit{tot}}(\overline {\boldsymbol{q}})=\overline {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}}), \\[-28pt] \nonumber \end{align}

(2.13b) \begin{align} \partial _t\boldsymbol{\widetilde {q}} +\boldsymbol{L}_{\textit{tot}} (\overline {\boldsymbol{q}},\widetilde {\boldsymbol{q}})=\widetilde {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}}), \\[0pt] \nonumber \end{align}

where, in the final form of the operators $\boldsymbol{N}_{\textit{tot}}(\overline {\boldsymbol{q}})$ and $\boldsymbol{L}_{\textit{tot}} (\overline {\boldsymbol{q}},\widetilde {\boldsymbol{q}})$ in (2.13), the modified eddy viscosity $\mu +\mu _t$ has been incorporated ( $\mu =\rho _0\nu$ is the dynamic viscosity). The TV is obtained from an incompressible Reynolds-averaged Navier–Stokes (RANS) with $k$ – $\omega$ -SST closure alongside the mean flow variable fields $\overline {\boldsymbol{q}}$ . These fields are then used to construct $\boldsymbol{L}_{\textit{tot}}$ .

2.3. Boundary conditions and mean flow computation

The walls $\varGamma _{\textit{wall}}$ are thermally insulated with no-slip velocity ( $\boldsymbol{u}=0$ ) condition. On the outlet boundary $\varGamma _{\textit{out}}$ , the stress-free condition $[(p+(2/3)\mu\boldsymbol{\nabla }\cdot\boldsymbol{u})\boldsymbol{I}-2\mu\boldsymbol{S}]\cdot\boldsymbol{\hat{x}}=0$ is imposed on the velocity field. To reduce computational expense, a zero shear stress condition $\boldsymbol{u}\boldsymbol{\cdot }\hat {\boldsymbol{y}}=0$ is set at the lower boundary $\varGamma _{\textit{slip}}$ , which does not affect the aeroacoustic dynamics of the side branch cavity on the upper side of the channel.

In order to replicate the experimental conditions of Bourquard et al. (Reference Bourquard, Faure-Beaulieu and Noiray2021) the inlet and outlet should not reflect waves when the equations for the fluctuations (2.14b ) are solved. These anechoic conditions are achieved by combining suitable boundary conditions at $\varGamma _{\textit{in}}$ and $\varGamma _{\textit{out}}$ and buffer zones damping wave propagation. Non-reflecting boundary conditions have been the topic of intense research for the simulation of compressible flows, e.g. Poinsot & Lele (Reference Poinsot and Lele1992), Colonius (Reference Colonius2004), Hu et al. (Reference Hu, Li and Lin2008) and Pirozzoli & Colonius (Reference Pirozzoli and Colonius2013). Given that the duct acts as an acoustic wave guide whose cut-on frequency for transverse mode propagation is significantly above the cavity whistling frequency and that the cavity is relatively far from the inlet and outlet boundaries, the one-dimensional plane wave approximation is used, with the following Riemann invariant wave definition: $R_\pm =u' + p'/{\rho _0(\pm 1+U_0/c)}$ , with $+$ denoting downstream and $-$ upstream propagating waves. Based on the work of Abom (Reference Abom2010), the acoustic boundary conditions are then set to $R_-|_{\varGamma _{\textit{out}}}=0, \ \ R_+|_{\varGamma _{\textit{in}}}=0$ . Acoustic waves are further damped by grid stretching and buffer zones of thickness $l_s$ that feature artificial viscosity and additional dissipative terms.

The mean flow is obtained by using a incompressible two-dimensional RANS solver with a $k$ – $\omega$ -SST turbulence model implemented in Ansys Fluent (Ansys Inc. 2024). Boundary conditions at the inlet $\varGamma _{\textit{in}}$ are uniform streamwise velocity $U_b$ , temperature $T=300$ K and turbulence intensity $5\,\%$ . Here, $\varGamma _{\textit{out}}$ is a pressure outlet and $\varGamma _{\textit{wall}}$ is the no-slip, adiabatic wall. The reader is referred to Appendix B for more details on mesh selection and boundary treatment.

2.4. Eigenvalue problem for a given mean flow

The compressible LNSE (2.13b ) are now turned into an eigenvalue problem. First, the nonlinear terms $\overline {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}})$ and $\widetilde {\boldsymbol{F}}_{\textit{co}}(\widetilde {\boldsymbol{q}})$ are neglected and a solution is sought in the form of an eigenmode expansion. Substituting a solution of the form $\widetilde {\boldsymbol{q}}=\widetilde {\boldsymbol{q}}(\boldsymbol{x})e^{\lambda t}+ \ \textrm{c.c.}$ to (2.13b ), neglecting nonlinear terms in $\widetilde {\boldsymbol{q}}$ and considering a given mean flow $\overline {\boldsymbol{q}}$ one gets

(2.14a) \begin{align} \boldsymbol{N}_{\textit{tot}}(\overline {\boldsymbol{q}})=0 , \\[-28pt] \nonumber \end{align}

(2.14b) \begin{align} \lambda \widetilde {\boldsymbol{q}} +\boldsymbol{L}_{\textit{tot}} (\overline {\boldsymbol{q}},\widetilde {\boldsymbol{q}})=0 . \\[0pt] \nonumber \end{align}

The above eigenvalue problem associated with the LNSE is solved with a finite element method based on the open-source software FreeFEM (Hecht Reference Hecht2012). The mesh of $N_E=42\,521$ triangular elements is constructed using FreeFEM meshing tools (see Appendix B for details and the mesh convergence study). A hybrid interpolation approach is employed, using a combination of P2 elements for velocity and P1 elements for state variables. After discretisation of the LNSE system, its eigenvalues and eigenmodes are computed by the shift-invert method. Solving the LNSE system produces a spectrum which typically exhibits a branch of stable eigenvalues, as shown in figure 2(a) for the case of the RANS mean flow obtained with $U_b=55$ m s $^{-1}$ . Their negative growth rate (imaginary part) is generally smaller when their frequency (real part) gets larger. In this spectrum, less stable eigenvalues can be found above the branch. These eigenvalues correspond to eigenmodes involving intense fluctuations along the shear layer developing at the cavity opening. In figure 2(a), the marginally stable one around $900$ Hz exhibits a fluctuating pressure field which corresponds to the $3/4$ wave eigenmode of the pure acoustic problem, and velocity fluctuations in the shear layer which correspond to the first Kelvin–Helmholtz (K-H) mode of the pure hydrodynamic, i.e. incompressible, problem. The axial component of the velocity field of this mode and its pressure field are shown in figures 3(a) and 3(d). In contrast to the results obtained with the incompressible LNSE by Boujo et al. (Reference Boujo, Bauerheim and Noiray2018), where K-H eigenmodes were always linearly stable, i.e. with negative growth rate, the compressible LNSE considered in this work can lead to linearly unstable, i.e. with positive growth rate, aeroacoustic eigenmodes for some ranges of cavity depth $H$ and bulk flow velocity $U_b$ , in agreement with experimental observations of Bourquard et al. (Reference Bourquard, Faure-Beaulieu and Noiray2021).

Figure 2.

Eigenvalue spectrum of the compressible LNSE problem. $(a)$ The case of $U_b = 55 \text{ m s}^{-1}$ and $H = 25$ cm with TV. All the eigenmodes are linearly stable including the 3/4 wave mode (blue triangle). $(b)$ Spectrum obtained with TV for an inlet turbulence intensity of $5\,\%$ (blue), superimposed on the spectrum calculated with $10\,\%$ turbulent intensity (red) for the RANS mean flow obtained with $U_b = 65 \text{ m s}^{-1}$ and $H = 25$ cm. The change of turbulent intensity does not visibly alter most of the eigenvalues, including the linearly unstable 3/4 wave mode.

Figure 3.

Dominant eigenmode of the compressible LNSE for the mean flow obtained with $U_b=65$ m s $^{-1}$ and $H=25$ cm. Coherent component of fluctuating velocity $(a{-}c)$ , with TV and turbulence intensity of $5\,\%$ $(a$ , $b)$ , with TV and turbulence intensity of $10\,\%$ $(c)$ , and pressure $(d)$ of the 3/4 mode for $5\,\%$ turbulence intensity.

Furthermore, the TV has a strong damping effect on the eigenvalues of the compressible LNSE problem. It substantially reduces the growth rate of all modes, particularly in the high-frequency regime. This is illustrated in figure 2 for the case of $U_b=65$ m s $^{-1}$ , $H=25$ cm treated with TV for the same RANS mean flow. As the RANS equations employ turbulence closure models, the self-consistent model is strictly speaking not parameter free, unlike in the case of Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014). As the turbulent intensity is often only known with high uncertainty, we demonstrate the sensitivity of the eigenvalue spectrum towards this parameter in figure 2(b) by comparing cases of $5\,\%$ and $10\,\%$ turbulence intensity. Notably, the eigenvalue corresponding to the $3/4$ wave mode remains unchanged, while some highly damped eigenvalues in the high-frequency regime are affected by turbulence intensity. We further identify modes at $1600$ and $2000$ Hz that remain stable for both inlet velocities. These modes are visibly separated from the remainder of eigenvalues with significantly lower real values, as they are related to K-H instabilities exciting cavity modes at higher frequencies, as in the case of the $5/4$ wave mode close to $2000$ Hz. The experiments of Bourquard et al. (Reference Bourquard, Faure-Beaulieu and Noiray2021) show that the supercritical Hopf bifurcation which separates cavity resonance condition from aeroacoustic limit cycles occurs close to this bulk flow velocity $U_b$ . Comparing figures 3(b) and 3(c) one can see that TV also modifies the mode shape by reducing the penetration depth of the fluctuation inside the cavity. All the results presented in the remainder of this paper are obtained with TV. Moreover, for all combinations of $U_b$ and $H$ considered in this work, one mode dominates the spectrum of eigenvalues of the compressible LNSE associated with the RANS mean flow. It corresponds to the 3/4 mode of the cavity with an acoustic pressure wavelength equal to approximately 4 $H$ /3.

2.5. Self-consistent saturation model

For cases where the compressible NS equations linearised around the RANS mean flow lead to a linearly unstable $3/4$ wave eigenmode, infinitesimal perturbations exponentially grow and an acoustic limit cycle is ultimately reached. The method to predict its amplitude is now explained. It is based on the self-consistent analysis of hydrodynamic limit cycles proposed by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014). They validated their method with the incompressible laminar self-oscillating wake behind a cylinder and developed it by building upon the work of Stuart (Reference Stuart1958). Assuming that the mode shape does not change with the amplitude, Stuart (Reference Stuart1958) solved an amplitude-dependent LNSE eigenvalue problem by subtracting a scaled mean component of the coherent Reynolds-stress divergence such that a marginally stable base flow solution is found. The scaling factor then provides an estimate of the stable limit-cycle amplitude. However, this approach often produces an incorrect estimate of the amplitude because the shape of the mode exerting the steady forcing on the mean flow usually depends on the oscillation amplitude. An elegant way to overcome this issue was proposed by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014). The saturated state is iteratively obtained by gradually increasing the value of the forcing amplitude and resolving the updated equations at each step. It artificially creates a feedback loop, where steady forcing from the unstable mode (see an example in figure 4) modifies the mean flow, which in turn leads to a modified mode shape and so on until the mode amplitude increment leads to a marginally stable mode. This development mimics the physical feedback mechanism that is responsible for the saturation process of an aeroacoustic limit cycle in the deep cavity by gradually varying the shape of the base flow (and the unstable mode), leading to accurate and robust convergence behaviour. As mentioned in the introduction, the self-consistent saturation analysis performed in this study accounts for both the effects of turbulence and compressibility. Turbulence effects are accounted for through the terms originating from the triple decomposition of the flow field and an eddy-viscosity closure model, as explained in the previous sections, while the effects of compressibility are resolved by the compressible equations for density, temperature and pressure in addition to the momentum equation. The resulting system of equations describing the self-consistent saturation model, including all nonlinear terms up to first order, can be found in Appendix C or in the work of Fani et al. (Reference Fani, Citro, Giannetti and Auteri2018). For the present low-Mach flow cases, one of the terms dominates the others by more than an order of magnitude and is the only one retained in the present model. Note that restricting the model to a subset of first-order terms yields a truncated version of the self-consistent formulation. However, for the present flow configuration, we refer to the work of Boujo et al. (Reference Boujo, Bauerheim and Noiray2018), where an LNSE approach applied to a mean flow obtained from LES accurately predicted coherent fluctuations. Boujo et al. (Reference Boujo, Bauerheim and Noiray2018) has further shown that the leading Reynolds stresses in this configuration project poorly onto the optimal forcing field, motivating the monochromatic approximation adopted in this work. Despite being only approximately self-consistent, we refer to the model as self-consistent, as it can capture the dominant saturation mechanism. The resulting system in solution-vector notation reads

(2.15a) \begin{align} \boldsymbol{N}_{\textit{tot}}\big(\overline {\boldsymbol{q}}^{(j+1)} \big)=-2A^2\textrm{Re} \big\{\overline {\rho }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)*}\boldsymbol{\cdot }\boldsymbol{\nabla }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)}\big\}\boldsymbol{\cdot }[0,\boldsymbol{1},0,0]^{\intercal }, \\[-28pt] \nonumber \end{align}

(2.15b) \begin{align} \lambda \widetilde {\boldsymbol{q}}^{(j)} +\boldsymbol{L}_{\textit{tot}} \big(\overline {\boldsymbol{q}}^{(j)} ,\widetilde {\boldsymbol{q}}^{(j)} \big)=0 , \\[0pt] \nonumber \end{align}

Figure 4.

$(a)$ Horizontal and vertical components of the steady forcing from the Reynolds stress of the coherent fluctuations for $U_b=65$ m s $^{-1}$ and $H=25$ cm at the first iteration. $(b)$ Evolution of the linear growth rate of the unstable aeroacoustic mode during the 6 steps of the self-consistent iterative process.

where $j\in \mathbb{N}^*$ is the iteration step, and $\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)}$ is the velocity component of the unstable eigenmode $\widetilde {\boldsymbol{q}}^{(j)}_{\textit{un}}$ , with $\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)*}$ its complex conjugate. Consistent with the single-harmonic truncation of the harmonic balance method, the term $\tilde {F}_{co}(\tilde {q})$ , which appears on the right-hand side of (2.13b ), is neglected on the right-hand side of (2.15b ). The first mean flow $\overline {\boldsymbol{q}}^{(1)}$ used to obtain $\widetilde {\boldsymbol{q}}^{(1)}_{\textit{un}}$ with (2.15b ) is the RANS solution. The forcing term on the right-hand side of (2.15a ) is presented in figure 4(a) for $U_b=65$ m s $^{-1}$ and $H=25$ cm. We have verified for the momentum equation of the mean flow that the restriction of coherent Reynolds-stress forcing to a single remaining term in (2.15a ) yields errors of $\approx 2\,\%$ based on a comparison of the $L_2$ norms of the forcing terms presented in Appendix C. In figure 4(b), one can see that the evolution of the linear growth rate of the unstable mode during the iteration process leads to marginal stability after six steps. The computation time for this result is of the order of $20$ min using parallelisation with $4$ cores in Ansys Fluent. The constant $A$ in (2.15a ) can be interpreted as the mode amplitude, and as proposed by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014) in their algorithm, its value is determined by requiring the value of $A_{\!f}^2=A^2 \Vert {2\textrm{Re} \{\overline {\rho }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)*}\boldsymbol{\cdot }\boldsymbol{\nabla }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)}}\Vert /\Vert {\overline {\rho }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)}}\Vert ^2$ to be fixed for ensuring convergence of the iteration process. Introducing $A_{\!f}$ is not strictly necessary for the iterative process. Its use explicitly shows the normalisation of the eigenvectors using the metric of the space over which they are defined. In the present discretised LNSE operator space the dot product of two vectors is defined using the mass matrix of the problem $\boldsymbol{B}(\overline {\boldsymbol{q}})$ with $\Vert {\tilde {\boldsymbol{q}}}\Vert ^2=\tilde {\boldsymbol{q}}^{\dagger }\boldsymbol{B}(\overline {\boldsymbol{q}})\tilde {\boldsymbol{q}}\simeq \iint _\varOmega \langle \overline {\rho }\tilde {u}_i^*\tilde {u}_i\rangle \text{d}\varOmega =1$ . An adequate tuning of $A_{\!f}$ is needed for ensuring that a marginally stable solution (vanishing growth rate) is reached and for minimising the number of steps (see an example in figure 5). Following Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2014), under-relaxation is also used to ensure robust convergence by updating the mean flow with $\overline {\boldsymbol{u}}^{(j+1)}=\gamma \overline {\boldsymbol{u}}^{(j-1)} +(1 - \gamma ) \overline {\boldsymbol{u}}^{(j)} , \ \text{for }j\geqslant 2$ , where $\gamma \in (0,1)$ is the under-relaxation factor. As shown by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2015), smaller values of $\gamma$ increase the robustness of convergence. Likewise, the number of iterations required to achieve a given relative change between successive iterations, $R_q^j = \|\overline {\boldsymbol{q}}^{(j+1)}-\overline {\boldsymbol{q}}^{(j)}\|/\|\overline {\boldsymbol{q}}^{(j)}\|$ , increases for smaller $\gamma$ . Typical under-relaxation factors in this work lie in the range $\gamma \in [0.1, 0.35]$ and remain constant throughout the iterations.

Figure 5.

$(a{-}c)$ Convergence of the growth rate during the self-consistent model iterations for different values of $A_{\!f}$ and for $U_b=68$ m s $^{-1}$ and $H=25$ cm. $(d)$ Final growth rate for monotonically increasing $A_{\!f}$ . The marginal stability, i.e. $\sigma =0$ rad s $^{-1}$ , is obtained for $A_{\!f}=0.3$ in this case.

Figure 5. Long description

The image contains four panels of line graphs. Panel A, B, and C show the convergence of growth rates during self-consistent model iterations. The x-axis represents iterations, and the y-axis represents growth rate in rad s^-1. Each panel shows a line with diamond markers indicating data points. Panel D shows the final growth rate for monotonically increasing values of A_f. The x-axis represents A_f, and the y-axis represents the final growth rate in rad s^-1. The graph shows a line with square markers indicating data points. The marginal stability is obtained for a specific value of A_f in this case.

Figure 6.

Iterative scheme for computation of the transient evolution of the saturation amplitude $A(T)$ .

2.6. Transient evolution of the limit-cycle amplitude

Instead of tuning the parameter $A_{\!f}$ manually, as shown in figure 5, in the following, we follow the approach first presented by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2015) to compute the transient evolution of the modal amplitude $A$ and eventually its saturation value. Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2015) have computed the transient evolution of a limit cycle for a marginally stable incompressible flow around a cylinder using a self-consistent model. The prediction of this model has been validated against DNSs, yielding good agreement in the final limit-cycle amplitude $A$ and the spatial–temporal evolution of the saturated flow field. We therefore investigate the applicability of the transient self-consistent saturation model of Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2015) to the low-Mach compressible turbulent flow in this work. The iterative scheme for computing the transient evolution of $A(T)$ is illustrated in figure 6. Note that, the inner loop structure is consistent with the saturation model in § 2.5 for a fixed value $A_{\!f}$ .

The evolution of the limit-cycle amplitude $A$ is described on a slow time scale $T$ governed by the modal growth rate $\sigma$ . The core idea is to iteratively adjust the limit-cycle amplitude $A(\sigma )$ in an outer loop, while performing the iterative updates of the mean flow variables $\overline {\boldsymbol{q}}$ for a fixed amplitude in each time step. Eventually, the amplitude $A(\sigma )$ reaches its saturated value, for which the growth rate vanishes, i.e. $\sigma = 0$ . Employing the quasi-steady assumption, the temporal evolution of the limit-cycle amplitude $A$ is given by the following ordinary differential equation (ODE) with an arbitrarily chosen initial amplitude $A_0$ :

(2.16) \begin{align} \frac {\mathrm{d}A(T)}{\mathrm{d}T} = A(T)\,\sigma (A)\quad A(0) = A_0\,. \end{align}

The ODE is discretised in time using a first-order implicit time-stepping scheme. For each $A(T)$ , the iterative procedure in (2.15) is employed until convergence is achieved. We define the convergence based on the relative change between consecutive eigenvalues, $R_\lambda ^j = \|\lambda ^{j+1}-\lambda ^{j}\|/\|\lambda ^j\| \lt \mathrm{tol}$ for $\mathrm{tol} = 2.5\times 10^{-4}$ . Once convergence is achieved, $A(T)$ is updated in time and a new value of $A_{\!f}(T)$ is obtained accordingly using the relation $A_{\!f}^2=A^2 \Vert {2\textrm{Re} \{\overline {\rho }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)*}\boldsymbol{\cdot }\boldsymbol{\nabla }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)}}\Vert /\Vert {\overline {\rho }\widetilde {\boldsymbol{u}}_{\textit{un}}^{(j)}}\Vert ^2$ , with $A_{\!f}$ being fixed within iterations of the inner loop in accordance with § 2.5.

Figure 7 illustrates the transient self-consistent saturation model over a duration of two time steps in the outer loop. Figure 7(a) shows the evolution of the modal growth rate $\sigma$ and the frequency $f$ . The respective markers indicate a new time step in the outer loop, which is consistent with the residual $R_\lambda^{j}$ in figure 7(c) falling below the specified tolerance, thereby defining convergence of the inner loop. Along with the decreasing $\sigma$ , a modulation towards higher frequencies is observed within each time step. The coherent fluctuation amplitude $A$ of the unstable mode increases continuously within each time step and between time steps, as shown in figure 7(b). The value of $A_{\!f}(T)$ (not shown) remains constant over a time step and is updated based on the converged values of $A(T)$ at the end of the inner iteration loops. It eventually converges to $A_{\!f} = 0.3$ , in accordance with figure 5, once the amplitude is saturated. Note that $A_{\!f}(0) = 0.01$ in the present case, but it can be selected arbitrarily as any value greater than zero to enable the unstable mode to grow as proposed by Mantič-Lugo et al. (Reference Mantič-Lugo, Arratia and Gallaire2015).

Figure 7.

Evolution of (a) growth rate $\sigma$ and frequency, (b) coherent fluctuation unstable mode amplitude $A$ , (c) residual $R_\lambda ^j$ with convergence criterion (dashed line) within two time steps in the outer loop $A(T)$ .

Figure 7. Long description

Panel A: A line graph shows the evolution of growth rate (in radians per second) and frequency (in hertz) over iterations. The growth rate is represented by a blue line with square markers, while the frequency is represented by an orange line with triangle markers. Both lines show fluctuations over the iterations. Panel B: A line graph depicts the evolution of the coherent fluctuation unstable mode amplitude (in meters per second) over iterations. The amplitude increases stepwise across the iterations. Panel C: A line graph illustrates the evolution of the residual with a convergence criterion (dashed line) within two time steps in the outer loop. The residual values fluctuate around the convergence criterion.

Figure 8.

Evolution of streamwise velocity of the unstable mode over the course of the iterative process for $U_b=65$ m s $^{-1}$ and $H=25$ cm.

Figure 9.

Evolution of streamwise Reynolds-stress divergence of the unstable mode (forcing term of the mean flow equations) over the course of the iterative process for $U_b=65$ m s $^{-1}$ and $H=25$ cm.