A direct numerical simulation of an oblique shock wave impinging on a turbulent boundary layer at Mach number 2.28 is carried out at moderate Reynolds number,

simulating flow conditions similar to those of the experiment by [Dupont et al., 2006].

The low-frequency shock unsteadiness, whose characteristics have been the focus of considerable
research efforts, is here investigated via the Morlet wavelet transform.
Owing to its compact support in both physical and Fourier spaces,
the wavelet transformation makes it possible to track the time evolution
of the various scales of the wall-pressure fluctuations.
This property also makes it possible
to define a local intermittency measure, representing a frequency-dependent flatness factor,
to pinpoint the bursts of energy that characterise the shock intermittency scale by scale.
As a major result, wavelet decomposition shows that the broadband shock movement is actually
the result of a collection of sparse events in time, each characterised by its own temporal scale.
This feature is hidden by the classical Fourier analysis, which can only show the time-averaged behaviour.

Then, we propose a procedure to process any relevant time series,
such as the time history of the wall-pressure or that of the
separation bubble extent,
in which we use a condition based on the local intermittency measure
to filter out the turbulent content
in the proximity of the shock foot and to isolate only the intermittent component of the signal.
In addition, wavelet analysis reveals the intermittent
behaviour also of the breathing motion of the recirculation bubble behind the reflected shock,
and allows us to detect a direct, partial correspondence between the most significant intermittent events
of the separation region and those of the wall-pressure at the foot of the shock.

MSC Codes   76N15, 75N99

§ INTRODUCTION

The interaction between shock waves and turbulent boundary layers represents a
major challenge for modern aerospace research, given its occurrence in a broad range of
engineering applications involving transonic, supersonic, and hypersonic systems.
Typical examples include high-speed intakes, over-expanded rocket nozzles,
transonic airfoil buffeting, and aerodynamics of high-speed vehicles or space launchers.
stbli must carefully be taken into account in the design process,
since they have the potential of harmfully impacting the performance of aerospace systems by
enhancing aerodynamic drag and heat transfer at the wall. The adverse pressure gradient imposed by the shock
can dramatically alter the structure of the boundary layer, inducing its thickening and, in the case of strong interactions,
large-scale separation. The latter is usually associated with concurrent complex phenomena as amplification of turbulence and noise,
formation of large-scale vortical structures, and unsteadiness involving a wide range of spatial and temporal scales.
Due to their relevance from the technological and scientific standpoint, stbli have been an active research field for more
than five decades, as demonstrated by the numerous review articles about the topic

[Green, 1970, Adamson & Messiter, 1980, Delery, 1985, Dolling, 2001, Gaitonde, 2015].

One of the most gripping aspects of stbli
is the low-frequency unsteadiness of the system,
which consists in the appearance of streamwise oscillations of the
separation shock characterised by a broadband motion at frequencies well below the typical
time scales associated with fine-grained turbulence in the
upstream boundary layer. It has been highlighted that
the shock motion exhibits unsteady features that are similar
among various stbli in different 2D and 3D canonical geometries,
such as compression ramps, reflected shocks, blunt and sharp fins.
Several experimental and numerical works have greatly contributed to shed
some light on the origin of this unsteadiness,
for which different mechanisms have been
proposed [Erengil & Dolling, 1991, Beresh et al., 2002, Humble et al., 2009],
mainly divided into two categories. On the one hand,
according to the upstream mechanisms, the source of the unsteadiness is
the presence of low-frequency flow structures in the incoming turbulent
boundary layer [Andreopoulos & Muck, 1987, Ganapathisubramani et al., 2009];
on the other hand, according to the
downstream mechanisms, the unsteadiness is related to the breathing
motion of the separation bubble behind the shock, which expands and contracts periodically [Touber & Sandham, 2009, Piponniau et al., 2009].
Although a general consensus has still not been reached,
it has been argued that both mechanisms always coexist
in all stbli: the downstream mechanism dominates for
strongly separated flows, whereas a combined mechanism
dominates for weakly separated flows [Clemens & Narayanaswamy, 2009].

In the last decade, several studies have exploited high-fidelity simulation data sets
to extract the relevant dynamical features of both compression ramps and oblique shock reflections
[Pirozzoli et al., 2010, Grilli et al., 2012, Nichols et al., 2017, Martelli et al., 2020, Bakulu et al., 2021].
For example, Priebe et al., 2016 applied dmd to the data obtained through
dns of a Mach 2.9 compression ramp stbli,

and found a strong similarity between dmd modes and
linear stability modes reported in the literature.
Furthermore, they revealed the existence of streamwise elongated regions of low and high momentum
extending from the shock foot in the downstream flow that the authors interpreted as reminiscent
of Görtler-like vortices, already observed in previous studies [Loginov et al., 2006].
Similar findings were reported by Pasquariello et al., 2017,
who performed a wall-resolved, long-time integrated les of
a high-Reynolds number stbli at Mach 3 with massive flow separation.
Their sparsity-promoting dmd identified two types of modes:
a low-frequency mode involving the main actors of the interaction (shock system,
separated shear-layer and separation bubble) associated with the
classical breathing motion of the recirculating flow, and medium-frequency modes
associated with the shear-layer vortices produced at the shock foot and convected downstream.

To better understand the unsteadiness of the stbli system,

time series extracted from the numerical or the experimental fluid domain
are usually analysed by means of the classical Fourier trasform.

This standard technique is typically used to
generate frequency spectra of the physical variables,
to determine the energy content of the signals at each frequency, and
has been applied extensively to characterise the motion of the reflected shock in stbli,
as reported in several works [Clemens & Narayanaswamy, 2014].
However, according to Bell & Lepicovsky, 1995,
such an approach is fundamentally justified only when stationary or periodic,
ergodic data are measured over a long time.
In fact, for highly unsteady and irregular time series,
Fourier transform may misrepresent localised features of the signal
in a time-averaged sense, which would distort the
representation of the actual physical phenomena taking place
[Huang et al., 1998].

The wavelet transform [Mallat, 1989, Daubechies, 1992] is instead an analysis
tool well-suited to the study of multi-scale and non-stationary processes. Indeed,
decomposing a time series into the time/temporal-scale space

makes it possible to detect the localised variations of the energy within the signal,
through a localised counterpart of the standard Fourier spectra.
As a result, wavelet transform allows the detection of intermittent or modulated features
in complex flows. On the basis of such properties, wavelet analysis
has provided valuable insight into fluid mechanics [Liandrat & Moret-Bailly, 1990],
as attested by its use for the study of turbulence [Farge, 1992, Farge et al., 1999],
aeroacoustics [Grizzi & Camussi, 2012, Camussi et al., 2017],
transonic buffet [Kouchi et al., 2016], and over-expanded nozzles [Martelli et al., 2017].
However, only very few works have used it to characterise the low-frequency unsteadiness
in stbli [Poggie & Smits, 1997, Guo et al., 2020],
and none of them performed such an analysis on a detailed dns database.

Moreover, given the fact that intermittency is defined by the presence of
localised bursts of high-energy activity,

wavelets can represent a suitable basis also to map and estimate the importance of rare but strong events:
a task that could be hardly fulfilled by the trigonometric functions used in the Fourier transform,
because of their infinite support in the time domain.

In this work, we present novel dns data of an oblique shock wave impinging on a turbulent boundary layer
inspired to the experiment performed at the IUSTI's hypo-turbulent supersonic wind tunnel
in Dupont et al., 2006, Dupont et al., 2008.
The flow we simulated is characterised by a Mach number equal to 2.28 and a
momentum-based Reynolds number $Re_\theta$ equal to 6900.

Several groups have previously studied the same flow case
either at the same conditions of the experiment, but using
les [Touber & Sandham, 2009, Morgan et al., 2013, Agostini et al., 2015],
or using dns, but at a reduced Reynolds number [Pirozzoli & Bernardini, 2011],
because of the hindering computational burden associated to the simulation.

Unlike stbli at low Reynolds numbers
[Dolling & Or, 1985, Adams, 2000, Pasquariello et al., 2014, Nichols et al., 2017],
in which the viscosity diffuses the separation shock into a milder compression fan,
the moderate Reynolds number considered allows the shock foot to penetrate
more deeply into the boundary layer. As a result, the inviscid pressure jump
of the oscillating shock leaves a direct footprint in the
wall-pressure, which exhibits a clear and definite intermittency.
For lower Reynolds number, instead, since the compression is less sharp,
the transition between the pressure states before and after the shock is smoother,
which results in an attenuated intermittency and a broader range of
frequencies covered by the shock unsteadiness [Ringuette et al., 2009].

Thanks to the possibility of considering the direct trace of the shock unsteadiness at the wall,
we are able to perform an exhaustive wavelet analysis
on the wall-pressure signal that demonstrates the power of a description of the stbli unsteadiness

in the time/time-scale domain.
The continuous wavelet transform is first directly applied to the unsteady signal
of the pressure at the wall, to identify the time evolution of the temporal scales characterising

the wall-pressure signature of the present stbli.

