2.1. Problem description

A depiction of the flow–acoustic interaction within tandem deep cavities is illustrated in figure 1(a): concurrent shear layers are formed when flow passes the entrances of two cavities; the internal acoustic forcing can originate from the excited acoustic eigenmodes. Similar to a single cavity, the shear layer perturbations can excite the quasi-trapped acoustic modes of each cavity and the acoustic oscillation amplifies the shear layers in a feedback way. Differently, more complex interactions appear as the two isolated cavities are connected by the global acoustic eigenmode. This complexity has been confirmed by acoustic measurements of Ziada & Bühlmann (Reference Ziada and Bühlmann1992), Okuyama et al. (Reference Okuyama, Tamura, Takahashi, Ohtsuka and Tsubaki2012) and Tonon, Willems & Hirschberg (Reference Tonon, Willems and Hirschberg2011b). The excited acoustic pulsations would not only be trapped inside the cavities, but also propagate outside along the main duct. More attention should be paid on the far upstream and downstream equipment. These clarifications can be confirmed by exploring the acoustic eigenmodes using acoustic modal analysis to solve the acoustic Helmholtz equation (Barbieri & Barbieri Reference Barbieri and Barbieri2006). In the acoustic modal analysis, tetrahedral elements were generated to discretize the acoustic domain, ensuring at least five grid nodes per acoustic wavelength. The inlet and outlet were treated as non-reflective acoustic soft boundaries, allowing acoustic waves to pass. Results are illustrated in figure 1(b), in which four configurations with ${\rm \pi} = L/{D_m} = 1,\; 4,\; 8,\; 16$ are selected for demonstration. Here, L is the distance between the vertical centrelines of the tandem deep cavities and ${D_m} = 40\ \textrm{mm}$ is the size of the square cross-section of the main duct. The length, depth and width of the two cavities are equally fixed at ${L_c} = 0.8{D_m}$, ${D_c} = 4{D_m}$ and ${W_c} = {D_m}$. Here, a depth-to-length ratio of 5 : 1 is utilized to effectively facilitate the excitation of the depth-oriented acoustic mode. A spanwise ratio of 1 : 1 between the main duct and tandem deep cavities promotes a quasi-two-dimensional flow feature around the centreplane of the main duct. With ${\rm \pi} = 1,4$ of close arrangement, their first-order eigenmodes demonstrate an antisymmetric standing-wave feature, as an out-of-phase acoustic pressure distribution is identified inside the tandem deep cavities. With ${\rm \pi} = 8,\; 16$, their first-order eigenmodes demonstrate a symmetric standing-wave feature, as an in-phase acoustic pressure distribution is identified inside the tandem cavities even in faraway arrangements, while an out-of-phase acoustic pressure is identified at the middle section of the main duct. In comparison, the second-order eigenmodes are harmonics of the first-order eigenmodes, but the acoustic pressure distribution is not a purely standing-wave feature due to propagation along the main duct.

Figure 1. (a) Schematic of flow–acoustic interaction within tandem deep cavities; (b) contour plots of the acoustic eigenmodes corresponding to the four configurations with ${\rm \pi} = 1,\; 4,\; 8,\; 16$; the first row represents the first-order mode and the second row represents the second-order mode.

2.2. Acoustic measurements

A schematic diagram of the experimental set-up for acoustic measurements of the ducted tandem deep cavities is illustrated in figure 2(a). This set-up includes the construction of an open-circuit wind channel. To attain a low background noise level and maintain a stable air temperature within the wind channel, the air-breathing running mode of the downstream fan was activated. As part of the set-up, a sizable silencer was constructed and integrated. To attain a consistent inflow characterized by low turbulence intensity, we employed a diffuser with a 6 : 1 contraction ratio, along with the incorporation of a honeycomb layer and two wire-mesh layers at the inlet section. Two acoustic liners, optimized for frequencies exceeding 200 Hz, were crafted and positioned both upstream and downstream of the test section to mitigate the reflection of acoustic waves within the wind channel. The total length of the experimental set-up extended to a measurement of 8 metres. Here, we focus on two variables, i.e. the mainstream Reynolds number $(R{e_{in}})$ and the distance ratio between the two deep cavities. The range of $R{e_{in}}$ varies from $0.2 \times {10^5}$ to $2 \times {10^5}$ considering the side length $({D_m})$ of the mainstream duct and the mainstream velocity (${U_{in}}$: 7.5–75 m s⁻¹). The maximum mainstream Mach number $(M{a_{in}})$ remains below 0.3, indicative of an incompressible air environment within the wind channel. A fully developed turbulent flow profile can be achieved as the leading edge of the upstream deep cavity located almost $50{D_m}$ downstream of the inlet diffuser. As a result, the inlet turbulence intensity was measured to be below 2 % using a hot-wire measurement; the maximum total sound pressure level was below 70 dB with high-speed running of the fan, far less than the acoustic oscillation intensity. A dynamic pressure transducer array with a total of 11 PCB-103B02 sensors was flush mounted at the backplate of the ducted tandem deep cavities to detect the frequency, amplitude and waveform of the excited acoustic modes. The pressure signals were simultaneously measured by an NI-cDAQ-9172 data acquisition equipment at a sampling frequency of 51.2 kHz. For demonstration purposes, the monitoring location was chosen to be at the endplate of the upstream cavity, which corresponds to the antinode of the acoustic eigenmode. Since this position is distant from the mainstream flow and acts as a flow dead zone, the measured pressure fluctuations can be considered as purely acoustic pulsations.

