relation (A18) and/or (A19). This introduces a cut in the complex frequency-plane which leads to a power law when a Laplace transform formalism is used for solving the initial value problem. For strong shocks in the Newtonian limit the acoustic flow in the shocked gas is smaller than that of the incompressible flow. The linear dynamics of the front is mainly governed by the vorticity wave leading to normal modes that are neutral to leading order in this limit. The linear equation for the wrinkles of the shock front then reduces to a simple wave equation [START_REF] Clavin | Analytical studies of the dynamics of gaseous detonations[END_REF][START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF]. Another property of wrinkled shock waves is the formation of Mach stems (triple points). For a long time, they have been systematically observed on the shock front in experiments [START_REF] Briscoe | Experimental and theoretical study of the stability of plane shock waves reflected normally from perturbed flat walls[END_REF][START_REF] Lapworth | An experimental investigation of the stability of planar shock waves[END_REF], even for a weak wrinkling, but no particular attention was payed to these nonlinear effects. However, such triple points propagating in the transverse direction of the shock front are basic ingredients of the cellular structure of unstable gaseous detonations [START_REF] Clavin | WilliamsMultidimensional stability analysis of overdriven gaseous detonations[END_REF][START_REF] Clavin | Analytical studies of the dynamics of gaseous detonations[END_REF]. A model equation for the dynamics of shock waves, including the Mach stem formation, has been recently derived for strong shocks in the Newtonian limit [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF]. In this limit the nonlinear effects are mainly governed by the vorticity wave leading to the formation of singularities of the slope of the front, representative of Mach stems. The compressible effects that are skipped to leading order in the analysis are however essential for the formation of the secondary shock in Mach stems. But, in the limit used in the theoretical analysis, they are not the nonlinear mechanism responsible for the formation of singularity of the front slope. They are a consequence of this singularity and do not affect neither the geometry not the dynamics of the shock front, see fig. D1. This is no longer the case for weak shock wave or when the departure from unity of the specific heat ratio is not sufficiently small. The objective of this work is to carry out a numerical analysis of the model equation and to compare the results with direct numerical simulations of shockvortex and turbulence-vortex interactions and also with a preliminary experiment.

For the sake of completeness and simplicity the derivation of the model equation is outlined in physical terms in section 2. The technical details can be found elsewhere [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF]. Numerical studies of the model equation are presented in section 3. In this section analogies with the model for the formation of large-scale structure of the Universe are pointed out. The results are compared with direct numerical simulations of shock-turbulence interaction in three dimension geometry performed at Stanford University [START_REF] Larsson | Reynolds-and Mach-number effects in canonical shockturbulence interaction[END_REF][START_REF] Larsson | Direct numerical simulation of canonical shock/turbulence interaction[END_REF] and also of shock-vortex interaction in two dimensional geometry performed at Normandie University, see section 4. In section 5 a preliminary experimental study of a plane shock wave reflected from a wavy wall, performed recently at Aix-Marseille University, is presented showing a qualitative agreement with the numerical study of the model equation, but also quantitative differences.

Model equation

The shock wave is considered as a discontinuity separating two flows of a polytropic gas. The conditions upstream and downstream from the shock front will be identified by the subscripts 1 and 2 respectively. The flow of the upstream gas is prescribed and is not pertured by the shock front since the propagation velocity is supersonic, M 1 > 1. The flow in the shocked gas extends to infinity downstream (the supporting piston is at infinity and is assumed to have no influence on the shock dynamics). The dynamics of the shock front is obtained by solving the Euler equations in the shocked gas, using the Rankine-Hugoniot conditions at the front and a boundedness condition at infinity downstream. Using the same configuration as in Appendix A, u and w denote the longitudinal and transverse flow velocity. A two-dimensional geometry is consider for simplicity, extension to 3D is straightforward. From now on we introduce the notation, u → u + u, ρ → ρ + ρ, p → p + p, where an overbar identify the unperturbed flow in the referential frame of the shock wave, α = w 1 = w 2 = 0. The analysis is performed in the distinguished limit for which the compressed gas flow is strongly subsonic relative to the shock front,

M 2 ≡ u 2 /a 2 ≪ 1, M 1 ≡ u 1 /a 1 ≫ 1, M 2 1 (γ -1) = O(1). (1)
According to the Rankine-Hugoniot condition (A16) in Appendix A, M

2 2 ≈ (γ - 1)/2+1/M 2 1 ≪ 1,
and a perturbation analysis is carried out by using the Mach number of the shocked flow (relative to the shock front) as a small parameter

ǫ ≡ M 2 , M 1 = O(1/ǫ), ǫ ≡ M 2 ≪ 1, u 2 /u 1 ≈ ǫ 2 , u 1 u 2 ≈ a 2 2 , a 2 /a 1 = O(1). (2)
2.1 Linear analysis for a polytropic gas

Consider the normal mode analysis outlined in Appendix A for a general material. For a polytropic gas, to leading order in the limit ǫ → 0, the dispersion relation (A18) reduces to S 2 + 1 = 0. This corresponds to neutral modes of frequency ω ≈ a 2 |k| solution to the wave equation

∂ 2 α ∂t 2 -a 2 2 ∂ 2 α ∂y 2 =0, (3)
where a 2 is the sound speed in the compressed medium. The physical insight is as follows. In the Fourier representation (A4) and according to (A6), the ratio l/|k| is of order ǫ so that the acoustic waves of the eigenmodes are propagating in the transverse direction quasi-parallel to the front. According to the second term in the bracket [ ] of (A7) and (A8), the order of magnitude of the velocity of the acoustic wave is,

δu (a) = O ǫ 2 δp 2f /ρ 2 u 2 ,δ w (a) = O (ǫδp 2f /ρ 2 u 2 ) , (4)
where, according to (A12), δp 2f /p 2 = O(ǫ 2 αt /u 1 ), see also the second equation (A4). Using (2) and the relation p 2 ≈ ρ 2 a 2 2 , equations (4) yield

δu (a) = O(ǫ 2 αt ),δ w (a) = O(ǫa 2 α ′ y ) (5)
where the notations αt ≡ ∂α/∂t and α ′ y ≡ ∂α/∂y have been used, x = α(y, t) denoting the position of the perturbed front. Therefore the velocity of the flow of the acoustic wave ( 5) is smaller by a factor ǫ 2 than the flow velocity just downstream to the shock, given by the linearized Rankine-Hugoniot relations (A13),

x =0: δu 2 = δu 2f =α t ,δ w = δw 2f = u 1 α ′ y . (6)
To leading order in the limit ǫ → 0, the linear dynamics is thus controlled by the vorticity wave

δu (i) ≈ αt (y, t -x/u 2 ),δ w (i) ≈ u 1 α ′ y (y, t -x/u 2 ). (7)
The incompressibility condition ∂δu (i) /∂x + ∂δw (i) /∂y =0t h e ny i e l d st h ew a v e equation [START_REF] D'yakov | The stability of shock waves: Investigation of the problem of the stability of shock waves in arbitrary media[END_REF] where the relation u 1 u 2 ≈ a 2 2 in (2) has been used. The vorticity wave ( 7) is a shear flow quasi-parallel to the front, |δu (i) /δw (i) | = O(ǫ). A sketch of the flow associated with an eigenmode propagating in one transverse direction (simple wave) is sketched in fig. D1 left.

