2.1. Flow configuration

The flow consists of the interaction between an incident shock wave and a laminar boundary layer developing on a flat plate. A 2-D sketch of the flow (2-D slice) is shown in figure 2. The presence of the impinging shock wave (also called the incident shock wave, noted as $(S_I)$ ) imposes a sharp adverse pressure gradient to the boundary layer of thickness $\delta$ ( $99 \,\%$ thickness) evaluated just upstream of the interaction region. As accurately described by Délery & Dussauge (Reference Délery and Dussauge2009), in the configuration presently studied, the adverse pressure gradient is strong enough to initiate a separation of the boundary layer that reattaches further downstream forming a closed separation bubble. The mean separation and reattachment points are respectively noted as $\overline {x}_S$ and $\overline {x}_R$ . The flow being subsonic in the wall vicinity, the pressure rise due to the incident shock wave is sensed upstream to the location, noted as $x_{sh}$ , where the incident shock wave would impact the wall (if the flow were non-viscous) because of slow acoustic waves, explaining the location of the separation point, $\overline {x}_S$ , upstream of the location of the incident shock wave impingement $x_{sh}$ . At the separation point, the deviation of the supersonic flow due to the separation leads to the formation of the so-called reflected shock wave, noted as $(S_R)$ . In fact, the incident shock is reflected at the apex of the separation bubble as expansion waves. At the reattachment point, the deviation of the supersonic flow due to the presence of the wall leads to compression waves that can also coalesce to form the so-called reattachment shock wave. For given Reynolds and Mach numbers, the extent of the separation bubble is driven by the intensity of the incident shock wave, let us say the pressure ratio from each side of the shock, prescribed by the shock angle $\alpha$ . In the following, we characterise this extent by the interaction length $L_{int}=x_{sh}-\overline {x}_S$ .

Figure 2. Sketch of a 2-D slice of the flow, where $(S_I)$ is the incident shock wave of angle $\alpha$ , $\delta$ is the boundary layer thickness just upstream of the interaction zone, $(S_R)$ is the reflected shock wave, $\overline {x}_S$ is the mean separation point, $\overline {x}_R$ is the mean reattachment point, $x_{sh}$ is the location of the impingement shock wave (if the flow were non-viscous) and $L_{int}=x_{sh}-\overline {x}_S$ is the interaction length.

The free-stream conditions are similar to the experimental and numerical test case documented in Degrez, Boccadoro & Wendt (Reference Degrez, Boccadoro and Wendt1987), with a free-stream Mach number $M=2.15$ and a free-stream Reynolds number $Re=1600$ , based on the boundary layer thickness just upstream of the interaction zone $\delta =1.6 \times 10^{-2}\,{\textrm {m}}$ . Compared with the configuration of Degrez et al. (Reference Degrez, Boccadoro and Wendt1987), where the whole flow is steady (fully laminar interaction), the incident shock wave angle $\alpha$ is increased up to $33.8^\circ$ . The interaction is thus strengthened, the separation bubble is enlarged and the dynamics of the separation bubble becomes unsteady with transition to turbulence triggered in the bubble (Ben Hassan Saidi Reference Ben Hassan Saidi2019). Table 1 summarises the physical parameters of the simulation.

Table 1. Flow parameters of the SWBLI.

2.2. Governing equations

To perform the DNS, we consider the dimensionless compressible Navier–Stokes equations expressed in a Cartesian coordinate system. The length and time scales for non-dimensionalisation are respectively $x_{sh}$ and $U_\infty / x_{sh}$ , with $U_\infty$ being the free-stream velocity. All the primitive variables are non-dimensionalised by the corresponding free-stream quantities except the thermodynamic pressure $P^*$ and the total energy per unit volume $\rho ^* E^*$ that are non-dimensionalised by the free-stream dynamic pressure $\rho _\infty U_\infty ^2$ . The dimensionless Navier–Stokes system is therefore written as

(2.1) \begin{equation} \frac {\partial \boldsymbol {U}}{\partial t} + \nabla \cdot \boldsymbol {F}(\boldsymbol {U}) - \nabla \cdot \boldsymbol {F_v}(\boldsymbol {U}, \nabla \boldsymbol {U}) = 0, \end{equation}

where $\boldsymbol {U}$ is the dimensionless vector of conservative variables, $\boldsymbol {F(U)}$ the dimensionless convective fluxes and $\boldsymbol {F_v(U, \nabla (\boldsymbol {U}))}$ the dimensionless diffusive fluxes that write respectively as

(2.2) \begin{equation} \begin{array}{lll} \boldsymbol {U}=\left [ \begin{array}{c} \rho \\ \rho \boldsymbol {u}\\ \rho E \end{array} \right ], \quad & \boldsymbol {F}=\left [ \begin{array}{c} \rho \boldsymbol {u}\\ \rho \boldsymbol {u} \otimes \boldsymbol {u} + P \mathbb {I}\\ (\rho E + P) \boldsymbol {u} \end{array} \right ] \quad \text {and} \quad \quad & \boldsymbol {F_v}=\left [ \begin{array}{c} 0\\ { \frac {1}{Re} \mathbf {\sigma }} \\ { \frac {1}{Re} \boldsymbol {u}. \mathbf {\sigma } + \frac {\mu }{ (\gamma - 1) Re \ Pr M^2} \nabla T } \end{array} \right ] . \end{array} \end{equation}

Here $\,\mathbb {I}$ stands for the identity matrix, $\rho$ is the dimensionless density, $\boldsymbol {u}=[u,v,w]^T$ is the dimensionless velocity vector, $E$ is the dimensionless specific total energy, $P$ is the dimensionless thermodynamic pressure and $\lambda$ is the dimensionless thermal conductivity; $P$ relates to the dimensionless conservative variables through the relationship

(2.3) \begin{equation} P = (\gamma - 1) \left ( \rho E - \frac {1}{2} \frac {(\rho \boldsymbol {u}) \cdot (\rho \boldsymbol {u})}{\rho } \right ), \end{equation}

and to the dimensionless static temperature $T$ through the following ideal gas dimensionless equation of state:

(2.4) \begin{equation} T = \gamma M^2 \frac {P}{\rho }. \end{equation}

The dimensionless viscous stress tensor is expressed as

(2.5) \begin{equation} \displaystyle {\mathbf {\sigma }= \mu ( \nabla \boldsymbol {u} + \nabla ^T \boldsymbol {u} ) - \frac {2}{3} \mu \left (\nabla \cdot \boldsymbol {u}\right )\boldsymbol {I} }. \end{equation}

We assume that the dynamic viscosity only depends on the temperature through Sutherland’s law. Given the Prandtl number, the dimensionless thermal conductivity is then deduced from the dimensionless dynamic viscosity $ \lambda = \mu Cp/\ Pr$ .

2.3. Numerical methods

The DNS of this flow is performed using an in-house finite-volume based DNS and/or large-eddy simulation solver parallelized using the MPI (message passing interface) library. The numerical methods employed are presented in detail in Ben Hassan Saïdi et al. (Reference Ben Hassan Saïdi, Fournier and Tenaud2020), where the ability of this code to compute high Reynolds compressible (turbulent and shocked) flows has already been demonstrated (including the boundary conditions presented in § 2.4). An one step (OS) procedure is employed that splits the resolution into the convective part and the viscous part of the problem. The convective fluxes are discretised by a high-order one step (OS) finite-volume scheme based on the Lax–Wendroff approach that is a seventh-order accurate coupled time and space approximation. It is coupled to an original monotonicity-preserving (MP) shock-capturing procedure to form the OSMP7 scheme (first introduced in Daru & Tenaud Reference Daru and Tenaud2004). A second-order centred finite difference scheme is used for the spatial discretisation of the diffusive fluxes. The temporal discretisation of the diffusive fluxes is obtained by means of a second-order Runge–Kutta time integration. In Ben Hassan Saïdi et al. (Reference Ben Hassan Saïdi, Fournier and Tenaud2020) we showed that the use of a higher-order scheme for computing diffusive fluxes has a negligible effect on the accuracy of the SWBLI calculation and is therefore unnecessary.

