B Bugeat

J.-Ch Robinet

J.-C Chassaing

P Sagaut

Low-frequency resolvent analysis of the laminar oblique shock wave/boundary layer interaction

Resolvent analysis is used to study the low-frequency behaviour of the laminar oblique shock wave/boundary layer interaction (SWBLI). It is shown that the computed optimal gain, which can be seen as a transfer function of the system, follows a first-order low-pass filter equation, recovering the results of Touber & Sandham (J. Fluid Mech., vol. 671, 2011, pp. 417-465). This behaviour is understood as proceeding from the excitation of a single stable, steady global mode whose damping rate sets the time scale of the filter. Different Mach and Reynolds numbers are studied, covering different recirculation lengths L. This damping rate is found to scale as 1/L, leading to a constant Strouhal number St L as observed in the literature. It is associated with a breathing motion of the recirculation bubble. This analysis furthermore supports the idea that the low-frequency dynamics of the SWBLI is a forced dynamics, in which background perturbations continuously excite the flow. The investigation is then carried out for three-dimensional perturbations for which two regimes are identified. At low wavenumbers of the order of L, a modal mechanism similar to that of two-dimensional perturbations is found and exhibits larger values of the optimal gain. At larger wavenumbers, of the order of the boundary layer thickness, the growth of streaks, which results from a non-modal mechanism, is detected. No interaction with the recirculation region is observed. Based on these results, the potential prevalence of three-dimensional effects in the low-frequency dynamics of the SWBLI is discussed.

Introduction

An accurate understanding of the dynamics of the shock wave/boundary layer interaction (SWBLI) is needed in many aerospace and aeronautical applications to predict, for example, flows around transonic airfoils, supersonic air intakes or deflected control surfaces of vehicles at transonic or supersonic speed. These interactions can lead to an increase of drag, to separation and to loss of performance. Moreover, shock wave/turbulent boundary layer interaction (SWTBLI) generally produces a low-frequency unsteadiness of the shock system. This unsteadiness can, for example, modify the thermal load or induce fatigue of the structure. It is well known that the flow separation may be at the origin of this unsteadiness in any regime, from the incompressible [START_REF] Weiss | Unsteady behavior of a pressure-induced turbulent separation bubble[END_REF][START_REF] Mohammed-Taifour | Unsteadiness in a large turbulent separation bubble[END_REF][START_REF] Le Floc'h | Measurements of pressure and velocity fluctuations in a family of turbulent separation bubbles[END_REF] to the hypersonic [START_REF] Gaitonde | Progress in shock wave/boundary layer interactions[END_REF][START_REF] Priebe | Turbulence in a hypersonic compression ramp flow[END_REF], taking a particular form for compressible flows when the separation is induced by a shock wave. A considerable amount of work has been carried out to investigate the steady and unsteady aspects (see [START_REF] Gaitonde | Progress in shock wave/boundary layer interactions[END_REF] for a review). Experimental research of the SWBLI started in the mid-1940s with the work of [START_REF] Ackeret | Investigation of compression shocks and boundary layers in gases moving at high speed[END_REF] and has shown a continuing interest since then. At this time, most of the experiments only measure steady aerodynamic quantities such as pressure distribution, skin friction and heat transfer rates. Detailed investigations of the phenomenon and its dependence on flow and boundary layer parameters have been already published by [START_REF] Adamson | Analysis of two-dimensional interactions between shockwaves and boundary layers[END_REF], [START_REF] Délery | Shock-wave boundary layer interactions[END_REF] and [START_REF] Stanewsky | Shock boundary layer interaction[END_REF].

It was only at the end of the 90s that the study of the dynamics of a SWTBLI was more systematically undertaken [START_REF] Dolling | Fifty years of shock-wave/boundary-layer interaction research: what next?[END_REF][START_REF] Knight | Advances in CFD prediction of shockwave turbulent boundary layer interactions[END_REF]. It has been shown that the dynamics is rich, spanning several frequency decades. Two main length scales are at play: the vorticity thickness of the shear layer coming from the separated region and linked to the disturbances developing along the shear layer and, on a more global scale, the separation length L that is related to lower frequencies and is associated with the dynamics of the separation bubble [START_REF] Dupont | Unsteadiness in compressible separated flows[END_REF][START_REF] Gaitonde | Progress in shock wave/boundary layer interactions[END_REF]. Three main frequency ranges can be identified. First, high frequencies, characterised by Strouhal numbers based on L typically in the range of St L ≥ 1, are associated with the transitional or turbulent boundary layer dynamics [START_REF] Bonne | Analysis of the two-dimensional dynamics of a mach 1.6 shock wave/transitional boundary layer interaction using a RANS based resolvent approach[END_REF]. Secondly, a vortex shedding resulting from the shear layer dynamics usually takes place at mid-frequency (0.1 ≤ St L ≤ 1). This can be observed either in incompressible recirculation bubble [START_REF] Weiss | Unsteady behavior of a pressure-induced turbulent separation bubble[END_REF][START_REF] Mohammed-Taifour | Unsteadiness in a large turbulent separation bubble[END_REF] or in the SWBLI [START_REF] Thomas | On the mechanism of unsteady shock oscillation in shock wave/turbulent boundary layer interactions[END_REF][START_REF] Agostini | Zones of influence and shock motion in a shock/boundary-layer interaction[END_REF]. Finally, a third dynamics is observed in the low-frequency range (St L < 0.1) [START_REF] Clemens | Low-frequency unsteadiness of shock wave/turbulent boundary layer interactions[END_REF]. This is the frequency domain we will mainly focus on from now on.