Weakly nonlinear analysis for a wrinkled shock in a quiescent medium

In the same conditions as before, acoustic waves are also negligible in the nonlinear dynamics of a weakly wrinkled shock propagating in a quiescent medium, as shown below. The Euler equations in the compressed gas can be written in the form,

∂u ∂t + u 2 ∂u ∂x = U - 1 ρ ∂p ∂x , ∂w ∂t + u 2 ∂w ∂x = W - 1 ρ ∂p ∂y , U ≡ δu ∂u ∂x + w ∂u ∂y ,W ≡ δu ∂w ∂x + w ∂w ∂y ,
where δu ≡ u -u 2 and U = W = 0 in the unperturbed flow. Mass conservation is added to these equations. A weakly nonlinear approximation is valid when the quadratic terms U and W are small compared to the unsteady terms. When expressed in terms of the vorticity wave [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF] and using the relation u 1 u 2 ≈ a 2 2 ,t h e quantities U and W take the form,

U ≈ 1 2 ∂H ∂x ,W ≈- 1 2 u 1 u 2 ∂H ∂y , (8)
where

H ≡ [-α2 t (y, t -x/u 2 )+a 2 2 α ′ 2 y (y, t -x/u 2 )].
The ratio of each quadratic terms in U and W to the linear terms of the Euler equations, evaluated from [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF], is of order

|α ′ y |/ǫ ≈|α t |/u 2 , (9)
where we have used [START_REF] Lapworth | An experimental investigation of the stability of planar shock waves[END_REF]. This quantity will be assumed to be small in the weakly nonlinear analysis, |α ′ y |/ǫ ≪ 1, so that the amplitude of the wrinkles are limited to small values such that |α ′ y |≈|α t |/a 2 ≪ ǫ. Then, a perturbation analysis is performed using the small parameter ε ≪ 1, defined such that |α ′ y |/ǫ = O(ε). To leading order, the analysis is simplified because the acoustic waves introduce terms that are even smaller than the quadratic terms, by a factor ǫ at least.

The source terms U and W in (8) satisfy the same equations as the linear vorticity wave, ∂U/∂t + u N ∂U/∂x = 0 and ∂W/∂t + u N ∂W/∂x = 0, so that they would introduce a secular contribution (a term growing linearly in time) to u and w if they were not balanced by a pressure term, as in the study of cellular detonations [START_REF] Clavin | Instability and nonlinear patterns of overdriven detonation in gases[END_REF][START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF]. However it is not worth carrying out such a complex calculation here. One only has to notice that the source terms U and W are zero for simple progressive waves, αt = ±a N α ′ y ⇒ H = 0. This corresponds to the fact that the shear flow associated with the progressive waves of normal modes is an exact solution of the incompressible Euler equations. More generally, for an initial disturbance of the front of finite size Λ, the source terms U and W cannot influence the front geometry after a short finite time, since the term H vanishes for t>Λ/a N , namely when the two progressive waves (propagating in opposite direction) no longer overlap. These terms do not play a significant role in the ultimate formation of cusps on wrinkles of small amplitude in the limit [START_REF] Briscoe | Experimental and theoretical study of the stability of plane shock waves reflected normally from perturbed flat walls[END_REF]. However this is no longer true for initial wrinkles whose amplitude is of the same order as the wavelength. This case is beyond the scope of the theoretical analysis. To leading order of the weakly nonlinear analysis, the Euler equations then reduce to their linear form and the flow takes the same form as in the linear analysis [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF], the pressure fluctuations being fully taken into account by the acoustic flow.

However, in the laboratory frame, the boundary conditions at the front (Rankine-Hugoniot relations) introduce quadratic terms leading to corrections of order ε to the dynamics [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF],

x = α : δu ≡ δu 2f (y, t) ≈ αt + u 1 α ′ 2 y ,w ≡ δw 2f (y, t) ≈ u 1 α ′ y , (10)
where the nonlinear terms that introduce corrections of a higher order (higher than ε) have been neglected. The shift of the front position, induced by the wrinkling, also introduces quadratic terms into the boundary values at the origin

x =0: δu ≡ u 2f (y, t) ≈ δu 2f -αu ′ x ,w ≡ w 2f (y, t) ≈ δw 2f -αw ′ x , (11)
where u ′ x (y, t) ≡ ∂δu (i) /∂x| x=0 , and w ′ x (y, t) ≡ ∂δw (i) /∂x| x=0 . Solving the vorticity wave with the boundary condition [START_REF] Clavin | Instability and nonlinear patterns of overdriven detonation in gases[END_REF], the incompressibility condition then leads to a nonlinear equation for the wrinkles in the form -(1/u 1 )∂u 2f /∂t+ ∂w 2f /∂y = 0. The nonlinear terms coming from the shift of origin in [START_REF] Clavin | Instability and nonlinear patterns of overdriven detonation in gases[END_REF] give no contribution. This can be shown by using the following relations

-u 2 ∂(αu ′ x )/∂t = -u 2 αt ∂δu (i) /∂x| x=0 + α∂ 2 δu (i) /∂x 2 | x=0 , ∂(αw ′ x )/∂y = α ′ y ∂δw (i) /∂x| x=0 + α∂ 2 δw (i) /∂x∂y| x=0 ,
since, according to the incompressibility condition, the sum of the two last terms is zero. When considering conditions for which H = 0, the sum of the two first terms is also zero since they give (1/2)∂H/∂t. Incompressibility then leads to a weakly nonlinear equation for the evolution of the front x = α(y, t),

|α ′ y |/ǫ ≪ 1: ∂ 2 α ∂t 2 -a 2 2 ∂ 2 α ∂y 2 + u 1 ∂ ∂t ∂α ∂y 2 =0, (12)
where the nonlinear term introduces a correction of order ε to the wave equation (3) coming from the Hugoniot condition in [START_REF] Larsson | Direct numerical simulation of canonical shock/turbulence interaction[END_REF]. By the way notice that, according to [START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF], the mean velocity of the wrinkled shock front is increased by a term u 1 < (∂α/∂y) 2 > where <.>denotes an average in the transverse direction and in time. This is not important here.

Introducing the non-dimensional quantities,

η ≡ y/Λ,τ = a 2 t/Λ,A (η, τ) ≡ α(y, t)/(ǫΛ) (13)
equation [START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF] takes the following non-dimensional, form free from parameter,

A ′ η ≪ 1: ∂ 2 A ∂τ 2 - ∂ 2 A ∂η 2 + ∂ ∂τ ∂A ∂η 2 =0, (14)
where (2) has been used. In this equation, A ′ η ≡ ∂A/∂η = α ′ y /ǫ is a small quantity, by definition of order ε, A ′ η = O(ε), see text below [START_REF] Larsson | Reynolds-and Mach-number effects in canonical shockturbulence interaction[END_REF]. Authors' name 2.3 Discussion of the dynamics of a weakly wrinkled shock front