The simulation is performed with a time step of $\Delta t \simeq 1.19 \times 10^{-4} \times L_{int} / U_\infty$ corresponding to a maximum Courant–Friedrichs–Lewy number of $0.5$ in the computational domain.

2.4. Computational domain and boundary conditions

The geometry of the 3-D domain is $\mathcal {D}=[0;250\delta ]\times [0;125\delta ]\times [0;62.5\delta ]$ . The domain is discretised using a Cartesian mesh with non-uniform spacing in the direction normal to the wall ( $z$ ). The mesh employs $800 \times 400 \times 202$ grid points in $(x\times y \times z)$ . In the normal-to-the-wall direction, the mesh is tightened close to the wall using a hyperbolic tangent law to obtain a minimum grid spacing over the plate of $\Delta z_{min}=0.0125 \delta$ . A grid convergence study (see Appendix B) has been performed in order to assess that the employed mesh is fine enough to compute precisely the flow in the interaction region.

A uniform flow is prescribed as the inlet of the domain as prescribed by the guidelines of the 4th HiOCFD workshop (Wang et al. Reference Wang2013; HiOCFD 2016) to compute the SWBLI in the configuration of Degrez et al. (Reference Degrez, Boccadoro and Wendt1987), which has been validated in Ben Hassan Saïdi et al. (Reference Ben Hassan Saïdi, Fournier and Tenaud2020) against reference experimental and numerical results. The incident shock wave is created by imposing the Rankine–Hugoniot relationships in the inlet plane at a height $z$ chosen so that the shock wave impinges the flat plate at the desired abscissa $x_{sh}=62.5 \delta$ . No-slip and adiabatic wall conditions are prescribed at the flat plate location $(z = 0)$ . Outlet time-dependent non-reflecting boundary conditions (Thompson Reference Thompson1987) are imposed at the top surface and at the downstream outlet boundary of the computational domain. Periodic conditions are used in the spanwise direction ( $y$ ).

The numerical strategy has already been validated for this type of simulation. Indeed, the same numerical schemes and boundary conditions have already been successfully employed to compute SWBLIs in the flat plate configuration. Please refer to Ben Hassan Saïdi et al. (Reference Ben Hassan Saïdi, Fournier and Tenaud2020) for more details.

2.5. Database

After the initial transient, corresponding to the initial propagation of the incident shock wave in the domain and the creation of the separation bubble and subsequent shock wave system, $3 \times 10^6$ iterations have been simulated with a time step of $\Delta t \simeq 1.19 \times 10^{-4} \times L_{int} / U_\infty$ . A snapshot has been recorded every $1000$ time steps so that a total of 3000 snapshots are available for analysis. Each snapshot consists in the state vector of primitive variables $\boldsymbol {q}(\boldsymbol {x}_i,t_l)=[\rho ,u,v,w,T]^T$ evaluated at each mesh point $\boldsymbol {x}_i$ (with $ 1 \leq i \leq M$ , $M$ being the number of mesh points considered) in the 3-D domain containing the interaction region and at each time of snapshot recording $t_l$ (with $1 \leq l \leq N_t$ , $N_t$ being the number of snapshots recorded). All data are gathered in the snapshot matrix

(2.6) \begin{equation} {\boldsymbol{Q}}=[\boldsymbol {q}_1\,\ldots \,\boldsymbol {q}_{N_t}]=\begin{pmatrix} \boldsymbol {q}(\boldsymbol {x}_1,t_1) & \cdots & \boldsymbol {q}(\boldsymbol {x}_1,t_{N_t}) \\ \vdots & & \vdots \\ \boldsymbol {q}(\boldsymbol {x}_{M},t_1) & \cdots & \boldsymbol {q}(\boldsymbol {x}_{M},t_{N_t}) \end{pmatrix}. \end{equation}

All the results presented in the rest of this paper are based on the exploitation of the data contained in this snapshot matrix.

All the power spectral densities (PSDs) of signals presented in § 2.6 have been computed using the periodogram method implemented in the scipy.signal package of python (Virtanen et al. Reference Virtanen2020) using Hann windowing. For all time signals studied in this section, the sampling characteristics are given in table 2. The total length of the signal considered is around 11 events at $St=0.02$ and 22 events at $St=0.04$ .

Table 2. Sampling characteristics of time signals.

For both SPOD and BSMD analysis algorithms introduced in §§ 3 and 4, the Welch’s method (Welch Reference Welch1967) is used to compute the discrete Fourier transform (DFT) of the state vector. The matrix of snapshots is split into $N_b= 10$ with an overlapping of $67 \,\%$ between the blocks. With this slicing, the total length of the signal considered for the analysis is around 30 events at $St=0.02$ and 60 events at $St=0.04$ .

In the following, for both SPOD and BSMD analyses, we consider the centred snapshot matrix

(2.7) \begin{equation} {\boldsymbol{Q}\,}^c=[\boldsymbol {q}^c_1\,\ldots \,\boldsymbol {q}^c_{N_t}]=\begin{pmatrix} \boldsymbol {q}^c(\boldsymbol {x}_1,t_1) & \cdots & \boldsymbol {q}^c(\boldsymbol {x}_1,t_{N_t}) \\ \vdots & & \vdots \\ \boldsymbol {q}^c(\boldsymbol {x}_{M},t_1) & \cdots & \boldsymbol {q}^c(\boldsymbol {x}_{M},t_{N_t}) \end{pmatrix}, \end{equation}

with $\boldsymbol {q}^c(\boldsymbol {x}_i,t_l)=\boldsymbol {q}(\boldsymbol {x}_i,t_l)-\boldsymbol {\overline {q}}(\boldsymbol {x}_i)$ , where $\boldsymbol {\overline {q}}(\boldsymbol {x}_i)$ is the vector of primitive variable averages in time and spanwise direction.

Note that, for the signal length, for the following SPOD and BSMD analyses, a convergence study has been performed. Indeed, these analyses have been performed with different numbers of snapshots (i.e. a different length of the signal and frequency resolution) leading to the same conclusions to those presented below.

2.6. Mean flow organisation and frequency content

The time mean longitudinal velocity averaged in the spanwise direction $\overline {u}$ is plotted in figure 3 in the interaction zone. The velocity is scaled by the free-stream velocity. The red line shows the isocontour $\overline {u} = 0$ . It allows us to clearly identify the separation bubble in which the mean longitudinal velocity is negative. The separation bubble also corresponds to the region in which the mean skin friction coefficient plotted in figure 4 is negative. The green line in figure 3 represents the $z/\delta =2.5$ plane in which we analyse the characteristic sizes of the flow structures in the transverse direction in figure 6(a). The green dots, referenced by numbers 1–4, mark the position of the probes placed in the mixing layer to analyse its characteristic frequencies by calculating the PSD of the velocity signal shown in figure 7.