Figure 2. (a) Experimental set-up of acoustic measurements for tandem deep cavities; (b) plots of r.m.s. of wall pressure fluctuation as a function of mainstream Reynolds number, the selected monitor is locating at the endplate of upstream deep cavity; (c) plots of SPSLs; (d) spectrogram for the four configurations with ${\rm \pi} = 1,\; 4,\; 8,\; 16$.

In figure 2(b), we plotted the root mean square (r.m.s.) of the measured acoustic pulsations, normalized by the inlet dynamic pressure head, as a function of the mainstream Reynolds number. Within specific ranges of $R{e_{in}}$, the flow passing through the tandem cavities can induce acoustic resonances owing to the aeroacoustic feedback occurring between concurrent shear layers and their associated natural acoustic eigenmodes. Comparatively, the configuration with ${\rm \pi} = 8$ gives rise to the strongest acoustic pulsations, then the configuration with ${\rm \pi} = 1$. As for ${\rm \pi} = 8$, an obvious range of $0.8 \times {10^5} \le R{e_{in}} \le 1.7 \times {10^5}$ can be identified with essentially intensified acoustic pulsations; the maximum amplitude reaches two times of the inlet dynamic pressure head at $R{e_{in}} = 1.17 \times {10^5}$. It is noteworthy that the resultant acoustic particle velocity ratio reaches ${u_a}/{U_{in}} = 0.76$, which corresponds to the high-amplitude regime according to Bruggeman et al. (Reference Bruggeman, Hirschberg, van Dongen, Wijnands and Gorter1991) and Tonon et al. (Reference Tonon, Hirschberg, Golliard and Ziada2011a). As for ${\rm \pi} = 1$, the Reynolds number range with intensified acoustic pulsation amplitudes remains similar to the configuration with ${\rm \pi} = 8$, the maximum amplitude reaches 1.5 times of the inlet dynamic pressure head at $R{e_{in}} = 1.03 \times {10^5}$ and the resultant acoustic particle velocity ratio reaches ${u_a}/{U_{in}} = 0.57$ within the high-amplitude regime. In this configuration, another peak value of 0.8 times the inlet dynamic pressure head, which has an acoustic particle velocity ratio of 0.38, can also be found at ${R{e_{in}} = 1.29 \times {10^5}}$. As for ${\rm \pi} = 4,\; 16$, the excited acoustic pressure pulsations were relatively attenuated, the maximum values are below 0.4 times the inlet dynamic pressure head and the acoustic particle velocity ratios drop to ${u_a}/{U_{in}} < {O}({10^{ - 1}})$, corresponding to a moderate-amplitude regime. The amplitude discrepancy can be explained in view of the acoustic damping coefficient $(\xi )$ of the excited acoustic eigenmodes. As the acoustic eigenmodes were determined by solving the eigenvalue problem of the acoustic Helmholtz equation, one can obtain the real part and imaginary parts of their eigenvalues $({\lambda _a} = a + bi)$. The real part (a) represents the eigenfrequency of an acoustic eigenmode, while the imaginary part can be used to calculate the acoustic damping coefficient $(\xi = b/\sqrt {{a^2} + {b^2}} )$. Here, a larger damping coefficient represents more potential dissipation within an acoustic system. In our study, the derived damping coefficients for a tandem deep cavity with close and half-wavelength arrangements are ${\xi _{{\rm \pi} = 1}} = 0.015$, ${\xi _{{\rm \pi} = 8}} = 0.0012$, respectively. Consequently, the half-wavelength arrangement is found to induce the most robust acoustic pulsations.

In terms of the frequency information, the sound pressure spectrum levels (SPSLs) were calculated through the fast Fourier transformation method to delineate the acoustic spectra at the pressure sensors. The spectrograms in figure 2(d) illustrate that the fundamental frequencies of the above-mentioned aeroacoustic instabilities are 480 Hz for ${\rm \pi} = 1$ and 520 Hz for ${\rm \pi} = 8$, which are close to the first-order acoustic eigenmodes in figure 1(b). Beyond this fundamental frequency, several high-order harmonic frequencies are also highlighted in the spectrograms. In accordance with Ziada & Lafon (Reference Ziada and Lafon2014), such a frequency lock-in phenomenon confirms the occurrence of flow–acoustic resonance by the natural acoustic eigenmodes. Four quantitative plots with respect to the intensified acoustic pulsations are plotted in figure 2(c). The discrepancies between the experimental measurements and acoustic modal analysis are around 2 %, which could be attributed to the turbulence interaction effect. Accounting for the flow convection effect tends to marginally elevate the predicted eigenfrequencies. In addition, the slight frequency changes of the measured pressure pulsations at different Reynolds numbers may be attributed to the transport speed of the shear layer vortex, according to Tonon et al. (Reference Tonon, Hirschberg, Golliard and Ziada2011a).