When limiting the attention to simple waves for the wrinkles propagating in a single transverse direction at velocity a 2 , to leading order in the weakly nonlinear analysis, equation [START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF] takes the form of a Burgers' equation (with zero dissipation) for the slope of the front that can be written, according to (B2), in the form ∂α ′ y /∂t ′ +u 1 α ′ y ∂α ′ y /∂y = 0, see Appendix B. This indicates that equation ( 12) leads to the formation of singularities (discontinuity of α ′ y ) in finite time, representative of triple points as sketched in fig. D1 right.

In this weakly nonlinear analysis, the nonlinear term in ( 14) comes from the first Rankine-Hugoniot jump condition in [START_REF] Larsson | Direct numerical simulation of canonical shock/turbulence interaction[END_REF]. Except for a different numerical coefficient (1 instead of 2), this term has the same form as the one that would be introduced by a simple Huygens construction of a weakly wrinkled front propagating at a constant velocity in the normal direction. However the wave equation in the transverse direction describing the linear dynamics of a wrinkled shock front, see [START_REF] D'yakov | The stability of shock waves: Investigation of the problem of the stability of shock waves in arbitrary media[END_REF], cannot be obtained by a Huygens construction. But, the relaxation mechanism of initial disturbances, in the form of a power law coming from compressible effects, is absent from the analysis. It should be obtained by pushing the perturbation analysis to a higher order. This is a difficult task that has not yet be performed. The problem is different from (and more difficult than) weakly unstable detonations for which the acoustic damping is exponential [START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF].

By the way, notice also that equation ( 14) could be relevant to describe microseisms, see equations ( 111) and ( 112) in [START_REF] Longuet-Higgins | At h e o r yo ft h eo r i g i no fm i c r o s e i s m s[END_REF].

Model equation for shock-vortex and shock turbulence interaction

Consider a cylindrical and subsonic vortex of diameter Λ and turnover velocity v e /a 1 <ǫ . With this limitation the pressure fluctuations have a negligible effect upon the shock dynamics. The axis of the vortex is parallel to the planar shock wave and its velocity relative to the upstream gas is neglected. The analysis is performed below for a vortex of finite size, Λ, and can be extended to real vortices. The analysis can be summarized as follows [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF]. A small distortion of the shock front is formed during the short lapse of time of shock-vortex interaction τ int ≡ Λ/u 1 (period of crossover). After this short lapse of time, t>τ int , the shock propagates in a quiescent medium and its evolution is described by [START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF]. This subsequent evolution involves a time scale much longer than the interaction time in the limit (1-2), Λ/u 2 ≫ τ int . The problem takes the form of a two-timescale problem in which the first stage (period of crossover) determines the initial condition for the second stage controlled by [START_REF] Lapworth | An experimental investigation of the stability of planar shock waves[END_REF].

During the crossover (first stage for v e ≪ a 1 ), a quasi-planar pressure pulse of transverse extension Λ (∆y Λ) propagating downward in the longitudinal direction x at the velocity a 2 +u 2 , is generated in the compressed gas. This pressure pulse is generated by a longitudinal impulsive source of lifetime τ int , and will not influence the subsequent front evolution. Simultaneously the transmitted vortex takes the form of a shear wave quasi parallel to the front [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF]. The shock-vortex interaction is sketched in fig. D2. Denoting u 1 (x/Λ,y/Λ) and w 1 (x/Λ,y/Λ) the longitudinal and transverse component of the flow velocity of the vortex in the upstream gas, the perturbation of the upstream flow velocity, felt by the shock wave during the crossover (u 1f , w 1f )i s

u 1f (t/τ int ,y/Λ) = u 1 (-u 1 t/Λ,y/Λ),w 1f (t/τ int ,y/Λ) = w 1 (-u 1 t/Λ,y/Λ)
where the origin of time corresponds to the beginning of cross-over. It is then shown in [START_REF] Clavin | Nonlinear analysis of shock-vortex interaction Mach stem formation[END_REF] that the local disturbance of the shock front is given by an integration with respect to time of the following expression

v e /a 1 ≪ ǫ, 0 <t<τ int : ∂α/∂t ≈ 2(a 2 /u 1 )u 1f , (15)
leading to a wrinkle of small amplitude

|δα|/Λ ≈ 2a 2 v e /u 2 1 ( 16
)
where v e denotes the turn-over velocity,

|u 1 |/v e = O(1), |w 1 |/v e = O(1)
. The transverse velocity w 1 of the vortex does not affect the shape of the front. Its role is limited to feed the vorticity flow propagating in the shocked gas, see fig. D2.

In the spirit of the two-time scale approximation the subsequent dynamics of the shock front is given by the solution to [START_REF] Leveque | Clawpack Software[END_REF] for the following initial condition,

t =0: ∂A/∂t =0, |A| =2v e /u 1 , (17)
which corresponds to a small reduced amplitude, smaller than ǫ 2 . A composite solution, valid at any time is then obtained by combining ( 12) and ( 15)

|u 1f | a 1 ≪ ǫ : ∂ 2 α ∂t 2 -a 2 2 ∂ 2 α ∂y 2 + u 1 ∂ ∂t ∂α ∂y 2 =2 a 2 u 1 ∂u 1f ∂t . ( 18
)
The extension to three-dimensional geometry is straightforward. By using the nondimensional quantities in [START_REF] Golda | APPSPACK (Asynchronous parallel pattern search package)[END_REF], equation ( 18) can be written for a two-dimensional shock surface in the following non-dimensional form,

∂ 2 A ∂τ 2 -△A + ∂|∇A| 2 ∂τ = ∂ψ ∂τ where ψ ≡ 2 u 1f a 2 . ( 19
)
Expressed in terms of the dimensional units of space and time, the non dimensional source term takes the form

∂ψ/∂τ = 2(Λ/u 2 )∂u 1 /∂x| x=-u1t ,
which is a term of small amplitude, of order 2v e /u 2 ≪ 1, rapidly fluctuating on the reduced time scale a 2 /u 1 = O(ǫ) since, according to (2), a 2 /a 1 = O(1). This is true under the condition v e /a 1 ≪ ǫ used by the weakly nonlinear analysis in the limit [START_REF] Briscoe | Experimental and theoretical study of the stability of plane shock waves reflected normally from perturbed flat walls[END_REF]. Equation ( 19) can be used as a model for the geometry and dynamics of a two-dimensional shock front propagating in a frozen turbulent flow u 1 (x, y), Λ and v e denoting the integral scale and the turbulence intensity respectively.

Numerical simulations of the model equation

In the following subsections we will use clawpack software [START_REF] Leveque | Clawpack Software[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] (a well known program for hyperbolic equations) to simulate the model equations ( 14) and (19).

Numerical simulations without forcing term