Then, we use a wavelet-based indicator to
evaluate the degree of intermittency of the wall-pressure, and to
extract from the signal its most energetic, intermittent part.
Inspired by Poggie & Smits, 1997 and by Camussi & Bogey, 2021,
we propose to use the wavelet transform to filter the wall-pressure and to separate the
contributions from the flow turbulence and the separation shock unsteadiness.
Our work tries thus to exploit the compact support of the wavelet basis to
highlight the intermittent nature of the large-scale unsteadiness of the system
and to define a practical procedure to extract relevant sparse features
from signals of interest for the system modelling.
The development of similar procedures to process the pressure – but also other quantities – is a

prerequisite for further investigations on stbli systems.

Finally, we consider the dynamics of the separation bubble behind the separation shock.
According to aforementioned theoretical, experimental, and numerical works
[Dupont et al., 2006, Piponniau et al., 2009, Pasquariello et al., 2017], the
motion of the shock system – especially for strong interactions – is strictly related
to the breathing motion of the separated region.
For this reason, we attempt to relate the intermittent features of the wall-pressure
at the foot of the separation shock to the dynamics of the separation bubble,
by replicating the wavelet intermittency analysis on a signal recording an estimate of the
recirculation region extent on a single vertical slice throughout the
duration of the simulation. Wavelet cross-spectra and wavelet coherence
are also considered to estimate the relationship between the wall-pressure at
the shock foot and the area of the separation bubble.

The paper is organised as follows: Section <ref> presents the
numerical setup of the simulations; Section <ref> describes the
database generated and the validation carried out; Section <ref>
presents the results of the analysis; finally, Section <ref>
reports some final comments.

§ NUMERICAL SETUP