In addition to frequency lock-in phenomenon at a single pressure monitor, we also identified the acoustic waveform from the dynamic pressure sensor array and explored their spatio-temporal evolution, which are illustrated in figure 3. The results further confirm the correspondence of the excited acoustic pulsations with the natural acoustic eigenmodes of the duct–cavity configurations. Here, the three flow conditions with ${\rm \pi} = 1$, $R{e_{in}} = 1.03 \times {10^5}$, ${\rm \pi} = 1$, $R{e_{in}} = 1.29 \times {10^5}$ and ${\rm \pi} = 8$, $R{e_{in}} = 1.18 \times {10^5}$ were analysed. Sinusoidal features of the identified acoustic waveform for the two configurations are highly consistent with the acoustic pressure distribution of their first-order acoustic eigenmodes. A plausible interpretation of the discrepancies in amplitudes will be discussed using the following numerical simulation results. A comparative view between the two configurations demonstrates that anti-phase propagation motions of the excited acoustic waves can be found in the tandem deep cavities in close proximity; while closely in-phase propagation motions are found for the tandem deep cavities in the latter configuration. Especially, the acoustic pressure pulsations in the lengthy main duct demonstrate anti-phase propagation motions with the above-mentioned synchronous tendencies between the faraway-arranged tandem deep cavities. Such an interesting special phenomenon was attributed the equivalence between the distance and sum of the tandem deep cavities, serving as a half-wavelength arrangement as their distance is equal to the half of the acoustic wavelength ($\lambda = {c_a}/{f_a} \approx 640\ \textrm{mm}$, with ${c_a}$ the sound speed and ${f_a}$ the eigenfrequency of the first-order acoustic eigenmode). These observations, in terms of the intensified pulsations, locked frequencies and in-phase or out-of-phase waveform of the measured acoustic pressure pulsations, strongly validate the occurrence of resonance between the ongoing concurrent shear layers and the natural acoustic eigenmodes.

Figure 3. Waveform identification (a), spatio-temporal evolution of the flow-excited acoustic pulsations along the ducted tandem deep cavities with ${\rm \pi} = 1$, $R{e_{in}} = 1.03 \times {10^5}$ (b) and $R{e_{in}} = 1.29 \times {10^5}$ (b); and ${\rm \pi} = 8$, $R{e_{in}} = 1.18 \times {10^5}$ (d,e). The white square symbols in (a) and (d) denote the sensor locations; in (b,c) and (e), the abscissa denotes the sensor number, starting at the endplate of the upstream cavity, the time intervals were sampled equally within one acoustic period corresponding to the fundamental frequency.

3.1. High-order spectral/hp element method

The flow–acoustic interaction problem can be described by the nonlinear Navier–Stokes equations, containing the continuity, momentum and energy equations. In the present study, both the incompressible and compressible governing equations are solved for comparison, which artificially highlights the acoustic forcing effect on the turbulent shear layers. A spectral/hp element discretization method, which is implemented in the opensource package of Nektar++ (Cantwell et al. Reference Cantwell2015; Moxey et al. Reference Moxey2020), is used to solve the governing equations. This discretization approach combines the inherent geometric flexibility of the finite-element method with the exponential convergence and superior accuracy of the spectral method. This is realized by employing high-order Lagrangian interpolations in the finite-element formulation's expansion base function. Consequently, resolution convergence can be achieved by refining the mesh distribution (h-type refinement), or by increasing the interpolation order, N_p, within the basis expansion (p-type refinement). Spectral/hp element discretization serves as the underpinning approximation for both continuous Galerkin (CG) and discontinuous Galerkin (DG) formulations. Fundamental information of the spectral/hp method can be found in Karniadakis & Sherwin (Reference Karniadakis and Sherwin2005), and an outline is provided herein. In general, the solution domain $\textrm{(}\varOmega \textrm{)}$ is segregated into ${N_e}$ non-overlapping elements $({\varOmega _e})$. Within each of these elements, we then seek approximate solutions, ${\boldsymbol{u}^\delta }$, belonging to the finite-element space

(3.1)\begin{equation}{\chi ^P} = \{ \phi \in {L^2}(\varOmega )|\phi {|_{{\varOmega _e}}} \in {\mathrm{{\mathcal{P}}}_P}({\varOmega _e})\} ,\end{equation}

where ${\mathrm{{\mathcal{P}}}_P}({\varOmega _e})$ is the space of polynomials of order N_p defined in the element ${\varOmega _e}$.

In this work, the incompressible Navier–Stokes equations are solved via a velocity correction scheme where velocity and pressure are typically decoupled. A weak pressure Poisson problem is formulated by taking the inner product over the domain ${\varOmega _e}$ with respect to the gradient of test function, $\boldsymbol{\nabla }v$, i.e.