Let us start with a numerical simulation of equation ( 14) in one dimension and without forcing term, using periodic boundary conditions in the transverse direction. A Authors' name typical result for the initial value problem with a sinusoidal initial shape and a zero initial derivative with respect to time is presented in fig. D3. The reduced length of the domain is 5 and the sinusoidal initial condition is τ = 0: A =0.3sin(2πη/5), ∂A/∂τ = 0. After a delay of a few reduced time a nonlinear pattern of standing wave is formed with the same wavelength as the initial condition but with crests, also called cusps (sharp transition of the slope) pointing toward the shocked gas and propagating in both transverse directions. The figure D3 corresponds to a final time of τ = 8 and a crest to crest amplitude of 0.5 compared to 0.6 at the initial condition. A stronger damping is observed for smaller wavelengths (not shown here), this could be influenced by the small viscous damping term introduced by the numerical code to treat the slope discontinuity. Further work has to be done in that direction. The result corresponds roughly to that of Burgers equation for the spatial derivative but with the difference that the transverse propagation described by the second order equation ( 14) is in both directions. The result of figure D3 is in direct correspondence with that obtained with the cellular detonation model [START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF] in which the wrinkles are triggered by a linear instability mechanism (absent here) in the vicinity of the threshold of a Hopf bifurcation at a finite wavelength. We now turn to simulation of [START_REF] Gilling | Tugen[END_REF] in two-dimensional geometry, still without forcing term, ψ = 0. We have already noted the analogy of ( 14) and [START_REF] Gilling | Tugen[END_REF] with the Burgers equation. The formation of large-scale structures in the Universe is currently analyzed via a multi-dimensional form of the Burgers equation [START_REF] Shandarin | The large-scale structure of the universe: Turbulence, intermittency, structures in a self-gravitating medium[END_REF][START_REF] Gurbatov | The Zeldovich approximation and the adhesion model[END_REF]. This analysis of mass distribution in the Universe is based on the Zeldovich approximation: each element of matter (including dark matter or anything else..) is moving in a straight line with a constant velocity which is a prescribed random field associated with small initial fluctuations in the early Universe, as they are "observed" from the cosmic microwave background. After the onset of singularities regions of high density form a disordered cellular structure called "foam" or "sponge" in the form of thin layers of high density delimiting large regions of low density. It is tempting to check that this property of coalescence, which is known for the solutions to the multidimensional Burgers equation with an initial random field, is also obtained with equation [START_REF] Leveque | Clawpack Software[END_REF]. The purpose here is not to perform a sophisticated and detailed analysis such as those extensively developed for a long time for the Burgers equation, see for example [START_REF] Vergalossa | Noullez Burgers'equation, Devil's staircases and the mass distribution for large-scale structures[END_REF]. Attention is limited here to point out a general property of the solution that can be relevant for the shock-turbulence interaction considered below. This property of solutions to [START_REF] Gilling | Tugen[END_REF] is illustrated in figures D4. The simulation is started with a white noise of small amplitude, |δA|≈0.01, on the 256*256 points of the computational square domain of size 5 with periodic boundary conditions. Figures D4 show the wrinkled front just after the beginning of the simulation and at time 4 respectively. Large cellular structures delimited by folds (also called crests or cusps), are observed at time 4. This illustrates the so-called "coalescence property" similar to that characterizing the solution to Burgers equation for random initial conditions.

Shock turbulence interaction

In this subsection a simulation of shock-turbulence interaction is presented using the model equation [START_REF] Gilling | Tugen[END_REF]. It would be expensive to simulate a three dimensional turbulence to be introduced into the model equation. This is not necessary here since our purpose is limited to point out a general geometrical characteristic of the patterns of folds that are formed on the front surface.

We simply use a synthetic turbulence generator [START_REF] Gilling | Tugen[END_REF], which produces an isotropic frozen velocity field with a k -5/3 cascade. The shock is then propagated at a con-stant velocity 1.5 through this frozen velocity field. Only the longitudinal component is used in equation [START_REF] Gilling | Tugen[END_REF]. The simulations are performed in a square box of length 5, using 256*256 points, with a turbulent kinetic energy 0.078, which correspond to ψ =0.396 and an integral scale of 0.15. A small integral scale is selected in order to point out the "coalescence property" already observed with a turbulent flame model [START_REF] Denet | Are small scale of turbulence able to wrinkle a premixed flame at large scale[END_REF]. The longitudinal component of the turbulent flow is shown in fig. D5a. Starting from a flat shock front, the two-dimensional field of the amplitude of the wrinkled front is shown at time 4 in fig. D5b. Strongly cusped folds defining a disordered pattern of cellular structures with a characteristic size larger than the integral scale are clearly observed. The front looks qualitatively similar to that observed in direct numerical simulations [START_REF] Larsson | Reynolds-and Mach-number effects in canonical shockturbulence interaction[END_REF][START_REF] Larsson | Direct numerical simulation of canonical shock/turbulence interaction[END_REF], except that no local quenching is observed, as it could have been anticipated since only small amplitudes are considered in [START_REF] Gilling | Tugen[END_REF]. However we were not able to find in the published DNS works [START_REF] Larsson | Reynolds-and Mach-number effects in canonical shockturbulence interaction[END_REF][START_REF] Larsson | Direct numerical simulation of canonical shock/turbulence interaction[END_REF] the data necessary to know whether or not the difference of length scale between the integral scale of the upstream turbulence and the typical size of the cellular structure of the wrinkled shock front is also observed in DNS.

In order to be more quantitative, the time evolution of the spectrum of the fluctuations of front position, relative to the mean front position (actually onedimensional spectrum averaged over samples at a fixed time), is shown in figure D6a. The spectra at 30 different times are plotted in this figure, showing the initial increase and the final saturation of the amplitude of the wrinkles. The energy spectrum of the upstream turbulence, E(k), is also plotted on the same figure. The comparison of the latter with the other spectra shows clearly that the characteristic size of the cellular structure of the shock front becomes larger that the integral scale of the turbulence as the time increases. An histogram of the front position of the last run (τ = 4) is shown in figure D6b. The skewness of the probability density function results from folds with crests pointing towards shocked gases. A similar skewness is observed in direct numerical simulations [START_REF] Larsson | Reynolds-and Mach-number effects in canonical shockturbulence interaction[END_REF].

DNS of shock vortex interaction