The results presented in this work are obtained from
Direct Numerical Simulations performed using
streams [Bernardini et al., 2021], a high-fidelity solver of the compressible Navier-Stokes equations,
targeted to canonical wall-bounded turbulent high-speed flows. The code, freely available online,[<https://github.com/matteobernardini/STREAmS>]

solves the complete set of Navier-Stokes equations for a perfect, heat conducting gas.

The equations are discretised on a Cartesian mesh and solved by means of the finite difference approach. The convective terms
are discretised using a hybrid energy-conservative shock-capturing scheme in locally conservative form [Pirozzoli, 2010].
In particular, a sixth-order, central, energy-preserving flux formulation is adopted in smooth regions of the flow, which ensures
a robust and accurate discretisation of the wide range of spatial and temporal scales typical of turbulence,
without relying on numerical (artificial) diffusivity. Shock-capturing capabilities are achieved through the Lax–Friedrichs flux vector splitting,
where the characteristic fluxes are reconstructed at the interfaces using a fifth-order,
weno reconstructions [Jiang & Shu, 1996].
The hybridisation between the central and the weno scheme is managed by the implementation of a shock sensor that computes the local smoothness
of the solution and identifies discontinuities into the flow. Sixth-order, central
finite-difference approximations are also applied for
the discretisation of the viscous terms. The resulting system of ordinary differential equations
is integrated in time by means of a third-order, low-storage, Runge-Kutta scheme.
The solver is written in Fortran 90 and presents an MPI parallelisation based on a classical
domain-decomposition. The current version is able to run on NVIDIA multi-gpu architectures through the CUDA Fortran paradigm.
Additional details on the numerical methods and the scalability performance of streams can be found in Bernardini et al., 2021.

§ DATABASE DESCRIPTION AND VALIDATION

The dns database is inspired to
the oblique shock-wave/turbulent boundary layer
interaction reported in Dupont et al., 2006, Dupont et al., 2008,
characterised by a free-stream Mach number $M_\infty = 2.28$,
an incidence angle of the shock generator $\phi = 8^\circ$,
and a Reynolds number $Re_\theta = 5100$.

In the simulation, the Reynolds number based on the
compressible momentum thickness of the incoming turbulent boundary layer
and on the freestream viscosity is equal to $Re_{\theta} \approx 6900$,
which corresponds to $Re_{\delta_2} \approx 3900$ if wall viscosity is considered.

The main properties of the upstream boundary layer are reported
in table <ref> for both the numerical and the experimental case.

Properties of the incoming boundary layer: experimental and numerical data.

–       $M$     $Re_{\theta}$     $Re_{\delta_2}$     $Re_{\tau}$     $C_f$     $\theta_{0,inc}/\delta_0$     $H$     $H_{inc}$

Experiment     2.28     5100     2915     –       $2.00 \cdot 10^{-3}$     0.11     –       –

DNS            2.28     6882     3900     1100     $1.98 \cdot 10^{-3}$     0.10     3.21     1.33

The simulation is conducted in a computational domain of size $L_x/\delta_0 \times L_y/\delta_0 \times L_z/\delta_0 = 70 \times 12 \times 6.5$,
the reference length $\delta_0$ being the thickness of the incoming boundary layer upstream of the interaction ($99 \%$ of the freestream velocity).
According to results of [Morgan et al., 2013],
the spanwise width of the computational domain that we use
is sufficient to avoid any confinement effect and thus
any artificial increase in the size of the separation bubble.
The computational mesh consists of $N_x \times N_y \times N_z = 8192 \times 1024 \times 1024$ grid nodes, equally spaced in the wall-parallel directions.
In terms of wall units evaluated in the upstream boundary layer, the spacings in the streamwise and spanwise directions
are $\Delta x^+ \approx 7$, $\Delta z^+ \approx 5$.
A stretching function is applied in the wall-normal direction to increase the resolution in the near-wall region and in the interaction zone.
At the wall, the spacing $\Delta y^+$ varies in the streamwise direction and it ranges between $\Delta y^+ \approx 0.6$
and $\Delta y^+ \approx 1.1$, upstream and downstream of the interaction respectively.

The boundary conditions are specified as follows. At the outflow, non-reflecting
conditions are imposed by performing a characteristic decomposition in the direction normal to the boundary [Poinsot & Lele, 1992].
A similar treatment is also applied at the top boundary away from the incoming shock, where instead the inviscid Rankine-Hugoniot
jump conditions are locally imposed to mimic the presence of a shock generator.
A characteristic wave decomposition is also employed at the bottom no-slip wall, where the wall temperature is set to
the recovery value for the upstream boundary layer.
A critical issue in the simulation of spatially evolving turbulent flows is the prescription of an inflow turbulence
generation method. In streams, velocity fluctuations at the inlet plane are imposed by means of a synthetic df
approach [Klein et al., 2003], extended to the compressible case by means of the strong Reynolds analogy [Touber & Sandham, 2009].
For stbli computations, this type of approach is preferable with respect to alternatives based on the recycling-rescaling procedure [Lund, T.S. et al., 1998],
because it does not introduce spurious frequencies that might interact with the dynamics of the reflected shock.
An efficient implementation of the method is here obtained using an optimised df procedure [Kempf et al., 2012],
whereby the filtering operation is decomposed in a sequence of fast one-dimensional convolutions.
The implementation requires the specification of the Reynolds stress tensor at the inflow plane, which is interpolated by
a data set of previous dns of supersonic boundary layers performed by the same group [Pirozzoli & Bernardini, 2011].
Finally, in the spanwise direction the flow is assumed to be statistically homogeneous, and periodic
boundary conditions are applied.
In order to characterise the incoming turbulent boundary layer,
figure <ref> shows a comparison of the van Driest–transformed mean velocity profiles
and of the density-scaled Reynolds stress components with the experimental data
from [Elena & Lacharme, 1988] and [Piponniau et al., 2009].
A satisfying agreement can be observed for all the distributions,
except for the wall-normal Reynolds stress component,
which is typically underestimated by measurements.

Comparison of (a) van Driest–transformed mean velocity profile and (b) density-scaled Reynolds
stress components for the incoming boundary layer with reference experimental data. Symbols denote
experiments by [Elena & Lacharme, 1988] (diamond, $M_\infty = 2.32$, $Re_\theta = 4700$) and
[Piponniau et al., 2009] (squares, $M_\infty = 2.28$, $Re_\theta = 5100$).

Various studies focused on the computational analysis of the confinement effects imposed by the presence of sidewalls
for stbli [Bermejo-Moreno et al., 2014, Poggie & Porter, 2019, Deshpande & Poggie, 2021].

In particular, using wall-modeled les with an equilibrium wall model
to replicate our same reference experiment, Bermejo-Moreno et al., 2014
found that the mean separation bubble behind the interaction

is characterised by a strong three-dimensionality imposed by the lateral confinement,
and observed that the size of the bubble is significantly larger than that predicted by spanwise-periodic simulations.
Therefore, to compare our results with the experimental measurements, in the following we use scaled interaction coordinates,
$x^* = (x-x_0)/L_{int}$ and $y^* = y/L_{int}$, where $L_{int}$ denotes the interaction length-scale, defined as the distance between the
nominal impingement location and the foot of the reflected shock, and $x_0$ denotes the mean position of the
ideal incident shock.[In the reference experimental work, $x_0$ is
equal to the mean position of the reflected shock.
As a result, our values of $x^*$

are quantitatively different from those in Dupont et al., 2006,
even if the actual positions are
For the present dns data, the ratio $L_{int}/\delta_0$ is equal to $3.30$,
a value larger than that previously obtained at reduced Reynolds number ($L_{int}/\delta_0 = 2.89$)
but smaller than the experimental measurement ($L_{int}/\delta_0 = 4.18$).
Table <ref> reports a comparison of the length scales obtained
with some references.

Comparison of the length scales of the stbli considered ($M_\infty \approx 2.28$, $\phi = 8^\circ$) with some references.
The superscript * indicates that the separation length has been estimated by the relation $L_{sep} = 0.8\,L_{int}$ [Clemens & Narayanaswamy, 2009].

–       $M$     $Re_{\theta}$     $L_{int}/\delta_0$     $L_{sep}/\delta_0$

Experiment [Dupont et al., 2006]     2.28     5100     4.18     3.40*

Present case                           2.28     6882     3.30     2.16*

[Morgan et al., 2013]                      2.28     4800     3.02     1.61*

[Aubard et al., 2013]                 2.25     3700     2.84     2.35*

[Pirozzoli & Bernardini, 2011]             2.28     2300     2.89     2.31*

[Touber & Sandham, 2009]                 2.30     5100     4.80     3.90*

[Pirozzoli & Grasso, 2006]             2.25     3700     2.17     1.18*

After an initial transient needed to develop the reflected shock and to achieve a statistically steady state,
the computation was advanced for a total time $T \approx 2000\,\delta_0/U_{\infty}$ which is much longer than
the previous dns simulations. It is worth highlighting that the time interval here considered allows us to cover several
low-frequency cycles (approximately 25 cycles), thus making possible an in-depth analysis of the wall-pressure signature unsteadiness.
The simulation time step is $\Delta t \approx 0.0008\,\delta_0/U_{\infty}$ and the wall-pressure history
is recorded with a sampling time $\Delta t \approx 0.04\,\delta_0/U_{\infty}$.

A comparison of the dns results with the piv data measured in the reference experiment is shown in figure <ref>,
which reports the contours of the mean streamwise velocity, mean vertical velocity, streamwise velocity fluctuations,
vertical velocity fluctuations, and Reynolds shear stress. In the figure, contour lines obtained from dns
are superposed onto heat maps derived from piv in a limited region close to the interaction zone.
The comparison highlights that the overall structure of the interaction is recovered effectively. The numerical simulation reproduces correctly the thickening of the upstream boundary
layer and the amplification of the turbulence across the interaction,
which are associated with the shedding of vortices in the shear layer
developing at the separation shock.
It is worth highlighting that, although a similar qualitative comparison was also previously observed at reduced Reynolds number,
the present dns results provide a better quantitative agreement with the experiment
with respect to previous works,
in particular with regard to the intensity of the turbulent fluctuations.
Figure <ref> shows an instantaneous contour of the non-dimensional temperature on a
vertical slice, in which the black line indicates the locus of the points characterised by the sonic Mach number.
The instantaneous separation region is highlighted by indicating the edge of the local separation bubble.
From the figure, we can observe the classical structure of the stbli.
The thickening of the boundary layer, induced by the adverse pressure gradient from the shock system, makes
the sweeps and ejections of the fluid towards and from the wall more evident.
Moreover, the reflected shock generates an upward deflection of the flow before the virtual
impingement point of the incident shock, thus promoting separation
between $x^* = -0.5$ and $x^* = 0.0$.
On top of the recirculation region, the innermost part of the incident shock interacts with the
vorticity of the supersonic boundary layer in a non-trivial way, inducing even subsonic regions
close to the intersection point between the incident and the reflected shocks.
Finally, the following expansion and recompression
that bring back the fluid to the original mean direction are also evident, especially in the top part of the

[$\bar{u}/U_\infty$  ]
[$\bar{v}/U_\infty$ ]
[$\overline{u'^2}/U_\infty^2$ ]
[$\overline{v'^2}/U_\infty^2$ ]
[$-\overline{u'v'}/U_\infty^2$ ]
Comparison of velocity and Reynolds stress components on the spanwise centre plane in the interaction region. Colour maps refer to experimental PIV data, whereas DNS data
are shown using contour lines. From top to bottom: mean streamwise velocity, mean vertical velocity, streamwise velocity fluctuations, vertical velocity fluctuations and Reynolds shear stress.

Instantaneous contour of the non-dimensional temperature on a vertical slice.
The black line indicates the points characterised by Mach number equal to 1,
whereas the white line highlights the edge of the instantaneous recirculation bubble.

[$\bar{p}_w/p_\infty$  ]

[$\sqrt{\overline{p'^2}}/p_\infty$  ]

Streamwise evolution of (a) mean wall-pressure and (b) wall-pressure fluctuation intensity across the interaction region.
$\bullet$ Experimental data [Dupont et al., 2006], – Present DNS,  –
Present DNS rms evaluated by integrating the frequency spectra up to the cut-off of the experimental sensors.

§ ANALYSIS OF THE RESULTS

In this section, we focus on the behaviour of the wall-pressure,
by investigating the frequency and the intermittent
content of the time signals collected around the shock region
by means of classical Fourier analysis and by means of wavelet decomposition.
Extracting the features of this signal
means understanding the behaviour of the shock unsteadiness,
which is the main responsible for the dynamics of the wall-pressure fluctuations
besides turbulence.
In addition, we also provide an analysis of the dynamics
of the separation bubble in the time/time-scale domain,
and we discuss the relationship between the unsteadiness of the shock
and the motion of the recirculating region.

§.§ Time-averaged statistics of the wall-pressure

At first, we focus on the behaviour of the mean pressure along the wall.
Figure <ref> shows the streamwise evolution of the mean wall-pressure,
obtained by averaging the signal in time and in the spanwise direction, together with the experimental
results from the reference work [Dupont et al., 2006].
Upstream of the interaction, the mean pressure is approximately constant and equal to the free-stream value.
Later on, the pressure increases smoothly because of the presence of the oscillating shock,
and after approximately $2\,L_{int}$ it reaches a plateau region.
The comparison with the experimental data is satisfactory.

Figure <ref> reports instead the comparison of the streamwise distribution of the
rms of the wall-pressure fluctuations.
The typical behaviour shown in other stbli is observable.

At first, downstream of the foot of the shock, the intensity of the fluctuations increases suddenly

because of the backward and forward movement of the reflected shock.
Then, in the separation and re-attachment region, the intensity decays,
although the level of the fluctuations remains always higher than
that of the attached boundary layer.
The numerical results (solid black line) seems to show an overestimation of the fluctuation
intensity, but if the rms is evaluated
integrating the frequency spectra up to the cut-off
of the experimental sensors (solid red line),

the results agree well with the measurements.

Friction coefficient and premultiplied power spectral density of the wall-pressure along the interaction. The horizontal  blue line
in the first plot indicates the null value of $C_f$, whereas
the vertical dashed lines indicate the position of the probes used for the sampling of the
wall-pressure signals.

§.§ Fourier analysis of the wall-pressure

The premultiplied psd of the wall-pressure signals has been computed for positions ranging
from the zone before the shock to the
relaxation zone, and the results are reported in figure <ref>.
This map reveals the frequencies that contribute the most to the global energy of the local
wall-pressure fluctuations,
and the characteristics that clearly emerge agree with what was shown in Dupont et al., 2006.
In order to better interpret the different regions of the interaction, we report on top of the
map also the distribution of the time- and spanwise-averaged friction coefficient $C_f = \tau_w / 1/2 \rho U_\infty^2$, where $\tau_w$ is the wall shear stress.

Consistently with the literature, the mean friction coefficient is characterised by
the presence of two negative peaks, where the upstream one has typically a larger amplitude in absolute value.
However, [Morgan et al., 2013] showed that as the Reynolds number increases, the absolute magnitude of the
first negative peak of the $C_f$ decreases, which thus explains the distribution observed in our case

leading to the appearance of two completely separate separation bubbles in the mean field.
However, the behaviour observed in the mean field is actually
made of the superposition of different instantaneous states,
with the separation bubble alternating between a smaller and a
significantly larger extent compared to what indicated by the negative
mean friction coefficient.
On the basis of these observations, we define as a separation length
the extent that goes from the first $C_f$ minimum at approximately
$x^* \approx -0.8$ to the point with null $C_f$ at approximately $x^* \approx -0.15$.
A separation length of $L_{sep}/L_{int} \approx 0.65$ or $L_{sep}/\delta_0 \approx 2.16$
is finally obtained, which is in line with the trend observed in [Morgan et al., 2013]
for increasing Reynolds number.

The frequencies in the contour of figure <ref> can be conveniently split into four separate zones, each involving its own characteristic temporal scales:

* the incoming turbulent boundary layer zone, $x^*< -0.8$, characterised by high-frequency contributions, with $S_L = f \, L_{int}/ U_\infty > 1$, and attached flow;

* the unsteady foot of the reflected shock, $x^* \approx -0.8$, characterised by low-frequency energy content.
The peak
is widespread along a broad range of frequencies, suggesting
an underlying shock unsteadiness more complicated than a simple periodic oscillation.
In this region, the adverse pressure gradient imposed by the moving shock reduces the friction coefficient
and forces the flow to start the separation;

* the interaction zone, $-0.8 < x^*< -0.2$, characterised by intermediate turbulent scales developing in the detached supersonic shear layer. In this region, the spectrum reaches also its overall maximum in the point of minimal friction,
which corresponds to the centre of the aft portion of the recirculation bubble;

* the relaxation zone, $x^* > -0.2$, characterised by medium-frequency fluctuations developing in the separation region and dominating the boundary layer even after the reattaching point.

One of the most interesting characteristics revealed by the spectral map is the substantial distinction between the
frequency content of the shock movement and
the frequency content of the attached boundary layer, the separated shear layer, and the relaxation region.
Figure <ref> shows the premultiplied psd in correspondence of the highlighted probes
at $x^* = -1.60$, $-0.79$, $-0.27$ and $0.97$

to better appreciate the spectral features of the fluctuations in points
representative of the above presented zones.
In particular, at $x^* = -0.79$, the broad and high-amplitude peak
at $S_L \approx 0.04$
reflects the low-frequency unsteadiness
of the shock motion, while the second wide bump with a peak at $S_L \approx 1.5$

is the trace of the incipient separated supersonic shear layer.

Premultiplied power spectral density of the wall-pressure signal at:

$\cdot \cdot \cdot$ $x^* = -1.60$, - - $x^* = -0.79$, -$\cdot$- $x^* = -0.27$,
– $x^* = 0.97$.
The vertical dashed line indicates the low-frequency peak at $S_L = 0.04$.

§.§ Wavelet analysis and detection of intermittency of the wall-pressure

Elements of the wavelet theory can be found in several texts [Kaiser & Hudgins, 1994, Mallat, 1999],
whereas Lewalle, 1994 and Farge, 1992 discussed its applications to fluid mechanics.
Here we report few elements to clarify the following discussion.

The wavelet transform of a continuous signal $g(t)$ is defined as:

\begin{equation}
G_{\Psi}(k,\tau) = k \int_{-\infty}^{+\infty} g(t)\Psi^{*}(k(t-\tau))dt \,,
\label{eq:def}
\end{equation}

where $\Psi$ is the wavelet mother function, $k$ is a dilatation parameter,
$\tau$ is the time-translation parameter and $*$ denotes the complex conjugate.
By varying the wavelet scale $k$ and translating along the time with the time shift $\tau$,
one can construct a picture showing the amplitude of any feature at a certain scale and
also its temporal evolution.
In this study, the Morlet wavelet has been chosen since higher
resolution in frequency can be achieved when compared with other mother functions.
An analytical expression for the complex Morlet wavelet is:

\begin{equation}
\Psi(t) = \pi^{-1/4} e^{i \omega_0 t} e^{-t^2/2}
\label{eq:mother}
\end{equation}

where $\omega_0$ is the non-dimensional frequency, which is equal to 6 here to
satisfy the admissibility condition [Farge, 1992].
If $\hat{g}(\omega)$ is the Fourier transform of $g(t)$:

\begin{equation}
\hat{g}(\omega) = \int_{-\infty}^{+\infty} g(t) e^{-i \omega t} d\omega\,,
\label{eq:Fourier}
\end{equation}

then it follows that:

\begin{equation}
G_{\Psi}(k,\tau) =\frac{1}{2\pi} \int_{-\infty}^{+\infty} \hat{g}(\omega) \hat{\Psi}^*\left(\frac{\omega}{k}\right) e^{i \omega t} d\omega\,.
\label{eq:invFourier}
\end{equation}

As a consequence, the wavelet transform at a given scale $k$ can be interpreted as
a band-pass filter in the Fourier space.
Finally, following the method of Meyers et al., 1993, the relationship between
the equivalent Fourier period and the wavelet scale can be derived
analytically by substituting a cosine wave of known frequency
into Equation <ref> and computing the scale $k$ at which the wavelet power spectrum reaches its maximum.
More details about the operative algorithms adopted to compute the wavelet transform in this work
can be found in Torrence & Compo, 1998.

In the following, we apply the wavelet analysis to the wall-pressure signal
at four relevant stations, respectively at
$x^* = -1.60$, $-0.79$, $-0.27$ and $0.97$
(see also figure <ref> for the position of the numerical probes).
In particular, we first focus on the regions
before and in the proximity of the foot of the reflected shock;

then, we address the recirculating region, and finally we consider the re-attachment location.

[$x^* = -1.60$ ]

[$x^* = -0.79$ ]

[$x^* = -0.27$ ]

[$x^* = 0.97$ ]

Wall-pressure signals (top) and
Morlet wavelet transform modulus $|G_{\Psi}(k,\tau)|$ (bottom) at

(a) $x^* = -1.60$,
(b) $x^* = -0.79$,
(c) $x^* = -0.27$,
(d) $x^* = 0.97$.
For each probe, only values of the coefficients above 30% of the relative maximum are shown.

On the top left, the green curve reports the original
wall-pressure signal at $x^* = -0.79$,
while the black curve reports its filtered-and-reconstructed part.
On the top right, the pdf of the original signal with
the corresponding fitted normal distribution in gray.
On the bottom left, contour of $LIM2$ of $p_w/p_\infty$.
Values below 3 are cut off, whereas the gray area indicates the scales used for the filtering procedure.
On the bottom right, the corresponding WFF.
The horizontal dashed lines indicate the peak low frequency.

Figure <ref> presents the modulus of the Morlet transform $G_{\Psi}(k,\tau)$,
the so-called scalogram, of the wall-pressure signals at
the four prescribed locations.

The abscissa reports the non-dimensional time
and the ordinate reports the temporal scale $k$ of the events in log-scale.
The blue color indicates a small magnitude of $|G_{\Psi}(k,\tau)|$,
while the red color indicates a larger wavelet modulus,
representing a more intense modulation of the amplitude
of the pressure fluctuations within a certain range of temporal scales $k$
and at the given time instant $\tau$. Indeed, the real part of the wavelet coefficient
is proportional to the amplitude of the fluctuation ($p'$), while the
imaginary part is proportional to its time variation ($\partial p'/\partial t$) [Abe et al., 1999].

Figure <ref> shows the modulus of the Morlet transform for the wall-pressure
recorded at $x^* = -1.60$,
which is inside the boundary layer and before the foot of the reflected shock.
It can be seen that, as expected, the
attached boundary layer is characterised by small-scale
eddies, which appear to be randomly distributed and fairly time-filling,
although the spotted appearance underlines the intermittent character of such turbulent structures.
Figure <ref> shows the wavelet modulus at $x^* = -0.79$,
immediately downstream of the foot of the reflected shock.
Given the position, the pressure signal from this probe is clearly modulated by the shock movement
in the amplitude and in the characteristic temporal scales,
or instantaneous frequencies if we consider the inverse of the temporal scales.
The contour confirms that the shock movement has a dominant content
at low frequencies/large time scales, as observed in the Fourier spectra.
However, it reveals also that

this motion is actually made of a collection
of events localised in time and characterised by different temporal scales,
larger than and well separated from the ones present in the upstream boundary layer.
The spectral content from the Fourier analysis is thus the result of
an average in time of such events, whose combined effect is the broad low-frequency peak
spotted in figure <ref>.
Figure <ref> characterises instead the pressure signal at $x^* = -0.27$,
in the small separated region surrounded by the detached supersonic shear layer.
Here, the wavelet map is representative of a shock-induced turbulent recirculating flow:
the shock footprint is well visible from the distribution in time of large scale events,
appearing and disappearing in time, but also the trace of the separated shear layer emerges
from the presence of small-scale, time-filling events,
similarly to the region of the upstream boundary layer.
Finally, figure <ref> shows the wavelet map at $x^*=0.97$,
the re-attachment location. Large-scale events are now more sparse, indicating that
the effect of the shock unsteadiness is weaker and a new attached boundary layer is developing.

The scalograms shown above reveal the presence of different events, or structures, at various scales,
and allows us to detect their degree of intermittency.

In this context, we define the intermittency as the
alternation of regimes with normal spectral content
and regimes with significant excess – or burst – of energy in a given range of scales.
Given the wavelet transform coefficients $G_{\Psi}(k,\tau)$, it is possible
to obtain the scale-time distribution of the energy density $|G_{\Psi}(k,\tau)|^2$ of
the wall-pressure signal $p(t)$ at specified scale $k$ and time $\tau$.

Taking into account this property, Meneveau, 1991 and Camussi & Bogey, 2021 suggested that an
effective indicator of the intermittency is the squared local intermittency measure, denoted as $LIM2$:

\begin{equation}
\label{eq:LIM2}
\end{equation}

where $\langle \bullet \rangle_\tau$ indicates the time average of the considered quantity.
The rationale for this statement lies in the fact that
$LIM2$ can be interpreted as a time-scale dependent measure of the flatness factor or kurtosis
of the input signal $g(t)$, which indicates the importance of rare events for the probability
distribution of the variable and is defined as the Pearson's index:

\begin{equation}
FF=\frac{\langle g^{'~4} \rangle}{\langle g^{'2}\rangle^2} \,,
\label{eq:FF}
\end{equation}

where $g'$ is the fluctuating part of $g$, and $\langle \bullet \rangle$ indicates the
expected value of the considered quantity.
Therefore, the $LIM2$ parameter will be equal to 3
when the probability distribution is Gaussian, while the condition $LIM2>3$ identifies
only those rare outliers contributing to the deviation of the wavelet coefficients
from a normal, Gaussian distribution.
According to Camussi & Bogey, 2021,
the energy increment at a certain scale $k$ and time $\tau$ is
associated with the passage of a coherent structure characterised by the scale $k$.

Figure <ref> (bottom left) shows the $LIM2$ field for the signal at
$x^*=-0.79$, which has a general kurtosis (see equation <ref>)
equal to 3.20.
Only the levels greater than the threshold value 3 are shown.
In this picture, one can observe the occurrence of the intermittent events
in the time-scale region characterising the shock unsteadiness (indicated in gray in the background),
which is the range of large temporal scales defining the broad low-frequency peak
in the Fourier spectrum (see figure <ref>).
Now the significant energy bursts appear more sparse in time
with respect to the representation offered
in figure <ref>.
In addition, we can also see the intermittency of the turbulent content
at smaller temporal scales.
Moreover, we notice that the rather long period considered for the simulation we carried out
is able to capture only few coherent intermittent events at large time scales. For this reason,
it is evident that, to characterise thoroughly the unsteadiness of stbli,
very long periods must be considered, in order to deconstruct the broadband low-frequency unsteadiness
reported in the literature into the actual scattered dynamics associated to the flow.

A sign of the presence of intermittent events is already evident from the pdf
of the signal (top right),
which shows a distribution slightly departing from a normal Gaussian function in
its external parts. However, although classical statistical indicators
suggest the relevance of outlier events with large deviation from the mean,
as indicated by the mildly non-gaussian pdf and kurtosis of the signal, they are not able to locate
and estimate precisely the single intermittent episodes,

which thus makes it hard to understand the physical mechanisms
behind the resulting dynamics.
The $LIM2$ measure instead allows one to spot the instants and the scales
of data that differ significantly from other observations, and thus makes it possible
to characterise more accurately the dynamics of the flow, which may be hidden in time-averaged analyses.

After the inspection of such a scenario, we propose
the following procedure to extract the characteristics of the signal
of interest [Consolini & De Michelis, 2005]:
first, we filter the original wall-pressure using equation <ref>,
considering only the wavelet coefficients of the large scales
associated with the low-frequency shock unsteadiness
(gray region in the $LIM2$ contour);
then, we reconstruct the intermittent signal $p_{w,I}$ only on the basis of the wavelet coefficients
related to the burst events
by means of the conditioned inverse continuous wavelet transform:

\begin{equation}
p_{w,I}(k_{min}, t)=\frac{1}{C_{\Psi}} \int_{k_{min}}^{\infty} dk \int_{-\infty}^{\infty} G_{\Psi,\,LIM2}(k,\tau)\Psi\left(\frac{\tau-t}{k}\right)d\tau,
\label{eq:invLIM}
\end{equation}

where $C_{\Psi}$ is a normalisation constant depending on the chosen wavelet,
$k_{min}$ is the smallest time scale considered by the large-scale-pass filtering,
and $G_{\Psi,\,LIM2}(k,\tau)$ are the wavelet coefficients related to the
intermittent events, extracted using the $LIM2$ condition ($LIM2(k,\tau) > 3$).

This procedure can be considered as an extension of the one adopted by Poggie & Smits, 1997,
in which the authors used the extremes in the global wavelet spectrum to filter the signal
in the space of the wavelet temporal scales.

Figure <ref> (top) shows the original wall-pressure in time at the foot of the reflected shock
(green line) together with its filtered-and-reconstructed intermittent part (black line).
The inspection of the reconstructed signal allows us
to appreciate that the small-scale turbulent fluctuations have been removed by the procedure, and
that only the pattern given by the events associated with the shock intermittency remains.

Therefore, it is possible to detect the signal shape
given by the movement of the shock, which is governed by
a predominantly two-dimensional mechanism (high span-wise coherence) according to
Sasaki et al., 2021.

The $LIM2$ map can be averaged in time to recover the flatness factor
as a function of the temporal scale [Camussi & Bogey, 2021],
and thus to reveal the most intermittent temporal scales.
This quantity is usually called wff and is defined as:

\begin{equation}
WFF(k)=\langle LIM2(k,\tau)\rangle_{\tau}=\frac{\langle|G_{\Psi}(k,\tau)|^4\rangle_{\tau}}{\langle|G_{\Psi}(k,\tau)|^2\rangle_{\tau}^2}\,.
\label{eq:WFF}
\end{equation}

Figure <ref> (bottom right) shows the wff for the corresponding $LIM2$ map
of the wall-pressure at the shock foot.
We can see that the wff is characterised by two peaks at large scales,
one at $k\,U_\infty/L_{int} \approx$ 15 and the other at $k\,U_\infty/L_{int} \approx$ 36.

However, the peak frequency computed by the Fourier and wavelet spectra
lays in between the scales indicated by the WWF, at $S_L = 0.04$ or $k\,U_\infty/L_{int} = 24.27$.
This result represents an indication that the most probable temporal scales of intermittency
are different from the most probable temporal scales of the overall unsteadiness.

§.§ Wavelet intermittency analysis of the separation bubble dynamics

In view of the many works sustaining that the low-frequency shock motion is strictly related to the dynamics of the separation bubble [Dupont et al., 2006, Piponniau et al., 2009, Pasquariello et al., 2017],
we then attempt to relate the behaviour of the wall-pressure at the foot of the shock
to the change in time of the separation region, to assess the eventual correspondence between these two features
and in the attempt of providing a characterisation of their link locally in time and time scale.
In fact, [Piponniau et al., 2009] highlighted that, far from weak interactions with incipient
separation, where there seems to be an actual correlation [Humble et al., 2007], and at least for the case of shock
reflections, there is no relevant coherence between the upstream superstructures
and the dynamic response of the system as suggested by [Ganapathisubramani et al., 2007].
Therefore, for the case under consideration, the dynamics of the separated bubble seems
clearly the main source of the low-frequency unsteadiness.

In order to detect the zone occupied by the recirculating flow,
we recorded the whole area occupied by negative streamwise
velocity at each instant,
assuming that this provides a rough estimate of the extent of the bubble.
The time signal of such area, indicated with $A_{u<0}/A_{ref}$, is reported in figure
<ref> (top) together with its corresponding filtered-and-reconstructed part.
On the right, we report also the pdf of the original signal which has a kurtosis equal to 3.68.
The extent of the separation bubble oscillates between two more probable values
but presents also many sudden variations at larger areas. The probability density function
reports thus a highly skewed, bimodal distribution with a long tail that reveals the presence of important
and intermittent fluctuations in the breathing motion of the separation area.
Contours of the Mach number for two different states during an intermittent event are also shown
in figure <ref>, in which the points with null streamwise velocity
are indicated with a white line that gives a hint of the area $A_{u<0}$.
In the first state, the flow is largely separated, and the deviation imposed on the
flow by the presence of the separation bubble is significant.
On the other hand, in the second state, the flow is almost
completely attached, and the deviation is practically negligible.[
A movie of the time evolution of the contour of the Mach number on
a vertical slice during the intermittent event around
$t\,U_\infty/\delta_0 \approx 1160$
is provided in the supplementary material.]
In order to analyse and compare the intermittency of the wall-pressure induced by the oscillation
of the shock and of the fluctuations of the separation extent,
we report a contour of the $LIM2$ measure also for $A_{u<0}$ (figure <ref>, bottom).
From the comparison with figure <ref>,
we can observe that some of the most important events of the $LIM2$ of the
wall-pressure at the shock foot correspond to the most important events of the corresponding measure
of the separation bubble (see for example $t U_\infty / \delta_0 \approx 200$, $1160$, and $1380$), even though the temporal scales are different.
The situation is even more evident when we look at Figure <ref>, which reports the time evolution of the $LIM2$ intermittency measure averaged in scale, considering only the large scale interval that characterises the low-frequency unsteadiness.
The intermittency of the pulsation of the separation bubble is thus transmitted to the
unsteadiness of the shock system, and thus to its wall-pressure signature at the wall, or vice versa.
However, although many of the intermittent events of the wall-pressure have a corresponding
trace in the $LIM2$ measure of the separation bubble area, other energy bursts

for the wall-pressure do not seem to have a direct relation
with the breathing mode of the system, e.g. $t U_\infty / \delta_0 \approx 485$, $887$, or $1270$.
This observation suggests that the low-frequency, intermittent unsteadiness of the shock
can only be in part related to the breathing motion of the separation bubble, and so that
other key actors may play a role in this phenomenon.

Finally, if we observe also the $WFF$ of the area signal (figure <ref>, bottom right),
we can see that the distribution of the flatness factor along the time scales follows closely the
behaviour of the $WFF$ of the wall-pressure (figure <ref>, bottom right).
As a result, also for the dynamics of the separation bubble extent,
the most probable intermittent scales are different from
the large temporal scales characterising the overall shock unsteadiness
(the broad low-frequency peak in the Fourier domain).

On the top left, the green curve reports the original
$A_{u<0}/A_{ref}$ signal,
while the black curve reports its filtered-and-reconstructed part.
On the top right, the pdf of the original signal.

The two peaks of the distribution are indicated with grey dotted lines.
On the bottom left, a contour of $LIM2$ of $A_{u<0}/A_{ref}$.
Values below 3 are cut off, whereas the gray area indicates the scales used for the filtering procedure.
On the bottom right, the corresponding WFF.

[$t \,  U_\infty / \delta_0 = 1158.1$]

[$t \,  U_\infty / \delta_0 = 1195.6$]

Instantaneous contours of the Mach number on a vertical slice for two different states during the intermittent event around $t  \,  U_\infty / \delta_0 \approx  1160$.
The white line indicates the points with null streamwise velocity.
On the left, the area occupied by the bubble is large; on the right, the bubble is almost absent.

Normalised scale-average of the $LIM2$ measure above threshold value 3,
considering only the large scales – low frequencies – used for the filtering procedure.

– Wall-pressure signal at the shock foot, - - Separation bubble area signal.

In order to further characterise the relationship between the signals of the wall-pressure at the
foot of the shock and the separation area, we consider also the information from the wavelet cross spectrum
and from the wavelet coherence [Torrence & Webster, 1999, Grinsted et al., 2004, Maraun et al., 2007].
Through such analyses, we are able to determine whether there are regions in
the time/time-scale domain with large common power and also if these have consistent phase relationship,
which could potentially suggest the causality between the two time series.
In particular, the wavelet cross spectrum of two time series $x$ and $y$,
with wavelet transform $G_\Psi^{x}$ and $G_\Psi^{y}$ respectively, is defined as:

\begin{equation}
G_{x,y} (k, \tau) = G_\Psi^{x} \left(G_\Psi^{y}\right)^*.
\end{equation}

with $*$ indicating the complex conjugation. The wavelet cross spectrum represents thus the covarying power of the two processes.
This means that an independent variance appearing only in one of the two single spectrum
does not appear in the cross spectrum, and that the cross spectrum
vanishes for independent processes. Moreover, for the related time and time scales, the complex argument
of the wavelet cross spectrum $\phi(G_{x,y})$ gives an indication of the local relative phase between the
original time series in the time/time-scale domain.
For example, to establish a simple cause and effect
relationship between two phenomena represented by the two time series, oscillations are expected to be locked in phase in the
regions with significant common power [Grinsted et al., 2004]. Moreover, starting from the definition of the wavelet cross spectrum,
it is possible to define a direct measure of the correlation between two signals throughout the time and the scales.
Torrence & Webster, 1999 defined the so-called wavelet coherence $R^2_{x,y}$ between two time series $x$ and $y$ as:

\begin{equation}
R^2_{x,y} (k, \tau) = \frac{|S(G_{x,y})|^2}{S(|G_\Psi^x|^2) S(|G_\Psi^y|^2)}\,,
\end{equation}

where $S$ is a smoothing operator in time and time-scale. The wavelet coherence has values between zero and one,
and has the relevant property of quantifying the linear relationship between two processes [Maraun et al., 2007].
Figure <ref> shows the wavelet coherence between the signal of the wall-pressure at the foot of the
shock and that of the separation area. In the first instance, we consider the original signals (figure <ref>).
At small time scales, the turbulent fluctuations results in the typical pattern with long,
intermittent, vertical strikes covering the high-frequency range.
The coherence is here mostly related to the properties of the small turbulent coherent structures.
In the large-scale – low-frequency – range, large areas of the contour show that the two signals present a significant
instantaneous correlation, although they are not phase locked.
However, the distribution of the peak regions seems to contradict the correspondence
of the intermittent events inferred from the comparison of the $LIM2$ measure of the two signals (figure <ref> - figure <ref>).
For this reason, we consider also the wavelet coherence of the two signals filtered and reconstructed according to
the procedure proposed above (figure <ref>). Given the large-scale-pass preliminary filtering
of the signal, the small time-scale region is not considered. From the figure, we can observe an important coherence
with strong events sparse in time, and most of all, we can recognise the trace of the most important intermittent
events pointed by the $LIM2$ comparison previously discussed. The signal processing isolates the intermittent content of the
signal, and the wavelet coherence analysis allows us to reveal that the intermittency of the wall-pressure at the
foot of the shock is strictly related to the pulsation of the separation area. However, if we look at the phase angle of
the cross wavelet between the filtered-and-reconstructed signals in figure <ref> for the
most coherent events, we can see that also for the intermittent signals the phase does not remain constant across the
time scales, indicating that a simple linear relationship between the intermittency of the breathing mode
and that of the unsteadiness of the shock can not be inferred.

[Original signals. ]

[Intermittent signals. ]

Wavelet coherence between the shock foot wall-pressure signal and the separation area.

Arrows indicate the relative phase from the wavelet cross spectrum for the most coherent regions.
Shaded region in the right contour is not significant given the preliminary large-scale-pass filtering.

Cross wavelet phase angle for times and time scales for which the wavelet coherence between
the filtered wall-pressure and separation area is larger than 0.6.

§ CONCLUSIONS

The present work analyses the wall-pressure fluctuations of the stbli
generated by an impinging oblique shock wave,
by leveraging the enhanced capabilities of the wavelet analysis
in the time/time-scale domain [Mallat, 1999].

The data are obtained numerically by means of a dns carried out at moderate Reynolds number, so that

the shock is able to penetrate into the
turbulent boundary layer, leaving a marked trace in the wall-pressure signal
that makes it possible to examine the shock unsteadiness
and its intermittency more adequately than in the cases with lower Reynolds
number [Pirozzoli, 2011].

Although the
classical Fourier analysis provides useful information
about the spectral content of the overall stbli unsteadiness,
other techniques may be better able to analyse the transient phenomena that are taking place.
Among them, given its compact support, the wavelet transform represents
a valuable tool that is still scarcely used in flows similar to the one considered here.

From the analysis of the wall-pressure signals and of the time evolution of the separation bubble,
we attest in this work

the ability of wavelet analysis to infer relevant information about the dynamics of stbli.

In the first place,
the validation of our results in terms of basic
statistics of the fluid field and of the wall-pressure
shows a notable improvement in the reproduction
of the reference experiment with respect to previous works.
Furthermore, the analysis of the wavelet transform modulus
reveals and clarifies
that the broadband low-frequency peak, usually observed in the Fourier spectrum
at the foot of the reflected shock in stbli [Clemens & Narayanaswamy, 2014],
is actually the result of a collection of sparse, intermittent events
with different temporal scales.
The Fourier spectral analysis shows thus only the time-averaged effect
of such behaviour and is not able to characterise the
local modulation of the pressure fluctuation, which makes it impossible to
understand precisely what is happening in the fluid flow correspondingly.

Then, we show that it is possible to characterise the energy and timing of
the intermittent events linked to the shock unsteadiness
on the basis of a wavelet-based local flatness factor,
called $LIM2$ [Meneveau, 1991, Camussi & Bogey, 2021].
The global kurtosis of the wall-pressure, related to the outliers in the shock signature,
is thus decomposed into its evolution in the time/time-scale domain, to localise
the most significant bursts of energy in the signal.
Inspired by the work of Poggie & Smits, 1997,
the identification of the frequency content characterising
the shock unsteadiness and the application of the
$LIM2$ condition have thus been used to propose a
signal processing procedure to, first,

separate the footprint of the shock motion from the turbulent
pressure fluctuations, and then, isolate its intermittent content.
The analysis of the obtained signal reveals that the intermittency of the wall-pressure
is mostly characterised by temporal scales that are different from those that
dominate the overall unsteadiness of the shock system according to Fourier analysis.

The proposed signal processing can represent a preliminary practice
for further investigations into the causality and the features of
the low-frequency unsteadiness of stbli systems,
by using for example conditional statistics based on the system intermittency.
In this kind of studies, an essential step is in fact the proper estimation of the
correlation between the motion of the reflected shock and the characteristic time scales
of the incoming turbulent boundary layer, the separated region and the re-attachment region:
isolating the low-frequency intermittent part of the associated time series

could help recognising relations that may be overshadowed
by the combination of the different components of the complete signals.

Motivated by the many studies in the literature attesting a relationship
between the unsteadiness of the shock and the pulsation of the separation bubble
[Piponniau et al., 2009],
we decided to study by means of the proposed procedure the
dynamics of the recirculation region behind the reflected shock.

The time signal recording the evolution of the separation bubble extent
presents a skewed probability distribution with
strong fluctuations between states with large and minimal separation respectively.
The analysis of the intermittency through wavelet transform
shows that also the motion of the bubble has strong, sparse energy bursts
at large time scales and reveals a connection between such events and
those highlighted by the $LIM2$ measure of the shock foot wall-pressure,
although it was not possible to estimate the statistics of such events
because of the insufficient number of occurrences.
In order to qualify and quantify the connection between the two
signals, wavelet cross spectra and wavelet coherence [Torrence & Webster, 1999]
have been considered for the signals processed by the proposed procedure.
The study of the coherence among the times and time scales
confirms that the intermittent contributions of the breathing bubble motion
and the shock foot wall-pressure are highly correlated during the
most important intermittent events. However, the corresponding cross wavelet phase angle
shows that there is not a simple linear relationship between the two signals,
which would have suggested a causality relation similarly to what recently considered by
[Sasaki et al., 2021].

In the end, our results confirm the hypothesis about the link between the
low-frequency shock unsteadiness and the breathing mode of the separation bubble,
even if our investigation points out the complexity of the breathing motion
and indicates that other phenomena may contribute to the dynamics of the overall system.

[Acknowledgements]We acknowledge CINECA for the availability of high performance computing resources and support during the pre-production stage of Marconi100 cluster. We also acknowledge Prof. S. Piponniau and Prof. P. Dupont
for having provided us with the experimental data used in the validation section.

[Supplementary data]A movie of the time evolution of a contour of the Mach number on a vertical slice, during the intermittent event around $t \, U_\infty / \delta_0 \approx 1160$ is available.
The white line indicates the locus of the points with null streamline velocity.

[Funding]This research received no specific grant from any funding agency, commercial or not-for-profit sectors.

[Declaration of interests]The authors report no conflict of interest.

[Data availability statement]The full data set of the DNS simulations is on the order of thousands of gigabytes. By contacting the authors, a smaller subset can be made available.

[Author ORCID]M. Bernardini, <https://orcid.org/
0000-0001-5975-3734>; G. Della Posta, <https://orcid.org/
0000-0001-5516-9338>; F. Salvadore, <https://orcid.org/
0000-0002-1829-3388>; E. Martelli, <https://orcid.org/

natexlab#1#1#1#1 #1#1 #1#1#1#1#1#1
#1#1#1#1 #1#1 #1#1

[Abe et al., 1999]
Abe, Y., Kiya, M. & Mochizuki, O. 1999
Low-frequency unsteadiness in turbulent elliptic wakes.  JSME Int.
J. Ser. B  42 (3),  349–354.

[Adams, 2000]
Adams, N. A. 2000  Direct simulation of the turbulent
boundary layer along a compression ramp at M = 3 and Retheta = 1685.
J. Fluid Mech.  420,  47–83.

[Adamson & Messiter, 1980]
Adamson, T. C. & Messiter, A. F. 1980  Analysis of
two-dimensional interactions between shock waves and boundary layers.
Annu. Rev. Fluid Mech.  12 (1),  103–138.

[Agostini et al., 2015]
Agostini, L., Larchevêque, L. & Dupont, P. 2015
Mechanism of shock unsteadiness in separated shock/boundary-layer
interactions.  Physics of Fluids  27 (12),  126103.

[Andreopoulos & Muck, 1987]
Andreopoulos, J. & Muck, K.C. 1987  Some new aspects
of the shock wave/boundary layer interaction in compression ramp flows.
J. Fluid Mech.  180,  405–428.

[Aubard et al., 2013]
Aubard, G., Gloerfelt, X. & Robinet, J.C. 2013
Large-eddy simulation of broadband unsteadiness in a shock/boundary-layer
interaction.  AIAA J.  51,  2395–2409.

[Bakulu et al., 2021]
Bakulu, F., Lehnasch, G., Jaunet, V., Goncalves da
Silva, E. & Girard, S. 2021  Jet resonance in truncated
ideally contoured nozzles.  J. Fluid Mech.  919,  A32.

[Bell & Lepicovsky, 1995]
Bell, W. & Lepicovsky, J. 1995 Analysis of flows between
rotor blades with the wavelet transform.  In 33rd ASME,
p. 18.

[Beresh et al., 2002]
Beresh, S.J., Clemens, N.T. & Dolling, D.S. 2002
Relationship between upstream turbulent boundary layer velocity
fluctuations and separation shock unsteadiness.  AIAA J.  40,

[Bermejo-Moreno et al., 2014]
Bermejo-Moreno, I., Campo, L., Larsson, J., Bodart,
J., Helmer, D. & Eaton, J. K. 2014  Confinement
effects in shock wave/turbulent boundary layer interactions through
wall-modelled large-eddy simulations.  J. Fluid Mech.  758,

[Bernardini et al., 2021]
Bernardini, M., Modesti, D., Salvadore, F. &
Pirozzoli, S. 2021  STREAmS: A high-fidelity accelerated
solver for direct numerical simulation of compressible turbulent flows.
Comput. Phys. Commun.  263,  107906.

[Camussi & Bogey, 2021]
Camussi, R. & Bogey, C. 2021  Intermittent statistics
of the 0-mode pressure fluctuations in the near field of Mach 0.9 circular
jets at low and high Reynolds numbers.  Theor. Comput. Fluid Dyn.
35 (2),  229–247.

[Camussi et al., 2017]
Camussi, R., Di Marco, A. & Castelain, T. 2017
Statistical analysis of the hydrodynamic pressure in the near field of
compressible jets.  Int. J. Heat Fluid Flow  64,  1–9.

[Clemens & Narayanaswamy, 2009]
Clemens, N. T. & Narayanaswamy, V. 2009 Shock/turbulent
boundary layer interactions: review of recent work on sources of
unsteadiness.  In 39th AIAA Fluid Dyn. Conf.,  p. 3710.

[Clemens & Narayanaswamy, 2014]
Clemens, N. T. & Narayanaswamy, V. 2014
Low-frequency unsteadiness of shock wave/turbulent boundary layer
interactions.  Annu. Rev. Fluid Mech.  46 (1),  469–492.

[Consolini & De Michelis, 2005]
Consolini, G. & De Michelis, P. 2005  Local
intermittency measure analysis of AE index: The directly driven and
unloading component.  Geophys. Res. Lett.  32 (5).

[Daubechies, 1992]
Daubechies, I. 1992 Ten Lectures on Wavelets.

[Delery, 1985]
Delery, J. M. 1985  Shock wave/turbulent boundary layer
interaction and its control.  Prog. Aerosp. Sci.  22 (4),

[Deshpande & Poggie, 2021]
Deshpande, A. S. & Poggie, J. 2021  Large-scale
unsteadiness in a compression ramp flow confined by sidewalls.  Phys.
Rev. Fluids  6,  024610.

[Dolling, 2001]
Dolling, D. S. 2001  Fifty years of
shock-wave/boundary-layer interaction research: What next?  AIAA J.
39 (8),  1517–1531.

[Dolling & Or, 1985]
Dolling, D. S. & Or, C. T. 1985  Unsteadiness of the
shock wave structure in attached and separated compression ramp flows.
Exp. Fluids  3 (1),  24–32.

[Dupont et al., 2006]
Dupont, P., Haddad, C. & Debiève, J. F. 2006
Space and time organization in a shock-induced separated boundary layer.
J. Fluid Mech.  559,  255–277.

[Dupont et al., 2008]
Dupont, P., Piponniau, S., Sidorenko, A. & Debiève,
J. F. 2008  Investigation by particle image velocimetry
measurements of oblique shock reflection with separation.  AIAA J.
46 (6),  1365–1370.

[Elena & Lacharme, 1988]
Elena, M. & Lacharme, J-P 1988  Experimental study of
a supersonic turbulent boundary layer using a laser doppler anemometer.
J. Mec. Theor. Appl.  7 (2),  175–190.

[Erengil & Dolling, 1991]
Erengil, M.E. & Dolling, D.S. 1991  Correlation of
separation shock motion with pressure fluctuations in the incoming boundary
layer.  AIAA J  29,  1868–1877.

[Farge, 1992]
Farge, M. 1992  Wavelet transforms and their applications
to turbulence.  Annu. Rev. Fluid Mech.  24 (1),  395–458.

[Farge et al., 1999]
Farge, M., Kevlahan, NKR, Perrier, V. & Schneider,
K. 1999  Turbulence analysis, modelling and computing using
wavelets.  Wavelets Phys.  pp. 117–200.

[Gaitonde, 2015]
Gaitonde, D. V. 2015  Progress in shock wave/boundary layer
interactions.  Prog. Aerosp. Sci.  72,  80–99, celebrating
60 Years of the Air Force Office of Scientific Research (AFOSR): A Review of
Hypersonic Aerothermodynamics.

[Ganapathisubramani et al., 2007]
Ganapathisubramani, B., Clemens, N. T. & Dolling, D. S.
2007  Effects of upstream boundary layer on the unsteadiness of
shock-induced separation.  JFM  585,  369–394.

[Ganapathisubramani et al., 2009]
Ganapathisubramani, B., Clemens, N. T. & Dolling, D. S.
2009  Low-frequency dynamics of shock-induced separation in a
compression ramp interaction.  J. Fluid Mech.  636,

[Green, 1970]
Green, J.E. 1970  Interactions between shock waves and
turbulent boundary layers.  Prog. Aerosp. Sci.  11,

[Grilli et al., 2012]
Grilli, M., Schmid, P. J., Hickel, S. & Adams, N. A.
2012  Analysis of unsteady behaviour in shockwave turbulent boundary
layer interaction.  J. Fluid Mech.  700,  16–28.

[Grinsted et al., 2004]
Grinsted, A., Moore, J. C. & Jevrejeva, S. 2004
Application of the cross wavelet transform and wavelet coherence to
geophysical time series.  Nonlinear Proc. Geoph.  11 (5/6),

[Grizzi & Camussi, 2012]
Grizzi, S. & Camussi, R. 2012  Wavelet analysis of
near-field pressure fluctuations generated by a subsonic jet.  J. Fluid
Mech.  698,  93–124.

[Guo et al., 2020]
Guo, S., Wu, Y. & Liang, H. 2020  Strong
interactions of incident shock wave with boundary layer along compression
corner.  Chin. J. Aeronaut.  33 (12),  3149–3157.

[Huang et al., 1998]
Huang, N. E., Shen, Z., Long, S. R., Wu, M. C.,
Shih, H. H., Zheng, Q., Yen, N.-C., Tung, C. C. &
Liu, H. H. 1998  The empirical mode decomposition and the
hilbert spectrum for nonlinear and non-stationary time series analysis.
Proc. R. Soc. A: Math. Phys. Eng. Sci.  454 (1971),

[Humble et al., 2007]
Humble, R., Elsinga, G., Scarano, F. & van Oudheusden,
B. 2007 Investigation of the instantaneous 3d flow organization of a
swtbli using tomographic PIV.  In 37th AIAA Fluid Dyn. Conf.
Ex.,  p. 4112.

[Humble et al., 2009]
Humble, R.A., Elsinga, G.E., Scarano, F. & van
Oudheusden, B.W. 2009  Unsteady aspects of an incident shock
wave/turbulent boundary layer interaction.  J. Fluid Mech.  635,

[Jiang & Shu, 1996]
Jiang, G.S. & Shu, C.W. 1996  Efficient
implementation of weighted ENO schemes.  J. Comput. Phys.
126,  202.

[Kaiser & Hudgins, 1994]
Kaiser, G. & Hudgins, L. H. 1994 A friendly guide to
wavelets, ,  vol. 300.  Springer.

[Kempf et al., 2012]
Kempf, A.M., Wysocki, S. & Pettit, M. 2012  An
efficient, parallel low-storage implementation of Klein's turbulence
generator for LES and DNS.  Comput. Fluids  60,

[Klein et al., 2003]
Klein, M., Sadiki, A. & Janicka, J. 2003  A
digital filter based generation of inflow data for spatially developing
direct numerical or large eddy simulations.  J. Comput. Phys.
186 (2),  652–665.

[Kouchi et al., 2016]
Kouchi, T., Yamaguchi, S., Koike, S., Nakajima, T.,
Sato, M., Kanda, H. & Yanase, S. 2016  Wavelet
analysis of transonic buffet on a two-dimensional airfoil with vortex
generators.  Exp. Fluids  57 (11),  1–15.

[Lewalle, 1994]
Lewalle, J. 1994  Wavelet transforms of some equations of
fluid mechanics.  Acta Mech.  104 (1),  1–25.

[Liandrat & Moret-Bailly, 1990]
Liandrat, J. & Moret-Bailly, F. 1990  The wavelet
transform - some applications to fluid dynamics and turbulence.  Eur. J.
Mech. B Fluids  9 (1),  1–19.

[Loginov et al., 2006]
Loginov, M. S., Adams, N. A. & Zheltovodov, A. A.
2006  Large-eddy simulation of shock-wave/turbulent-boundary-layer
interaction.  J. Fluid Mech.  565,  135–169.

[Lund, T.S. et al., 1998]
Lund, T.S., Wu, X. & Squires, K.D. 1998
Generation of inflow data for spatially-developing boundary layer
simulations.  J. Comput. Phys.  140 (2),  233–258.

[Mallat, 1989]
Mallat, S.G. 1989  A theory for multiresolution signal
decomposition: the wavelet representation.  IEEE PAMI  11 (7),

[Mallat, 1999]
Mallat, S. 1999 A wavelet tour of signal processing.

[Maraun et al., 2007]
Maraun, D., Kurths, J. & Holschneider, M. 2007
Nonstationary gaussian processes in wavelet domain: Synthesis,
estimation, and significance testing.  Phys. Rev. E  75,

[Martelli et al., 2017]
Martelli, E., Ciottoli, P.P., Bernardini, M., Nasuti,
F. & Valorani, M. 2017  Detached eddy simulation of shock
unsteadiness in an over-expanded planar nozzle.  AIAA J.
55 (6),  2016–2028.

[Martelli et al., 2020]
Martelli, E., Saccoccio, L., Ciottoli, P. P., Tinney,
C. E., Baars, W. J. & Bernardini, M. 2020  Flow
dynamics and wall-pressure signatures in a high-Reynolds-number
overexpanded nozzle with free shock separation.  J. Fluid Mech.
895,  A29.

[Meneveau, 1991]
Meneveau, C. 1991  Analysis of turbulence in the
orthonormal wavelet representation.  J. Fluid Mech.  232,

[Meyers et al., 1993]
Meyers, S. D., Kelly, B. G. & O'Brien, J. J. 1993
An introduction to wavelet analysis in oceanography and meteorology: With
application to the dispersion of yanai waves.  MWR  121 (10),
2858 – 2866.

[Morgan et al., 2013]
Morgan, B., Duraisamy, K., Nguyen, N., Kawai, S. &
Lele, S.K. 2013  Flow physics and RANS modelling of oblique
shock/turbulent boundary layer interaction.  J. Fluid Mech.
729,  231–284.

[Nichols et al., 2017]
Nichols, J. W., Larsson, J., Bernardini, M. &
Pirozzoli, S. 2017  Stability and modal analysis of
shock/boundary layer interactions.  Theor. Comput. Fluid Dyn.
31 (1),  33–50.

[Pasquariello et al., 2014]
Pasquariello, V., Grilli, M., Hickel, S. & Adams,
N. A. 2014  Large-eddy simulation of passive
shock-wave/boundary-layer interaction control.  Int. J. Heat Fluid Flow
49,  116–127, tSFP8.

[Pasquariello et al., 2017]
Pasquariello, V., Hickel, S. & Adams, N. A. 2017
Unsteady effects of strong shock-wave/boundary-layer interaction at high
Reynolds number.  J. Fluid Mech.  823,  617–657.

[Piponniau et al., 2009]
Piponniau, S., Dussauge, J. P., Debiève, J. F. &
Dupont, P. 2009  A simple model for low-frequency unsteadiness
in shock-induced separation.  J. Fluid Mech.  629,

[Pirozzoli, 2010]
Pirozzoli, S. 2010  Generalized conservative approximations
of split convective derivative operators.  Journal of Computational
Physics  229 (19),  7180–7190.

[Pirozzoli, 2011]
Pirozzoli, S. 2011  Numerical methods for high-speed
flows.  Annu. Rev. Fluid Mech.  43,  163–194.

[Pirozzoli & Bernardini, 2011]
Pirozzoli, S. & Bernardini, M. 2011a
Direct numerical simulation database for impinging shock wave/turbulent
boundary-layer interaction.  AIAA J.  49 (6),  1307–1312.

[Pirozzoli & Bernardini, 2011]
Pirozzoli, S. & Bernardini, M. 2011b
Supersonic turbulent boundary layers - DNS database.

[Pirozzoli & Grasso, 2006]
Pirozzoli, S. & Grasso, F. 2006  Direct numerical
simulation of impinging shock wave/turbulent boundary layer interaction at
M=2.25.  Phys. Fluids  18 (6),  065113.

[Pirozzoli et al., 2010]
Pirozzoli, S., Larsson, J., Nichols, J.W., Bernardini,
M., Morgan, B. & Lele, S.K. 2010  Analysis of unsteady
effects in shock/boundary layer interactions. Proceedings of the 2010 CTR
summer program.  Stanford University.

[Poggie & Porter, 2019]
Poggie, J. & Porter, K. M. 2019  Flow structure and
unsteadiness in a highly confined shock-wave–boundary-layer interaction.
Phys. Rev. Fluids  4,  024602.

[Poggie & Smits, 1997]
Poggie, J. & Smits, A. J. 1997  Wavelet analysis of
wall-pressure fluctuations in a supersonic blunt-fin flow.  AIAA J.
35 (10),  1597–1603.

[Poinsot & Lele, 1992]
Poinsot, T.S. & Lele, S.K. 1992  Boundary conditions
for direct simulations of compressible viscous flows.  J. Comput. Phys.
101,  104–129.

[Priebe et al., 2016]
Priebe, S., Tu, J. H., Rowley, C. W. & Martín,
M. P. 2016  Low-frequency dynamics in a shock-induced separated
flow.  J. Fluid Mech.  807,  441–477.

[Ringuette et al., 2009]
Ringuette, M. J., Bookey, P., Wyckham, C. & Smits,
A. J. 2009  Experimental study of a mach 3 compression ramp
interaction at retheta = 2400.  AIAA J.  47 (2),

[Sasaki et al., 2021]
Sasaki, K., Barros, D. C., Cavalieri, A. V. G. &
Larchevêque, L. 2021  Causality in the shock wave/turbulent
boundary layer interaction.  Phys. Rev. Fluids  6,  064609.

[Torrence & Compo, 1998]
Torrence, C. & Compo, G. P. 1998  A practical guide
to wavelet analysis.  Bull. Amer. Meteor.  79 (1),  61 –

[Torrence & Webster, 1999]
Torrence, C. & Webster, P. J. 1999  Interdecadal
changes in the enso–monsoon system.  J. Clim.  12 (8),
2679 – 2690.

[Touber & Sandham, 2009]
Touber, E. & Sandham, N. D. 2009  Large-eddy
simulation of low-frequency unsteadiness in a turbulent shock-induced
separation bubble.  Theor. Comput. Fluid Dyn.  23 (2),