Lionel Agostini

Lionel Larchevêque
email: lionel.larcheveque@polytech.univ-mrs.fr.

Pierre Dupont

Jean-François Debiève

Jean-Paul Dussauge

Zones Of Influence And Shock Motion In A Shock Boundary Layer Interaction

This paper aims at describing the main features of a shock reflection on a turbulent boundary layer. The data used for this analysis are the results of Large Eddy Simulations of the interaction carried out with three different shock intensities, from incipient to fully separated cases. Computational results are validated vs experiments obtained for the same interaction geometries. The main space-time properties of the leading shock motions are described together with their links with the other regions of the flow. In particular, information about the origin of the shock motion is derived from the correlations between shock motion and unsteady pressure field. It is shown that the shock motion reveals the flow unsteadiness found in the interaction region. Nomenclature α Angle of the characteristics (see Eq. 3) α ′ Angle of the Mach wave (see Eq. 6) ∆x + Streamwise cell dimension in wall unit ∆z + Spanwise cell dimension in wall unit δ Boundary layer thickness

St LP Strouhal number of the disturbances along the leading shock.

U Flow velocity

U 1 Velocity downstream the incident shock U η Propagation velocity along characteristics (see Eq. 4)

U ζ

Propagation velocity along Mach waves produced by convected sources (see Eq. 7)

U c Convection velocity

V P Propagation velocity along the shock V PC Propagation velocity along the shock for disturbances associated with characteristics (Eq. 5).

V PMW Propagation velocity along the shock for disturbances associated with Mach waves (Eq. 8).

X * Longitudinal dimensionless coordinate (x -X 0 )/L X 0 Mean position of the foot of the leading shock

I. Introduction