As the incident boundary layer is perfectly laminar, it is particularly prone to separation when subjected to the adverse pressure gradient imposed by the incident shock wave, compared with turbulent boundary layers. The resulting mean separation length is $L=\overline {x}_R-\overline {x}_S\simeq 74.68 \delta$ , with $\overline {x}_S\simeq 13.44 \delta$ and $\overline {x}_R\simeq 88.12 \delta$ respectively the time mean separation and reattachment points averaged in the spanwise direction. The interaction length is $L_{int}=x_{sh}-\overline {x}_S=49.06 \delta$ . The time mean height of the bubble, averaged in the spanwise direction, is $\overline {h}=3.5 \delta$ . Moreover, the mean separation bubble is quasi-symmetric with respect to the vertical axis passing through the apex of the bubble. This length and aspect of the mean bubble is consistent with comparable simulations (Robinet Reference Robinet2007; Song & Hao Reference Song and Hao2023). For similar Mach and Reynolds numbers, the length of the separation bubbles for SWBLIs with laminar incident boundary layers is around one order of magnitude bigger than their counterpart for SWBLIs involving turbulent incident boundary layers; see, for example, Adler & Gaitonde (Reference Adler and Gaitonde2018) where $\overline {L}=4.0 \delta$ or Ben Hassan Saïdi et al. (Reference Ben Hassan Saïdi, Fournier and Tenaud2020) where $\overline {L}=4.92 \delta$ . For turbulent SWBLIs, the symmetry of the separation bubble is also broken, the second part of the bubble being less long than the first part, which is due to the increased mixing rate in the shear layer. The same symmetry breaking (less marked) is observed for SWBLIs with a sufficiently strongly forced laminar incident boundary layer (Sansica Reference Sansica2015; Diop et al. Reference Diop, Piponniau and Dupont2019; Mauriello et al. Reference Mauriello, Larcheveque and Dupont2022).

The difference of topology of the mean separation bubble of the present flow configuration with respect to turbulent SWBLIs calls for a comment about the characteristic length used to build non-dimensionalised frequencies. Indeed, the main frequencies involved in the interactions have been introduced in the literature for turbulent interactions as Strouhal numbers based on the separation length $St_L=fL/U_e$ . For the present flow configuration, as the bubble is symmetric, the interaction length $L_{int}$ is more suitable to be consistent with the characteristic frequencies highlighted in turbulent interactions. In the following, the frequencies will therefore be expressed as Strouhal numbers based on the interaction length $St_{L_{int}}=fL_{int}/U_e$ .

Figure 3. Time mean longitudinal velocity field averaged in the spanwise direction $\overline {u}$ . The red line shows the isocontour $\overline {u} = 0$ . The green line is parallel to the wall at a height of $z/\delta = 2.5$ . Coordinates of probes in the $(x,z)$ basis, from 1–4: $(13.44 \delta ,0.3\delta )$ , $(43.75,3.5\delta )$ , $(68.75 \delta , 2.5 \delta )$ , $(92.50 \delta ,0.3 \delta )$ .

Figure 4. Distribution of the spanwise averaged skin friction along the flat plate (). Blasius laminar boundary layer solution (). Turbulent correlation from Cousteix (Reference Cousteix1989) ().

Figure 5. Vortical structures highlighted by isosurfaces of the discriminant criterion coloured by the magnitude of the longitudinal component of the velocity. Shock waves are highlighted by isosurfaces of $\lVert \nabla P \rVert$ . The upper left insert shows a zoomed view of the vortical structures around the separation bubble. The upper right insert shows an overhead view of the vortical structures in the separation bubble and in the downstream boundary layer. On the upper right insert, the vertical black lines indicate, from left to right, the mean separation line (), the mean line of incident shock impingement on the apex of the separation bubble (), the mean line of shock impingement on the flat plate () and the mean reattachment line ().

A snapshot of the flow is shown in figure 5. The discriminant criterion, introduced in Chong et al. (Reference Chong, Soria, Perry, Chacin, Cantwell and Na1998), already used by Pirozzoli & Grasso (Reference Pirozzoli and Grasso2006) in the context of SWBLI simulations, is used to identify the vortex structures present in the flow. Shock waves are highlighted by isosurfaces of $\lVert \nabla P \rVert$ coloured in black. The upstream boundary layer not being artificially forced by modes resulting from a stability problem (Sansica et al. Reference Sansica, Sandham and Hu2014, Reference Sansica, Sandham and Hu2016) or a free-stream turbulence method, the only fluctuations supplied to the flow come from shock waves (leading edge, separated or the incident shock waves) that disrupt the flow mainly two-dimensionally, although the natural most unstable normal modes are 3-D. Figure 5 shows clearly this characteristic. In what we seek to show, the existence of a nonlinear interaction between medium-frequency modes generating low-frequency modes, this has no consequence whether the convective modes exited are 2-D or 3- D. Large spanwise vortices, corresponding to Kelvin-Helmholtz rolls that progressively develop, are clearly visible in the shear layer edging the separation bubble between the reflected and the incident shocks. After the incident shock impingement, the shear layer is populated by elongated structures in the streamwise direction. Inside the separation bubble, 3-D structures are visible in the downstream part of the separation zone. For the particular value of the discriminant criterion shown in figure 5, no coherent structure is visible in the early part of the interaction. The dynamical activity inside the separation bubble therefore seems to be mainly concentrated in its second part. The dynamic activity of the separation bubble causes the rapid transition of the boundary layer after the reattachment point. The turbulent nature of the boundary layer downstream of the interaction zone is confirmed by the downstream distribution of the mean skin friction coefficient, as visible in figure 4.

Figure 6. Spanwise length scales in the interaction zone. (a) Time averaged spanwise wavenumber of the longitudinal component of velocity $u$ in a plane parallel to the flat plate at a height of $z/\delta = 2.5$ . The vertical dashed lines indicate the $\beta =0.25$ and $\beta =2$ wavenumbers (). The horizontal dashed lines indicate, from bottom to top, the separation point (), the raising shear layer (), the descending shear layer () and the reattachment point (). (b) Isocontours of the longitudinal component of velocity $u=-0.041$ coloured by $u$ . Slice ( $y$ plane) of the density $\rho$ .

In order to physically characterise the longitudinal structures identified in figure 5, the spectral content of the velocity signal in the transverse direction is studied. To this end, the time averaged one-dimensional spatial Fourier transform in the spanwise direction of the longitudinal component of the velocity ( $u$ ) was calculated in different planes parallel to the wall, defined by $z=cste$ ( $z/\delta =0.3,0.5,1,1.5,2,2.5,3,3.5,4,4.5$ ). Figure 6(a) shows the evolution of these spectra as a function of the abscissa along the flat plate in the $z/\delta = 2.5$ plane marked by the green line in figure 3. The wavenumber is scaled using the same scaling as in Bugeat et al. (Reference Bugeat, Robinet, Chassaing and Sagaut2022) for comparison. The abscissa is normalised by using the mean interaction length $L_{int}$ and the mean separation point $\overline {x}_S$ . In this figure, the cyan and green horizontal dashed lines indicate the intersection between the plane under consideration and the shear layer bounding the separation bubble. The line at $(x-\overline {x}_S)/L_{int}=0.445$ indicates the shear layer upstream of the impact of the incident shock, which we call the ‘rising shear layer’. The line at $(x-\overline {x}_S)/L_{int}=1.2$ indicates the shear layer downstream of the impact of the incident shock, which we call the ‘descending shear layer’. The vertical dashed lines indicate the $\beta =0.25$ and $\beta =2$ wavenumbers that were highlighted as the two maxima of the $G(\beta )$ curve, where $G$ is the optimal gain of a low-frequency resolvent analysis, and that were respectively associated with Görtler-type structures and streaks in Bugeat et al. (Reference Bugeat, Robinet, Chassaing and Sagaut2022). Just upstream of the descending shear layer (i.e. inside the bubble) we can clearly see a peak at $\beta =0.25$ . In the descending shear layer, we see a plateau containing the value $\beta =2$ . These same peaks are found in the spectra of the other planes (not shown here), indicating the presence of Görtler-type structures inside the bubble and streaks in the descending shear layer. This location of Görtler-type structures and streaks is in agreement with the results of Bugeat et al. (Reference Bugeat, Robinet, Chassaing and Sagaut2022). The optimal modes of low-frequency resolvent analysis are therefore found in the DNS data of the present unforced transitional SWBLI. Analysis of the instantaneous 3-D fields highlights these longitudinal structures and illustrates their location. For example, figure 6(b) shows the isocontour of the longitudinal velocity field $u=-0.041$ next to a 2-D slice coloured by $\rho$ to locate the separation bubble. The wavelength of the Görtler-type structures and the streaks clearly manifest themselves in the isocontour of velocity in the zones identified by the spectral analysis.