(3.2)\begin{align}\int_{{\varOmega _e}} {\boldsymbol{\nabla }v\boldsymbol{\cdot }\frac{{\partial {\boldsymbol{u}^\delta }}}{{\partial t}}\,\textrm{d}\kern0.06em \boldsymbol{x}} + \int_{{\varOmega _e}} {\boldsymbol{\nabla }v\boldsymbol{\cdot }\mathrm{{\mathcal{N}}}({\boldsymbol{u}^\delta })\,\textrm{d}\kern0.06em \boldsymbol{x}} =- \int_{{\varOmega _e}} {\boldsymbol{\nabla }v\boldsymbol{\cdot }\boldsymbol{\nabla }p\,\textrm{d}\kern0.06em \boldsymbol{x}} + \int_{{\varOmega _e}} {\boldsymbol{\nabla }v\boldsymbol{\cdot }\nu {\nabla ^2}{\boldsymbol{u}^\delta }\,\textrm{d}\kern0.06em \boldsymbol{x}} ,\end{align}

where $\mathrm{{\mathcal{N}}}({\boldsymbol{u}^\delta }) = {u^\delta }\boldsymbol{\cdot }\boldsymbol{\nabla }{u^\delta }$. Using the identity ${\nabla ^2}{u^\delta } =- \boldsymbol{\nabla } \times \boldsymbol{\nabla } \times {u^\delta } + \boldsymbol{\nabla }(\boldsymbol{\nabla }\boldsymbol{\cdot }{u^\delta })$, we can enforce the divergence to be zero by setting the last term to zero. By integrating the first, second and last terms in (3.2) by parts, we obtain the weak pressure equation

(3.3)\begin{equation}\begin{array}{ccccc} & \int_{{\varOmega _e}} {\boldsymbol{\nabla }v\boldsymbol{\cdot }\boldsymbol{\nabla }{p^{n + 1}}\,\textrm{d}\kern0.06em x} = \int_{{\varOmega _e}} {v\boldsymbol{\nabla }\boldsymbol{\cdot }\left( {{{\dfrac{{\partial {u^\delta }}}{{\partial t}}}^{n + 1}} + \mathrm{{\mathcal{N}}}{{(u)}^{n + 1}}} \right)\textrm{d}\kern0.06em x} \\ & \quad - \int_{\partial {\varOmega _e}} {v\left( {{{\dfrac{{\partial u}}{{\partial t}}}^{n + 1}} + \mathrm{{\mathcal{N}}}{{(u)}^{n + 1}} + \nu \boldsymbol{\nabla } \times \boldsymbol{\nabla } \times {u^{n + 1}})} \right)\boldsymbol{\cdot }\boldsymbol{n}\,\textrm{d}\kern0.06em x} , \end{array}\end{equation}

where $\partial {\varOmega _e}$ is the boundary of the elements, $\boldsymbol{n}$ indicates the normal to the element and the superscript $n + 1$ denotes the time step we aim to compute. By using a backward approximation of the time derivative, along with a consistent extrapolation for the nonlinear term $\mathrm{{\mathcal{N}}}{({u^\delta })^{n + 1}}$, the pressure at time step $n + 1$ can be calculated. Then, the velocity at step $n + 1$ can be solved by using the pressure just calculated at $n + 1$. For compressible simulation, we substitute ${\boldsymbol{u}^\delta } \in {\chi ^P}$ into the compressible Navier–Stokes equations, multiply the resulting equation by a test function $\boldsymbol{\nu } \in {\chi ^P}$ and apply the divergence theorem to obtain the weak form

(3.4) \begin{align}&\int_{{\varOmega _e}} {\frac{{\partial {\boldsymbol{u}^\delta }}}{{\partial t}}\boldsymbol{\cdot }\boldsymbol{\nu }\,\textrm{d}\kern0.06em \boldsymbol{x}} + \int_{\partial {\varOmega _e}} {(\boldsymbol{f}_n^c({\boldsymbol{u}^\delta }) - \boldsymbol{f}_n^v({\boldsymbol{u}^\delta },\boldsymbol{\nabla }{\boldsymbol{u}^\delta }))\boldsymbol{\cdot }\boldsymbol{\nu }\,\textrm{d}s}\notag\\ &\quad - \int_{\partial {\varOmega _e}} {(\boldsymbol{f}_i^c({\boldsymbol{u}^\delta }) - \boldsymbol{f}_i^v({\boldsymbol{u}^\delta },\boldsymbol{\nabla }{\boldsymbol{u}^\delta }))\boldsymbol{\cdot }\frac{{\partial \boldsymbol{\nu }}}{{\partial {x_i}}}\,\textrm{d}\kern0.06em \boldsymbol{x}} = 0.\end{align}

The subscript n denotes the normal flux on the boundary of element $(\partial {\varOmega _e})$. Equation (3.4) currently encompasses operations exclusive to local elements. However, to keep the solution within a space of bounded variation and to allow information to propagate between elements requires some elemental coupling. In a CG setting, the approximate solution adheres to the C₀-continuous requirement by ensuring the polynomial expansion is continuous across elemental boundaries. In a DG setting, the normal fluxes appearing in the boundary integral are replaced by unique numerical fluxes to ensure appropriate flux quantities are continuous between elements.

3.2. Implicit large eddy simulation