The interaction between a two-dimensional inviscid vortex and a relatively strong normal shock is analyzed by numerical simulation using the the Spectral Difference solver described in Appendix C. Compared to previous studies using (W)ENO-type reconstruction techniques and spatial discretization with 3rd/4th-order accuracy on structure Cartesian grids [21,[START_REF] Rault | Shock-wave-vortex interactions at high Mach numbers[END_REF], the numerical solver used here allows higherorder (7th-order) computations over unstructured meshes and features a sub-cell shock capturing approach, see Appendix C. The rectangular computational domain of dimensions 2H ×H consists of 200×80 quadrilateral elements, see fig. D7. Setting the order n equal to 7, the total number of degrees of freedom (DoF) is 784 000. In the region where the shock is located, namely, in the middle of the domain, the mesh is refined in the horizontal direction by means of the following stretching function (mirrored about the centerline):

x(η) H = tanh(C str η) tanh(C str ) -1, (20)
with η ∈ [0 : 1] and the stretching coefficient C str =1.5. Inflow (prescribed density and velocity) and outflow (prescribed pressure) conditions are set at the left and the right boundaries, respectively, whereas periodic conditions are used at the top and bottom boundaries. The choice of periodic boundaries is justified by the fact that the present analysis is only concerned with the short time interval during Authors' name which the interaction between the vortex and the shock takes place. Within this time-frame, the perturbations emanating from the shock-vortex interaction region and traveling in the vertical direction do not have enough time to cross the periodic planes and perturb the region of interest.

The initial conditions and the prescribed boundary values are determined according to the Rankine-Hugoniot jump conditions (A17) such as to obtain a Mach 2 stationary shock, M 1 = 2. The static pressure p 2 , namely, the outflow pressure, and M 1 = 2 are the only prescribed quantities. In order to establish a well resolved shock for the mesh used, a precursor computation without the vortex is initially performed and let run until any initial perturbation issuing from the shock region has left the domain. The established shock region is depicted, greatly magnified within a region of width 0.04H, in fig. D8. The comparison between the contours of the Mach number (top picture) and the artificial viscosity (bottom picture) clearly shows that subcell resolution is attained, as expected for such a high-order computation. Also note the good level of localization provided by the sensor, which is clearly switching on the artificial viscosity only inside the element where the shock is located (the linear smoothing algorithm is obviously spreading the artificial viscosity to the closest neighbour elements).

Once the stationary shock is established, a homentropic Taylor vortex [START_REF] Taylor | On the dissipation of eddies,T h eS c i e n t i fi cP a pe r so fS i rG e o ff r e yI n g r a mT a y l o r[END_REF] is placed in the middle of the left region at five vortex radii upstream of the shock using the following relations:

ρ(x) ρ 1 γ-1 = p(x) p 1 γ-1 γ =1- γ -1 2 M 2 v exp(1 -ξ 2 v -η 2 v ), (21)
u 1 -u(x) u 1 η v = w(x) u 1 ξ v = M v M 1 exp 1 -ξ 2 v -η 2 v 2 , (22)
where M v represents the vortex Mach number whereas ξ v =(x -x v0 )/R v and η v = (y -y v0 )/R v are the spatial coordinates in the vortex reference frame normalized with respect to the vortex radius R v =0.1H (x v0 and y v0 are the coordinates of the vortex center at the beginning of the simulation). Two different values of the vortex Mach number are analyzed, namely, M v =0 .4 and 0.8. These values are outside the validity domain of ( 14) and (18-19) (strong shock and weak vortex). However as seen thereafter, for M v =0.4, a good agreement is observed between the shock surface evolution modeled by [START_REF] Leveque | Clawpack Software[END_REF] and DNS. The time history of the contours of the density in the region of interest is depicted in fig. D9 for the M v =0 .8 test. According to the classification presented in [21], although both values of the selected shock Mach numbers fall within the region of strong interactions (with a Mach reflection pattern), the former is relatively close to the weak interaction region, whereas the latter is well within the strong interaction region. Indeed, by observing the computed Schlieren plots depicted in figures D10 and D11, a very marked Mach reflection pattern can be seen in the M v =0.8 test case, whereas, for the lower value of M v , the two cusps of the front are still present but are extremely weak and barely visible in fig. D10d. In order to make distinction between the two test cases, the computations performed with M v =0.4 and 0.8 are referred in the following to as moderate and strong shock-vortex interaction tests, respectively. The time evolution of the shock profiles for the moderate and strong interaction tests are depicted in fig. D12, along with the numerical solutions of the model equation ( 14), obtained with a sixth order PADE scheme [START_REF] Lele | Compact finite difference schemes with spectral like resolution[END_REF] adding a small amount of second order artificial viscosity to handle the cusps. The initial conditions used for solving [START_REF] Leveque | Clawpack Software[END_REF] are the profiles extracted from the two-dimensional simulations using the maximum pressure gradient as a marker of the shock location. Each profile is translated along the x-axis according to the corresponding normalized time. Note that the discontinuity in the three profiles between a time of 0.2 and 0.4 for M v =0 .8 (fig. D12b) is due to the particular post-processing algorithm which identifies the shock at the location of the maximum pressure gradient (in the norm). This extremum-one for each value of y/H-is detected in this case on the relatively strong reflected shock, see fig. D11. The two cusps at the front where the Mach stem is formed can be clearly seen in the plots. The velocity in the laboratory frame of the upward moving cusp is approximatively 0.72 a 2 for M v =0.4 and 0.66 a 2 for M v =0 .8. This should be analyzed by taking into account the downward component of the velocity of the transmitted vortex which is not negligible in the numerical simulations in contrast to the limit used in the theoretical description. The two cusps, one propagating upward and the other downward, are observed. The agreement of equation ( 14) with the full simulation is satisfactory for the M v =0 .4. The M v =0 .8 case is too far from the conditions used in the theoretical analysis for providing a satisfactory agreement.

Experimental study of a plane shock wave reflected from a wavy wall

In order to test the model equation ( 14) and more particularly the solution for an initial condition given by sinusoidal shape of the front, an experiment of shock reflexion on a wavy wall was recently performed in a shock tube by the group at IUSTI directed by Lazhar Houas and Georges Jourdan. A similar experiment was performed a long time ago for studying the linear relaxation power law of the amplitude of the front [START_REF] Briscoe | Experimental and theoretical study of the stability of plane shock waves reflected normally from perturbed flat walls[END_REF]. The attention is focused in the IUSTI experiment on the trajectory of the triple points. The facility is described elsewhere [START_REF] Jourdan | An e wv a r i a b l ei n c l i n a t i o n shock tube for multiple investigations[END_REF]. An incident planar shock wave propagating at 514 m/s in air (Mach number 1.5) is reflected from a wavy wall having a one dimensional sinusoidal shape (amplitude 1 mm, wavelength 2 cm), see the three pictures of fig. D13 left. The propagation Mach number of the reflected shock is 1.43 corresponding to 325 m/s for the velocity in the laboratory frame and to 536 m/s for the propagation velocity relative to the upstream gas put in motion at 239 m/s by the incident shock before reflexion. A high speed camera (125 000 frames/second) is used to record the Schlieren pictures of the reflected shock. The three pictures of fig. D13 left are taken at 120 µs, 200 µs and 280 µs respectively after reflexion. The structure of the flow and of the shock front just after reflexion (not shown here) is relatively complex, as the flow features that arise from the impact of shock wave on a concave cavity [START_REF] Skews | KleinFlow features resulting from shock wave impact on cylindrical cavity[END_REF]. Few time after reflexion, 20 µs, typically the front shape takes a strongly cusped and periodic cell structure with crests (cusps) pointing toward the concave part of the wavy wall. At approximately 50 µs that is about 1.7 cm from the wall, these cusps split into two triple points, one propagating upward in the transverse direction and the other downward. A quite similar feature is also observed concerning the singularities of slope at the early times of the numerical study of the model equation ( 14) starting with a sinusoidal shape and a zero derivative with respect to time. After less than 100 µs, at a distance about 3.4 cm (1.7 times the wavelength) from the wall, the front is quasi-planar with two arrays of regularly spaced triple points moving transversally at constant velocity in the opposite directions. The trajectories of these triple points form a diamond-shaped pattern shown in fig. D13 right. This result is qualitatively similar to the above numerical results obtained with the model equation, see fig. D3 and also to markings left by cellular detonations on soot-coated foils at the walls [START_REF] Clavin | WilliamsMultidimensional stability analysis of overdriven gaseous detonations[END_REF][START_REF] Clavin | Analytical studies of the dynamics of gaseous detonations[END_REF][START_REF] Clavin | Diamond patterns in the cellular front of an overdriven detonation[END_REF].

However there is a quantitative discrepancy, the transverse velocity of the triple points in the experiment is about 310 m/s that is smaller than the sound speed in the shocked gas of the reflected shock wave, evaluated to be 446 m/s by using the Rankine-Hugoniot conditions. Assuming that the shocked gases are quiescent (at rest), this subsonic velocity is not compatible with the propagation of the secondary shock of Mach stems. Further experimental investigations are required to decipher the nature of the triple points.

Conclusions

The model equations ( 14) and ( 19) derived for strong shocks in the Newtonian limit (1) and for very subsonic disturbances of the upstream flow capture important features of the dynamics of shock waves, notably concerning the formation of singularities of the slope of the shock front and the two-dimensional network formed by the corresponding folds of the shock surface. However the front disturbances are weak for such conditions. Quantitative discrepancies with experiments and DNS, performed in more realistic conditions, exist and require further investigations.

Among the questions to be addressed the two following ones are worth mentioning. Is it true in real situations, that the typical length scale of the wrinkled front of a shock wave propagating in a turbulent flow becomes larger than the integral scale, as it is predicted by equation [START_REF] Gilling | Tugen[END_REF]? What is the nature of the triple points observed in the two-dimensional experiments of wrinkled shock fronts ?

According to (A16) for γ 5/3, equation (A18) has two pairs of imaginary roots (neutral normal modes), S = ±iΩ, with Ω > 1, so that the right hand side of (A18) is negative. Therefore the ± sign in the expression ±i √ Ω 2 -1=± √ 1+S 2 in the left hand side of (A18) should be the same as that of S = ±iΩ, to give

-2M 2 Ω Ω 2 -1=-Ω 2 1+1/M 2 1 +1< 0. (A19)
Then, according to (A6), each Ω corresponds to a single value l>0, il ± = ±il

l |k| =   M 2 Ω+ √ Ω 2 -1 1 -M 2 2   > 0, Ω ≡ ω a 2 |k| 1 1 -M 2 2 > 1.
(A20)

A.3 Non-radiating condition

The reduced frequency ω/a 2 |k| > 0 can also be expressed in terms of l/|k| and M 2 by solving the quadratic equation for Ω 2 obtained from the first equation (A20). This yields two solutions Ω 2 (1 -

M 2 2 )=(1+M 2
2 )l 2 /k 2 +1 ± 2(l/k)M 2 l 2 /k 2 + 1, that can be written ω 2 /(a 2 2 k 2 )= l 2 /k 2 +1± M 2 l/k 2
, corresponding simply to the solutions of the wave equation (A5). From these solutions, only the one with the minus sign has to be retained

ω = a 2 l 2 + k 2 -u 2 l, (A21)
This is because the solution to (A19) should satisfy the inequality

ω<a 2 √ l 2 + k 2 yielding (1 -M 2 2 )Ω 2 -1 <l 2 /k 2 . The latter is checked for Ω > 1b yi n t r o d u c i n g the expression (A19) for 2M 2 Ω
√ Ω 2 -1i n t ol 2 /k 2 obtained from (A20). Introducing the wave vector κ ≡ le x +ke y where e x and e y are the unitary vectors along the axes Ox and Oy respectively, the frequency in (A21) takes the form ω = κ.[a 2 κ/ √ l 2 + k 2 -u 2 e x ], so that, in the reference frame of the unperturbed (planar) shock, the component along the x-axis of the propagation velocity of the sound wave is

u 2 -a 2 l/ √ l 2 + k 2 . The non-radiating condition, u 2 -a 2 l/ √ l 2 + k 2 < 0, can be written M 2 < (l/|k|)/ (l/k) 2 + 1 that is (l/|k|) > M 2 /(1 -M 2 2 ) 1/2 ,
a condition which is satisfied by the expression of l/|k| in (A20) since Ω > 1. The fronts with non-radiating neutral modes are stable [START_REF] Kontorovich | Conserning the stability of shock waves[END_REF], in contrast to those characterized by sound emission since they correspond to an unbounded reflexion coefficient of acoustic waves emitted in the shock gases.

Appendix B. Burgers equation

Using reduced amplitudes of order unity, Ã ≡ A/ε, equation ( 14) can be written

à = O(1) : ∂ 2 à ∂τ 2 - ∂ 2 à ∂η 2 + ε ∂ ∂τ ∂ à ∂η 2 =0, (B1)
showing the two-time-scale nature of the problem, a wave equation with both frequency and wavelength of order unity and a slowly varying amplitude on the time scale τ ′ ≡ ετ . When limiting the attention to simple waves for the wrinkles propagating in a single transverse direction at velocity, for example in the direction of negative η, the solution to ( 14) is looked for in the form Ã(η, τ, ετ)=A(η +τ, τ ′ ) and ∂/∂τ → ∂/∂η + ε∂/∂τ ′ . When the terms of order ε 2 are neglected in front of those of order ε, equation (B1) reads 2ε∂ 2 A/∂η∂τ ′ + ε∂(∂A/∂η) 2 /∂η = 0, yielding

∂A ′ ∂τ ′ + A ′ ∂A ′ ∂η =0, (B2)
where A ′ ≡ ∂A/∂η.

Appendix C. The numerical scheme for shock-vortex interaction

The shock-vortex interaction numerical tests are conducted by solving the unsteady two-dimensional Euler equations for gas dynamics:

∂U ∂t + ∂F i ∂x i = 0, (i =1, 2), ( C1
)
where t represents the time, x =(x 1 ,x 2 ) the position vector, U = ρ ρu 1 ρu 2 ρe T is the vector of conservative variables, and F i is for the flux vector

F i =     ρu i ρu 1 u i + δ 1i p ρu 2 u i + δ 2i p (ρe + p)u i     . (C2)
In the above equations, ρ is the fluid's density, u =(u 1 ,u 2 ) is the velocity vector, δ ij is the Kronecker delta, p is the pressure and e is the total energy (internal + kinetic), which is related to the pressure assuming the validity of the ideal gas law,

ρe = p γ -1 + 1 2 ρu i u i , (C3)
where γ = c p /c v is the ratio between specific heat capacities at constant pressure and volume.

The above equations are solved using the high-order Spectral Difference (SD) method for unstructured spatial discretizations. The formulation of the equations on quadrilateral grids is similar to the one described in [START_REF] Sun | High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids[END_REF]. In the SD method, in particular, the solution in each element is reconstructed from n solution points using a degree n -1 polynomial for each coordinate direction. Likewise, degree n polynomials are used to reconstruct the fluxes from n+1 flux points. In the present study, the solution points in 1D are chosen to be the Gauss-Legendre quadrature points, whereas the flux points are selected to be the Gauss-Legendre quadrature points of order n-1 plus the two end points -1 and 1. This particular combination of solution and flux points, in fact, can be proved to be linearly stable for all orders of accuracy and optimal in reducing aliasing errors and providing good conditioning for the non-linear case [START_REF] Jameson | A proof of the stability of the spectral difference method for all orders of accuracy[END_REF][START_REF] Jameson | On the non-linear stability of flux reconstruction schemes[END_REF]. The reconstructed fluxes are only element-wise continuous, but discontinuous across cell interfaces. Therefore a Riemann solver is employed to compute a common flux at cell interfaces to ensure conservation and stability. In the current implementation, the Roe solver with entropy fix [START_REF] Ro | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] is used. The left and right states here represent the solution on both sides of the shared edge flux point.