We now want to document the frequency content of the dynamics of the interaction zone. The dynamic structures mainly visible in figure 5 seem to be concentrated in the shear layer. We therefore begin by examining the frequency content of the mixing layer. To do this, we examine the PSD of the longitudinal velocity signal measured at different points in the shear layer. These points are shown in a 2-D plane in figure 3, where their coordinates in the $(x,z)$ plane are specified. In practice, for these coordinates in the $(x,z)$ plane, a probe was placed at each mesh point in the $y$ direction to measure the velocity signal and compute the PSD. The PSDs were then averaged over span.

We can see that in the zone of the separation point (probe 1), the energy level of the oscillations in the shear layer is very low and localised in the low- and medium-frequency ranges characterising the breathing and flapping phenomena. There are clear peaks at frequencies $St_{L_{int}}=0.04,0.0667$ and $0.1$ . This frequency content is amplified in the rising mixing layer. Indeed, the signal measured by probe 2 has approximately the same frequency content as probe 1, but at a higher energy level.

We saw above that the passage through the incident shock is associated with a distortion of the vortex structures, which are mainly oriented in the spanwise direction upstream of the shock and preferentially in the streamwise direction downstream of the shock. These structures are rapidly amplified downstream of the incident shock. Just downstream of the incident shock (probe 3), the energy of the fluctuations is amplified compared with the level observed at probe 2. The peaks measured at low and medium frequencies are still present, but the spectrum is filling out and these frequencies are less predominant. In addition, higher frequencies appear. This signal is massively amplified in the descending shear layer. In the attachment zone (probe 4), the energy of the signal is almost quadrupled compared with probe 3. This amplification is also selective. The frequencies already highlighted, $St_{L_{int}}=0.04,0.0667$ and $0.1$ , are particularly peaked in the reattachment zone. There is also a further filling in of the spectrum at high frequencies concomitant with the generation of fine structures associated with the transition to turbulence in this region.

These measurements in the shear layer seem to indicate that the mixing layer is forced from the separation point by the characteristic SWBLI frequencies at low and medium frequencies. The rising shear layer amplifies these frequencies by a factor of around 10. However, the energy levels in the rising shear layer are low relative to those in the descending shear layer. The latter shows a massive amplification of fluctuations in the low- and medium-frequency ranges, and in particular, of the frequencies characteristic of the SWBLI unsteadiness $St_{L_{int}}=0.04,0.0667$ and $0.1$ . These observations are consistent with various previous works, in particular, those of Sansica et al. (Reference Sansica, Sandham and Hu2014) and Mauriello et al. (Reference Mauriello, Larcheveque and Dupont2022) who have shown that the low-frequency content of the interaction zone is created by the dynamics of the shear layer and seems to result from the nonlinear evolution of the latter in the downstream part of the recirculation bubble. These authors also showed that disturbances at these frequencies then travel upstream inside the recirculation bubble.

Figure 7. Power spectral densities of the spanwise averaged longitudinal velocity $u$ signal recorded on probes shown in figure 3. (a) Probe 1. (b) Probe 2. (c) Probe 3. (d) Probe 4.

In order to study the frequency content of this information feedback in the recirculation bubble, we plot in figure 8(a) the evolution of the spanwise averaged premultiplied PSD of the longitudinal velocity signal ( $u$ ) measured in a plane parallel to the flat plate, in very close proximity to the wall $(z/\delta =0.3)$ . For each abscissa along the flat plate, the PSD at this abscissa is normalised and premultiplied. This makes it possible to determine the predominant frequencies for each abscissa, but the relative energy differences between abscissas are erased by this representation. (The relative difference in energy between abscissa is however shown in figure 8(b), which is commented on below.) In this figure, the horizontal dashed lines indicate, from bottom to top, the positions of the separation point, the rising shear layer, the impact of the incident shock on the flat plate and the reattachment. We can see that the spectral content in the attachment zone is present throughout the bubble in the near wall, but with a progressive narrowing of the significant frequency range towards the low- and medium-frequency range as we move upstream in the bubble. Thus, the disturbances reaching the foot of the rising shear layer have an exclusively low- and medium-frequency spectral content with a cutoff frequency of approximately $St_{L_{int}}\simeq 0.1$ . We can therefore see that the rise of the disturbance inside the bubble selects the low and medium frequencies. This result had already been shown in Mauriello et al. (Reference Mauriello, Larcheveque and Dupont2022). In order to determine whether this feedback into the bubble is damped, we integrate the PSD between the $St_{L_{int}}=0.02$ and $0.1$ for each abscissa along the flat plate to obtain the streamwise evolution of the expected power in this frequency range $E(x)=\int _{St=0.02}^{St=0.1} \text{PSD}(u(x)) {\mathrm d}St$ . This power distribution is shown in figure 8(b). This figure clearly shows a large peak of power in the reattachment region, in agreement with the analysis of the spectra shown in figure 7. We can clearly see an exponential decrease in power inside the bubble. The upwelling of disturbances is therefore damped inside the bubble. The signal reaches the detachment zone with low but finite energy levels. The shear layer is therefore forced into the separation zone by a signal of finite energy and quasi-continuous spectrum in the frequency range $St_{L_{int}}\in [0.02,0.1]$ . However, figure 8(a) shows a high selectivity of the mixing layer. Indeed, the quasi-continuous spectrum becomes strongly peaked in the shear layer at the frequencies already highlighted in figure 7(a), i.e. $St_{L_{int}}=0.04,0.0667$ and $0.1$ .

Figure 8. Power spectral density and power of the longitudinal velocity signal close to the wall inside the separation bubble. (a) Distribution along the flat plate of the spanwise averaged premultiplied and normalised PSD of the longitudinal component of velocity $u$ measured in a plane parallel to the flat plate at a height of $z/\delta = 0.3$ . The vertical dashed lines indicate significant frequency peaks in the low- and medium-frequency range (). The horizontal dashed lines indicate, from bottom to top, the reflected shock foot (), the separation point (), the crossing of the rising shear layer (), the incident shock impingement location () and the reattachment point (). (b) Distribution along the flat plate of the expected power in the range $St_{L_{int}}\in [0.01,0.1]$ of the longitudinal component of velocity $u$ measured in a plane parallel to the flat plate at a height of $z/\delta = 0.3$ . The vertical dashed lines indicate, from left to right, the separation point (), the crossing of the rising shear layer (), the incident shock impingement location (), the crossing of the descending shear layer () and the reattachment point ().