Currently, employing direct numerical simulation of the Navier–Stokes equations for high Reynolds number flows remains unaffordable. An alternative approach is LES, in which the larger-scale motions are resolved, whereas the smaller scales are modelled. The underlying rationale for LES is underpinned by the fact that the larger scales carry the majority of the flow energy, owing to their magnitude and strength, and are responsible for most transport mechanisms; hence, they should be precisely simulated. Conversely, the smaller scales, which have a comparatively minor influence on the mean flow, can be approximated. Furthermore, these scales tend to be more uniform and isotropic, making them more amenable to modelling than the larger scales. In this study, the implicit-LES (iLES) methodology will be adopted, interpreting the unresolved eddies through numerical dissipation. We employ a spectral vanishing viscosity (SVV) in the iLES to carry out a reasonable numerical representation of under-resolved turbulence. The SVV technique was first introduced into the spectral Fourier method by Tadmor (Reference Tadmor1989) with the primary objective of suppressing high-frequency oscillations without compromising the physical characteristics of the flow at lower frequencies. This ensures stability in the simulations while retaining solution accuracy. Kirby & Sherwin (Reference Kirby and Sherwin2006) introduced a multi-dimensional SVV operator within the solution domain. In this approach, the following term is added to the governing equations to achieve the artificial SVV dissipation:

(3.5)\begin{equation}{S_{VV}} = {v_{svv}}\mathop \sum \limits_{i = 1}^d \frac{\partial }{{\partial {x_i}}}\left[ {{Q^d}\ast \frac{{\partial \boldsymbol{u}}}{{\partial {x_i}}}} \right]\!,\end{equation}

where ${v_{svv}}$ represents the effective dissipation strength and ${Q^d}$ denotes the SVV kernel that determines the manner in which dissipation impacts the elementwise polynomial modes within the spectral/hp element discretization. The numerical study here is carried out based on the SVV-iLES method with a DG kernel in the spectral/hp space. The dissipative attributes of this methodology align closely with those of the conventional DG method, which is considered a viable strategy for achieving a harmonious balance between dissipation and precision. Especially for the compressible iLES, Roe's approximate Riemann solver (Moura, Sherwin & Peiró Reference Moura, Sherwin and Peiró2015) was combined with the DG discretization to deal with the unfiltered Navier–Stokes equations without a sub-grid scale model. This methodology, credited for its commendable numerical dissipation characteristics and numerical stability, has garnered widespread adoption in the spectral/hp method for under-resolved turbulence simulations (Moura et al. Reference Moura, Aman, Peiró and Sherwin2020; Mengaldo et al. Reference Mengaldo, Moxey, Turner, Moura, Jassim, Taylor, Peiro and Sherwin2021). In addition, in high-order compressible simulations, numerical aliasing can induce numerical instabilities, particularly in flow regimes dominated by convection. These instabilities might manifest as oscillations in the numerical solution, which do not genuinely represent physical phenomena but are artefacts of the numerical approach. Generally, there are two kinds of nonlinearities: partial differential equation (PDE) nonlinearities, associated with nonlinear and quasi-linear fluxes, and geometrical nonlinearities, linked to deformed or curved meshes. Given that Nektar++ addresses nonlinearities of PDE and geometric aliasing separately during the projection and solution of the equations, the selection of quadrature points should be predicated on the maximum order of the nonlinearities (Mengaldo et al. Reference Mengaldo, De Grazia, Moxey, Vincent and Sherwin2015).

3.3. Numerical set-up

The primary geometric parameters of the computational domain align with those detailed in § 1.1. Based on the experimental results from figure 2(b), two configurations that exhibit the most pronounced acoustic disturbances have been chosen. One of these configurations corresponds to ${\rm \pi} = L/{D_m} = 1$, and it is referred to as a ‘close arrangement’. The other configuration, with ${\rm \pi} = 8$, is termed a ‘half-wavelength arrangement’. For the close arrangement, the mainstream Reynolds number is $R{e_{in}} = 1.03 \times {10^5}$, and this corresponds to an inlet Mach number of $M{a_{in}} = 0.12$; and for the half-wavelength arrangement, $R{e_{in}} = 1.17 \times {10^5}$ and $M{a_{in}} = 0.13$. As for the computational domain, the total streamwise length is approximately ${L_{tot}} = 30{D_c}$, with an inlet distance approximately ${L_{in}} = 8{D_c}$ away from the upstream cavity; and the spanwise size is extended to ${D_z} = 0.8{D_m}$. Before conducting the three-dimensional (3-D) simulation, a series of 2-D simulations were conducted to see the independence of the h-type mesh and sensitivity of the p-type order. The entire computational domain is meshed with structured hexahedral and quadrilateral grids. The specific grid configurations and polynomial set-ups are listed in table 1. Additionally, the nodal basis is evaluated for a number of Gauss–Lobatto–Legendre quadrature points, which should be greater than N_p. This implementation facilitates polynomial de-aliasing of the nonlinear terms. Regarding the boundary conditions, the bottom surface and all the wall surfaces of the tandem deep cavities were assigned a no-slip condition. Simultaneously, non-reflective boundary conditions using a Riemann-invariant treatment were set at the inlet and outlet to eliminate the influence of reflected acoustic waves on the simulations. The outflow boundary is located at a relatively distant downstream position from the cavity, ensuring ample distance for the vortex dissipation. In incompressible iLESs, SVV was implemented to deal with small-scale turbulence structures, a second-order-accurate stiffly stable implicit–explicit (IMEX) scheme was used for time marching method and the CG method was used for spatial discretization. Besides, the advection term is explicitly integrated, whereas the viscous term is implicitly integrated. In compressible iLESs, the DG method was implemented as the spatial discretization method, along with Roe's Riemann solver to compute the convective flux $\boldsymbol{f}_n^c$, and the symmetric interior penalty method (Hartmann & Houston Reference Hartmann and Houston2006) for the computation of viscous flux $\boldsymbol{f}_n^v$. An implicit time advancing method is carried out using the second-order, singly diagonally implicit Runge–Kutta (DIRK) method. This allows us to take large time steps, and hence the time step is chosen only based on physical time resolution considerations and not on numerical stability, the resultant Courant–Friedrichs–Lewy (CFL) number was around 1.