The low-frequency range has been intensely studied since the early 2000s in order to unveil its physical origin [START_REF] Dussauge | Unsteadiness in shock wave boundary layer interactions with separation[END_REF]. [START_REF] Dupont | Space and time organisation of a shock wave/turbulent boundary layer interaction[END_REF]; [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF] have shown that, for an interaction strong enough, low-frequency unsteadiness is observed around St L 0.03 and is characterised by oscillations of the separated shock [START_REF] Dussauge | Unsteadiness in shock wave boundary layer interactions with separation[END_REF][START_REF] Humble | Unsteady aspects of an incident shock wave/turbulent boundary layer interaction[END_REF]. Scenarios attempting to model this low-frequency dynamics have been devised. The first of them is related to the interaction of large structures upstream of the boundary layer with the shock leading to low-frequency response of the shock [START_REF] Ganapathisubramani | Effects of upstream boundary layer on the unsteadiness of shock-induced separation[END_REF]). However, [START_REF] Dussauge | Shock/boundary-layer interactions: possible sources of unsteadiness[END_REF] have shown that the influence of downstream conditions, especially in the recirculation zone, is more significant than the upstream conditions with respect to low frequencies when the interaction is strong. [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF] proposed a model based on a mass balance of the system 'shear layer -separated zone' where coherent structures in the shear layer feed the recirculation zone, which increases up to a critical size beyond which it empties, causing a breathing of the bubble and consequently the movement of separation of the shock [START_REF] Agostini | Mechanism of shock unsteadiness in separated shock/boundary-layer interactions[END_REF]. [START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shock-wave/boundary-layer interactions[END_REF] furthermore characterised the system using a combination of theoretical and numerical model. They showed that the separated shock foot acts as a low-pass filter with respect to white noise. This result is in agreement with the linear interaction analysis [START_REF] Ribner | Convection of pattern of vorticity through a shock wave[END_REF][START_REF] Robinet | Critical interaction of a shock wave with an acoustic wave[END_REF] of the shock response to a harmonic infinitesimal perturbation.

The three-dimensional (3-D) dynamics resulting from the development of 3-D disturbances, independently of aspects strictly related to the effects of turbulence, has been studied more recently. Mainly three types of phenomena have been identified. When the upstream boundary layer is turbulent and the SWBLI weak (incipient separated zone), streaks from the upstream boundary layer and associated with longitudinal vortices can interact with the shock system [START_REF] Ganapathisubramani | Effects of upstream boundary layer on the unsteadiness of shock-induced separation[END_REF][START_REF] Ganapathisubramani | Low-frequency dynamics of shock-induced separation in a compression ramp interaction[END_REF][START_REF] Di Renzo | Crossflow effects on shock wave/turbulent boundary layer interactions[END_REF]. However, when the SWBLI is strong enough to generate a significant separated zone, the latter can induce a streamline curvature of the boundary layer and initiate a centrifugal Görtler-like instability [START_REF] Pasquariello | Unsteady effects of strong shock-wave/boundary-layer interaction at high Reynolds number[END_REF][START_REF] Zhuang | Görtler-like vortices in an impinging shock wave/turbulent boundary layer interaction flow[END_REF]. When the SWBLI is sufficiently strong, the separated zone becomes three-dimensional and the emergence of a global non-oscillating instability can be observed. This characteristic is very general and can be observed from subsonic [START_REF] Theofilis | On the origin of unsteadiness and three-dimensionality in a laminar separation bubble[END_REF][START_REF] Rodriguez | Structural changes of laminar separation bubbles induced by global linear instability[END_REF][START_REF] Rodríguez | Self-excited primary and secondary instability of laminar separation bubbles[END_REF] to supersonic regime [START_REF] Robinet | Bifurcations in shock wave / laminar boundary layer interaction: global instability approach[END_REF][START_REF] Hildebrand | Simulation and stability analysis of oblique shock-wave/boundary-layer interactions at mach 5.92[END_REF]. The role that these 3-D effects may play in the low-frequency unsteadiness remains, however, unclear.

From a numerical point of view, a substantial body of work has been done over the years, including direct numerical simulations and large eddy simulations. The works of Touber & Sandham (2009a,b) and [START_REF] Pirozzoli | Direct numerical simulation database for impinging shock wave/turbulent boundary-layer interaction[END_REF] made it possible for the first time to characterise the dynamics with reasonable statistical convergence, showing a very good agreement with Institut Universitaire des Systèmes Thermiques Industriels experiments [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF]. In particular, these numerical studies also showed the broadband nature of the low-frequency dynamics. [START_REF] Priebe | Direct numerical simulation of a reflected-shockwave/turbulent-boundary-layer interaction[END_REF], for a reflected shock, and [START_REF] Wu | Analysis of shock motion in shockwave and turbulent boundary layer interaction using direct numerical simulation data[END_REF], for a compression ramp, detailed the low-frequency dynamics and recovered the scaling proposed in the literature. They also observed a large amplification of the turbulence through the SWTBLI. [START_REF] Priebe | Low-frequency unsteadiness in shock wave/turbulent boundary layer interaction[END_REF] for the compression ramp have studied the physical mechanism that drives the shock motion. In their simulations, the flow undergoes low-frequency changes of topology in the interaction region, including the breaking up of the recirculation bubble and the shedding of vortical structures. In addition, the growth of energetic turbulent structures in the shear layer was found to be modulated at low frequency. This could imply a modulation of the shear layer entrainment rate which is consistent with the scenario of [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF]. These authors thus suggest that the low-frequency dynamics is related to the dynamics of the separated zone. [START_REF] Aubard | Large-eddy simulation of broadband unsteadiness in a shock/boundary-layer interaction[END_REF] and [START_REF] Priebe | Low-frequency unsteadiness in shock wave/turbulent boundary layer interaction[END_REF] then showed that the dynamics of the separated zone occurs at medium frequency and is linked to Kelvin-Helmholtz instabilities in the shear layer and that the low-frequency dynamics is only a modulation of the latter and corresponds to the breathing of the separated zone.

Resolvent analysis, or resolvent-based modelling, has been intensively used in the fluid dynamics community in recent years. This approach is based on the singular value decomposition (SVD) of the resolvent operator that provides an optimal orthonormal basis for an external forcing (input) of the system and its associated linear response (output). The optimality is here defined as the way for a given energy input, at a given frequency, to force the system in order to trigger the largest possible energy growth. This is therefore particularly suited to the study of convective instabilities in noise amplifier flows. In such flows, ranges of frequency of external perturbations are amplified in space and time through linear mechanisms, while others are damped [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF]. Resolvent analysis was successfully applied, for example, to boundary layer [START_REF] Åkervik | Global two-dimensional stability measures of the flat plate boundary-layer flow[END_REF][START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF] and jet [START_REF] Nichols | Global modes and transient response of a cold supersonic jet[END_REF][START_REF] Garnaud | The preferred mode of incompressible jets: linear frequency response analysis[END_REF] flows. Because non-modal effects are taken into account [START_REF] Schmid | Nonmodal stability theory[END_REF], algebraic instabilities like streaks or non-modal mechanisms such as the Orr mechanism can furthermore be modelled via this approach [START_REF] Monokrousos | Global three-dimensional optimal disturbances in the blasius boundary-layer flow using time-steppers[END_REF]. Resolvent analysis has also been intensively used in the modelling of turbulent flow since [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF] noticed it could be used to identify coherent, large scale structures. While turbulence is a highly nonlinear phenomenon, the (linear) resolvent analysis yet captures the growth of these structures. The interpretation, and consequently the relevance of the resolvent analysis, is that nonlinear forcing terms, whose knowledge is not ultimately required here, provide background perturbations that grow into energetic coherent structures through linear mechanisms. Several authors later discussed the potential of the resolvent approach in reduced-order modelling of turbulent flows [START_REF] Semeraro | Stochastic and harmonic optimal forcing in subsonic jets[END_REF], noticing, in particular, the link between resolvent and spectral-proper orthogonal decomposition modes [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF].

Resolvent analysis has been recently used to study the dynamics of the SWBLI. [START_REF] Sartor | Unsteadiness in transonic shock-wave/boundary-layer interactions: experimental investigation and global stability analysis[END_REF] used a turbulent mean flow of a transonic SWBLI over a bump and compared it with their experimental results. They obtained similar features at low and medium frequencies, supporting the idea of a forced dynamics as their mean flow was found to be globally stable. [START_REF] Bonne | Analysis of the two-dimensional dynamics of a mach 1.6 shock wave/transitional boundary layer interaction using a RANS based resolvent approach[END_REF] studied a transitional case, in which they identified three ranges of frequencies and proposed a scenario for the low-frequency dynamics. In our paper, we focus exclusively on the low-frequency dynamics of the SWBLI, that we analyse by means of the resolvent analysis. Unlike the two papers cited above, a laminar base flow will be considered, for which similar features of the low-frequency dynamics can be observed [START_REF] Yao | Direct simulation of transional shock boundary-layer interactions at mach 2[END_REF][START_REF] Boin | 3d steady and unsteady bifurcations in a shock-wave/laminar boundary layer interaction: a numerical study[END_REF][START_REF] Sandham | Transitional shock-wave/boundary-layer interactions in hypersonic flow[END_REF][START_REF] Sansica | Forced response of a laminar shock-induced separation bubble[END_REF], 2016). Our objective is first to understand to what extent a resolvent-based approach can model the behaviour of the SWBLI for this range of frequency, where the dynamics does not result from convective or global instabilities. This will also provide insights into the physics of the system that will be discussed. In particular, by conducting our analysis for different Mach and Reynolds numbers, we aim at characterising the scales involved in this dynamics. Besides, despite considering a laminar flow, this analysis could, to some extent, support the idea that the turbulent SWBLI has a forced dynamics at low frequency. Indeed, if the intrinsic forcing emanating from nonlinear interactions in the turbulent case reasonably projects onto the optimal forcing found in resolvent analysis, it is then possible to assume that a similar response will be found. All of the above will be carried out considering 2-D perturbations only, similarly to the two previously cited papers. The second objective is to explore the resolvent analysis of the SWBLI for 3-D perturbations for which, to our best knowledge, no work has been published. By 3-D perturbations, we mean that Fourier modes will be considered in the homogeneous spanwise direction, characterised by their wavenumber β. The aim is to investigate the role that streaks or any 3-D effects may play in the low-frequency dynamics of the system.

The paper is structured as follows. The governing equations and the theoretical and numerical approaches are described in § 2. The methodology to generate the set of base flows under study is then detailed in § 3. In § 4, the resolvent analysis is carried out for 2-D perturbations on different base flows. A link with global stability analysis will be proposed, leading to a model of the optimal and sub-optimal gain obtained from the SVD of the resolvent operator. The scaling of key quantities will be shown, before discussing these results with respect to the literature. In § 5, 3-D perturbations will be considered in the resolvent analysis. Only one base flow will be studied, for which the newly obtained features of the dynamics will be classified into two regimes depending on their wavenumbers. This section ends with a discussion on the 3-D effects. A general conclusion can be found in § 6.

Theoretical and numerical approach

Governing equations

The compressible nonlinear Navier-Stokes equations are considered for the conservative variables q = (ρ, ρu, ρE) T , where ρ, u and E are the density, the velocity vector and the specific total energy, respectively. The components of the velocity vector are noted u, v and w in the streamwise direction x, the wall-normal direction y and the spanwise direction z, respectively. The variables T, p, μ, κ and γ , standing for temperature, pressure, dynamic and kinematic viscosity, thermal conductivity and heat capacity ratio, are also introduced. I is the identity tensor. All variables are non-dimensional using their value at infinity, except the pressure that is normalised by ρ ∞ u 2 ∞ . The reference length scale and time scale are δ * 0 and δ * 0 /u ∞ , where δ * 0 is the compressible displacement thickness at a chosen streamwise location x, defined as

δ * 0 = ∞ 0 {1 -ρu/(ρ ∞ u ∞ )} dy.
More details on the choice of δ * 0 for our problem will be given in § 3.2. The governing equations read ∂ρ ∂t

+ ∇ • (ρu) = 0, (2.1a) ∂ ∂t (ρu) + ∇ • ρu ⊗ u + pI - 1 Re τ = 0, (2.1b) ∂ ∂t (ρE) + ∇ • (ρE + p) u - 1 Re τ u - κ PrRe(γ -1)M 2 ∇T = 0. (2.1c)
The fluid is assumed to be a perfect gas whose equation of states and total energy read

p = 1 γ M 2 ∞ ρT, E = p ρ(γ -1) + 1 2 u • u. (2.2a,b)
Furthermore assuming the fluid to be Newtonian, the viscous stress tensor verifies