We are also interested in the pressure forces exerted on the wall in the interaction zone. To this end, we analyse the evolution in the streamwise direction of the premultiplied PSD normalised to the wall pressure shown in figure 9(a). A similar evolution of the spectrum can be observed between the points of separation and reattachment in the sense that we move from the low/medium-frequency range to a wider range, extended to higher frequencies. However, we can clearly see that in the first part of the separation bubble, the wall pressure spectrum predominantly highlights the spectral signature of the shear layer, which is peaked at frequencies $St_{L_{int}}=0.04,0.0667$ and $0.1$ . In both respects, the wall pressure spectrum is very similar to the spectra presented in Mauriello et al. (Reference Mauriello, Larcheveque and Dupont2022). In this simulation, we therefore record pressure oscillations at the foot of the reflected shock dominated by the low- and medium-frequency range $St_{L_{int}}\in [0.02,0.1]$ that characterises the SWBLI unsteadiness. The evolution in the streamwise direction of the power of the parietal pressure signal in this frequency range is shown in figure 9(b). As with the near-wall velocity signal, there is a quasi-exponential attenuation of the wall pressure perturbations as we move upwards from the reattachment point. We then see an increase in the first quarter of the bubble as we get closer to the separation point and the shear layer moves closer to the wall. The signal power in the separation zone is an order of magnitude lower than the power in the reattachment zone.

Figure 9. Spanwise averaged PSD and power of the wall pressure along the flat plate. (a) Distribution along the flat plate of the spanwise averaged premultiplied and normalised PSD of the wall pressure. The vertical dashed lines indicate significant frequency peaks in the low- and medium-frequency range. The horizontal dashed lines indicate, from bottom to top, the reflected shock foot (), the separation point (), the incident shock wave impingement location () and the reattachment point (). (b) Distribution along the flat plate of the expected power in the range $St_{L_{int}}\in [0.02,0.1]$ of the wall pressure. The vertical dashed lines indicate, from left to right, the separation point (), the incident shock impingement location () and the reattachment point ().

Figure 10. History of the spanwise averaged abscissa along the flat plate of (a) the separation point $x_S$ and (b) the reattachment point $x_R$ .

The complex dynamics described above are also reflected in the breathing and flapping cycles of enlargement and shrinkage of the separation bubble, which will be described in § 3 documenting the SPOD analysis we carried out. This low- and medium-frequency dynamics of the bubble can also be seen by analysing the temporal evolution of the position of the separation and reattachment points, respectively, $x_S$ and $x_R$ , shown in figure 10. It can be seen that the reattachment point is subject to oscillations whose maximum extrusions correspond to about $7\delta$ while the separation point is subject to very weak oscillations with maximum amplitudes of the order of $0.4 \delta$ . This order-of-magnitude difference in the amplitudes of the oscillations of $x_S$ and $x_R$ is consistent with the order-of-magnitude difference observed in the power of the oscillations of the parietal pressure and velocity signals in these two regions. Similarly, in agreement with the wall pressure and velocity signals in these regions of the flow, the spectral content of $x_S(t)$ oscillations is mainly located in the low- and medium-frequency range with significant peaks at frequencies $St_{L_{int}}=0.04,0.0667$ and $0.1$ , while the oscillations of $x_R(t)$ present a wider band spectrum with peaks also localised at these frequencies, as shown in figure 11.

Figure 11. Spanwise averaged PSDs of (a) the separation $x_S$ and (b) the reattachment point locations $x_R$ .

In the present flow configuration, the low energy level of perturbations coming from the reattachment point reaching the separation point explains the low amplitude of the separation point and associated reflected shock foot (not shown here) oscillations. This point is a salient feature of the flow studied here, which makes it different from configurations with a forced laminar or fully turbulent boundary layer. As already stated in the introduction of this paper, in this respect, the present case study represents a limit case in which the mechanisms at the origin of the low-frequency dynamics appear to be present at low intensity and are not disturbed by the natural frequencies of the incident turbulence that, we believe, makes it a good candidate for identifying the intrinsic underlying mechanisms of this dynamics. In the cases of forced or turbulent incoming boundary layers, the more intense dynamics of the mixing layer (due to upstream forcing) has two consequences: (i) the energy of the perturbations at reattachment, and therefore, the energy of the signal emitted from reattachment is much higher; (ii) the separation length, i.e. the distance between the zone of signal emission and the separation point, is smaller. The combination of these two characteristics means that the disturbances received at the separation point are more energetic than in the case studied in this paper. In the forced and turbulent cases, disturbances are then more likely (for a given strength of interaction) sufficiently energetic to induce higher amplitude oscillations at the separation point.

3.1. Results

The SPOD spectrum of eigenvalues $\lambda _m$ is shown in figure 12. We can see that the SPOD spectrum describes a marked dynamic in the low- and medium-frequency range, with a clear drop-off beyond that. We can clearly see that the energy decreases with the number of the mode, so that the first mode contains about twice as much energy as the third. In the low- and medium-frequency range, almost $60 \,\%$ of the flow energy is contained in the first two modes. Interestingly, the three main frequencies for mode 1 are the frequencies found in local DNS measurements (see figure 7): $St_{L_{int}}=0.04,0.0667$ and $0.1$ . In particular, we have already shown that these three frequencies, which are present throughout the shear layer, are largely dominant in the rising shear layer. We can also see that mode 2 signs the frequency $St_{L_{int}}=0.04$ and its first subharmonic $St_{L_{int}}=0.02$ as well as the sum of the two ( $St_{L_{int}}=0.06$ ) and another frequency at $St_{L_{int}}=0.073$ . The three frequencies highlighted by the second mode (i.e. $St_{L_{int}}\simeq 0.02,0.06$ and $0.073$ ) are well present in the spectra of the DNS signals measured locally in the shear layer shown in figure 7. However, unlike the frequencies highlighted by the first mode, these frequencies reach energy levels comparable to the former in the descending shear layer, downstream of the impact of the incident shock, and are comparatively less energetic in the upstream part of the separation bubble. The first mode therefore seems to bring out the dominant frequencies in the whole of the interaction zone, while the second mode brings out frequencies that are mainly present in the downstream part of the interaction zone.

Figure 12. The SPOD spectrum of the state vector $\boldsymbol {\hat {q}}$ for each mode ( $m_j,\,\,j=1,\ldots ,10$ ) expressed as a relative energy given in percent.

We are now interested in the shape of the dominant modes in the entire bubble. As an illustration, figure 13 shows the modulus of the first mode at frequency $St=0.04$ (characteristic frequency of the separation bubble breathing) for the three velocity components $u$ , $v$ and $w$ . For the longitudinal component $u$ , we can see that the activity of the mode is mainly localised in the shear layer. In particular, large amplitudes of the mode are observed in the rising shear layer. In the descending shear layer, the amplitude increases the closer we get to the abscissa of the reattachment point. In the attachment zone, strong amplitudes are also observed close to the wall and amplify downstream of the interaction. We also observe activity in the longitudinal velocity mode localised on the reflected shock as well as the reattachment shock, materialising the oscillations of these shocks at this low frequency, characteristic of SWBLI unsteadiness. These oscillations of the reflected and reattachment shocks are even more visible in the shape of the vertical velocity mode $w$ , whose activity is localised in these zones as well as in the shear layer, with less intensification than for the longitudinal velocity in the reattachment zone. The activity of the mode for the spanwise velocity component $v$ is completely localised inside the separation bubble.