Table 1. Flow conditions for numerical simulations using high-order spectral/hp element method. Here, DOF denotes degree of freedom; 2-/3-D denotes two-/three-dimensional; Com./Incom. denotes compressible/incompressible simulation; DIRK denotes diagonally implicit Runge–Kutta method and IMEX stable implicit–explicit (IMEX) scheme for time marching; DG and CG denote discontinuous and continuous Galerkin method, respectively.

5.1. Numerical set-up

The decoupled aeroacoustic simulation is classified into three major steps, the first obtains the steady-state or time-averaged flow field from flow simulation, the second maps the flow information from the fluid mesh to the acoustic mesh and the third numerically solves the compressible LNSEs. Here, the shear stress transport (SST) turbulence model is used to determine the steady-state turbulent flow within the tandem deep cavities, as much less numerical diffusion and better prediction of the eddy viscosity can be achieved. Note that the time-averaged flow fields from the previous LES results were not used as a base flow, as they have been already excited by the acoustic resonance oscillations and superimposed with the turbulence fluctuations. The SST-simulated mean flow quantities containing the turbulent viscosity $({\mu _\tau })$ were interpolated from the flow mesh onto the acoustic mesh during the mapping step. The LNSEs around the mapped base flow were numerically solved with a generalized minimal residual solver (Saad & Schultz Reference Saad and Schultz1986) and a Galerkin/least-square finite element (FE) stabilization scheme. A gradient term stabilization was introduced to suppress the Kelvin–Helmholtz instabilities by the reactive terms of the LNSEs. The numerical set-up is depicted in figure 16(a), where two configurations, one with a close arrangement and the other with a half-wavelength arrangement, were established. The boundary conditions for their flow simulations were kept consistent with the previous LESs. To discretize the fluid domain and acoustic domain, two sets of meshes were generated. Triangular elements and boundary layer elements were adopted within the computational domains. For flow meshes, the thickness of the first cell was determined such that y + value was almost close to 1, and approximately 0.2 and 0.5 million elements were finally generated within the two configurations. For the acoustic mesh, approximately10 nodes per acoustic wavelength were set to describe the acoustic wave propagation without numerical dissipation. The frequency range of the incident acoustic waves was set to 200–3500 Hz with a step of $\Delta f = 10\ \textrm{Hz}$, yielding a maximum element size of 0.01 mm with additional mesh refinement around the cavity entrances. The number of acoustic meshes reached 0.06 million for one configuration and 0.08 million for the other configuration. Note that the incident acoustic waves were planar waves, as their frequencies were below the cutoff frequency of the main duct. During acoustic simulations, the acoustic source regions were added at the two cavity endplates and harmonic pressure forcing was set as a source term to LNSEs. This was implemented by writing an expression into the LNSEs to define the direction, amplitude, frequency, phase and wavenumber information. For the close arrangement, a phase difference of $180\mathrm{^\circ }$ was applied between the two endplate acoustic sources, therefore out-of-phase acoustic forcing was created and acted on the concurrent shear layers; for the half-wavelength arrangement, in-phase acoustic forcing was created at the two cavity endplates. Besides, regions with a perfectly matched layer (PML) were also set at the inlet, outlet and cavity endplates to dissipate the incoming acoustic waves and prevent acoustic wave reflections. Here, PML imposes a complex-valued coordinate transformation (frequency domain) or inverse Laplace transformation (time domain) to satisfy the Sommerfeld radiation condition that can effectively maintain its wave impedance and eliminate wave reflections (Johnson Reference Johnson2021). An adiabatic condition was applied to the wall surface to make the temperature perturbation at the walls consistent with the pressure perturbations using the Neumann condition. To save computational cost, the acoustic simulations were conducted in the frequency domain rather than the time domain.

Figure 16. Numerical set-up (a) for the decoupled aeroacoustic simulation using linearized Navier–Stokes equations with out-of-phase (b) and in-phase (c) external acoustic forcing imposed on the closely arranged and half-wavelength-arranged tandem deep cavities, respectively. The enlarged panels highlight the acoustic scattering around the cavity entrances.