τ = 2μ S - 1 3 Tr (S) I , with S = 1 2 ∇ ⊗ u + (∇ ⊗ u) T (2.3)
Sutherland's law is used for the dynamic viscosity

μ(T) = T 3/2 1 + T s /T ∞ T + T s /T ∞ , (2.4)
where T s = 110.4 K and T ∞ = 288 K. The thermal conductivity κ is also assumed to follow Sutherland's law [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]. Three non-dimensional numbers have been introduced that control the flow: Re, M and Pr, the Reynolds, Mach and Prandtl numbers, respectively, defined as

Re = ρ ∞ u ∞ δ * 0 μ ∞ , M = u ∞ c ∞ , Pr = μ ∞ c p κ ∞ , (2.5a-c)
where c ∞ is the speed of sound at infinity and c p the heat capacity. While Pr is set to Pr = 0.72 in order to model the fluid as air, the influence of Re and M will be studied.

Note that a geometrical parameter φ will also be introduced. It is defined as the angle between the streamwise direction and the shock wave impinging the boundary layer. Finally, (2.1) will be recast in the dynamical system form as

∂q ∂t = N (q), (2.6)
where N is the nonlinear differential operator of the Navier-Stokes equations.

Global stability

In this paper, the SWBLI is studied via two types of linear analysis. The linearised Navier-Stokes equations, considered in their semi-discrete form, are written as

∂q ∂t = Jq , (2.7)
where J and q are the Jacobian matrix and the state-vector perturbations around the base flow q, respectively. The base flow q(x, y) corresponds to a steady 2-D laminar solution of the nonlinear equations (2.1). A Fourier transform in both time and space (in the spanwise direction z) can then be performed, allowing the perturbations to be studied as the following global modes

q (x, y, z, t) = q(x, y) exp(i(βz -ωt)) + c.c., (2.8)
with ω the angular frequency and β and the spanwise wavenumber. Considering ω ∈ C and β ∈ R, the global stability (Theofilis 2011) of the system can be studied by solving the eigenvalue problem iω q = J q.

(2.9)

If an eigenvalue is found such that the imaginary part of ω is positive, then the flow is globally unstable. In this paper, only globally stable flows will be considered: if not externally sustained, any linear perturbations will eventually decay.

Resolvent analysis

Even though globally stable, flows such as boundary layer or mixing layer flow can exhibit convective instabilities [START_REF] Huerre | Local and global instabilities in spatially developing flows[END_REF]). In such flows, termed noise amplifier flows, studying the linear forced dynamics is relevant to understand and model their behaviour. This can be achieved using resolvent analysis [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: a linearized approach[END_REF], which has been shown to provide a good model for convective instabilities over ranges of frequency where large energy growths take place [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF]. Besides, this approach gives the location and structure of the optimal harmonic forcing of the flow, i.e. the forcing field that triggers the largest energy growth. More generally, resolvent analysis provides a transfer function of the system, characterising the response to a forcing for each frequency [START_REF] Schmid | Nonmodal stability theory[END_REF].

The linear forced dynamics can be studied by introducing a harmonic forcing term f in (2.9) as

-iωq ˆ = Jq ˆ + f ˆ,
(2.10) with ω ∈ R, which now represents the forcing frequency that one wants to study.

Introducing the resolvent matrix R = (-ωI -J) -1 , the relationship between any forcing and its linear response is then given by q = R f .

(2.11)

One can now look for the forcing field that triggers the most energetic response by computing the optimal gain σ defined as

σ 2 (ω, β) = max f q 2 E f 2 F .
(2.12)

Here, two norms must be introduced to measure the energy of the forcing and the response fields. The norm introduced by [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part I)[END_REF] is used for the response while the Euclidean norm is used for the forcing. More details can be found in [START_REF] Bugeat | 3D global optimal forcing and response of the supersonic boundary layer[END_REF]. It can be shown that the optimal gain σ corresponds to the largest singular value of the resolvent matrix modified by the norm matrices [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF]. The associated optimal forcing and response are then the first left and right singular vectors, respectively. From (2.12), σ can be seen as a transfer function as it measures the energy of the response with respect to that of an input forcing. In the following, the largest (or first) singular value will be called the optimal gain and will be referred to as σ or σ 1 without any differences. The subsequent singular values σ 2 , σ 3 , . . . such that σ 1 > σ 2 > σ 3 > • • • are called sub-optimal gains and are associated with the subsequent singular vectors which form an orthonormal basis.

Numerical methods

Base flow computation

In order to get a 2-D steady base flow, the nonlinear Navier-Stokes equations (2.1) are solved using a finite volume computational fluid dynamics solver. Spatial discretisation of convective fluxes is performed using the advection upstream splitting method scheme [START_REF] Liou | A sequel to AUSM: AUSM+[END_REF]) associated with a fifth-order monotone upstream-centered schemes for conservation laws extrapolation [START_REF] Boin | 3d steady and unsteady bifurcations in a shock-wave/laminar boundary layer interaction: a numerical study[END_REF]. Viscous fluxes at cell interfaces are obtained by a second-order centred finite difference scheme. The unsteady equations are marched in time until a steady state is reached. An implicit dual time stepping method with local time step is used [START_REF] Jameson | Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings[END_REF]. More details about this solver and its validation can be found in [START_REF] Boin | 3d steady and unsteady bifurcations in a shock-wave/laminar boundary layer interaction: a numerical study[END_REF]. A rectangular numerical domain is defined as in figure 1. The streamwise length L x is set such that the distance between the impinging shock x imp and the outlet of the domain x out is 70δ * 0 . The height L y is set such that the reflecting shock leaves the domain at the outlet boundary. This value depends on the physical case considered. In any case, the base flow eventually used in the stability and resolvent calculations is cropped to L y = 50δ * 0 . A Cartesian mesh is set with a geometrical progression from the wall. The number of points is 600 × 270 in the x and y directions, respectively. The independence of the results from the aforementioned numerical parameters is checked in Appendix A. The following boundary conditions are used. At the inlet, oblique shock conditions are imposed: for y < y shock , a parallel flow is set using the Mach number at infinity while the corresponding downstream jump conditions are used for y > y shock . The latter condition is also set at top boundary of the domain. Note that y shock is calculated after choosing the angle of the shock φ and the location of the impinging shock x imp . At the outlet, variables are extrapolated from the inside of the domain. At the bottom, an adiabatic flat plate is considered with no-slip conditions imposed on the velocity.

I n c i d e n t s h o c k E x p a n s i o n f a n F o o t s h o c k R e a t t a c h m e n t s h o c k L x L y x in x imp x out y shock φ M ∞ Recirculation L Figure 1.
Sketch of the numerical domain and typical structure of the oblique shock wave/boundary layer (OSWBLI).

Linear solvers

The global stability and resolvent analyses are based on the linearised Navier-Stokes equations, which revolve around the computation of the Jacobian matrix J. In order to compute it, the discretised Navier-Stokes equation are linearised via a finite difference approximation [START_REF] Mettot | Computation of eigenvalue sensitivity to base flow modifications in a discrete framework: application to open-loop control[END_REF]. This is achieved by using the discrete residual N (2.6) of the nonlinear Navier-Stokes used in the base flow solver, ensuring numerical consistency between the base flow calculation and the linear solvers. The adaptation of this framework proposed by [START_REF] Bugeat | 3D global optimal forcing and response of the supersonic boundary layer[END_REF] for 3-D perturbations in base flows featuring one homogeneous direction is used. Exhaustive details of the resolvent calculation can also be found in [START_REF] Bugeat | 3D global optimal forcing and response of the supersonic boundary layer[END_REF] where, in particular, the boundary conditions, Chu's norm matrix for conservative variables and the Krylov subspace algorithm used to compute the singular values are fully described. A sponge zone of length 10δ * 0 is also implemented at the top and downstream boundaries of the domain using the approach described by [START_REF] Agarwal | Calculation of sound propagation in nonuniform flows: suppression of instability waves[END_REF]. Finally, the PETSc [START_REF] Balay | PETSc users manual[END_REF], SLEPc [START_REF] Hernandez | SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems[END_REF] and MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF] open libraries are called to compute the large scale linear systems and eigenvalues problems featuring sparse matrices.

Base flow

Validation

In order to validate the base flow calculation, the flow conditions are matched to those of the experimental investigation of [START_REF] Degrez | The interaction of an oblique shock wave with a laminar boundary layer revisited. An experimental and numerical study[END_REF]. The free-stream Mach number is M = 2.15 and the angle of the incident shock measured from the horizontal axis is φ = 30.8 • (corresponding to a flow deflection angle of 3.81 • ). The Reynolds number, based on the streamwise length from the leading edge at which the shock impinges the flat plate, is Re x = 0.96 × 10 5 . Comparisons with the experimental results reported by [START_REF] Degrez | The interaction of an oblique shock wave with a laminar boundary layer revisited. An experimental and numerical study[END_REF] are displayed in figure 2. This shows that our simulation is in good agreement with experimental data, both in terms of wall pressure measurements and location of the separation bubble.