Figure 13. First spatial SPOD mode $(m_1)$ averaged in the spanwise direction at $S_t=0.04$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

The spatial modes at frequencies $St_{L_{int}}=0.0667$ and $0.0934$ are qualitatively almost identical to those at frequency $St_{L_{int}}=0.04$ . For the sake of brevity, they are shown in Appendix A, in figures 29 and 32 respectively for $St_{L_{int}}=0.0667$ and $St_{L_{int}}=0.0934$ .

In figure 14 we plot the real part of the first spatio-temporal mode $\boldsymbol {\unicode{x1D6F7}}_1(\boldsymbol {x}_i,t)$ of the three velocity components for the Strouhal number $St = 0.04$ , characteristic of the breathing of the recirculating bubble. We plot regularly spaced snapshots in a period $T=1/St$ . For each speed component, a supplementary movie of the animation of the spatio-temporal evolution of the associated mode is available at https://doi.org/.

Figure 14. Visualisation of $\boldsymbol {\unicode{x1D6F7}}_1(\boldsymbol {x},t)$ velocity $u$ at $St = 0.04$ , averaged in the spanwise direction. (a) Velocity component $u$ . (b) Velocity component $v$ . (c) Velocity component $w$ .

Figure 15. Representation of $\boldsymbol {\unicode{x1D6F7}}_1(x,z=\text {cste},t)$ for $z/\delta = 1.48$ . (a) Velocity component $u$ . (b) Velocity component $v$ . The vertical black lines indicate, from left to right, the separation point, the crossing of the rising shear layer, the crossing of the incident shock wave, the crossing of the descending shear layer and the reattachment point.

In the longitudinal velocity mode $u$ , we can clearly see that the oscillations of the reflected shock in the longitudinal direction are in phase with the oscillations of the shear layer in this direction, and therefore, in phase with the oscillations of the separation point. We can also see that the breathing cycles follow the following sequence: the value of the mode changes in the downstream part of the bubble before being propagated upstream in the first part of the bubble, eventually contaminating the entire shear layer. The upstream propagation of information in the bubble is shown in figure 15, where we plot versus time (for a single breathing period) the value of $\boldsymbol {\unicode{x1D6F7}}^u_1$ along the flat plate measured at a constant height above it $z/\delta =1.48$ . In the first part of the bubble, we clearly distinguish regions of the same sign oriented with a negative slope in the $(x,t)$ frame, denoting an upstream propagation of this sign. This behaviour is qualitatively observed at all altitudes in the separation bubble. The upstream propagation from the latter part of the separation bubble is also visible in the spanwise velocity mode as shown in figures 14(b) and 15(b). This result seems to confirm that the downstream part of the separation plays a decisive role in the breathing dynamics of the separation bubble. More precisely, it confirms the analysis of § 2.6: the low-frequency oscillations of the separation bubble are driven by oscillations at this frequency in its latter part (materialised by the change in sign of the mode in the reattachment zone) propagated upstream in the separation bubble.

A similar qualitative behaviour is observed in the results obtained for the frequencies $St_{L_{int}} = 0.0667$ and $St_{L_{int}}=0.0934$ , the latter being characteristic of the flapping of the separation bubble. The figures spatial evolution of the first mode at these frequencies can be found in Appendix A, in figures 30, 31 and 33, 34 respectively for $St_{L_{int}} = 0.0667$ and $St_{L_{int}}=0.0934$ . The similarity of the spatio-temporal modes for breathing and flapping may suggest a link to be determined between these two oscillation modes.

4.1. The BSMD methodology

The BSMD algorithm is a modal decomposition that reveals the presence of triadic nonlinear interactions from multidimensional data. To this end, an integral measure of the pointwise bispectral density between frequency components of the vector of primitive variables is introduced as

(4.2) \begin{equation} b(f_k,f_l)=E\left [ \int _{\Omega }^{} \hat {\boldsymbol {q}}_k^* \circ \hat {\boldsymbol {q}}_l^* \circ \hat {\boldsymbol {q}}_{k+l} \, {\mathrm d}\boldsymbol {x} \right]=E \left[\hat {\boldsymbol {q}}_{k \circ l}^H \boldsymbol {\mathcal {W}} \hat {\boldsymbol {q}}_{k + l} \right] = E\left [ \langle \hat {\boldsymbol {q}}_{(k \circ l)} , \hat {\boldsymbol {q}}_{(k + l)} \rangle \right ], \end{equation}

where $\hat {\boldsymbol {q}}_{(k \circ l)}=\hat {\boldsymbol {q}}(\boldsymbol {x},f_k) \circ \hat {\boldsymbol {q}}(\boldsymbol {x},f_l)$ (with $\circ$ is the Hadamard product), $\hat {\boldsymbol {q}}_{(k + l)} = \hat {\boldsymbol {q}}(\boldsymbol {x},f_{k+l})$ and $(.)^H$ is the conjugate transpose; $\boldsymbol {\mathcal {W}}$ is the diagonal matrix of spatial quadrature weights and $\Omega$ is the spatial domain over which the flow is defined.

Using the Welch’s slicing of the snapshot matrix, two modal projections are constructed: the cross-frequency fields

(4.3) \begin{equation} \boldsymbol {\Phi }^{[i]}_{k \circ l} (\boldsymbol {x},f_k, f_l) = \sum _{j=1}^{N_b} a_{i}^{[j]}(f_{k+l}) \hat {\boldsymbol {q}}_{(k \circ l)}^{[j]} \end{equation}

and the bispectral modes

(4.4) \begin{equation} \boldsymbol {\Phi }^{[i]}_{k + l} (\boldsymbol {x},f_{k+l}) = \sum _{j=1}^{N_b} a_{i}^{[j]}(f_{k+l}) \hat {\boldsymbol {q}}_{(k + l)}^{[j]}, \end{equation}

which share a common set of expansion coefficients $a_{i}^{[j]}$ with $ 1 \leqslant j \leqslant N_b$ .

The goal of BSMD is to compute modes that optimally represent the data in terms of the integral bispectral density defined in (4.2). That is, we seek the set of expansion coefficients $a_1$ that maximises the absolute value of $b(f_k,f_l)$ as defined in (4.2) with $\hat {\boldsymbol {q}}_{(k \circ l)}$ and $\hat {\boldsymbol {q}}_{(k + l)}$ replaced by their modal expansion. We therefore target the set of expansion coefficients that verifies

(4.5) \begin{equation} \boldsymbol {a}_1 = \text {arg} \max _{\parallel \boldsymbol {a}_1 \parallel = 1} \mid E \left[ \boldsymbol {\Phi }^{[1]}_{k \circ l} \boldsymbol {\mathcal {W}} \boldsymbol {\Phi }^{[i]}_{k + l} \right] \mid , \end{equation}

where the coefficient vector is required to be a unit vector in order to guarantee boundedness of the expansion.

More details about the BSMD methodology can be found in Schmidt (Reference Schmidt2020), especially about the construction and solving of this optimisation problem.

As a result of the BSMD analysis, we obtain optimal cross-frequency fields and bispectral modes, respectively $\mathbf {\Phi }^{[1]}_{k \circ l}$ and $\mathbf {\Phi }^{[1]}_{k + l}$ , as well as the mode bispectrum

(4.6) \begin{equation} \lambda _1(f_k,f_l) = E\left[\mathbf {\Phi }^{[1]H}_{k \circ l} \mathbf {W} \mathbf {\Phi }^{[1]}_{k + l} \right] \in \mathbb {C}. \end{equation}