Figures 16(b) and 16(c) illustrate the successfully applied acoustic forcing fields to the tandem deep cavities, their frequencies and intensities were kept in consistency with the abovementioned experimental results. After the incident acoustic waves were generated at the cavity endplates, they propagate across the turbulent shear layers, enter the main duct and radiate outside through the PMLs. Four snapshots with different phases were utilized for demonstration; out-of-phase acoustic forcing was observed for the closely arranged configuration and in-phase acoustic forcing was for the half-wavelength configuration, confirming that the present LNSE model can mimic the internal acoustic forcing excited by the acoustic eigenmodes. Furthermore, upon closer examination of the isoline patterns within the highlighted panels, we observed a sequence of acoustic pulsations surrounding the entrances of the cavities and the downstream wall surfaces. These observations suggest the presence of acoustic scattering phenomena. These locations were nearly coincident with the vortices within the shear layers, providing confirmation that the current decoupled aeroacoustic simulations successfully captured the interactions between the stable turbulent shear layers and the external acoustic forcing.

To further confirm the present results, two measures, i.e. coherent kinetic energy $({\hat{G}_k})$ and coherent vorticity convection $({\hat{G}_\omega })$, are used to delineate the linear response of the turbulent shear layers to the external acoustic forcing by referring to Boujo et al. (Reference Boujo, Bauerheim and Noiray2018)

(5.1a,b)\begin{equation}{\hat{G}_k} = \frac{1}{{u_{ac}^2}}\left( {\iint_{{\Omega_c}} {({{\tilde{u}}^2} + {{\tilde{v}}^2})\,\textrm{d}\kern0.06em x\,\textrm{d}y} } \right)\!,\quad {\hat{G}_\omega } = \frac{1}{{{u_{ac}}}}\left( {\iint_{{\Omega_c}} {({{\tilde{\omega }}_z}\bar{U})\,\textrm{d}\kern0.06em x\,\textrm{d}y} } \right)\!,\end{equation}

where ${\varOmega _c}$ is the integrated domain spanning the streamwise extension of each cavity, ${\tilde{\omega }_z} = {\partial _x}\tilde{v} - {\partial _y}\tilde{u}$ is the spanwise vorticity of the coherent response. The first measure is useful to evaluate the kinetic energy of the acoustically excited fluid perturbation, the second is associated with the Coriolis force and represents the acoustic power production. Globally, both measures exhibit substantial amplification followed by dissipation as the frequency of the applied acoustic forcing increases, as shown in figure 17. For the closely arranged configuration, we observe two distinct peaks near 550 Hz and 1100 Hz, whereas for the half-wavelength-arranged configuration, two peaks are identified in close proximity to 620 Hz and 1270 Hz. Obviously, the first peaks can be related to the resonated acoustic eigenmodes but are overpredicted, and the second peaks are harmonics of the first peaks. The overprediction can be attributed to the influence of convective flow effects. Typically, turbulence interaction effects tend to shift resonance towards lower frequencies, while convective flow effects have the opposite effect, pushing it towards higher frequencies (Meliga Reference Meliga2017). This phenomenon also accounts for the frequency difference observed between the two tandem cavities. The spatial patterns of coherent perturbations within the shear layers will be examined in the upcoming section, focusing on these two dominant frequencies as indicated by the LNSE results.

Figure 17. Linear response of turbulent shear layers to the harmonic acoustic forcing: (a) coherent kinetic energy by out-of-phase acoustic forcing in closely arranged tandem deep cavities; (b) coherent kinetic energy by in-phase acoustic forcing in half-wavelength-arranged tandem deep cavities; (c) coherent vorticity convection by out-of-phase acoustic forcing; and (d) coherent vorticity convection by in-phase acoustic forcing. The frequency range of incident acoustic wave varied between 100 and 3000 Hz with a frequency interval of 10 Hz. The broad frequency bands from the compressible iLES simulations are also annotated for comparison.

5.2. Coherent perturbations of turbulent shear layers

The linear response of the steady-state turbulent shear layers to harmonic acoustic forcing is characterized in terms of the phase-dependent evolution of streamwise velocity perturbation and the superimposed streamlines, as shown in figure 18. Noting that the phase variation gives rise to different acoustic pressure gradients acting on the shear layers. Globally, the response is localized within the shear layer and downstream boundary layer, which is reasonable with the large flow gradient and turbulent viscosity here. Distinctive wavepacket structures originate from the cavity entrances, and they exhibit a streamwise expansion along the main duct while simultaneously contracting vertically as they approach the cavity's end. Simultaneously, weak velocity perturbations are generated at the upstream cavity corner and amplified when transported towards the downstream corner by the streamwise convecting shear layers; this process is a typical response of shear layer instability.

Figure 18. Phase-dependent evolution of coherent streamwise velocity perturbation $(\tilde{u})$ and coherent energy density $({G_k})$ by (a,c) out-of-phase acoustic forcing of closely arranged tandem deep cavities, and (b,d) in-phase acoustic forcing of half-wavelength-arranged tandem deep cavities.