High-order methods are highly susceptible to overshoots and undershoots in the neighborhood of discontinuities. These spurious oscillations can be damped using artificial viscosity. Note however that, in order to avoid detrimental effects on the solution accuracy, the effects of artificial viscosity should be localized to a small part of each of the cells through which the shock passes. In the present study, the method proposed in [START_REF] Persson | Sub-cell shock capturing for discontinuous Galerkin methods,A I A AP[END_REF]33] is adopted. This approach, which was originally developed for the Discontinuous Galerkin (DG) scheme, combines a highly selective spectral sensor, based on the modal decomposition via orthogonal polynomials, with a consistently discretized artificial viscosity. The method smooths out the discontinuities in the solution to a width that is appropriately resolved by the mesh and the polynomial approximations, thus obtaining subcell resolution for high-order discretizations. Moreover, if p = n -1 is the polynomial order of the solution, the amount of artificial viscosity is set to scale as 1/p and, for high order polynomials, shock profiles which are much thinner than the element size can be obtained. The method is briefly summarized below.

The solved equation set is augmented with a Laplacian diffusion term,

∂U ∂t + ∂F i ∂x i = ∂ ∂x i ǫ ∂U ∂x i , ( C4
)
where ǫ represents the added artificial viscosity to be introduced in regions where the flow is under-resolved, such as across shocks and contact discontinuities. In order to sense these regions accurately, the level of smoothness of a selected variable ψ, hereafter referred to as the sensor variable, is determined by checking the rate of decay of the expansion coefficients, or modal coefficients, of the solution in an orthogonal basis. In the present case, normalized orthogonal Legendre polynomials are adopted and the sensor variable is written as