The bispectrum mode is the fundamental result of the BSMD analysis. Indeed, significant triadic interactions are indicated by local maxima of the modulus of this quantity, also called bispectrum mode amplitude. Bispectral modes $\mathbf {\Phi }^{[1]}_{k + l}$ are linear combinations of Fourier modes and can be interpreted as observable physical structures resulting from quadratic interactions between frequencies $f_k$ and $f_l$ . The multiplicative cross-frequency fields $\mathbf {\Phi }^{[1]}_{k \circ l}$ , on the contrary, are maps of phase alignment between two frequency components that may not directly be observed. Once significant triadic interactions are identified in the map of bispectrum mode amplitude, interaction maps can be computed as

(4.7) \begin{equation} \varPsi _{k,l} (\mathbf {x},f_k,f_l)= | \mathbf {\Phi }^{[1]}_{k \circ l} \circ \mathbf {\Phi }^{[1]}_{k+l} |. \end{equation}

This field indicates the location of the triadic interaction involving the triad $(f_k,f_l,f_{k+l})$ in the flow field, as it quantifies the average local bicorrelation between the frequencies in the domain.

To sum up, the BSMD analysis process involves calculating the map of mode bispectrum amplitude $\lVert \lambda _1(f_k,f_l) \rVert$ to identify significant quadratic interactions. These are the interactions that contribute most to the flow dynamics. Once these interactions have been identified, the associated interaction map $\varPsi _{k,l} (\mathbf {x},f_k,f_l)$ is studied. This allows us to identify the regions of the flow where this triadic interaction occurs. Finally, the associated bispectral mode $\mathbf {\Phi }^{[1]}_{k+l}$ , which is the physical mode resulting from this interaction, is visualised and interpreted.

In the following we present the BSMD analysis of the database introduced in § 2.5. The interpretation of the results will be performed following the process that has just been outlined above, focusing on the low/medium-frequency dynamics of the interaction. In particular, we will be interested in identifying and documenting nonlinear links between medium frequencies and low frequencies.

4.2. Mode bispectrum amplitude

The modulus of the mode bispectrum $\lVert \lambda _1(St_k,St_l) \rVert$ obtained when applying BSMD to the snapshot matrix (2.7) is shown in figure 16. Only the region of non-redundant information, as demonstrated by Schmidt (Reference Schmidt2020), is shown. Moreover, the interactions between under-resolved frequencies for which two or less periods are resolved through the Welch Fourier transform estimation have been greyed out on the graph. The region above the horizontal axis of the graph represents sum interactions $St_1+St_2=St_3$ , the region below the horizontal axis represents difference interactions $St_1-St_2=St_3$ . The interactions located on the horizontal axis correspond to interactions for which $St_2=0$ resulting in $St_3=St_1$ and, therefore, correspond to the linear evolution of the flow around its mean. In our analysis, we do not consider local maxima located on the line $St_1=-St_2$ as they result in modes of frequency $St_3=0$ corresponding to the mean flow distortion.

We can clearly identify strong quadratic couplings (highlighted by the high value of $\lVert \lambda _1(St_k,St_l) \rVert$ ) for $(St_1,St_2) \in [0,0.20]^2$ , namely in the low- and medium-frequency ranges. The intensity of interactions decreases sharply at higher frequencies. In the low- and medium-frequency range, the following three types of interesting couplings (highlighted by the high value of $\lVert \lambda _1(St_k,St_l) \rVert$ ) can be identified, which will be analysed in more details in the following paragraphs.

Figure 16. Modulus of the complex mode bispectrum of the centred state vector $\boldsymbol {q}^c$ in the $(St_1 , St_2 )$ plane. Relevant types of interactions are circled or indicated by arrows: interactions between medium frequencies creating low frequencies (arrow); interactions between low frequencies creating medium frequencies (dashed arrow); interactions reflecting the linear evolution of the flow around its mean. Interactions between unresolved frequencies are greyed.

1. (i) The interactions in the difference interaction region correspond to a cascade of interactions generating low frequencies from medium frequencies. This cascade is indicated by the solid black arrow.

2. (ii) The interactions in the sum interaction region correspond to a cascade of sum interactions generating medium frequencies from low frequencies. This cascade is indicated by the dashed arrow.

3. (iii) The interactions circled correspond to interactions associated with the linear evolution of the flow around the mean flow that contribute significantly to the dynamics of the flow.

In the following, in order to further interpret the physical role of the triadic interactions that we have highlighted, for each type of interaction, we analyse in more detail the involved frequencies, as well as their characteristic interaction maps and the shape of the associated bispectral modes.

4.3. Cascade of interactions generating low frequencies from medium frequencies (type I)

In order to get more insight in the low- and medium-frequency range, the map of mode bispectrum amplitude is reproduced in figure 17 by zooming in these ranges.

In the region of the map below the horizontal axis (difference interaction region), there is a large continuous zone of interaction of significant amplitude defined by $(St_1,St_2)\in ]0,0.18]\times ]0,-0.18]$ . These interactions are difference interactions, each interaction therefore results in the creation of a lower frequency $St_3$ from two interacting frequencies $St_1$ and $St_2$ . As this zone is continuous, a multitude of cascades of successive interactions exist, creating low frequencies from interactions in the mid-frequency range characteristic of the well-documented flapping of the shear layer $St \in [0.1,0.18]$ .

As an illustration, an example of such a cascade is shown in figure 18. Each interaction is given a number (between brackets) that is reported on the mode bispectrum amplitude map (figure 17). We would stress that the cascade considered here is just one example of a multitude of other possible paths. The first column (interaction of numbers (1), (2), (3) and (4)) contains interactions between medium frequencies present in the flow (typical frequencies of the well-documented flapping mode) that are written in red. The presence of these medium frequencies in the flow is sufficient to initiate a cascade of interactions that eventually creates low frequencies. Indeed the interactions of the first column generate lower frequencies. In column two, we report interactions of these generated lower frequencies with medium frequencies (interaction of numbers (5), (6) and (7)). These interactions generate low frequencies typical of the breathing mode, written in orange. We show in the third column that interactions between frequencies generated from interactions of the first and second columns also contribute to the generation of low frequencies characteristic of the breathing (interaction of numbers (8), (9), (10) and (11)).

Figure 17. Modulus of the complex mode bispectrum of the centred state vector $\boldsymbol {q}^c$ in the $(St_1 , St_2 )$ plane: zoom on medium/low frequencies. Sequences of interactions of each type are highlighted and numbered. Interactions between unresolved frequencies are greyed.

Figure 18. Example of triadic interactions sequence between medium frequencies eventually creating low frequencies. The frequencies are expressed as $St_{int}$ . Medium frequencies at the origin of the sequence are written in red. The resulting low frequencies are written in orange. Each interaction has a number between parentheses reported in figure 17.

This example of cascade shows how the presence in the flow of the well-documented medium frequencies characteristic of the flapping generates low frequencies through a multitude of cascades of difference interactions.

The interaction map of an emblematic triadic interaction of type I is shown in figure 19 for the three components of velocity. Other interactions of type I are not shown here for the sake of brevity, but they exhibit the same qualitative behaviour. We clearly identify that the nonlinear interactions constituting the cascade creating energy at low frequencies from medium frequencies are mainly localised in the second part of the separation bubble and especially in the region of reattachment. This point is consistent with the analyses and interpretations exposed in previous sections, namely the direct analysis of the DNS database (§ 2.5) and the SPOD analysis (§ 3), in which the low-frequency breathing dynamics of the separation bubble has been shown to be driven by dynamical activity taking place in the downstream part of the separation bubble. Moreover, this result confirms and refines the results of Sansica et al. (Reference Sansica, Sandham and Hu2014) and Mauriello et al. (Reference Mauriello, Larcheveque and Dupont2022). Indeed, both works suspected a role of nonlinear interactions in generating low frequencies in the downstream part of the interaction zone, although the precise determination of the location of significant triadic interactions was difficult due to the local nature of the methods employed.