Set of base flows studied

A set of base flows is computed in order to carry out stability and resolvent analyses for different separation lengths L. Different Mach numbers are considered, ranging from 2.00 (incipient separation) to 2.35 (large separation length compared with the boundary layer length scale). First, the Reynolds number Re, based on the compressible boundary layer thickness δ * 0 at which the shock impinges the plate, is kept constant at Re = 1100 using the following protocol. For each Mach number, a boundary layer base flow is first computed without the impinging shock. At each x-location, δ * is then calculated by integrating the quantity 1 -(ρu)/(ρ ∞ u ∞ ) in the normal direction. This allows us to detect the x-location where Re = 1100 and to use it as an input in a second base flow computation with a shock impinging at this point. Thus, the reference length scale δ * 0 corresponds to the boundary layer thickness at the location of the impinging shock but prior to the shock interaction. In other words, it is not the effective thickness resulting from the shock interaction. The velocity field obtained at M = 2.20 and Re = 1100, which will be used as the reference case, is shown in figure 3. By increasing the Mach number, the pressure gradient felt by the boundary layer is increased. The separation length L/δ * 0 is readily computed by detecting the separation and reattachment points along the plate, which are the locations where ∂u/∂y at the wall becomes negative and positive again, respectively. The separation length L/δ * 0 is found very sensitive to M as its value increases over more than one order of magnitude for the Mach numbers considered in table 1. In order to assess the role of viscous effects in the low-frequency receptivity of the SWBLI, different Reynolds numbers are also considered, ranging from 600 to 2200. For one given Mach number, this corresponds to different x-locations at which the shock impinges. Note that because δ * ∼ √ x, this location varies over a wide range of values, e.g. by a factor 10 at M = 2.20 for which the maximum value of Re considered is 1900. As shown in table 1, increasing Re leads to increasing the recirculation bubble as there is less streamwise momentum near the wall to resist the pressure gradient generated by the shock. Finally, note that the angle of the shock is kept constant at φ = 30.8 • for all the computations presented in this paper.

Linear dynamics of 2-D perturbation at low frequency

Resolvent analysis: optimal forcing and response

Results