ψ(ξ)= n i=1 ψi P * i-1 (ξ), (C5)
where P * i (ξ) is the normalized Legendre polynomial of degree i

P * i (ξ)=P i (ξ) 2i +1 2 ⇒ 1 -1 P * i (ξ)P * j (ξ)dξ = δ ij , (C6)
and ψi are the modes of the sensor variable, which can be obtained from its nodal values through a simple matrix multiplication with the inverse Vandermonde matrix constructed using the selected polynomial basis [START_REF] Hesthaven | Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications,S p r i n g e rS c i e n c e + B u s i n e s sM e d i a[END_REF].

After introducing the truncated expansion

ψ(ξ)= n-1 i=1 ψi P * i-1 (ξ), (C7)
a resolution indicator [START_REF] Persson | Sub-cell shock capturing for discontinuous Galerkin methods,A I A AP[END_REF]33], or sensor, inside each element is defined as

s e = log 10 (ψ -ψ, ψ -ψ) e (ψ, ψ) e , (C8)
where (•, •) e is the standard inner product in the element. From Eqs. (C5) and (C7), the orthonormality of the Legendre polynomial basis (cf. Eq. (C6)), in particular, leads to the following simple expressions for the inner products in Eq. (C8):

(ψ -ψ, ψ -ψ) e = 1 -1 ψ2 n P * n-1 (ξ)P * n-1 (ξ)dξ = ψ2 n , (C9) (ψ, ψ) e = 1 -1 n i=1 n j=1 ψi ψj P * i-1 (ξ)P * j-1 (ξ)dξ (C10) = n i=1 n j=1 ψi ψj 1 -1 P * i-1 (ξ)P * j-1 (ξ)dξ = n i=1 n j=1 δ ij ψi ψj = n i=1 ψ2 i . (C11)
Hence, the artificial viscosity is evaluated as

ǫ e =      0 for s e <s 0 -κ, ǫ0 2 1+sin π(se-s0) 2κ for s 0 -κ ≤ s e ≤ s 0 + κ, ǫ 0 for s e >s 0 + κ, (C12)
where ǫ 0 , s 0 and κ are additional parameters to be chosen empirically. In the present implementation, ǫ 0 is computed from the spectral radius of the flux Jacobian and the element size h as

ǫ 0 = λ max h/(n -1), (C13)
where λ max = max( √ u i u i + a) is the maximum wave speed in the whole domain and a is the speed of sound [START_REF] Barter | Shock capturing with PDE-based artificial viscosity for DGFEM: Part I. formulation[END_REF].

With regards to the threshold s 0 and sensor tolerance κ, these are computed via an auto-calibration algorithm, which, through a manufactured solution, is able to determine the optimal values for the shock-sensor parameters for an arbitrarily selected value of n. In particular, s 0 is computed from the value of the sensor when a relatively sharp hyperbolic tangent profile is assumed for the solution. The relevant nodal values take the form

ψ m.f. (ξ i )= 1 2 1 + tanh tr(i -1 -n/2) + 1/2 (n -1)/5 , (i =1,...,n), (C14)
where tr(•) represents truncation to the lower integer and is used to prevent the manufactured signal from being perfectly centered in the element of odd order n, which is particularly undesired for low orders (for n = 3, for example, this would make the manufactured signal linear and cause the threshold s 0 to be too low and the sensor to trigger artificial viscosity in smooth regions of the flow). If s m.f. is the value of the sensor for the manufactured solution, computed from Eq. (C8), then s 0 = s m.f. -3, and κ =1. (C15)

Although the described calibration algorithm yet involves some level of empiricism on the determination of the artificial viscosity parameters-namely, the correction to s m.f. in Eq. (C15) and the value of κ itself-the resulting values were found to work well in a variety of different test cases with the order n ranging from 3 to 9.

A particularly desired feature of this strategy effectively appears to be the removal of the order-dependency of the optimal parameters.