Figure 19. Spanwise averaged interaction maps for a triadic interaction (interaction (7) in figures 17 and 18) illustrating the sequence of interactions creating low frequencies from medium frequencies (interactions of type I. Frequencies involved: $St_1=0.12$ , $St_2=-0.08$ , $St_3=0.04$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

The low-frequency bispectral mode resulting from the triadic interactions of type I evoked in figure 19 is shown in figure 20 for the three components of velocity. For all components of velocity, we clearly see that the bispectral mode is very similar to the first SPOD mode at the same frequency ( $St=0.04$ , characteristic of the breathing) shown in figure 13. This result strongly suggests that the low-frequency breathing dynamics is fuelled by interactions between medium-frequency modes characteristic of flapping.

Figure 20. Spanwise averaged bispectral mode for a triadic interaction (interaction (7) in figures 17 and 18) illustrating the sequence of interactions creating low frequencies from medium frequencies (interactions of type I. Frequencies involved: $St_1=0.12$ , $St_2=-0.08$ , $St_3=0.04$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

4.4. Cascade of interactions creating medium-frequency modes from low frequencies (type II)

By analysing figure 17, we identify a region of high amplitudes of the mode bispectrum above the horizontal axis (sum interactions), defined by $(St_1,St_2)\in [0,0.15]\times ]0,0.1]$ . This region is symmetrical to the region of type I interactions analysed in the previous paragraph. It has an analog effect on flow dynamics. However, as it is populated by sum interactions, it consists of a multitude of cascades of interactions generating medium frequencies from interactions between low frequencies instead of the opposite.

In the same way as for interactions of type I, we show in figure 21 a particular path for a cascade of such interactions. Each interaction is given a number (between brackets) that is reported on the mode bispectrum amplitude map (figure 17). To form this cascade, we only assume the existence of two low frequencies marked in red, i.e. here $St=0.02$ and $St=0.0267$ . The first stage of the cascade consists of the self-interactions of these frequencies and the interactions of these frequencies with each other, creating higher frequencies. We then unfold a cascade of interactions between the frequencies created and show how medium frequencies (marked in orange) can be created in one or two stages. Again, as the region of high amplitude sum interactions is continuous, a multitude of cascades exist that has the same effect on the flow: creating medium frequencies, characteristic of the flapping, from low frequencies, characteristic of the breathing of the separation bubble.

Figure 21. Example of a cascade of triadic interactions creating medium frequencies from low frequencies. The frequencies are expressed as $St_{int}$ . Low frequencies at the origin of the cascade are written in red. The resulting medium frequencies are written in orange. Each interaction has a number between parentheses reported in figure 17.

Figure 22. Spanwise averaged interaction map for a triadic interaction (interaction (18) in figures 17 and 21) illustrating the cascade of interactions creating medium-frequency modes from low frequencies (interactions of type II. Frequencies involved: $St_1=0.04667$ , $St_2=0.04667$ , $St_3=0.0933$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

The interaction map of an emblematic triadic interaction of type II is shown in figure 22 for the three components of velocity. Other interactions of type II are not shown here for the sake of brevity, but they exhibit the same qualitative behaviour. Again, these interaction maps clearly shows that the cascade of interactions creating medium-frequency modes from low frequencies mainly take place in the second part of the separation bubble and especially in the region of reattachment. This result is consistent with the SPOD analysis of the dominant dynamics at medium frequency (mode 1 at $St=0.0934$ ) that shows a driving role of the reattachment region’s dynamics in the flapping’s dynamics. The direct analysis of the DNS database also revealed a strong intensification of the medium-frequency dynamics in the second part of the separation bubble.

The medium-frequency bispectral mode resulting from the triadic interactions of type II evoked in figure 22 is shown in figure 23 for the three components of velocity. The medium-frequency bispectral modes resulting from the cascade of interactions of type II are very similar to the medium-frequency flapping modes highlighted by the SPOD analysis and shown in Appendix A. This result suggests a feedback loop between the medium and low frequencies. Indeed, the analysis of the triadic interactions of type I (undertaken in § 4.3) suggest that low-frequency modes (breathing) are fuelled by triadic interactions between medium frequencies characteristic of flapping. In turn, the current analysis of interactions of type II suggest that the medium-frequency flapping dynamics is fuelled by a cascade of triadic interactions between low-frequency modes characteristic of the breathing of the separation bubble.

Figure 23. Spanwise averaged bispectral mode for a triadic interaction (interaction (18) in figures 17 and 21) illustrating the cascade of interactions creating medium-frequency modes from low frequencies (interactions of type II. Frequencies involved: $St_1=0.04667$ , $St_2=0.04667$ , $St_3=0.0933$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

4.5. Modes arising from the linear evolution of the flow (type III)

We clearly see a range of high amplitudes along the horizontal axis of the map of the mode bispectrum amplitude. These points, circled in figure 16, correspond to interactions for which $St_1 \in [0.02,0.18]$ and $St_2=0$ . Hence, $St_3=St_1$ and the corresponding bispectral modes can be interpreted as the linear evolution of the flow around its mean, according to Schmidt (Reference Schmidt2020). Interestingly, the range of frequencies of these modes corresponds to the range of low and medium frequencies characteristic of the SWBLI unsteadiness, including breathing dynamics and flapping dynamics. The flow therefore shows a strong linear tendency of the flow to oscillate around it’s mean field at medium and low frequencies characterising respectively the flapping and the breathing of the bubble. In particular, we observe high intensities at the frequencies identified by the direct analysis of the DNS data and the SPOD analysis (peak frequencies of the first mode), namely $St_1=0.04$ , $St_1=0.0667$ and $St_1=0.0934$ . These three modes are reported in figure 17 as numbers $(27)$ , $(28)$ and $(29)$ .

We show below the bispectral modes for the modes oscillating at $St_1=0.04$ and $St_1=0.0934$ in figures 24 and 26, respectively. These modes are qualitatively very similar to the low-frequency breathing and medium-frequency flapping modes obtained by SPOD analysis. It shows that these linear modes are contributing to the breathing and flapping modes of the separation bubble.

The interaction maps corresponding to these modes are shown in figures 25 and 27. The amplitude of these maps are mainly significant in the second part of the separation bubble and especially in the reattachment region. It clearly shows that these medium- and low-frequency modes highlighted are driven by the dynamical activity in the downstream part of the interaction, similarly to the modes resulting from interactions of type I and II.

Figure 24. Spanwise averaged bispectral mode for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1=0.04$ , $St_2=0$ , $St_3=0.04$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

Figure 25. Spanwise averaged interaction map for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1=0.04$ , $St_2=0$ , $St_3=0.04$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Longitudinal velocity $v$ . (c) Longitudinal velocity $w$ .

Figure 26. Spanwise averaged bispectral mode for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1=0.0934$ , $St_2=0$ , $St_3=0.0934$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Spanwise velocity $v$ . (c) Vertical velocity $w$ .

Figure 27. Spanwise averaged interaction map for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1=0.0934$ , $St_2=0$ , $St_3=0.0934$ . The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$ . (b) Longitudinal velocity $v$ . (c) Longitudinal velocity $w$ .