The optimal gain is computed for different frequencies at M = 2.20 and Re = 1100 (figure 4). At low frequency (St L < 10 -1 ), the optimal gain monotonically decreases. The resolvent modes at St L = 10 -4 are shown in figure 5. The optimal forcing is spread over a large portion of the domain, admitting maxima near the separation point and along the right part of the recirculation bubble. The optimal response is concentrated above the bubble, following its shape. Moreover, even though a drop of optimal gain is observed at low frequency, the optimal forcing and response vectors are found independent of the forcing frequency: the same velocity and pressure fields are observed at any low Strouhal numbers. Noticeable differences start to appear for frequencies above St L 2 × 10 -2 . An interpretation will be proposed in § 4.1.3. At higher frequencies, a maximum of optimal gain is detected at St L 2 (figure 4). It is associated with a convective instability that starts developing in the recirculation region and continues to grow further downstream in the boundary layer (figure 6). This is reminiscent of the first mode instability of the supersonic boundary layer, which is the compressible counterpart of the Tollmien-Schlichting instability in incompressible boundary layers [START_REF] Mack | Boundary-layer linear stability theory[END_REF]. In order to efficiently trigger this instability, the optimal forcing is located upstream. Its tilted structure shows the concurrent action of the non-modal Orr mechanism, as usually observed in boundary layer flows (Ehrenstein & Gallaire 2008;[START_REF] Bugeat | 3D global optimal forcing and response of the supersonic boundary layer[END_REF].

Global stability analysis

A global stability computation is performed for the set of base flows described in § 3.2. The eigenvalue spectrum at M = 2.20 and Re = 1100 is presented in figure 7. Each eigenvalue has a negative growth rate ω i , meaning that the system is globally stable. This has also been verified for all the base flows considered in this study. The two least stable modes are found to be steady (S1 and S2, St L = 0) while the third and fourth modes are unsteady (T1 and T2, St L / = 0). Subsequent modes at higher frequency feature even larger damping rate. The velocity field of the mode S1 is mostly localised around the recirculation bubble, following its shape (figure 8a). The mode T1 shares common characteristics but features a phase opposition between the upstream and downstream region of the bubble (figure 8b).

The striking outcome of this global stability analysis is that the least stable global mode S1 is very similar to the optimal response found in the resolvent analysis in the previous section (figure 5b), which was observed for any low frequencies. This means that, at low frequency, the optimal response results from the excitation of the least stable global mode. These results suggest that the receptivity at low frequency is a modal phenomenon, as opposed to the non-modal mechanisms usually observed in convective instabilities [START_REF] Cossu | Global measures of local convective instability[END_REF]. Note that, because the system is stable, a continuous forcing is required to excite the flow around its base state.

It may appear somewhat unsettling that the receptivity, which pertains to the unsteady behaviour of the system, seems to be driven by a steady mode. This observation is analysed in the next section.

A low-pass filter model

To clarify how a steady mode could play a role in the low-frequency dynamics of the system, a model based on global stability analysis is proposed to describe the behaviour of the optimal gain. Two sources can cause the singular values of the resolvent operator to increase [START_REF] Schmid | Nonmodal stability theory[END_REF]). The first one appears when the forcing frequency is close to the eigenvalues of the Jacobian matrix (as found by a stability analysis). The second source is related to the non-normality of the eigenvectors of the Jacobian matrix (i.e. the global modes). It is now well known that the linear superposition of stable modes can lead to energy growth [START_REF] Farrell | Optimal excitation of perturbations in viscous shear flow[END_REF], termed non-modal growth [START_REF] Schmid | Stability and Transition in Shear Flows[END_REF]. This idea is translated, in the resolvent framework, into the concept of pseudo-resonance through which increases of optimal gain can occur away from the frequencies of the eigenmodes (see for example [START_REF] Bucci | Roughness-induced transition by quasi-resonance of a varicose global mode[END_REF].

From the insight obtained from the previous stability analysis ( § 4.1.2), we would like to test if the optimal gain proceeds from purely modal effects. Thus, non-modal phenomena are here discarded from the model. Furthermore, if one and only one global mode drives the optimal gain, then the singular value of the resolvent should vary as σ ∼ 1/d, where d is the distance in the complex plane between the forcing frequency St L and the eigenvalue ω/2π = (ω r + iω i )/2π of this global mode. This distance is readily obtained as

d = (St L -ω r /2π) 2 + (ω i /2π) 2 . (4.1)
Because of the resemblance previously observed between the optimal response and the global mode S1, this mode is a natural candidate to test the purely modal receptivity model. Its eigenvalue is then plugged into (4.1). Since this mode is steady (ω r = 0), the low-frequency model of the optimal gain can eventually be recast as

σ (St L ) = σ 0 1 + St L ω (S1) i /2π 2 , ( 4.2)
where σ 0 = lim St L →0 σ (St L ) is the value of the optimal gain when St L goes zero. Here, σ 0 is not predicted a priori by the model, but is obtained by one computation of the optimal gain at zero frequency (or any small frequency where the optimal gain experiences a plateau). But since the absolute value of the optimal gain ultimately carries little significance [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF], this value is not essential to this model. The aim of this model is rather to detect the drop of optimal gain in order to analyse low-frequency receptivity of the system. It can be recognised that (4.2) is that of a first-order low-pass filter. The optimal gain being by definition analogous to a transfer function (see (2.12)), describing its behaviour as a filter is consistent. In (4.2), the cutoff frequency is given by the damping rate ω (S1) i of the global mode S1, which is the least stable global mode. Even though this mode is steady, a harmonic forcing of non-zero frequency is still able to excite it. When the forcing frequency St L goes to zero, the gain does not depend on the frequency as the distance in the complex plane d tends to the constant value |ω (S1) i /2π|. However, increasing St L until the order of ω (S1) i /2π and above affects the value of the optimal gain that then decreases as the forcing frequency is pushed away from the eigenvalue of the global mode it excites. Thus, the time scale appearing in this low-pass filter model is not set by the frequency of the global mode (which would be zero here), but by its damping rate.

Test of the model for different Mach and Reynolds numbers

The optimal gain model proposed in (4.2) is tested for different Mach and Reynolds numbers. Setting Re = 1100, three Mach numbers are considered: M = 2.10, M = 2.20 and M = 2.30, with separation lengths ranging from 22δ * 0 to 88δ * 0 (table 1). Figure 9 shows that the low-pass filter model accurately detects the drop of optimal gain that appears between St L = 10 -3 and 10 -2 . In other words, the cutoff frequency is correctly predicted by the damping rate of the global mode S1. This is observed for every Mach numbers. Disagreements appear above St L 2 × 10 -2 , where the low-pass filter underestimates the actual gain. At these frequencies, the structure of the optimal response starts to change. This means that the global mode S1 does not drive the dynamics of the flow alone any longer, but that other modes are getting involved. Other Reynolds numbers, below and above the previous value Re = 1100, are considered while keeping M = 2.2. Below Re = 1100, the same agreement as previously described is observed (figure 10). As the Reynolds number increases (Re = 1600), discrepancies start to appear for smaller frequencies than St L = 10 -2 but the location of the drop of optimal gain is still correctly captured. At Re = 1900, the model becomes irrelevant. The growing mismatch between the low-pass filter model and the optimal gain as Re is increased can be caused by the action of other modes, previously damped by viscous effects. As a result, the model based on the excitation of the global mode S1 alone cannot predict the optimal gain anymore, even though this mode does not directly depend on viscous effect (in the next section, it will be shown that the global mode S1 is independent of the viscous length scale δ * 0 ).

σ 2.0 ×10 2 0 M = 2.10 4.0 ×10 2 6.0×10 2 St L 10 -1 10 -2 10 -3 σ 5.0×10 2 5.0 ×10 2 1.0 ×10 3 1.5 ×10 3 2.0 ×10 3 2.5 ×10 3 0 0 M = 2.20 1.0×10 3 1.5×10 3 St L 10 -1 10 -2 10 -3 M = 2.30 St L 10 -1 10 -2 10 -3 (a) (b) ( c)

Scaling

The low-pass filter model in (4.2) contains the damping rate ω (S1) i

of the mode S1. This parameter sets the unique time scale of the model. Its scaling relative to a reference is presented in figure 11 as a function of the non-dimensional length of the recirculation bubble L/δ * 0 . Note that each point is associated with a different pair of Reynolds and Mach numbers (see table 1). It is found that ω

(S1) i δ * 0 /u ∞ ∼ (L/δ * 0 ) -1 , meaning that ω (S1) i ∼ u ∞ /L.
= fL/u ∞ = (ω (S1) δ * 0 /2πu ∞ )(L/δ * 0 )
, leading to values around St L 1 × 10 -2 . While of the same order of magnitude as the data of the literature, the cutoff frequency is slightly underestimated. This motivates to push the resolvent analysis a step further by looking at the sub-optimal modes.

Resolvent analysis: sub-optimal forcing and response

Results

So far, only the optimal gain σ 1 and its associated forcing and response modes have been considered. Subsequent singular values of the resolvent operator σ 2 and σ 3 , with σ 1 > σ 2 > σ 3 (i.e. sub-optimal gains), are now computed. Figure 12 compares these three gains as a function of the forcing frequency. At low frequency, a factor of approximately three is observed between σ 1 and σ 2 , and approximately two between σ 2 and σ 3 . A drop of sub-optimal gains occurs at higher frequencies that those noted for σ 1 . In particular, the decrease of σ 2 is observed at Strouhal numbers between 10 -2 and 10 -1 . These values lie in the range of the experimental data of the low-frequency unsteadiness more closely than those of σ 1 . Our attention will then be focused on the behaviour of σ 2 and its associated modes at low frequencies. The forcing mode follows the shape of the recirculation bubble but is mostly concentrated around the separation point as shown in figure 13(a). Unlike the optimal forcing in figure 5(a), the phase varies in the normal direction. This is also the case of the associated response (figure 13b), whose phase varies in the streamwise direction. Finally, as previously observed for σ 1 , the sub-optimal modes associated with σ 2 do not depend on St L in the low-frequency domain. This could suggest that the sub-optimal gain is also based on a modal excitation. In the next section, this possibility is discussed by reconsidering the first-order low-pass filter model.

ω i (S1) δ 0 * /u ∞ ∼ (L/δ 0 * ) -1 L/δ 0 * 10 -

Sub-optimal gain: another low-pass filter

The response mode associated with σ 2 (figure 13b) shares some features with the unsteady global mode T1 (figure 8b) computed in § 4.1.2. Besides, the low-frequency behaviour of σ 2 as a function of St L is similar to that of σ 1 , the drop of gain being shifted to higher frequencies. It is then tempting to build a model of σ 2 based on the resonance of the global mode T1, derived from the distance in the complex plane between the real forcing frequency and its eigenvalue. However, because this mode is unsteady, a resonance model would lead to a peak of gain at the frequency of the mode T1 and a decrease of the gain as the frequency goes to zero. These two features are not observed in σ 2 (figure 12). Such a resonance model is then doomed to fail. Instead, (4.2), which is that of a first-order low-pass filter, satisfies the observation of σ 2 at low frequency. Therefore, we suggest using this equation to model σ 2 by plugging in the damping rate of the mode T1 instead of the mode S1 given its resemblance with the response mode associated with σ 2 . This equation becomes

σ 2 (St L ) = σ 2,0 1 + St L ω (T1) i /2π 2 . (4.3)
As previously mentioned, this approach is different from an approach based on the pure resonance of the global mode T1; instead, we only assume that the damping rate of this mode is the relevant time scale at play. The 'purely modal dynamics' interpretation elaborated for σ 1 (see § 4.1.3) then no longer holds. Other modes, either by normal or non-normal effect, are understood to also contribute to the low-frequency dynamics. This is furthermore expected as the response mode of σ 2 is actually not strictly identical to the mode T1. This model is tested for the same set of base flow as previously investigated. At Re = 1100, the location of the drop of the sub-optimal gain σ 2 , observed shortly after St L = 10 -2 , is captured for the three Mach numbers considered as seen in figure 14. This means that the cutoff frequency of the low-pass filter, based on the damping rate of the global mode T1, is indeed relevant. The slope of the drop is also predicted by the first-order filter until St L = 10 -1 , except for M = 2.10 for which the behaviour of σ 2 changes at lower frequencies. Now setting M = 2.20, the agreement of the cutoff frequency is also correctly captured for Reynolds numbers lower than or equal to Re = 1100 (figure 15). At higher Reynolds number, the model over-estimates this frequency while still reasonably characterising the behaviour at Re = 1600. For Re = 1900, the model completely fails (figure 15c,d), as already noted for σ 1 .

Scaling

The damping rate ω (T1) i of the global mode T1 has been shown to capture reasonably well the low-frequency drop of σ 2 . Its scaling is considered by plotting its value against the non-dimensional separation length for the set of Reynolds and Mach numbers studied (figure 16). Similarly to § 4.1.5, the damping rate is found to scale with u ∞ /L since ω

(T1) i δ * 0 /u ∞ ∼ (L/δ * 0 ) -1
. This again means that the Strouhal number St L associated with the cutoff frequency is constant for any Mach and Reynolds numbers. Its value is found around St L 4 × 10 -2 , in good agreement with the range of experimental data of the low-frequency unsteadiness. Before discussing the implication of this result, the bubble dynamics is studied in the next subsection.

Bubble dynamics

The first and second resolvent modes are now compared from the perspective of the bubble dynamics. In order to observe its motion, one response mode at one chosen frequency is linearly added to the base flow. The frequency considered here is St L = 10 -3 but its value is actually not important since, as previously mentioned, the resolvent modes do not depend on St L at low frequency. Snapshots at six different time steps of a periodic cycle, i.e. six different phases, are shown in figure 17. In each case, the distorted recirculation bubble is shown. It corresponds to the location in space where the total streamwise velocity is zero, calculated as ū + Aû = 0. The scalar A is an arbitrary amplitude chosen as A = 0.3 to ease visualisation. The undistorted bubble of the base flow (ū = 0) is 1), as a function of the separation length. Here, δ * 0 is the boundary layer thickness at Re = 1100 rather than that at each Reynolds number. also plotted as a reference. The optimal resolvent mode, associated with σ 1 , produces a breathing motion: the bubble periodically contracts and expands relative to the basic state (figure 17a,c,e,g,i,k). This dynamics is similar to the bubble motion found in the model of the low-frequency unsteadiness developed by [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF]. The dynamics produced by the sub-optimal mode, shown in figure 17(b,d,f,h,j,l), is different.

(c) ( d) (e) ( f ) (g) ( h) (i) ( j) (k) ( l)
Figure 17. Bubble dynamics at M = 2.20 and Re = 1100 resulting from the addition of the optimal (a,c,e,g,i,k) and sub-optimal (b,d, f,h,j,l) low-frequency resolvent modes to the base flow. The blue dashed line is the distorted recirculation bubble corresponding to the location where ū + Aû = 0 (using A = 0.3). Its evolution at 6 different phases of a periodic cycle is shown from top to bottom. The red line is the recirculation bubble of the base flow (ū = 0).

The bubble is now periodically displaced upstream and downstream from the reference state. No displacement is observed in the wall-normal direction.

Discussion