u 2 a 2 + u 2 l a 2 Figure D2 .

 2D2 Figure D2. Sketch of the shock-vortex interaction for a strong shock and a very subsonic vortex, in the Newtonian limit, see equation[START_REF] Briscoe | Experimental and theoretical study of the stability of plane shock waves reflected normally from perturbed flat walls[END_REF].

Figure D3 .

 D3 Figure D3. One dimensional simulation in a domain of length 5, starting from an initial condition 0.3 * sin(2 * π * x/5). Note the periodic diamond-shaped pattern formed by the trajectory of cusps (crests) pointing towards shocked gases (a single cell is shown). The shock propagates downwards. The downward propagation velocity is arbitrary in this picture.

Figure D5 .

 D5 Figure D5. Shock-turbulence interaction described by the model equation (19). (a) Fluctuations of the longitudinal component of the upstream frozen velocity field in a cross section perpendicular to the direction of propagation of the shock. (b) Fluctuations of the front position at time τ =4 .N o t et h a tt h et y p i c a l length scale of the size of the wrinkles of the shock front is larger than the integral scale of the upstream turbulence.

Figure D6 .

 D6 Figure D6. Statistical properties of the solution presented in figure D5. (a) Spectra of the fluctuations of the front position at different times. The wavenumber is plotted on the horizontal axis and the amplitude on the vertical axis. The amplitude increases with time. The curve with the largest amplitude corresponds to time τ =4f o rw h i c has a t u r a t i o ni so b s e r v e d .T h e s es p e c t r aa r ec o m p a r e dt ot h es p e c t r u mo ft h e upstream turbulent flow E(k)( l o w e rc u r v e )s h o w i n gt h ed i ff e r e n c eo fl e n g t hs c a l e sb e t w e e nt h ew r i n k l e d front and the upstream turbulent flow. (b) Histogram at time τ =4o ft h ea m p l i t u d eo ft h ew r i n k l e so f the front, showing a skewness which is due to the existence of lines of crest (folds, cusps).

Figure D7 .Figure D8 .

 D7D8 Figure D7. Computational grid: 200 × 80 quadrilateral elements.

8 Figure

 8 Figure D12. Time evolution of shock profiles. Bold line: DNS. Line: Solution of Eq. (14).

  Figure D13. Shock tube experiment of shock reflexion from a wavy wall (one-dimensional wrinles). Propagation in air with a Mach number 1.5 and 1.43 for the incident shock and the reflected shock respectively. The reflected shock propagates from left to right. Left: Schlieren pictures at 3 different times after reflexion: 120 µs, 200 µs, 280 µs. Triple points moving in opposite transverse directions are clearly seen on the reflected shock. Right: Diamond shaped pattern of the trajectories of the triple points

  Figure D1. Left: Flow of shocked gas associated with a progressive normal mode of a wrinkled shock front propagating into a quiescent gas (2-D geometry). The picture shows the upward propagation of the progressive wave associated with a harmonic disturbance of the front. The shear flow of the vorticity wave is in red. The non-radiative acoustic wave is in blue. To leading order the acoustic wave propagates in ad i r e c t i o nq u a s i -p a r a l l e lt ot h es h o c kf r o n ta n dt h ep h a s ev e l o c i t y ,ω/k,o ft h ew r i n k l e so nt h ef r o n ti s equal to the sound speed in the shocked gas a 2 .R i g h t :S k e t c ho ft h eM a c hs t e m sf o r m e db yt h en o n l i n e a r evolution of the progressive disturbance propagating upwards.

	!"#$%&'$"()
			u 2
			a2
	u 1	u 2 u 2	% + " ; -.1 + '7 % 8 ) ! " 3 :# % ( 1% .1 8 ) /	-) + " 3 ( % #6 / -, " + <
				=&1 2/ &13 )
		u 2	*%+,'-.)$/
	5" #.1+ 1.6 '7 %8 )/ 0-, )% #'9 "7 4	0.#12&)'2"13.4
	y
		weak acoustic pulse
	t a 2 vortex initial u 1	s h e a r w a v e	burst acoustic strong vortex transmitted	x
		u 2 t	a 2 t

Acknowledgement

Dr Lazhar Houas and Georges Jourdan are acknowledged for their guidance in the shock tube experiment.

Appendix A. Dispersion relation

Consider a planar shock wave propagating in the negative x direction at a constant (supersonic) velocity u 1 in a uniform medium, M 1 ≡ u 1 /a 1 > 1, where a denotes the sound speed. This requires an external device to trigger the shock, for example a piston at constant velocity in the shocked gas. We assume that such a piston is at infinity and does not generate disturbances. In the reference frame of the unperturbed planar solution, the unperturbed front stands perpendicular to the x-axis at x = 0 and the shocked material flows at a constant (subsonic, ) velocity

represent the perturbed position of the shock front at transverse position y and at time t. In order to save the notations we used a 2D geometry, the extension to 3D is straightforward. The upstream medium is unperturbed and is assumed to be uniform. The Euler equations are solved in the shocked material. The boundary conditions are given at x = 0 by the linearized Hugoniot relations and a boundedness condition is enforced downstream in the shocked material for x →∞. For any physical quantity f we introduce the decomposition f = f + δf,wheref represents the unperturbed solution. Vorticity and entropy are generated in the compressed material at the shock and convected downstream in the x-direction at the velocity u 2 . The flow is thus decomposed into a vorticity wave (superscrit i) and an acoustic wave (superscript a)

where u and w denote the longitudinal and transverse component respectively and where the pressure fluctuations are fully taken into account by the acoustic flow (δu (a) , δw (a) ). The flow (δu (i) , δw (i) ) is incompressible, ∂u (i) /∂x + ∂w (i) /∂y = 0. For an harmonic perturbation of the front position, a normal-mode decomposition is used,

where k is the transverse wave vector (a real quantity), and σ = s ± iω (ω>0) is a complex number whose real part s is the linear growth rate (or damping rate if s<0) and the imaginary part ω is the frequency of oscillation in time.

A.1 Analysis for a general material

The acoustic waves propagating in the shocked material, solution to the equation of linear acoustic take the form

where the quantities l ± are solutions to the second-order algebraic equation,

.

(A6)

The constant of integration p2f is given by the linearized Rankine-Hugoniot condition at x = 0. Its expression in (A4) is given for simplicity for a polytropic gas where γ is the ratio of specific heats, see the general relation below in (A12). Boundedness condition at x →∞implies Re(il ± ) < 0. The vorticity wave is then computed by using the boundary condition at x = 0 obtained by subtracting the acoustic flow from the flow given by the Rankine-Hugoniot condition (ũ 2f ,w 2f )

The incompressibility condition of the entropy wave

then leads to an equation for σ/(a 2 |k|)

where S and the ± sign are the same as in (A6) and where the following relation

has been used. Introducing the non-dimensional parameters r and n characterizing the material and the strength of the shock wave,

where dp/dρ -1 | 2f is the slope of the Hugoniot curve at the state of the shocked material in the plane (p,1 /ρ), the Rankine-Hugoniot conditions yield