Shock/boundary layer interactions are found in many aeronautical applications such as air intakes or over expanded nozzle flows. Fully separated interactions have been widely studied in the past [START_REF] Dolling | Fifty years of shock-wave/boundary-layer interaction research: what next[END_REF][START_REF] Beresh | Relationship between upstream turbulent boundary layer velocity fluctuations and separation shock unsteadiness[END_REF][START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Ganapathisubramani | Effects of upstream coherent structures on lowfrequency motion of shock-induced turbulent separation[END_REF][START_REF] Wu | Analysis of shock motion in shockwave and turbulent boundary layer interaction using direct numerical simulation data[END_REF]. Whatever the flow configuration (compressions ramps, shock reflection, blunt fin or over-expanded nozzles), low frequency shock motions are observed when the flow is separated, although their origin is not always clearly established [START_REF] Ganapathisubramani | Effects of upstream coherent structures on lowfrequency motion of shock-induced turbulent separation[END_REF][START_REF] Wu | Analysis of shock motion in shockwave and turbulent boundary layer interaction using direct numerical simulation data[END_REF][START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Touber | Oblique shock impinging on a turbulent boundary layer: low-frequency mechanisms[END_REF][START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shockwave/boundary-layer interactions[END_REF]. In the case of the reflection of a shock wave on a turbulent boundary layer at moderate Reynolds number, it has been recently proposed to relate the low frequency motion of the leading shock (denoted reflected shock in see Fig. 1) to the dynamics of the separated region developped downstream [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF]. The link between the leading shock motion and the low frequency unsteadiness was not fully specified although some elements have been examined. For example, the relation between the low frequency breathing of the separated region and the dynamics of the large coherent scales formed in the mixing layer originating at the foot of the leading shock has been discussed. [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF]. These structures produce frequencies at least one order of magnitude higher than the low frequency shock motions. Their influence on the shock unsteadiness has still to be explored. Moreover recent experimental investigations of various shock reflections have shown that incipient separation presents several similarities with separated cases, with evidence of low frequency shock unsteadiness as in separated cases [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Souverein | Application of a dual-plane particle image velocimetry (dual-PIV) technique for the unsteadiness characterization of a shock wave turbulent boundary layer interaction[END_REF].

Recent work by Touber and Sandham [START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shockwave/boundary-layer interactions[END_REF] has proposed an analysis in which the response of the shock system to external perturbations is examined, providing an estimate of the low frequency range of the interaction and justifying a former heuristic theory. This approach however was not well adapted to the determination of the origin of the low frequencies. The results presented here provide a quite different viewpoint. To some extent, it complements the results quoted here above. From simulation of comparable parameters, the details of the physical mechanism are explored, giving rise to the upstream/downstream influence on the shock motion; it also identifies spatially the zones controlling the motion of the different part of the shock system, along with the frequency ranges which are involved. All these analyzes are applied to a shock reflection configuration described in the following.

Fig. 1: Schlieren visualization of the fully separated interaction. Adapted from [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF].

A Mach 2.3 shock reflection was set up in the IUSTI supersonic wind tunnel. The flow deflection angle θ can be set from θ = 5.5°to θ = 9.5°causing flow separation, from incipient to fully separated.

The organization of the flow is illustrated by a Schlieren visualization of the interaction in Fig. 1. Results include (time-resolved) hot wire measurements along the leading shock and velocity field measurements from Particle Image Velocimetry (not resolved in time). Details on the experimental set up and the results can be found in Refs. [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF][START_REF] Souverein | Effect of interaction strength on unsteadiness in shock-wave-induced separations[END_REF][START_REF]Unsteady effects of shock wave induced separation[END_REF]. The Reynolds number based on the momentum thickness of the upstream boundary layer is about 5000, a convenient value to perform Large Eddy Simulation of such a flow.

The present work complements the experimental results by bringing time resolved data obtained from Large Eddy Simulations of the same configuration.

The numerical method and aerodynamic parameters are described in section II. The unsteadiness of the leading shock is characterized in section III. As will be shown, it involves two distinct frequency bands.

Each one will be related with phenomena occuring in particular regions of the flow and the main features of the communication paths from these regions to the shock are described in section IV.

II. Numerical simulation of the interaction

A. Computational model and parameters

Since the flow under study is fully turbulent with a moderate Reynolds number and since it develops low-frequency unsteadiness, a turbulent modeling relying on the Large-Eddy Simulation method appears to be a good candidate to perform time-accurate computations. As a matter of fact, LES has been proved to be a suitable modelling to resolve reflecting shock/turbulent boundary layer interactions of the type described above. The LES have produced accurate predictions of the mean features of the flow [START_REF] Garnier | Large-eddy simulation of the shock/boundary layer interaction[END_REF] and of the lowfrequency unsteadiness, in good agreements with experiments [START_REF] Touber | Comparison of three large-eddy simulations of shock-induced turbulent separation bubbles[END_REF].

The various computations are performed using ONERA's FLU3M solver that has been extensively used in the recent years to analyze successfully compressible flows either by DES, LES and DNS [START_REF] Deck | Unsteadiness of an axisymmetric separating-reattaching flow: Numerical investigation[END_REF][START_REF] Simon | Numerical simulation of the compressible mixing layer past an axisymmetric trailing edge[END_REF][START_REF] Riou | Compressibility effects on the vortical flowover a 65°sweep delta wing[END_REF]. The numerical scheme is designed to be able to capture the shock while meeting the LES requirement of very low dissipation in the turbulent region [START_REF] Garnier | Stimulated detached eddy simulation of three-dimensional shock / boundary layer interaction[END_REF]. This is achieved by adding the dissipative part of the Roe scheme [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF], modulated by Ducros' sensor [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF], to a second order centered scheme. The subgrid filtering is implicitly provided by the mesh and the subgrid modeling relies on the selective mixed-scale subgrid model, well suited for compressible wall bounded flows [START_REF] Lenormand | Subgrid-scale models for Large-Eddy Simulation of compressible wall bounded flows[END_REF]. Time integration is achieved by means of a secondorder accurate implicit Gear scheme [START_REF] Gear | Numerical Initial Value Problems in Ordinary Differential Equations[END_REF] with a timestep of 2.5 × 10 -7 s. yielding maximal CFL numbers of 12 so as to make the implicit time filtering negligible with respect to the implicit grid filtering. The resulting non-linear system is solved iteratively at every timestep with 6 sub-iterations, yielding a reduction of the residuals of 1.5 order of magnitude at worst and more than 2.5 order of magnitude on average.

Following [START_REF] Garnier | Large-eddy simulation of the shock/boundary layer interaction[END_REF], the computations aim at reproducing only the mid-span part of the wind-tunnel, excluding the influence of the side walls. Periodic boundary conditions are consequently used in the spanwise direction, with a spanwise extent of the computational domain equal to 1.6 boundary layer thicknesses δ, corresponding at least to 3 times the expected height of the separation bubble. The domain inflow and outflow boundaries are located 10δ away of the interaction region and rely on characteristics boundary condition.

The mesh is designed to match the resolution required to solve wall-bounded flows by LES with values ∆x + ≃ 40, y + 1 ≃ 0.9 and ∆z + ≃ 16. The grid stretching in the vertical direction is mild in order to be able to accurately capture the recirculation bubble, yielding a number of cells located in the initial boundary layer thickness of 95. The mesh is stretched in the streamwise direction starting from 5δ downstream of the interaction region in order to progressively damp the convected vortices. The total cell count is equal to 5.5 millions. Such a moderate grid size allows a very long time integration ranging from 0.5 million to 2 millions timesteps, corresponding to 0.125 to 0.5 seconds and being large enough to resolve 50 to 150 periods of the low-frequency oscillations. These numbers of periods are comparable to or larger than values found in previous studies [START_REF] Wu | Analysis of shock motion in shockwave and turbulent boundary layer interaction using direct numerical simulation data[END_REF][START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shockwave/boundary-layer interactions[END_REF][START_REF] Garnier | Stimulated detached eddy simulation of three-dimensional shock / boundary layer interaction[END_REF][START_REF] Kawai | A parametric investigation of oblique shock / turbulent boundary layer interaction using LES[END_REF].

Two additional LES of the case with the largest separation have been carried out for validation purposes.

the initial boundary layer thickness. The timestep was reduced to 2 × 10 -7 s. for this simulation. The second computation is based on the same grid resolution and the same timestep as the original one but with the span of the computational domain extended to 11δ. Both computations were run for 35 ms., a time large enough to obtain converged first and second order statistics and allowing the spectral analysis of about 10 periods of low-frequency oscillations.

Velocity, density and pressure data over various horizontal, vertical and transverse planes are stored on disk with a 200 kHz sampling rate, except for the computation on the refined mesh for which a 500 kHz sampling rate was set. The 200 kHz rate is not high enough to sample adequately the data in the incoming boundary layer region for which the maximum in the pressure spectrum is located at about 70 kHz. However comparisons between data sampled at 200 kHz and 500 kHz demonstrate that the former rate is high enough to ensure a negligible aliasing in the interaction region. In this region, less than 5% of the total density power spuriously originates from the aliasing for every frequency in the ranges of interest, namely the low frequency unsteadiness and the intermediate frequencies corresponding to the large scales of the mixing layer/shedding regions.

B. Boundary conditions

An adiabatic boundary condition is enforced at the wall and characteristic-based boundary conditions are used at the outflow and upper boundaries. Turbulent boundary conditions at the inflow are set using a variant of the Synthetic Eddy Method [START_REF] Jarrin | A synthetic-eddy method for generating inflow conditions for les[END_REF][START_REF] Pamiès | Generation of synthetic inflow data for large-eddy simulation of spatially evolving wall-bounded flows[END_REF] suited to supersonic boundary layers: the method described in [START_REF] Pamiès | Generation of synthetic inflow data for large-eddy simulation of spatially evolving wall-bounded flows[END_REF] was modified in the following way:

• Increase by a factor of ten of the density of the synthetic eddies yielding a more Gaussian-like probability distribution of the velocity fluctuations. It makes it possible to damp the tails of the distribution associated with fluctuations of huge kinetic energy resulting in non-physical high values of the local Mach number.

• Addition of density fluctuations computed from the streamwise velocity fluctuations using the Strong Reynolds Analogy coupled with the linearized ideal gas law and zero pressure fluctuations so as to achieve zero mean.

• Renormalization of the strength of the eddy locally in order to set at the inflow exactly the wanted Reynolds stress profiles even in regions of overlapping modes.

• Use of temporal characteristics to achieve non-reflectiveness of the boundary condition in the subsonic region of the boundary layer.

Following [START_REF] Pamiès | Generation of synthetic inflow data for large-eddy simulation of spatially evolving wall-bounded flows[END_REF], 5 modes were distributed over the boundary layer height, using the proposed loci, lengthand timescales but with Gaussian shape functions only.

It has been verified by comparisons with the experiments and with an auxiliary LES computation of a freely transitioning boundary layer that this method is able to recover accurate turbulence statistics and the expected organization within a development length smaller than 10δ 0 . At the reference station x = 260 mm, located 9.5δ 0 downstream of the inflow and at least δ 0 upstream of the interaction region, the values of the friction velocity computed either from the derivative at the wall or using Clauser's method differ by less than 4% and than 1.5% respectively from the experimental values [START_REF] Laurent | Turbulence d'une interaction onde de choc-couche limite sur paroi plane adiabatique ou chauffée[END_REF][START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF] estimated from Clauser's method.

Moreover the differences between the computations and the experiments regarding the displacement and the momentum thicknesses are lower than 2%.

These concordances can be assessed by looking at the mean and fluctuating velocity profiles plotted in wall unit in Fig. 2. Note that the slight differences seen in the van Driest velocity profiles can be almost suppressed by using Clauser's method to estimate the friction velocity from the LES data. The most noticeable discrepancy is found at the end of the log region where the streamwise fluctuations are underestimated with respect to the experiments. However this discrepancy was also found, though slightly weakened, in the LES computation of a freely transitioning boundary layer mentioned above. Consequently it has not to be associated with the SEM method. Lastly, the method has been found to be free of low to medium frequency forcing found especially for the pressure in preliminary tests [START_REF] Larchevêque | Experimental and numerical study of unsteadiness in boundary layer / shock wave interaction[END_REF] for which the compressible variant [START_REF] Sagaut | Turbulent inflow conditions for large-eddy simulation of supersonic and subsonic wall flows[END_REF] of the Lund recycling method was used.

C. Flow parameters

Three LES were made to reproduce flow cases in incipient, mild and full separation conditions, respectively corresponding to flow deviation of 5.5°, 8.0°and 9.5°in the experiments performed at IUSTI.

However, these shock reflection experiments were sensitive to lateral effects due to the finite span of the x = 260 mm, 9.5δ0 downstream of the inflow plane: standard LES (solid), refined LES (dashed), experiments [START_REF] Laurent | Turbulence d'une interaction onde de choc-couche limite sur paroi plane adiabatique ou chauffée[END_REF] using Pitot tube for the mean velocity and LDA for the fluctuations (circles). More recent experiments [START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF] using LDA measurement at the location x = 240 mm, 7.4δ0 downstream of the inflow plane, have been added (triangles). The friction velocity is computed using the Clauser method for all the experiments. wind tunnel, as recently shown by the RANS computations carried out for the European UFAST project [START_REF]Unsteady effects of shock wave induced separation[END_REF] and the DES by Garnier [START_REF] Garnier | Stimulated detached eddy simulation of three-dimensional shock / boundary layer interaction[END_REF]. This leads to experimental sizes of interactions being up to 20 % larger than the lengths found in the computations not taking into account the side walls. Such differences may be of low influence on the unsteadiness for the θ = 8.0°and θ = 9.5°cases since the change in size of the interaction region is not large enough to induce a change in the separation state. Therefore computations associated with these cases are carried out for the same flow deviation angle as the experiments.

The incipient case should on the contrary be very sensitive to any reduction in the size of the interaction that possibly could result in changing the onset of separation. The flow deviation angle of the incipient computation has consequently been adjusted iteratively to θ = 6.3°so as to result in an interaction length L (see Fig. 1) defined as the distance from the foot of the leading shock to the extrapolation down to the wall of the incident shock, equal to the experimental value of L = 25 mm found in the experiments. The flow computed from this setup exhibits a similar separation state as the experiment, as demonstrated by the distribution of the reverse flow probability plotted in Fig. 3.

Note that the interaction lengths L and some other geometrical parameters useful for further analyzes are listed for the three interactions in Tab. 1.

D. Validation

The accuracy of the LES can first be evaluated by analyzing the sensitivity of the computations to the definition of the computational grid for the 9.5°case. The divergence between the standard computation, the LES carried out on a refined mesh and the simulation with a wider span have first to be analyzed for the incoming boundary layer. No difference at all is found between the standard and the span-extended computation in that region whereas the comparison between the standard and refined computations shows that metrics of the boundary layer such as the friction velocity, momentum thickness and shape factor differ by less than 1.5% between the two simulations. When focusing on the interaction region, the standard and refined LES exhibit the very same properties:

separation length as well as interaction length differ by less than 1%, as seen in Fig. 4. Moreover almost identical geometry of the separated region and streamwise evolution of the friction coefficient are seen in this figure. The computation with an enlarged span also shows very similar features but with lengths reduced by 2% (interaction) and 6% (separation).

This has to be associated with the rather large asymptotic values of the spanwise correlations coefficient found in some part of the interaction for the LES with the 1.6δ span: c uu (ξ = 0.8δ) = -0.26 close to the wall at the beginning of the separation and c uu (ξ = 0.8δ) = -0.17 in the shedding region. The narrow span computation of the 9.5°case appear therefore to be slightly confined. However no other difference beyond the confinement-induced enlargement has been observed between the physics found in the narrow and wide computations. This is consistent with the analysis of one of the computations described in [START_REF] Touber | Low-order stochastic modelling of low-frequency motions in reflected shockwave/boundary-layer interactions[END_REF] subject to a much higher confinement. Particularly the unsteady behavior of the interaction, as deduced from the wall pressure spectra, appears to be unaffected once normalized using L. Note that the short duration of the enlarged computation results in a noticeable level of statistical uncertainty for the spectra but this can be counterbalanced for the most part by a space averaging in the spanwise direction.

Maximum absolute values of the streamwise correlation coefficient of velocity were found lower than 0.1 at every location of the 6.3°and 8.0°computed flowfields. It is therefore very unlikely that these computations are subject to spanwise confinement and no extra computation with an extended span has been done.

A grid independence study has been carried out for 4.8°computations not described in this paper by using the same standard and refined grids as the 9.5°case. Grid convergence was fulfilled for the θ = 4.8°c ase similarly to the 9.5°case. Consequently it is believed that refined grids should not bring significant improvements for the 6.3°and 8.0°computations.

The ability of the computations to reproduce the main features of the flow unsteady behavior can be evaluated by looking at the streamwise evolution of the wall pressure spectra plotted in Fig. 5. All three computations as well as the θ = 8.0°and θ = 9.5°experiments exhibit a similar development of the interaction, as described in [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF]: the first part of the interaction is dominated by low frequency component, followed by a transition region of increasing characteristic frequencies related to the development of the mixing layer that eventually result in a shedding process with an almost constant frequency, lower than the typical ones in the incoming boundary layer.

Assessment of the fidelity of the LES can be performed in a more quantitative way by normalizing both lengths and frequencies to take into account the differences between experiments and computations regarding the size of the interaction. The dimensionless coordinates X * = (x -X 0 )/L and Y * = y/L are used, where X 0 is the mean position of the foot of the leading shock. In this representation, the interaction extends from X * = 0 to X * = 1. In separated cases, previous works [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Erengil | Unsteady wave structure near separation in a Mach 5 compression ramp interaction[END_REF] have used a non dimensional shock frequency, or Strouhal number St L , defined as :

St L = f.L U 1 ( 1
)
where f is the frequency, L the interaction length and U 1 is the velocity downstream the incident shock. It shock motions around St L ≃ 0.03 (see [START_REF] Dussauge | Unsteadiness in shock wave boundary layer interactions with separation[END_REF]). Such a value is clearly seen in Fig. 5(d No such well defined low frequency content was found in the experimental spectra related to the incipient separation despite evidences of energy content in this range [START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF][START_REF] Souverein | Application of a dual-plane particle image velocimetry (dual-PIV) technique for the unsteadiness characterization of a shock wave turbulent boundary layer interaction[END_REF]. This appears when looking at the premultiplied power spectrum of Fig. 6 computed by the Burg method from the time series of a hot wire located on the mean location of the leading shock at y = 1.8δ. The streamwise momentum premuliplied spectrum from the incipient LES computation at the same location nonetheless seems to result in a slightly noticeable bump roughly centered at St L = 0.03 but this difference between experiments and computations may be again due to the difference in the temporal length of the signal. When the experimental signal is restricted to the same duration as the LES computation, spectra from the experiments and the LES show very similar features with the exception of the highest frequencies, because of the strain gauge effect of the hot-wire.

It is therefore concluded from Fig. 5-6 that the present LES are able to accurately reproduce the dominant unsteady features found in the experiments while giving access to time-resolved data encompassing the whole flowfield. In the next part the main results obtained from the LES computations will be presented, both for the mean fields as well as for the unsteady aspects. Here the height H i has been used to normalize the ordinate y. When changing the deviation angle θ, the shock angle varies, so that the triangles formed by the crossing of the incident and leading shocks cannot be similar. However, in practice, as seen in Table 1, the differences are small, the pattern is almost unchanged, and normalizing by H i or by L would lead to the same evolutions. The two frequency ranges

St L ≃ 0.03 and St L ≃ 0.5 are observed with different importance. Moreover, the spectra normalized by L 2 appear to collapse at low frequencies (St L ≃ 0.03). This is in agreement with [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF] in which it is shown that the measured length of leading shock excursion varies roughly like the length of interaction. As this amplitude of the shock motion is closely related to the low frequencies, the observed collapse of the spectra when normalized by L 2 is not surprising. This suggests still more firmly that L is the pertinent length scale for these low frequency motions, leading to self similar behavior in this frequency range.

Low frequency unsteadiness is expected in the separated case. However, it is more surprising to find it also in the incipient case at the same Strouhal number. In previous work, the same characteristic low and intermediate Strouhal numbers were observed in the interactions for separated cases [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF]. They have been respectively related to some breathing of the separation bubble and to the development of vortical structures of Kelvin Helmholtz type in the mixing layer downstream of the separation shock [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Piponniau | A simple model for low-frequency unsteadiness in shock-induced separation[END_REF]. The former physics was energy-dominant. In the incipient case, it seems that identical sources of unsteadiness are found, but in the normalizing of the displacement at low frequencies by the interaction length yields a same amplitude as in the separated cases. The origin of both low-and intermediate-frequency shock motions will be discussed in the next sections.

B. Correlation fields in the interaction

The pressure fields are now considered to investigate where the characteristic frequencies St L ≃ 0.03 and St L ≃ 0.5 are present in the field and how they can influence the leading shock motions. Two pressure correlation fields have been computed, associated with the two characteristic frequencies. They are derived from the following procedure :

• Firstly, the pressure fluctuations of the whole flowfield and the time series of the streamwise location of the shock at a given elevation are band-passed around the considered frequencies, with a bandwidth of 0 ≤ St L ≤ 0.1 for the low frequency range and of 0.3 ≤ St L ≤ 0.8 for the intermediate frequency range.

• Then, the cross-correlations between pressures and the shock streamwise location at the selected elevation are computed.

This procedure has been applied for the three interactions under consideration.

Figure 8(c) presents the fields of the resulting cross-correlation coefficients in the low frequency range for the mildly separated interaction (θ = 8.0°). The coefficients have been computed by taking as reference the shock motion at elevation y/H i = 1.4 for practical reasons. Note however that the results were found insensitive to the location along the shock. Figure 8(c) puts in evidence a high level of correlation between the shock and the initial part of the interaction, whereas the expansion wave, the second part of the interaction and the downstream flow are in anti-correlation with the shock position. Similar results were already obtained for shock reflection as well as for compression ramp from experimental [START_REF] Dolling | Fifty years of shock-wave/boundary-layer interaction research: what next[END_REF][START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF][START_REF] Thomas | On the mechanism of unsteady shock oscillation in shock wave/turbulent boundary layer interaction[END_REF][START_REF] Debiève | Dependence between shock and separation bubble in a shock wave / boundary layer interaction[END_REF] as well as numerical [START_REF] Touber | Oblique shock impinging on a turbulent boundary layer: low-frequency mechanisms[END_REF][START_REF] Larchevêque | Experimental and numerical study of unsteadiness in boundary layer / shock wave interaction[END_REF] unsteady wall pressure data. They are now generalized in the whole field of the interaction. It can be seen This is due to both an increase in the Strouhal numbers and a decrease in the convective velocity of the shed vortices in region 0.8 X * 2 when θ is increased. Furthermore, the pattern extents further upstream into the mixing layer region for the fully separated case. The wavelength is slightly lower in that region. This is in rather good agreement with the value determined experimentally from two-point wall pressure measurements that were associated with the mixing layer structures [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF].

It consequently appears that the cross-correlation pattern found at the end of the mixing layer and in the shedding region is directly related to convected vortices. The nature of the pattern found along the leading shock is however less clear and further analyzes are required to clarify that point.

C. Propagation and convection velocities in the interaction

Convection velocity of vortical structures in the mixing layer

In previous experimental works, two point measurements of unsteady wall pressure have been performed. Evidences of links were found between the leading shock and the other regions of the interaction, depending on the frequency range [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF]. Typical convection velocities were determined and associated with each frequency domain.

At low frequency, no convective behavior was found in the interaction, except for the 9.5°case, within a small region in the vicinity of the reattachment point [START_REF] Debiève | Dependence between shock and separation bubble in a shock wave / boundary layer interaction[END_REF].

On the opposite, the intermediate frequencies were observed to be convective. In the region 0.2 X * 0.5, phase velocities deduced from wall pressure measurements were found to be rather insensitive to the value of the deviation angle θ, with a typical value of 170 m/s. This value was associated with the convective velocities U c of the structures in the mixing layer [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF].

These experimental results have been confirmed by considering the unsteady wall pressure data computed by LES. Moreover, LES give an access to unsteady data over the whole flowfield, allowing a more direct evaluation of U c by using phase velocities deduced from velocity measurement along the center line of the mixing layer. Convective velocities ranging roughly from 160 m/s for the fully separated case to 200 m/s for the incipient case are obtained, in good agreement with previous values.

LES data were also be used to compute the propagation velocity V P associated with the corrugation of the shock. Time delays τ opt optimizing the cross-correlation of the shock location between two points taken along the leading shock separated by a distance (ξ) have been determined. The reference point is located just above the crossing point of the two shocks, whereas the second point is placed at higher altitudes, along the leading shock. Times delays are then converted into propagation velocity V P by the relation :

V P = ξ τ opt (2)
The optimum delay times for intermediate frequencies (0.3 < St L < 0.8) along the leading shock have been evaluated for the three interaction intensities. They are reported in Fig. 9. The local velocity of propagation along the shock is equal to the local slope of the curve ξ(τ opt ) . Several behaviors can be observed depending on the state of separation and on the distance from the reference point:

• For the incipient separation case (θ = 6.3°) the expansion fan is merged with the leading shock (see [START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF]). The propagation velocity along the shock is found approximately equal to 400 m/s.

• For the strongest shock intensity (θ = 9.5°) the expansion wave intersects the leading shock outside of the computational domain. The region below the intersection point H e of the leading shock and the expansion fan has a nearly constant convection velocity of about 250 m/s.

• For the shock of intermediate intensity (θ = 8.0°) the computational domain is such that it is possible to evaluate the propagation velocity in both regions: below and above the intersection of the expansion fan with the leading shock. In the first region, between the shock crossing and the intersection with the expansion, the convection velocity is the same as in the θ = 9.5°case (250 m/s), while above this region, where the leading shock and the expansion collapse, the same convection velocity as in the incipient case is obtained (400 m/s).

To summarize, the shock kinematics at intermediate frequencies can be split into two parts on each side of the about 400m/s. The possibility of a link between these two different velocities and the kinematics of the vortical structures found in the mixing layer and in the shedding region will be analyzed in the next section.

IV. Pressure disturbance propagation within the flow

A. Reminders about characteristics and disturbance propagation

An attempt will be made to interpret some of the previous results in terms of characteristic directions. This is often possible in flows with pressure gradients, since in turbulent flows, even in supersonic conditions, outside of the viscous sublayer, the pressure force is larger than friction. Therefore, outside of this thin viscous zone, the equations of motion reduce to Euler equations.

Under the assumption of an inviscid and supersonic flow, the hyperbolic equations governing compressible flows can be recast into the characteristic form. In this case, pressure fluctuations propagate along the characteristic lines, which are known as Mach waves. Shock waves are located at the coalescence of such lines. Therefore the directions and the speed of transmission of information described in the previous sections will be compared to the properties of characteristic directions.

Such analyses are carried out by considering the mean flowfield and are consequently invalid for the high frequency range. The upper limit of validity can be estimated to St L ≪ 1 as shown in appendix.

Interpretation of the data based on this model will therefore be meaningful for the St L ≃ 0.03 frequency range. It is also expected it will help analyzing the flow kinematics at intermediate frequencies considered, with their properties at time 3∆t. After a time ∆t, the source has moved along a distance ∆t U c .

St L ≃ 0.5 ∆t U. ∆t U.2 ∆t U.3 ∆ t C .1 ∆ t C .2 ∆ t C .3 C h a r a c te
As we consider the situation at time 3∆t, the wave front has traveled over a distance 2∆t U , and its radius is 2∆t C. For a perturbation emitted at time 2∆t, the reasoning is the same and corresponds to the third circle, of radius ∆t C, in Fig. 10(b). Simple geometrical similarity considerations show that the envelope of the wave fronts (Mach wave) is a straight line of slope α ′ defined in Eq. 6.

In this case of moving sources, disturbances propagate along the Mach waves with a speed U ζ given by the companion relation of formula 4:

U ζ = (U -U c ) 2 -C 2 . ( 7
)
Moreover, when the Mach wave produced by such a convected structure intersects the shock wave of angle σ with respect to the direction of convection, the propagation speed V PMW of the disturbance along the shock is equal to :

V PMW = U c sin(α ′ ) sin(α ′ -σ) (8)
Note that the two previous relations are obtained for U and U c being collinear. It has been checked that this condition condition holds approximately for the present flow cases.

Geometrical and kinematical properties of the disturbance propagations described by Eqs. marized in Fig. 11. It is obvious from this figure that steady and convective sources radiate in different directions. Consequently this general framework will be used in next sections to seek for the sources located within the interaction region that could be responsible for the propagation of perturbation along the leading shock.

B. Regions of influence for the shock unsteadiness at low frequencies

As already mentioned, no convective properties were found over the low frequency range for the most part of the interaction region [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF]. Consequently the domain of influence of the low frequency unsteadiness has rather to be searched for by considering the characteristics. This domain is visualized in Fig. 8(a) for the 8.0°interaction by plotting the field of mean characteristic. As seen on this figure, the isolevels of crosscorrelation follow closely the characteristic lines. The regions of origin for the disturbance along the shock at low frequency can therefore be identified as follow :

• The region of the leading shock located below H e is affected by disturbances coming from the first part of the separated zone 0 < X * < 0.5.

• The kinematics of region above H e is governed by disturbance sources located at the foot of the expansion fan with X * ≃ 0.6 -0.7.

It appears consequently that the low frequency motion of the whole leading shock has to be associated with the low frequency motion of the separation bubble and/or of the interaction zone. This result seems to hold for the fully separated case θ = 9.5°in Fig. 8(e) but also, more surprisingly, for the incipient case θ = 6.3°, although only in the y ≥ H e region because of the evanescence of the y < H e region in this case, as seen in

Flow regions influencing the leading shock above the intersection with the expansion fan

Performing the same analyzes as in the previous section leads to the conclusion that Mach waves radiated from the mixing layer region are not able to cross the leading shock above its intersection H e with the expansion fan. Other possible sources for Mach waves can be looked for, for example sources located further downstream in the vortex shedding region X * ≥ 0.7. A difficulty arises: the convection velocity of the vortical structures is far from being constant in that region, for instance there is a 20% increase between X * = 1 and X * = 3 for the θ = 6.3°incipient case. This results in a noticeable level of uncertainty in the definition of the paths followed by the Mach waves possibly radiated from the shed vortices.

St LP = V P L λ P U 1 (9)
where λ P is the wavelength of the disturbances. The value of λ P above H e can be accurately inferred only for the incipient case because of the restricted size of the sampling domain. It reads λ P = 47 mm, yielding

St LP = 0.41 through Eq. 9. This is in acceptable agreement with the St L ≃ 0.5 value associated with the shedding process [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF] when taking into account the rather large uncertainty on the St L measurement in the shedding region because of the turbulent nature of the interaction. The postulated scenario appears therefore to be plausible.

V. Conclusions

The unsteady behavior of a shock reflection on a turbulent boundary layer at Mach number 2.3 has been investigated using results of Large Eddy Simulations. Several cases have been computed, ranging from incipient to full separations. The computations have been carefully compared with experiments with a special emphasis put on the unsteadiness of the interaction region. The long-time computation led to wall pressure spectra with good statistical convergence from which a low and an intermediate frequency bands have been put in evidence whatever the state of separation. The corresponding Strouhal numbers of about 0.03 and 0.50 are similar to the experimental values [START_REF] Dupont | Space and time organization in a shock induced boundary layer[END_REF].

These two frequency bands are also found when the spectra of the shock location are considered. The respective proportion of the energy content associated with these two bands varies with the state of separation.

It is however demonstrated that a normalization based on the interaction length L results in a collapse of all the spectra in the low frequency range for all separation states from incipient to fully separated. In conclusion, the shock kinematics seems to be only the mirror of the physical phenomena localized in the separated zone, whatever the frequency range.

∂ρ ′ ∂t + u j ∂ρ ′ ∂x j + ρ ′ ∂ u j ∂x j + ∂ρu ′′ j ∂x j + ∂ρ ′ u ′′ j ∂x j = 0 (10)
We want to compare the first two terms of Eq. 10, using a decomposition of the form ρ ′ = ρ exp i k. x -ωt . The first term becomes -iωρ ′ , while the second reads i u j k j ρ ′ . The time derivative term can therefore be neglected if:

|ω| ≪ |k j u j | , (11)
or, by recasting Eq. 11 with frequency f and using U 1 and L as typical velocity and length scales in the interaction region:

f L U 1 ≪ 1 ⇔ St L ≪ 1 ( 12
)
The conclusion is that for St L ≪ 1, the equations of motion reduce to the steady Euler equations, and therefore pressure perturbations follow the characteristics defined by the mean field.

1 δ 0

 10 Boundary layer thickness at the inflow λ P Wavelength of the disturbances along the leading shock σ Angle between the shockwave and streamlines τ opt Optimal delay time of the cross-correlation function θ Elevation of the crossing point between the leading shock and the expansion fan H i Elevation of the crossing point between the incident and leading shocks L Length of the interaction M c Convective Mach number St L Strouhal number (see Eq. 1)

Fig. 2 :

 2 Fig. 2: Van Driest transformed mean velocity profile (a) and fluctuating velocity profiles (b) at the reference locationx = 260 mm, 9.5δ0 downstream of the inflow plane: standard LES (solid), refined LES (dashed), experiments[START_REF] Laurent | Turbulence d'une interaction onde de choc-couche limite sur paroi plane adiabatique ou chauffée[END_REF] using Pitot tube for the mean velocity and LDA for the fluctuations (circles). More recent experiments[START_REF] Dupont | Investigation by particle image velocimetry measurements of oblique shock reflection with separation[END_REF] using LDA measurement at the location x = 240 mm, 7.4δ0 downstream of the inflow plane, have been added (triangles). The friction velocity is computed using the Clauser method for all the experiments.

Fig. 3 :

 3 Fig. 3: Percentage of reverse flow from PIV measurements (a) and LES data (b).

  θ

Figure 2 (

 2 Fig. 2(b) also collapse correctly.

Fig. 4 :

 4 Fig.4: Separation line and geometry of the lower part of the incident and leading shocks (see Fig.1): standard LES (solid), refined LES (dashed) and standard LES with increased span (dotted).

  (a) LES, θ = 6.3°(b) LES, θ = 8.0°( c) LES, θ = 9.5°(d) Experiments, θ = 8.0°F ig. 5: Frequency-streamwise location map of the magnitude of the wall pressure premultiplied spectra normalized by the local value of the variance (arbitrary scale). Abscissa: frequency in Hz, with Strouhal number StL of values 0.03 and 0.5 indicated. Ordinate: streamwise coordinate, with dimensionless streamwise coordinate X * values of 0 and 1 indicated.

Fig. 6 :

 6 Fig. 6: Premultiplied power spectra of streamwise momentum on the mean shock at elevation y/δ = 1.8: experiments, θ = 5.5°with the full-duration signal (10 s, bold solid), experiments, θ = 5.5°(dot) and LES, θ = 6.3°(dash) with signals of duration 125 ms.

  ) for the experiment at θ = 8°. It is recovered for the LES of interaction with separation at θ = 8°and θ = 9.5°in Figs.5(b) and 5(c), although the former case exhibits peaks spread on each sides of the St L = 0.03 line rather than a smooth bump centered at St L = 0.03 . This difference is due to the rather short duration of the θ = 8°c omputation with about 50 periods of the low frequency shock motion, as checked by computing spectra from experimental data restricted to a similar duration that exhibit similar uneven patterns.

III.

  Shock motions and regions of influenceA. Power spectrum of the shock positionsThe time-resolved data are used in first place to characterize the kinematics of the leading shock. The LES data are considered to track the position of the shock at each time. Its location is determined at every timestep and for every elevation by seeking for the maximum of the pressure gradient in the direction normal to the shock. The premultiplied spectra of shock positions at the normalized elevation y/H i = 2, where H i is the elevation of the crossing point of the incident and leading shocks, are plotted versus the Strouhal number in Fig.7.

2 Fig. 7 :

 27 Fig.7: Premultiplied power spectra of the shock streamwise location at elevation y/Hi = 2 normalized by the interaction length L; dash-dot: incipient separation (θ = 6.3°); solid: mild separation (θ = 8.0°); dash: full separation (θ = 9.5°).

in Figs. 8 (

 8 a) and 8(e) that similar properties are also found for the θ = 6.3°and θ = 9.5°flow deviations, respectively.Results obtained for intermediate frequencies are shown in figures 8(b), (d) and (f). It was tried to choose the same y/H i value for the reference point in all cases, large enough to have signals not directly polluted by the fluctuations of the shear layer. This led to the value of y/H i = 3 for the incipient case and for the 8.0°d eviation. In the 9.5°case, it was not possible to reach such high values because of the smaller relative extend of the computational domain; a compromise at y/H i = 1.8 was adopted to minimize the contribution of the turbulent layers. Finally, it was checked that the results are just weakly sensitive to choice of y/H i . When focusing on the mild separation case θ = 8.0°in Fig.8(d), two distinct periodic patterns are observable. The first pattern is found along the leading shock while the second pattern is located downstream of the foot of the expansion wave and spans the vortex shedding region (X * 0.7). More detailed information on these patterns can be gained by looking at the incipient and fully separated cases.Considering the θ = 6.3°incipient case, Fig.8(b) provides the visualization of a large part of the leading shock located above the crossing point with the incident shock located at elevation y/H i = 1. Consequently the alternated correlated and anti-correlated regions found along the leading shock in Fig. 8(d) are more clearly seen here. They are on the contrary barely seen in Fig. 8(f) corresponding to the θ = 9.5°case since they are mostly located above the domain on which the unsteady data have been sampled. The second correlation pattern, located in the shedding region downstream of X * ≃ 0.7 is also found for both the incipient and fully separated cases in Figs. 8(b) and 8(f), respectively. Note that the wavelength does not scale with L, ranging from λ ≃ 1.7L for the incipient case to λ ≃ 1.1L for the fully separated case. Low frequencies (St L ≤ 0.1), incipient interaction (θ = 6.3°) with the field of mean characteristics superimposed (white lines). Medium frequencies (0.3 ≤ St L ≤ 0.80), incipient separation (θ = 6.3°). The leading shock location is outlined with a dashed-dotted line. Low frequencies (St L ≤ 0.1), mildly separated interaction (θ = 8.0°) with the field of mean characteristics superimposed (white lines). Medium frequencies (0.3 ≤ St L ≤ 0.8), mildly separated interaction (θ = 8.0°). The leading shock location is outlined with a dashed-dotted line. Low frequencies (St L ≤ 0.1), fully separated interaction (θ = 9.5°) with the field of mean characteristics superimposed (white lines). (f) Medium frequencies (0.3 ≤ St L ≤ 0.8), fully separated interaction (θ = 9.5°) with the field of the Mach wave for Uc = 160 m/s superimposed in white lines. The leading shock location is outlined by a dashed-dotted line.

Fig. 8 :

 8 Fig. 8: Maps of cross-correlation between the band-passed pressure and the band-passed time series of the shock streamwise location at the elevation denoted by a star. The dashed line corresponds to the sonic line.

Fig. 9 :

 9 Fig. 9: Separation distance ξ vs optimal time delay τopt along the shock; incipient separation (θ = 6.3°): •; mild separation (θ = 8.0°): +; full separation (θ = 9.5°): ⋄. Black solid line and gray dotted line correspond respectively to velocities equal to 400 m/s and to 250 m/s.

α

  C h a r a c te r is ti c s (b) Radiation by moving sources: Mach waves.

Fig. 10 :

 10 Fig. 10: Propagation of pressure disturbances.

Fig. 11 :

 11 Fig. 11: A general scheme for characteristics and Mach waves.

Fig. 8 (C. Regions of influence for shock unsteadiness at intermediate frequencies 1 .

 81 Fig. 8(a).

  Theses paths nonetheless agree well with the correlation patterns of the shedding region at intermediate frequencies seen in Fig.8(b)-8(f), demonstrating the relevancy of the computation. However none of such radiated Mach wave appear to be able to induce disturbances on the leading shock. Moreover, the resulting propagation velocity along the leading shock above H e would be estimated at V PMW ≃ 290 ± 15 m/s for the θ = 6.3°case and V PMW ≃ 260 ± 10 m/s for the θ = 8.0°case. Both values are far from the value of the sole value V P ≃ 400 found in Sec. III C 2 for both cases. Consequently Mach wave radiation associated with shed vortices is probably not responsible for disturbance propagation.One has therefore to consider disturbance propagating from steady sources despite the fact that the intermediate frequency range is rather convective in nature. The propagation velocity induced by steady sources is estimated by Eq. 5 to V PC = 425 m/s and V PC = 415 m/s for the θ = 6.3°and θ = 8.0°cases, respectively. Theses values agree reasonably with the value of V P measured above H e but an admissible process based on steady sources taking origin from a mostly convective physics remains to be described.All characteristic lines that could induce disturbances on the upper part of the leading shock originate from the foot of the expansion fan, as seen in Figs.8(a), (c) and (e). It has been observed in movies computed from LES data pass-banded around St = 0.5 that this narrow region is very strongly modulated by the shedding of the vortical structures at the end of the mixing layer. It is then inferred that these almost periodic modulations may act as a stationary source of pressure fluctuations which then propagate along the characteristics of the expansion fan. This scenario implies that the frequency/Strouhal number of the shedding and the frequency/Strouhal number of the disturbances found along the leading shock are equal. The value of the Strouhal number St LP associated with the shock disturbances can be computed by:

  For the intermediate frequency range, perturbation of the shock location propagates along the shock with propagation velocities depending on the region under consideration. Joint analysis of the shock location-pressure cross-correlation maps and of the fields of characteristics and Mach waves leads to an identification of the sources inducing the motion of the leading shock. The displacement of the shock at low frequency is associated with the fluctuations of the interaction bubble whereas the motions at intermediate frequency are produced by the large vortical structures of the mixing layer. The analysis of the propagation of the perturbations along the characteristics and Mach waves has lead to an analytical determination of velocities in agreement with the different propagation velocities obtained along the leading shock at intermediate frequency.The kinematics of the leading shock is split into two parts depending on the location with respect to the crossing point between the shock wave and the expansion fan:• The first part is affected at low frequency by the motion of the beginning of the interaction region and at intermediate frequency by Mach wave radiation from the vortical structures produced by the mixing layer.• The second part of the shock is affected by disturbances issued from the foot of expansion fan for both the low and intermediate frequency bands. These disturbances are associated with the breathing of the interaction region for the low frequency range. At intermediate frequencies, the disturbances are created by the crossing of vortical structures shed from the mixing layer.

Table 1 :

 1 Geometrical parameters of the interactions.

	δ, mm	L, mm	H i , mm	Hi δ	H e , mm	He δ
	6.3°11	25	7.5	0.68	7.5	0.68
	8.0°11	38	12	1.09	23.5	2.14
	9.5°11	63	21	1.91	56	5.09

Acknowledgments

This work was granted access to the HPC resources of IDRIS under the allocation 2009-021877 made by GENCI (Grand Equipement National de Calcul Intensif).

though in that range extra care will have to be taken to validate the results. Different paths of propagation for pressure disturbances are possible. Simple cases are considered here.

If the source of perturbation is stationary in a flow at velocity U , the wave front grows at a speed equal to the sound speed C, and is convected at the speed U . It is known that it remains inside of a Mach cone, or dihedral in two dimensions, and the disturbance follows a characteristic line of slope α given by : α = sin -1 (C/U ) .

(

where U and C denote the flow velocity and the speed of sound, respectively. Note that the propagation of the disturbance results from the motion of the tangency point of the wave front represented in simple cases by the sphere or the circle convected at velocity U , and growing at a rate equal to sound speed, as sketched in Fig. 10(a). In such conditions the propagation velocity U η along this characteristics is equal to:

Taking into account the fact that for the present cases the angle between the characteristic lines and the shock is rather small, the propagation speed V PC of the disturbance along the shock can be approximated by:

When the source of perturbation moves with a velocity U c , radiation should be considered in a frame of reference moving at velocity U c . Such sources consequently radiate Mach waves of slope α ′ :

where M c is the convective Mach number, see for example [START_REF] Bailly | Subsonic and supersonic jet noise predictions from statistical source models[END_REF]. This is illustrated by Fig. 10(b).

Some comments can be made on this figure. At time t = 0, a perturbation is produced at the origin. The situation is examined at time 3∆t in the example given in Fig. 10(b), for which the wave front has traveled over a distance 3∆t U , and its radius is 3∆t C. Perturbations produced at intermediate times ∆t and 2∆t are

Appendix

It is possible to underline particular properties of low frequencies and to relate them to the mean field. Let us consider the continuity equation written for density and velocity fluctuation using Favre's averaging (the analysis could also be made on momentum equations and would lead to the same conclusions).