In their review paper, [START_REF] Clemens | Low-frequency unsteadiness of shock wave/turbulent boundary layer interactions[END_REF] ended their discussion on whether the SWBLI behaved as a forced dynamical system at low frequency. Our resolvent analysis supports the idea that, indeed, the low-frequency behaviour results from a forced dynamics, as opposed to an intrinsic dynamics such as that of a vortex street behind cylinder. From the resemblance between their dynamic mode decomposition and stable global stability modes, [START_REF] Nichols | Stability and modal analysis of shock/boundary layer interactions[END_REF] suggested that the system could then be seen as a damped oscillator. Our study brings some new lights to this idea by showing that one stable global mode drives the resolvent analysis at low frequency. The fact that the optimal gain follows a first-order low-pass filter equation is furthermore consistent with the work of [START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shock-wave/boundary-layer interactions[END_REF] who were the first to model the low-frequency dynamics of the SWBLI as a filter of this kind. From this perspective, the dynamics requires a range of low frequencies that force the system and are filtered by the transfer function which, in the resolvent framework, is related to the singular values of the resolvent operator. Therefore, while carried out in the laminar regime, our analysis provides insights that the SWBLI in the turbulent regime could possess a forced dynamics at low frequency, in which fluctuations over a wide range of scales would force the flow (McKeon & Sharma 2010) and sustain its low-frequency behaviour. This would explain why low-frequency unsteadiness is observed in the turbulent regime rather that the laminar regime unless the latter is externally forced as in [START_REF] Sansica | Forced response of a laminar shock-induced separation bubble[END_REF]. In the present study, because resolvent analysis assumes by nature the presence of a forcing field, the low-frequency behaviour of the system can be qualitatively recovered even though the flow is laminar.

For both the optimal and sub-optimal gain, the damping rate of one stable global mode provides the time scale of the low-pass filter (the cutoff frequency). These frequencies scale as (L sep /u ∞ ) -1 , which was expected as these modes are related to the bubble dynamic. This results in a constant cutoff Strouhal number St L for any Mach and Reynolds numbers. This scaling is indeed observed both experimentally and numerically for different configurations of SWBLI, with frequencies reported from St L 3 × 10 -2 to St L 6 × 10 -2 [START_REF] Dussauge | Unsteadiness in shock wave boundary layer interactions with separation[END_REF]. With values presently found around St (2) L 4 × 10 -2 for the optimal and sub-optimal gains, respectively, St

(2) L compares more favourably with the literature. However, the dynamics of the bubble of the optimal mode is, on the one hand, related to a breathing motion that strongly resembles that at play in the low-frequency unsteadiness of the SWBLI [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Clemens | Low-frequency unsteadiness of shock wave/turbulent boundary layer interactions[END_REF]. It is also strikingly similar to that observed experimentally in a turbulent incompressible recirculation bubble [START_REF] Weiss | Unsteady behavior of a pressure-induced turbulent separation bubble[END_REF], also occurring at Strouhal numbers of the order of 10 -2 . On the other hand, the sub-optimal mode that produces a very different type of motion, which is expected given the orthogonality of the resolvent modes. The resulting bubble motion is similar to that obtained by [START_REF] Nichols | Stability and modal analysis of shock/boundary layer interactions[END_REF] in their global stability analysis; this is indeed the same unsteady stable mode. However, no conclusive element allows us to link it to the low-frequency behaviour reported in the literature. This leads us to hypothesize that the optimal (rather than sub-optimal) gain and modes, associated with the bubble breathing, indeed describe the low-frequency dynamics of the SWBLI. The fact that St

(1) L is underestimated in the present study may be due to the laminar nature of the SWBLI understudied while results in the literature usually deal with turbulent SWBLI. The shape of the bubbles being different between laminar and turbulent regimes this could account for the quantitative mismatch. Finally, while the range of values of separation length studied in the present work is large enough to characterise the scaling of the low-frequency time scale, the Mach number is restricted to a small range around M = 2.20. Similar analyses carried out at much different Mach numbers are then still needed to assess the universality of these findings.

The larger the gain separation σ 1 /σ 2 between the optimal and sub-optimal gain, the more likely the optimal response obtained through resolvent analysis is to model the system. More rigorously, this also depends on the projection of the effective forcing of the system (whether it is from external perturbations or intrinsic, turbulent nonlinear interactions) onto the forcing resolvent modes [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. Here, σ 1 /σ 2 3 is found at low frequency. This value not being significantly greater than 1, if we assume that the optimal response drives the system (which we do based on the previous discussions on the bubble dynamics), then the actual forcing can be assumed to project reasonably well onto the optimal forcing mode. This gives us some insight into the actual forcing in the flow. The optimal forcing is indeed maximum around the separation point and along the separation line on the downstream part of the bubble (figure 5). Thus, some forcing is required in these regions in order to maintain the low-frequency bubble dynamics. The limit of the resolvent analysis is eventually reached when seeking for the origin of the forcing. Because of its linear nature, this framework does not allow us to discuss how the low-frequency structures can be produced by nonlinear interactions. Thus, both boundary layer structures convected to the interaction region [START_REF] Ganapathisubramani | Effects of upstream boundary layer on the unsteadiness of shock-induced separation[END_REF]) exciting the separation point and upstream travelling mechanisms [START_REF] Sansica | Instability and low-frequency unsteadiness in a shock-induced laminar separation bubble[END_REF][START_REF] Bonne | Analysis of the two-dimensional dynamics of a mach 1.6 shock wave/transitional boundary layer interaction using a RANS based resolvent approach[END_REF] forcing the downstream part of the recirculation bubble could be valid candidates to sustain the low-frequency dynamics of the SWBLI.

Linear dynamics of 3-D perturbation at low frequency

Resolvent analysis

The low-frequency dynamics of 3-D perturbations, characterised by their non-zero spanwise wavenumber β, is studied in this section. The analysis is here limited to the base flow calculated at M = 2.20 and Re = 1100. The influence of β on the optimal gain at St L = 10 -4 is shown in figure 18. A first peak is observed for low wavenumbers around β = 0.25. The corresponding wavelength is λ z = 25δ * 0 , which is of the order of the recirculation length (approximately half of it, see table 1). The associated optimal forcing and response can be visualised by plotting the isosurface of a selected amplitude in the 3-D physical space (figure 19). The optimal response is located in the recirculation region and features a phase shift in the streamwise direction. A strong resemblance can be noted with the 3-D global mode found by [START_REF] Robinet | Bifurcations in shock wave / laminar boundary layer interaction: global instability approach[END_REF]. The possibility of a modal resonance to model the optimal gain, as developed previously for 2-D perturbation, will be explored in the next section. The optimal forcing field is made of elongated structures that lay upstream from the recirculation bubble until the separation point. The energy density dE(x) of the response can be defined at each streamwise location as dE(x) = 0.5ρ(uu * + vv * + ww * ) dy. An analogous quantity can be defined for the forcing field. The profiles of these energy densities confirm that the response is localised in the separation region as it experiences a rapid decay downstream from it (figure 20). The forcing energy is shown to spread over a larger portion of the domain and reaches a peak in the vicinity of the separation point.

In addition of the local maximum at β = 0.25, the optimal gain features a peak of larger magnitude around β = 2. The corresponding wavelength is λ z = 3δ * 0 and is now of the order of the boundary layer thickness. Elongated streamwise vortices located upstream of the recirculation bubble are found in the optimal forcing field (figure 21). A further inspection of the wall-normal organisation of its components (available in Appendix B, figure 28) shows that they correspond to streamwise vortices. Downstream of the recirculation bubble, streaks of streamwise velocities are observed in the boundary layer, for which profiles reveal that u v and u w (figure 28). These features indicate the action of the lift-up mechanism that is ubiquitous in shear flows at low frequency and non-zero spanwise wavenumbers [START_REF] Landahl | A note on an algebraic instability of inviscid parallel shear flows[END_REF]. Besides, these results are very similar to the 3-D optimal response in a supersonic boundary layer without impinging shock [START_REF] Bugeat | 3D global optimal forcing and response of the supersonic boundary layer[END_REF]. The energy profiles reveal the convective nature of this instability as perturbations keep growing until the downstream end of the domain (figure 20). Most of the energy is contained downstream from the recirculation area whereas most of the forcing energy is localised upstream from it. These results hence show that the growth of streaks does not interact significantly with the recirculation region.

Link with global stability

Eigenvalue spectra of the global stability problem are given in figure 22 for the wavenumbers β = 0.25 and β = 2 that correspond to maxima of the optimal gain. In both cases, the flow is found to be globally stable. Similarly to the analysis developed for 2-D perturbations in § 4, the relation between the optimal gain behaviour at low frequency and the global stability features is explored for these two spanwise wavenumbers.

Low β: excitation of a steady global mode

At β = 0.25, the least stable global mode is steady (mode S1) as shown in figure 22(a). Its structure is shown in figure 23. It strongly resembles that of the optimal response previously observed. For 2-D perturbations, these observations led us to formulate the hypothesis of a purely modal dynamics at low frequency, driven by a single mode. The first-order low-pass filter model (4.2) that results from this assumption and that is based on the damping rate of the 3-D global mode S1 is tested in figure 24. Once again, the model correctly predicts the cutoff frequency of the optimal gain. The slope of the gain drop is also reasonably captured until St L 10 -2 , where additional modes participate in the forced dynamics. The low-frequency dynamics of 3-D perturbations at low wavenumbers thus behaves the same way as that of 2-D perturbations, which was expected as it is the limit when β → 0. Finally, another peak of optimal gain can be observed at higher frequency in figure 24. This maximum is associated with the convective instabilities that are known to exhibit larger growth rates for β / = 0 in supersonic flows [START_REF] Mack | Boundary-layer linear stability theory[END_REF], resulting in larger values of the optimal gain [START_REF] Bugeat | 3D global optimal forcing and response of the supersonic boundary layer[END_REF].

Larger β: pseudo-resonance

The lift-up mechanism leading to the formation of streaks is a linear algebraic instability, as theorised by [START_REF] Ellingsen | Stability of linear flow[END_REF]. As such, linear stability theories based on a modal approach cannot capture it. Conversely, the resolvent analysis is particularly suited for their calculation as non-modal effects are taken into account [START_REF] Schmid | Nonmodal stability theory[END_REF]. We now suggest showing the failure of the first-order low-pass filter model (based the excitation of a single mode) as a way to point out the non-modal character, or pseudo-resonance, of the optimal gain at wavenumbers of the order of β = 2. The least damped mode of the spectrum is the unsteady modes T1, followed by the steady mode S1 (figure 22b). Their structure is shown in figure 25. They share similarities with the optimal response already presented in figure 21(b) as their growth takes place in the downstream region of the domain. Nevertheless, noticeable differences remain, which suggests that the optimal response does not proceed from the excitation of a unique global mode. Two low-pass filters model based on the modes S1 and T1 are tested in figure 26. Both models underestimate the cutoff frequency and underpredict the optimal gain, confirming the non-modal nature of the optimal response at these wavenumbers.

Discussion

Depending on their spanwise wavenumber, the 3-D optimal response at low frequency can result from two different mechanisms, underpinned by a modal or non-modal nature. At low wavenumbers, whose length scales are of the order of L sep , optimal perturbations interact with the recirculation bubble. Their energy is localised in this region similarly to what was observed for 2-D perturbations. Additionally, the low-pass filter model based on the excitation of a unique steady global mode reasonably describes the low-frequency behaviour of the optimal gain σ . Low wavenumbers hence behave in continuity to the 2-D results (β = 0) presented in § 4. Noticeably larger values of σ are, however, reached for 3-D perturbations. This raises the question of whether the dynamics of the system at low frequency and low wavenumber is dominated by 2-D or 3-D perturbations. On the one hand, larger values of the gain indeed indicates that more energetic 3-D perturbations are expected in the flow, assuming the same forcing energy and the same projection of this forcing onto the resolvent modes for each wavenumber. These two assumptions are, at this point, impossible to verify. For example, if the forcing generated by nonlinear interactions in a turbulent flow is preferentially two-dimensional, then one may still expect the 2-D response to be dominant, even though it is associated with a lower gain value. Thus, no definitive conclusion can be drawn. On the other hand, the cutoff frequency based on the 3-D global response is St L = 6 × 10 -3 . This underestimates the usual value of St L = 3 × 10 -2 even more than the 2-D perturbations, for which a value of St L = 1 × 10 -2 was found. This is in favour of the prevalence of 2-D perturbations at low frequency, even though a perfect quantitative agreement between our laminar analysis and the experimental data on turbulent flow is not expected. This is also backed up by previous work which restricted their analysis to 2-D perturbations [START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shock-wave/boundary-layer interactions[END_REF][START_REF] Sansica | Forced response of a laminar shock-induced separation bubble[END_REF] The possibility that the low-frequency behaviour of the SWBLI may be a 3-D phenomenon remains a promising lead to explore further given some arguments found in our resolvent analysis. It could be interesting to investigate whether the optimal spanwise wavenumber would then scale with the recirculation bubble. A space-Fourier transform along the spanwise direction of the time-Fourier modes from experimental or numerical data would be required to explore this phenomenon. This, however, may not be manageable experimentally if the spanwise width of the facility is not wide enough compare with recirculation length scale. Recent studies have besides pointed out the influence of the sidewalls in SWBLI experiments [START_REF] Rabey | Two-dimensional unsteadiness map of oblique shock wave/boundary layer interaction with sidewalls[END_REF][START_REF] Xiang | Corner effects for oblique shock wave/turbulent boundary layer interactions in rectangular channels[END_REF]. From a computational point of view, the same challenge of ensuring sufficiently wide domain is expected, additionally coupled with the costly computational time necessary to statistically converge the low-frequency time-Fourier transform. Such a study, if achievable, could help to characterise the 3-D organisation of the forcing field and eventually understand its role in the low-frequency dynamics.

While low wavenumbers exhibit a modal and localised (around the bubble) dynamics, optimal responses at larger wavenumbers have been shown to ignore the recirculation region as they are driven by a convective instability growing downstream in the boundary layer, as a result of non-modal effects. These findings are consistent with those of [START_REF] Dwivedi | Transient growth analysis of oblique shock-wave/boundary-layer interactions at mach 5.92[END_REF], who recently carried out a transient growth analysis in a SWBLI at M = 5.92. Such analysis also embeds the non-modal effect and can in fact be seen as the temporal counterpart of the resolvent analysis [START_REF] Sipp | Dynamics and control of global instabilities in open-flows: a linearized approach[END_REF]. The authors detected streaky structures evolving in time and convected in space, similar to those obtained in the present study. Our resolvent analysis showed that streaks are associated with larger optimal gain values than the low wavenumber perturbations. This is not surprising as it results from a convective instability which efficiently extract energy from the base flow. The role played by the streaks in the low-frequency dynamics (St L 3 × 10 -2 ) is not obvious. Contrary to low wavenumber perturbations, streaks are not directly involved in this dynamics since they do not interact with the recirculation region and are expected to scale with δ * 0 rather than L. At first, it then seems that they could only be of interest regarding the transition to turbulence of the system. However, the large energy growth they trigger could quickly lead to nonlinear interactions that may play a role in the forcing of the bubble through upstream feedback [START_REF] Bonne | Analysis of the two-dimensional dynamics of a mach 1.6 shock wave/transitional boundary layer interaction using a RANS based resolvent approach[END_REF]. Studying the development of streaks in a turbulent SWBLI, seen as coherent structures, would then be of interest. Following the work of [START_REF] Sartor | Unsteadiness in transonic shock-wave/boundary-layer interactions: experimental investigation and global stability analysis[END_REF], this could be achieved by performing this 3-D resolvent analysis on a turbulent mean flow.

Conclusion

A resolvent analysis was carried out on a laminar SWBLI in order to study the low-frequency behaviour of the system for both 2-D and 3-D perturbations. The optimal gain, that can be seen as a transfer function of the system, was computed over several decades of Strouhal numbers. A model of optimal gain based on the excitation of a single global mode was developed, leading to the equation of a first-order low pass-filter, in agreement with the work of [START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shock-wave/boundary-layer interactions[END_REF]. Our analysis was shown to be valid for 2-D perturbations over different base flows calculated from different pairs of Reynolds and Mach numbers, but was not relevant at the largest Reynolds numbers tested. The time scale of the model, i.e. the cutoff frequency of the filter, was found to be determined by the damping rate of the least damped global mode. This damping rate scales as u ∞ /L. Thus, the low-frequency dynamics is associated with a constant Strouhal number based on the recirculation length, as observed experimentally and numerically. The value obtained here is St L 1 × 10 -2 , which slightly underestimates the classical values of the literature. This mismatch could be attributed to the laminar nature of the flow studied in this paper. Ultimately, this work carried out in the laminar regime suggests that the low-frequency dynamics of the SWBLI in the turbulent regime could be a forced dynamics, sustained by the background turbulent fluctuations. The resulting bubble dynamics, which would be sustained by a background forcing produced by nonlinear interactions in a turbulent flow, is then found to exhibit a breathing motion.

Assuming the spanwise direction as homogeneous, a 3-D resolvent analysis showed that two regimes of spanwise wavenumbers β could be identified. At low β, the same modal mechanism as that of 2-D perturbations, leading to a first-order low-pass filter, was observed. Larger values of optimal gain were obtained for wavenumbers of the order of the recirculation length. The possibility of a prevailing 3-D dynamics was then discussed. At larger wavenumbers, streaks generated by the non-modal lift-up effect were detected. Their growth does not interact with the recirculation region, suggesting that streaks are not directly involved in the low-frequency dynamics of the SWBLI. Finally, from the perspective of transition to turbulence, 3-D resolvent analysis remains an efficient tool to characterise the energy growth of streaks and provides their optimal wavenumber and forcing location.

Future work investigating the forcing terms (both two-and three-dimensional) produced by nonlinear interactions in a turbulent SWBLI, either through direct numerical simulation or experimental observations, would be a decisive follow-up to this work. While the low-frequency behaviour of the SWBLI was suggested to be understood as a forced dynamics, the physical origin of the forcing is missing. Ultimately, this could also provide some insights to assess the prevalence of 3-D effects in the dynamics of the system.

Appendix A. Mesh convergence

The influence of the main numerical parameters on the eigenvalues computation of the global stability problem is tested. Eigenvalues associated with the first steady (S1) and unsteady (T1) global modes, that appear in the low-pass filter models in § § 4.1.3 and 4.2.2, are considered. The reference set of parameters used throughout the paper is given on the first line of table 2. Independency with respect to the number of grid points N x and N y is verified in the cases A and B. The influence of the height of the domain L y /δ * 0 and the length of the domain downstream of the impinging shock L DS x /δ * 0 are assessed in the cases C and D. Overall, the unsteady mode is the most easily converged mode: relative differences of eigenvalues between the reference case and the other ones are equal to or below than 0.1 % for the range of tested parameters. As for the steady mode, relative difference of eigenvalues lower than 1 % are found, which is deemed acceptable.

Appendix B. Profiles of the 3-D resolvent modes

Profiles of the forcing and response of the 3-D low-frequency resolvent modes at β = 0.25 and β = 2 are shown in figures 27 and 28 respectively. The x-location of each profile corresponds to that of the maximum of energy density as defined in § 5.1. Note that this location is different for the forcing and the response (see figure 20). As a result, the boundary layer thickness is different on each plot.

Figure 2 .

 2 Figure 2. Validation against the experimental data from Degrez et al. (1987). (a) Wall pressure normalised by the minimum pressure p 0 upstream from the shock. (b) Skin friction coefficient c f . The blue squares correspond to the experimental separation and reattachment points.

Figure 3 .

 3 Figure 3. Local Mach number of the base flow at M = 2.20, Re = 1100.

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Optimal gain at M = 2.20, Re = 1100.

Figure 7 .Figure 8 .

 78 Figure 7. Global stability spectrum at M = 2.20, Re = 1100. Steady and unsteady modes are referred to with the letters S and T, respectively.

Figure 9 .

 9 Figure 9. Comparison between the optimal gain (red squares) and the low-pass filter model from (4.2) (blue line) for different Mach numbers (Re = 1100).

Figure 10 .

 10 Figure 10. Comparison between the optimal gain (red squares) and the low-pass filter model from (4.2) (blue line) for different Reynolds numbers (M = 2.20).

Figure 12 .

 12 Figure 12. Optimal and sub-optimal gains at M = 2.20 and Re = 1100.

Figure 13 .

 13 Figure 13. Resolvent mode 2 at St L = 10 -4 , real part of the streamwise velocity. (a) Forcing. (b) Response.

Figure 14 .

 14 Figure 14. Comparison between the sub-optimal gain σ 2 (red squares) and the low-pass filter model in (4.3) (blue line) for different Mach numbers (Re = 1100).

Figure 15 .

 15 Figure 15. Comparison between the sub-optimal gain σ 2 (red squares) and the low-pass filter model in (4.3) (blue line) for different Reynolds numbers (M = 2.20).

Figure 16 .

 16 Figure 16. Damping rate of the global mode T1 for different Mach and Reynolds numbers (see table1), as a function of the separation length. Here, δ * 0 is the boundary layer thickness at Re = 1100 rather than that at each Reynolds number.

Figure 18 .

 18 Figure 18. Optimal gain as a function of the spanwise wavenumber at St = 10 -4 (Re = 1100, M = 2.20).

Figure 19 .

 19 Figure 19. Optimal resolvent mode at β = 0.25 and St L = 10 -4 , real part of w. Iso-surfaces at 40 % and -40 % of the maximum absolute value are plotted in red and blue. The dark grey surface is the recirculation bubble of the base flow (homogeneous in z).

Figure 20 .

 20 Figure 20. Energy density profiles of the optimal resolvent modes at β = 0.25 (a) and β = 2 (b), St L = 10 -4 . The profiles are normalised by their maximum value.

Figure 21 .Figure 22 .

 2122 Figure 21. Optimal resolvent mode at β = 2 and St L = 10 -4 , real part of w of the forcing (a) and real part of u of the response (b). Isosurfaces at 20 % and -20 % of the maximum absolute value are plotted.

Figure 23 .Figure 24 .

 2324 Figure 23. Global mode S1 at β = 0.25, M = 2.20 and Re = 1100. Isosurfaces of 40 % and -40 % of the maximum absolute value of w are plotted in red and blue.

Figure 25 .Figure 26 .

 2526 Figure 25. Global modes S1 (a) and T1 (b) at β = 2 and St L = 10 -4 , real part of u. Isosurfaces of 20 % and -20 % of the maximum absolute value amplitude are plotted in red and blue.

Figure 27 .Figure 28 .

 2728 Figure 27. Profiles of the optimal resolvent mode at β = 0.25, M = 2.20, Re = 1100, St L = 10 -4 .

Table 1 .

 1 Set of base flows studied with the length of the separation region L obtained in each case.

			* 0
	2.00	1100	0.00
	2.05	1100	2.83
	2.10	1100	22.0
	2.15	1100	43.8
	2.20	1100	63.8
	2.25	1100	78.2
	2.30	1100	88.4
	2.35	1100	94.5
	2.10	1600	45.4
	2.10	2200	77.6
	2.20	600	23.1
	2.20	850	42.7
	2.20	1300	80.3
	2.20	1600	104
	2.20	1900	128

  Instead, the recirculation length L is the relevant length scale. The aforementioned scaling furthermore shows that the associated Strouhal number St L , based on the the separation length, is then constant. Previous works investigating the low-frequency behaviour of the SWBLI have consistently noted a constant Strouhal number from St L 3 × 10 -2 to St L 6 × 10 -2[START_REF] Dussauge | Unsteadiness in shock wave boundary layer interactions with separation[END_REF]. To compare our results with this value, ω

	Because 0 , as explained in § 3.2) have been tested, several Reynolds numbers (based on several δ * this result ensures that ω (S1) i is independent of δ * 0 . (S1) i can be translated into a Strouhal number
	using St L

  Figure 11. Damping rate of the global mode S1 for different Mach and Reynolds numbers (see table 1), as a function of the separation length. Here, δ * 0 is the boundary layer thickness at Re = 1100 rather than that at each Reynolds number.

	10 -2
		M = 2.10
		M = 2.15
		M = 2.20
		M = 2.25
		M = 2.30
	10 -3	M = 2.35
	4
	10 1	10 2

  Declaration of interests. The authors report no conflict of interest. ORCIDs. B. Bugeat https://orcid.org/0000-0001-5420-7531; J.-Ch. Robinet https://orcid.org/0000-0002-3529-6003; J.-C. Chassaing https://orcid.org/0000-0002-6221-2303; P. Sagaut https://orcid.org/0000-0002-3785-120X.

	Author

Table 2 .

 2 Mesh convergence of the most unstable mode at M = 2.30, Re = 1100. The reference set of numerical parameters, used in all the calculation of the paper, is the first one of the table. ω ref corresponds to the difference between the eigenvalue of the reference case and that obtained in each other case.

	Name N x	N y	L DS x /δ * 0	L y /δ * 0	ω (S1) × 10 4	| ω (S1) ref | |ω (S1) ref |	ω (T1) × 10 3	| ω (T1) ref | |ω (T1) ref |
	Ref	600 270	70	50	6.673i	-	3.272 + 2.589i	-
	A	900 270	70	50	6.731i	0.87 %	3.271 + 2.589i	0.02 %
	B	600 400	70	50	6.651i	0.33 %	3.272 + 2.585i	0.10 %
	C	680 270	90	50	6.737i	0.96 %	3.271 + 2.589i	0.02 %
	D	600 343	70	75	6.631i	0.63 %	3.272 + 2.589i <0.01 %


domain: sde shs

text:
Pauline Bremond
email: pauline.bremond@inrae.fr

Pierre Balzergue

Frédéric Grelot

Valériane Marry

Maxime Modjeska

David Nortes Martinez

Contribution of a systematic and long-term observation method to improve the understanding of flood impacts

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Floods are the natural hazard that generates the most damage in the world, but many aspects of how this damage occurs are still poorly understood. Indeed, the impacts of floods can concern all types of stakes (dwellings, economic or agricultural activities, public infrastructures, networks, etc.), whether they are directly impacted by the flood or indirectly by the propagation of disturbances through society, or even in a temporal delay compared to the original event. Moreover, these impacts may differ according to the nature of the flood (e.g. overflow, runoff, salt water flooding). Finally, the impacts are the result of complex biophysical processes but also of protection, repair and adaptation actions implemented by affected people.

Our objective is to highlight the diversity of flood impacts, including the issue of their persistence. To this end, we propose an original method that combines post-flood observation surveys with longterm monitoring of individuals exposed to floods. This methodology is implemented in the framework of the system of observations of the impacts of floods (so-ii, http://so-ii.org) in the Greater Montpellier area coordinated by our team. For the quantitative surveys, based on our past observation and modelling experiences, we have developed questionnaires that aim to collect the impacts of flooding, following an event, in an exhaustive and cumulative manner. These questionnaires were applied following a flood event that occurred on September 20th 2020 in the north of Montpellier. For long-term monitoring, we set up networks of impact observers (inhabitants and farmers for the moment), with whom we agreed to work over the long term (about fifteen years). With these groups, we have defined an observation protocol, combining individual interviews to better specify each person's situation (e.g. individual vulnerability, characteristics of buildings) with participatory workshops to discuss the results, orient future work, and pool each person's experiences. In this paper we present the intermediate results of this ongoing research. We also discuss the relevance and limitations of such an approach to improve the understanding and assessment of flood damage.

Powered by TCPDF (www.tcpdf.org)