The present study can be discerned through a comparative analysis between the two tandem cavities in different arrangements. In the case of close arrangement, application of out-of-phase acoustic forcing resulted in opposing yet notably intensified velocity perturbations at the downstream cavity. Additionally, distinctive vortex structures with varying rotation directions were observed. This phenomenon supported the iLES finding with a more expanded shear layer within the downstream cavity. Such amplification can be confirmed in figure 18(c), in which the coherent energy density was calculated and plotted along the cavity entrances, almost one level higher energy density can be found along the downstream cavity than the upstream cavity. As for the half-wavelength arrangement, the applied in-phase acoustic forcing yielded similar but slightly intensified velocity perturbations at the downstream cavity, and the rotation directions of vortex structures were also same. Here, the discrepancy between the LNSE result and iLES result is attributed to the downstream acoustic intensity; as in iLES, the resonance intensity is quite attenuated compared with the upstream cavity. When coming to the phase-dependent evolution process of the coherent vortex structures, the number and spatial scale of the coherent vortex structures within the cavity entrances are vary periodically.

Under the influence of acoustic forcing, a large-scale vortex structure initially occupies the cavity entrance but evolves into two smaller-scale vortex structures in the subsequent phase. Subsequently, they reconfigure into a single large-scale vortex structure, characterized by streamwise displacement and an opposite rotation direction. This process can induce a streamwise growth of energy density from the upstream corner but with fluctuations when approaching the downstream corner. Besides, alternating variations of the coherent vortex structures and energy density can be identified between the two cavities in close arrangement, but synchronous variation of the half-wavelength-arranged configuration, consolidating the coherent perturbations of the turbulent shear layers.

Figures 19 and 20 illustrate the phase-dependent evolution of the synthetic streamwise velocity, synthetic spanwise vorticity and shear layer momentum thickness in response to the external acoustic forcing acting on the tandem deep cavities. These three quantities took the time-averaged flow fields into consideration, and can demonstrate the dynamic responses of mean flow to the periodic acoustic forcing. In figure 19(a) for the tandem deep cavities in close arrangement, snapshots of streamwise velocity, spanwise vorticity and superimposed velocity streamlines depict the formation, convection and collapse processes of the coherent vortex structures within the two cavity entrances. A phase difference between the two cavities demonstrates the occurrence of alternating motion. This can be further explored by the variations of shear layer momentum thickness, as plotted in figure 20(a). Periodic oscillations of the shear layer can be observed by tracing the phase-dependent locations of the maximum thickness, indicating a periodic flapping motion of the turbulent concurrent shear layers. During the first half-acoustic period, the maximum thickness of the upstream shear layer gradually increased and the corresponding locations gradually transported and approached the downstream wall surfaces; during the second half-period, after the shear layer impinging on the downstream wall surface, the maximum thickness gradually decreased and their locations gradually approached the upstream trailing edge. When looking into the downstream cavity, the shear layer momentum thickness was significantly amplified with larger values. However, the evolution process demonstrates a phase lag with the upstream shear layer. As the shear layer has already impinged on the trailing edge within the first half of an acoustic period; while in the second half, the locations of maximum thickness appeared almost at the centre of the cavity entrance. Such phase-differentiated shear layer perturbations by the out-of-phase acoustic forcing are in consistency with the alternating resonance behaviour.

Figure 19. Phase-dependent evolution of coherent shear layer perturbations inside tandem deep cavities, presented by synthetic streamwise velocity $(\tilde{u} + \bar{U})$ and synthetic spanwise vorticity $({\tilde{\omega }_z} + {\bar{\varOmega }_z})$: (a) with out-of-phase acoustic forcing of the close arrangement; (b) with in-phase acoustic forcing of the half-wavelength arrangement.

Figure 20. Phase-dependent evolution of coherent shear layer perturbations, represented by momentum thickness: (a) closely arranged and (b) half-wavelength-arranged tandem deep cavities.

In figure 20(b), phase-similar shear layer perturbations by the in-phase acoustic forcing are found in the half-wavelength arrangement, which are consistent with the abovementioned synchronous resonance behaviour of the turbulent shear layers. Even the momentum thickness demonstrates a slight intensification at the downstream cavity; the locations of the maximum thicknesses are close in each cavity at all of the four phases. Two kinds of flow structures are also found from LNSE results, i.e. the accompanying recirculation zone below the shear layer vortex structures and the separation bubble along the downstream main duct wall surface. The recirculation zone is formed by the turbulent entrainment effect and its oscillation is driven by the shear layer vortices. Regarding the separation bubble, it originated when the downstream shear layer vortex structures impinged upon the downstream corner and subsequently amplified, convecting along the main duct wall surface. Similar flow structures, characterized by recirculation zones and separation bubbles, can also be observed in Wang et al. (Reference Wang, Deng, Mao, He and Liu2020) and Ho & Kim (Reference Ho and Kim2021). This separation bubble was also attributed to the divergence of Reynolds stress, which can lead to an increase in the thickness of the shear layer (Mantič-Lugo & Gallaire Reference Mantič-Lugo and Gallaire2016). Now, focusing on the half-wavelength arrangement, as shown in figure 19(b), the flow structures resulting from in-phase acoustic forcing bear a resemblance to the previously mentioned findings. These include coherent vortex structures along the cavity entrances accompanied by recirculation zones inside the cavity volumes and separation bubbles along the downstream main duct wall surface. Nonetheless, two notable differences can be observed. First, there is a synchronous evolution process between the two tandem cavities, which is consistent with the iLES results. Second, there is the appearance of an additional separation bubble originating from the upstream cavity corner. This is logical as these two cavities are spaced further apart.