B Maxwell
email: brian.maxwell@case.edu.

J Melguizo-Gavilanes

Origins of instabilities in turbulent mixing layers behind detonation propagation into reactive-inert gas interfaces Origins of instabilities in turbulent mixing layers behind detonation propagation into reactive-inert gas interfaces

published or not. The documents may come

Origins of instabilities in turbulent mixing layers behind detonation propagation into reactive-inert gas interfaces

Brian Maxwell, Josué Melguizo-Gavilanes

I. INTRODUCTION

Studying detonation-interface interactions is relevant to industrial safety where upon accidental leaks/releases, layers of fuel and inert gas may form due to differences in molecular weight of the gases and/or insufficient mixing. In detonationbased propulsion applications, such as rotating detonation engines this flow configuration is also present. For example see Refs. 1 and 2. Detonation propagation in uniform reactive mixtures has been studied extensively. Gaseous detonations exhibit a characteristic cellular structure, which effectively describes the motion and collisions of transverse waves passing along the wave front forming triple points. Traditionally, detonations are characterized as having either a regular or an irregular structure. Regular detonations, typical of fuel mixtures involving high levels of argon dilution, have very structured patterns with cell widths (λ ) that can be unambiguously determined. [START_REF] Strehlow | Gas phase detonations: Recent developments[END_REF][START_REF] Austin | The role of instability in gaseous detonation[END_REF][START_REF] Pintgen | Direct observations of reaction zone structure in propagating detonations[END_REF] Irregular detonations, typical of fuel mixtures involving low levels of argon dilution, tend to exhibit much more stochastic-looking structures where various length scales are present. In general, such irregular detonations tend to exhibit flow fields that appear turbulent. [START_REF] Austin | The role of instability in gaseous detonation[END_REF][START_REF] Subbotin | Collision of transverse detonation waves[END_REF][START_REF] Oran | A study of detonation structure: The formation of unreacted gas pockets[END_REF][START_REF] Radulescu | The ignition mechanism in irregular structure gaseous detonations[END_REF][START_REF] Kiyanda | Photographic investigation into the mechanism of combustion in irregular detonation waves[END_REF][START_REF] Maxwell | Influence of turbulent fluctuations on detonation propagation[END_REF] In these cases, the presence of rapidly spreading shear layers, vortex rollup, and unburned pockets of reactive gas contribute towards the observed irregularity, arguably resulting in detonation cell widths that do not remain constant. Detonation propagation into non-uniform mixtures, however, has been studied much less extensively.

Up to now, re-initiation of detonation through a gap of inert gas, [START_REF] Bone | A photographic investigation of flame movements in gaseous explosions VII -the phenomenon of spin in detonations[END_REF][START_REF] Bull | Initiation of unconfined gas detonations in hydrocarbon-air mixtures by a sympathetic mechanism[END_REF][START_REF] Lean | Experiments on the transmission of explosions across air gaps[END_REF][START_REF] Teodorczyk | Interaction of detonation with inert gas zone[END_REF] and detonation propagation into fuel concentra-tion gradients aligned with the detonation wave propagation direction [START_REF] Strehlow | Transient studies of detonation waves[END_REF][START_REF] Thomas | The behavior of detonation waves at concentration gradients[END_REF][START_REF] Kuznetsov | Detonation propagation, decay, and reinitiation in nonuniform gaseous mixtures[END_REF] have received the most attention. Nonetheless, situations where fuel concentration gradients lie perpendicular, or oblique, to the direction of the propagating detonation wave are not very well understood. In the early works of Oran et al. [START_REF] Oran | Numerical simulations of detonation transmission[END_REF] and Tonello et al., [START_REF] Tonello | Mechanisms of detonation transmission in layered H2-O2 mixtures[END_REF] detonation propagation into a discrete array of reactive/inert gas were considered. Specifically, the simulations of Oran et al. [START_REF] Oran | Numerical simulations of detonation transmission[END_REF] were found useful to observe what happens when an unsteady detonation travels perpendicular to an interface that separates two different reactive mixtures. It was found that upon interaction, a detonation failure followed by a transient series of shock reflections were necessary to re-establish the complex detonation structure coupled to both reactive mediums. While complex large-scale gas-dynamic structures were observed as the detonation traversed from one layer to another, the structure of the thin mixing layer, on orders much smaller than the detonation reaction length, was ultimately believed to be inconsequential to the global behavior of the detonation. [START_REF] Kessler | Gas-phase detonation propagation in mixture composition gradients[END_REF] To explore the influence of much larger mixing layers, on orders larger than the reaction length scales of the detonation, experiments of Ishii and Kojima 21 considered the case of detonation propagation into comparatively large diffuse mixing layers containing reactive and inert gases. In these experiments, it was found that diffusion of the inert gas into the reactive mixture resulted in a curved detonation front as it passed through the mixing layer. Lieberman and Shepherd [START_REF] Lieberman | Detonation interaction with a diffuse interface and subsequent chemical reaction[END_REF][START_REF] Lieberman | Detonation interaction with an interface[END_REF] also considered such arrangements, but isolated diffuse from sharp interfaces, where the mixing layer, δ , is of the order of, or much smaller than λ , respectively. For diffuse interfaces, δ ∼ λ , the authors determined that the observed wave curvature resulted

Accepted to Phys. Fluids 10.1063/5.0113073

Origins of instabilities in turbulent mixing layers behind detonation propagation into reactive-inert gas interfaces 2 from a decrease in available heat release, and therefore a decrease in leading shock velocity as it traveled in the diluted mixture. The subsequent increase in induction delay lead to the creation of a decoupled turbulent mixing zone (TMZ) behind the transmitted shock, separated by a gap of shocked inert gas. Such features were also observed for detonation propagation into sharp interfaces of inert-reactive gas (δ ≪ λ ) but in the absence of strong wave curvature effects. Their analysis suggested that such turbulent mixing was mainly due to the velocity difference across the TMZ, separating shocked and burned gases behind the wave interaction with the inertreactive gas interface. In both cases, diffuse and sharp interfaces, secondary combustion in the TMZ was observed when the reactive mixture was separated by pure oxygen. The role and influence of this secondary combustion in supporting the transmitted shock wave remains unclear.

More recently, detonation propagation in the presence concentration gradients and confinement has also been investigated. [START_REF] Rudy | Critical conditions of hydrogen-air detonation in partially confined geometry[END_REF][START_REF] Boeck | Detonation propagation in hydrogen-air mixtures with transverse concentration gradients[END_REF] While a curved detonation front and decoupling of the reaction zone from the leading shock wave were observed upon interaction with the mixing layer in all the previously cited studies, little focus was placed on experimentally examining the TMZ and the source and influence of secondary combustion as the detonation propagated into the inert gas, except, for the work of Lieberman. [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] Investigation of the source of turbulent mixing behind detonation waves upon interaction with gas interfaces requires further attention.

There have been several attempts to model detonation wave interactions with interfaces that lie perpendicular to the initial wave flow direction. [START_REF] Oran | Numerical simulations of detonation transmission[END_REF][START_REF] Kessler | Gas-phase detonation propagation in mixture composition gradients[END_REF][START_REF] Calhoon | Detonation wave propagation in concentration gradient[END_REF][START_REF] Houim | The influence of acoustic impedance on gaseous layered detonations bounded by an inert gas[END_REF][START_REF] Gaathaug | Detonation failure in stratified layers -the influence of detonation regularity[END_REF][START_REF] Reynaud | A computational study of the interaction of gaseous detonations with a compressible layer[END_REF][START_REF] Han | Role of transversal concentration gradient in detonation propagation[END_REF][START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF][START_REF] Melguizo-Gavilanes | Detonation-diffuse interface interactions: failure, re-initiation and propagation limits[END_REF] Kessler et al., [START_REF] Kessler | Gas-phase detonation propagation in mixture composition gradients[END_REF] who investigated the effects of layers containing gradients in mixture equivalence ratio, found that the structure of detonations propagating through the mixture depended on the width of the stratified reactive layer relative to the detonation cell size. Han et al. [START_REF] Han | Role of transversal concentration gradient in detonation propagation[END_REF] later found that steeper gradients in equivalence ratio had a destabilizing effect on the detonation wave, tending towards marginal behaviour. For sharp interfaces of reactive mixtures with inert gas, Houim and Fievisohn [START_REF] Houim | The influence of acoustic impedance on gaseous layered detonations bounded by an inert gas[END_REF] found that differences in acoustic impedance between the mixtures influenced the detonation wave. For equal acoustic impedance across the interface, unstable detonations were observed, while inert gases with much greater or much lower acoustic impedance permitted stable detonation propagation in the reactive gas layer. The work of Reynaud et al. [START_REF] Reynaud | A computational study of the interaction of gaseous detonations with a compressible layer[END_REF] found that the reactive mixture sensitivity, i.e., its activation energy, was found to influence the minimum reactive layer height under confinement by an inert gas that permitted detonation propagation. The effect of chemistry modelling on the prediction of quenching limits was assessed in Taileb et al. [START_REF] Taileb | Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer[END_REF] It was found that a detailed description of the chemistry was in good agreement with experimental observations and that simplified models fitted to reproduce realistic ignition delay times overpredicted the quenching limits significantly. Note that most of the studies cited used inviscid formulations. While we do not dismiss such an approach, subgrid mixing unavoidably occurs through the introduction of artificial dissipation; these solutions may therefore be sensitive to and influenced by changes in resolution. This shortcoming was re-cently addressed in Melguizo-Gavilanes et al. [START_REF] Melguizo-Gavilanes | Detonation-diffuse interface interactions: failure, re-initiation and propagation limits[END_REF] , where the CLEM-LES approach [START_REF] Maxwell | Influence of turbulent fluctuations on detonation propagation[END_REF] was used to resolve cellular structures to study failure, re-initiation and propagation limits for detonation wave interaction with a diffuse interface created by a gravity current. The aforementioned approach has been used successfully in many configurations including two-and threedimensional planar detonation waves, [START_REF] Maxwell | Influence of turbulent fluctuations on detonation propagation[END_REF][START_REF] Maxwell | The Kolmogorov number in large eddy simulations of calorically perfect gas detonations: Tuning parameter or a constant?[END_REF] detonation quenching and re-initiation through porous media, [START_REF] Maxwell | Turbulent Combustion Modelling of Fast-Flames and Detonations Using Compressible LEM-LES[END_REF][START_REF] Maxwell | Modelling of the transition of a turbulent shock-flame complex to detonation using the linear eddy model[END_REF] and pulse compression detonation engines. [START_REF] Maxwell | Numerical simulations of compression and detonation strokes in a pulse compression detonation system[END_REF] The present work applies the same CLEM-LES methodology to gain insight on the role of turbulent mixing on the observed wave structure for interactions of gaseous detonation waves with sharp reactiveinert interfaces. The main objective of this study is thus to thoroughly examine the source of instabilities in the TMZ that forms upon the interaction of a detonation wave with a sharp interface following past experiments of Lieberman and Shepherd. [START_REF] Lieberman | Detonation interaction with an interface[END_REF]

II. NUMERICAL METHODOLOGY

A. The compressible linear eddy model for large eddy simulation (CLEM-LES)

For the highly compressible and transient flow at hand, the governing Navier-Stokes equations are filtered using the large eddy simulation (LES) methodology. For a calorically perfect gas, the conservation equations for mass, momentum, energy, and subgrid kinetic energy are:

∂ ρ ∂t + ∇ • ( ρ ũ) = 0, (1)
∂ ρ ũ ∂t +∇ • ( ρ ũ ⊗ ũ) + ∇ p -∇ • ρ(ν + ν t ) ∇ ũ + (∇ ũ) T - 2 3 (∇ • ũ) Î = 0, (2)
∂ ρ ẽ ∂t +∇ • ( ρ ẽ + p) ũ -ũ • τ - γ γ -1 ∇ • ρ( ν Pr + ν t Pr t )∇ T = -Q ω, (3)
∂ ρk sgs ∂t +∇ • ( ρ ũk sgs ) -∇ • ρν t Pr t ∇k sgs = ρν t ∇ ũ + (∇ ũ) T - 2 3 (∇ • ũ) Î • (∇ ũ) -ρε,( 4
)
where ρ, p, e, T , u, and k sgs refer to density, pressure, specific sensible + kinetic energy, temperature, velocity vector, and subgrid kinetic energy, all of which are normalized by the quiescent reactive mixture properties (see Ref. 35). The total energy (internal + kinetic) and the equations of state are given by ẽ = p ρ(γ -1)

+ 1 2 ũ • ũ + 1 2 k sgs , (5)
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073

p = ρ T . (6)
Favre-average (LES) filtering is achieved through f = ρ f / ρ, where f represents one of the many state variables (ρ, p, e, T , u, and k sgs ). Other usual properties to note are the heat release, Q, the ratio of specific heats, γ, the kinematic viscosity, ν, the filtered viscous tensor, τ , and the identity matrix, Î.

The turbulent viscosity and dissipation are modelled according to

ν t = 1 π 2 3C κ 3/2 √ k sgs ∆, (7)
ε = π 2k sgs 3C κ 3/2 / ∆, (8)
whose derivations in these forms are detailed in Ref. 35.

Here, ∆ is the minimum grid spacing (which corresponds to the LES filter size), and C κ is the Kolmogorov number. Finally, like in any other reactive LES approach, the chemical reaction source term, ω, requires closure. This is achieved using the CLEM sub-grid modeling strategy. [START_REF] Maxwell | Turbulent Combustion Modelling of Fast-Flames and Detonations Using Compressible LEM-LES[END_REF] The micro-scale mixing and chemical reaction are handled entirely on the subgrid, through a supplementary simulation of a 1D sample of the flow field within each fully refined LES cell. The system of unfiltered equations that was solved on the sub-grid included the conservation of energy and reactant mass fraction (Y ):

ρ DT Dt - γ -1 γ ṗ -ρ ∂ ∂ m ρ 2 ν Pr ∂ T ∂ m = - γ -1 γ Q ω + ḞT , (9)
ρ DY Dt -ρ ∂ ∂ m ρ 2 ν LePr ∂Y ∂ m = ω + ḞY . ( 10
)
The source terms, ḞT and ḞY , account for the effect of turbulence on the sub-grid in the form of random "stirring" events; [START_REF] Menon | The linear-eddy model[END_REF] ṗ accounts for temporal rates of change in pressure, which are obtained entirely from the large-scale simulation, Eqs. (1) to (4). Note that the partial derivatives are taken along m, which is a one-dimensional mass weighted coordinate whose mapping to Cartesian spatial coordinates is given by

m(x,t) = x x 0 ρ(x,t)dx. (11)
A first order one-step irreversible Arrhenius chemistry, Y → P, closes the governing equations, whose rate is given by

ω = -ρAY exp(-E a /T ). ( 12
)
where E a and A are the activation energy and pre-exponential factor, respectively. For simplicity, differences in molecular weight between the inert gas, reactive mixture and detonation products were neglected. As a result, we consider only the transport of a single reactant species (Y ), whose value of Y = 1 (or Ỹ = 1 on the filtered LES-scale) represents the unburned reactive mixture while Y = 0 (or Ỹ = 0) is used to denote the inert gas or detonation products.

B. Numerical implementation

To solve the governing equations, Eqs. ( 1)-( 4), a second order exact Godunov method [START_REF] Gottlieb | Assessment of Riemann solvers for unsteady one-dimensional inviscid flows of perfect gases[END_REF] is used. To keep the computation tractable at the scales considered, Adaptive Mesh Refinement (AMR) [START_REF] Falle | Body capturing using adaptive cartesian grids[END_REF] is also used to enable sufficient resolution in regions where strong gradients exist, such as shocks, reaction zones, and contact surfaces. For this study, the refinement criteria was defined as follows, a computational cell is refined if Ỹ > 1 × 10 -4 and ρ > 1.1ρ 0 , where ρ 0 is the density of the fresh mixture ahead of the detonation wave. Cells were also flagged as needing refinement when density changes of more than 1% occurred between grid levels; special care is taken to ensure smooth transitions across fine-coarse cell boundaries. For the LEM-subgrid Eqs. ( 9) and (10) were solved using implicit integration. For additional details on the implementation of the LEM and its coupling to the LES, the reader is referred to Refs. 10, 35, and 36.

C. Numerical domains, initial and boundary conditions

We studied detonation propagation into sharp interfaces separating reactive mixture (2.5C 2 H 4 +3O 2 ) from inert gas (N 2 ) at interface angles ranging from α = 0 • to 60 • , at 15 • intervals. The initial distribution of inert and reactive gas for a general case, and the special case of α = 0 • , where the interface horizon has been located midway through the height of the domain, are shown in Figs. 1a and1b, respectively. The total domain size of each simulation was 8100 ∆ 1/2 long (x = 0.81 m) by 1000 ∆ 1/2 high (y = 100 mm); nearly to scale with the experiments of Lieberman and Shepherd, [START_REF] Lieberman | Detonation interaction with an interface[END_REF] which had a test section height of y = 150 mm. In each case, a corresponding ZND profile was initialized in the first 100 ∆ 1/2 with D = 1.1D CJ in order to avoid startup errors. The subsequent detonation front was permitted to propagate up to x = 2000∆ 1/2 before encountering the interface, such that its velocity and cellular structure would be established. Initially, k sgs is null everywhere. Background turbulence in the quiescent fluid was thus neglected, whereas the subgrid turbulent kinetic energy was entirely produced by the detonation-sharp interface interaction. Finally, symmetric boundary conditions were specified on the top and bottom boundaries, while a physical wall boundary was specified at x = -100∆ 1/2 , that is, no-slip for velocity (Dirichlet), and zero gradient (Neumann) for the remaining variables.

D. Model calibration and grid convergence study

For the one-step kinetic scheme adopted, the model parameters, Q = 116 and E a = 27.8, were tuned to reproduce the correct post-shock ignition delay times, τ ind , for a fuel rich ethylene-oxygen mixture (2.5C 2 H 4 +3O 2 ) at an initial pressure of p 0 = 11 kPa. The ratio of specific heats (γ = 1.126) was determined at the detonation post-shock Von-Neumann state, for a shock travelling at the CJ-condition. Also, the pre-exponential factor, A = 40.7, and diffusion coefficients (ν = 0.00786, Pr = 0.77, Pr t = 1, and Le = 1) were chosen such that the one-step model reproduced the correct half reaction length (∆ 1/2 = 0.1 mm), and the laminar flame speed at post-shock conditions, for a shock traveling at 50% of the theoretical CJ speed for the given quiescent mixture (S L = 6.64 m/s at M D = 0.5M CJ ). This post-shock state was chosen in order to have a measurable flame speed prior to auto-ignition of the gas. The reference values for τ ind , ∆1/2 , and S L were determined with Cantera 41 using the chemical kinetic mechanism of Wang et al., [START_REF] Wang | A comprehensive mechanism of C2Hx and C3Hx fuel combustion[END_REF] which was specifically designed for C2 and C3 chemistry. The values of S L shown in Fig. 2 were determined from separate unsteady 1D simulations. As can be seen, grid convergence is demonstrated, and a minimum resolution of ∆ = ∆ 1/2 /128 was required to capture the correct S L using the one-step model, to within 1% of the steady solution obtained using Cantera. The performance and limitations of this model in properly capturing S L and flame structure as a function of resolution for the LEM-subgrid formulation is well documented in Ref. 36. Finally, ∆ (the LES filter size) was specified so that the average experimental cell size corresponding to the fuel mixture considered was reproduced. Figure 3 shows the influence of changing the LES filter size ∆, for a fixed C κ = 1.5 (the standard Kolmogorov constant value) and minimum number of subgrid elements per LES cell (N), on the cellular structure for a detonation wave propagating into uniform reactive mixture (no inert-gas interface present). The average cell structure obtained did not change beyond ∆ = ∆ 1/2 /8 with N = 16. Note that the same effective resolution ( ∆eff = ∆ 1/2 /128) was required to resolve S L on the subgrid, see Fig. 2. Increasing N, while holding ∆ constant, did not significantly influence the overall cell sizes; compare frames b) and c) in Fig. 3. Figure 4 shows a numerical soot foil for a detonation propagating into uniform mixture compared, to scale, with an experimental soot foil. [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] A mean cell size of λ sim = 2.8 mm was found by postprocessing the numerical data using an autocorrelation procedure, [START_REF] Sharpe | Statistical analysis of cellular detonation dynamics from numerical simulations: one-step chemistry[END_REF] which agrees well with the experimentally reported values, [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] where 2 mm < λ exp < 4 mm. A resolution of ∆ = ∆ 1/2 /8 with N = 16 subgrid elements per LES cell is thus used for the simulation results shown in the following sections as it was shown to resolve features of interest. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073

III. RESULTS

Before attempting a detailed comparison with the work of Lieberman and Shepherd, [START_REF] Lieberman | Detonation interaction with an interface[END_REF] and carrying out the analysis for detonation-sharp interface interactions, it is important to verify that the CLEM-LES approach is capable of reproducing the growth rates reported in classical experiments of compressible shear layers. [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] In Appendix A we include our verification and validation efforts, and we show that our methodology is indeed capable of capturing the correct shear layers growth rates for a wide range of velocity ratios. 5a) for α = 45 • , a transmitted shock-TMZ complex formed as the wave travelled through the interface. Since there is no reactive gas present, the TMZ served only to mix detonation products with the inert gas. As a chemical activity indicator, the reaction rate was superimposed onto the density fields in Fig. 5. In this particular case, all of the reactive mixture converts to products in the vicinity of the detonation wave, and therefore the evolution of the TMZ, for the most part, occurs between two inert gases. The transmitted shock wave angles from the experiment and simulation were found to be in agreement, with values of β exp = 70 • and β = 68.8 • , respectively. Note that, the angles α and β were measured relative to the horizontal (x-axis) and taken as positive in the clockwise direction. Additional details on how shock angles were measured numerically are given in the next section. The numerical simulation also qualitatively captured the formation of a Mach stem upon the transmitted shock reflection with the upper wall, as indicated in Fig. 5, although the details of the Mach stem size were not analysed here. In Ref. 23, replacing the inert gas with oxygen resulted in a much more pronounced Mach stem. The authors attributed this observation to a postulated secondary combustion, where shock-heated oxygen was able to mix with the unburned reactants present in the TMZ, hence changes in γ. [START_REF] Ben-Dor | Shock Wave Reflection Phenomena[END_REF] The flow structure for α = 0 • is shown in Fig. 5b. As mentioned above, the lack of reactivity in the inert gas also lead to the decoupling of the transmitted shock wave and TMZ in this case. The transmitted shock angle was found to be β = 26.33 • . Notably, for both cases, the incident shock angle relative to the initial interface angle (i.e. βα) were in good agreement with values of 23.8 • and 26.33 • , respectively. Although much less pronounced, the presence of a Mach stem was also observed for α = 0 • .

The CLEM-LES methodology captures quite well the experimentally observed qualitative features, which include the formation of a TMZ separated by the incident shock by a gap of dense gas. Additionally, the incident shock angle for α = 45 • was also quantitatively captured, with only a 1.2 • difference between experiment and simulation. This overall agreement was considered satisfactory for investigating the influence of interface orientation on the flow field dynamics, and also to gain further insight onto the main mechanisms contributing to the TMZ development. As the interface angle, α, was varied, the same features were always observed, with only minor differences in relative quantitative measurement of the shock angle (β -α), and qualitative size of the TMZ (always within a few degrees).

B. Self-Similarity of the wave evolution and growth of the turbulent mixing layer

From Fig 5, the interaction of the detonation wave with the sharp interface appears to yield a self-similar structure in its gas-dynamic evolution (i.e. a structure similar to itself at different times), with the exception of turbulence in the TMZ. To demonstrate the self-similar nature of the gas-dynamic evolution, and to determine the statistical self-similarity of the TMZ evolution, flow fields from several instances in time were ensemble-averaged, separately for each simulation, using the locus of the leading shock wave and reactive-inert gas interface as a spatial reference point. The averaged density gradient fields for the five different α cases considered, are shown in Fig. 6. For α = 0, 15, 30 and 45 • , 540 to 750 time instances were ensemble-averaged to produce smooth profiles; for α = 60 • , only 250 instances were available (limited by the domain height and frequency at which the data was collected). The averaging procedure included only times sufficiently before the detonation front encountered the bottom wall, such that the detonation cell structure was not influenced by restrictions in the number of cells within the local thickness of the re-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. action zone. For all cases, a clear incident shock location and collapse of the shock onto the reference node, where the leading shock wave meets the reactive-inert gas interface, were observed. Also, for all cases except α = 60 • , the TMZ thickness and gap size between the transmitted shock and TMZ were well defined. Note that a TMZ and gap were also present for α = 60 • through averaging of the temperature and density fields, although with an increased degree of measurement uncertainty/deviation due to less data available for the averaging process. Next, we focus our discussion on the α = 0 • case to demonstrate the self-similar nature of the wave dynamics and TMZ evolution.

PLEASE CITE THIS ARTICLE AS

As was done for the density gradient field, shown in Fig. 6, the ensemble-averaged density and temperature fields were also obtained for each simulation, using as a reference the location where the leading shock wave meets the reactiveinert gas interface. Figure 7 shows the vertical profiles of density and temperature extracted at various locations along the x-coordinate for α = 0 • (the remaining values of α studied were found to exhibit similar behavior). In this figure, the y-coordinate, was re-scaled by the size of the gap and TMZ measured along y for a given value of x, ∆y (TMZ+gap) .

The features of interest are the leading shock, and the start and end of the TMZ, which are used to determine the size of the gap and TMZ, and also the incident shock angle. Here, the average position of the transmitted shock wave is determined at the (x, y) coordinate where ρ > 2. The start and end of the TMZ were determined from the temperature field where the post-shock temperature increased by 4%, and where the temperature reached within 4% of the fully burned temperature, respectively. We chose this particular threshold, as it was found to best collapse our data in Fig. 7, where the measured quantities in the burned products were found to contain some residual fluctuations, even after ensemble-averaging up to 750 instances. The burned reference temperatures, for each case, were determined using the gas-dynamic analysis summarized later on in Section IV A. Upon normalizing by ∆y (TMZ+gap) in Fig. 7, for each density and temperature profile extracted along a given x value, they collapsed to a single profile providing evidence of the self-similar nature of the resulting structure. Note, however, that these self-similar curves were constructed sufficiently far from the node, at least 100∆ 1/2 (10 mm) such that initial curvature effects of the transmitted shock were inconsequential. Also noted are zoomed in portions of Fig. 7 that reveal the degree of the low frequency residual fluctuations not removed through the averaging process.

Owing to the self-similarity of the wave dynamics, it is possible to extract various angles of interest from each simulation; these are given in table I. β is the incident shock angle as indicated in Fig. 5. θ TMZ and θ gap are the angles associated with the TMZ size and gap size, as indicated in Fig. 6. Also shown, are the corresponding angles experimentally measured by Lieberman [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] for α = 45 • . The angles were determined from the ensemble-averaged temperature and den-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073 sity profiles, similar to the ensemble-averaged density gradient profiles shown in Fig. 6, using a linear regression procedure with the locus of the leading shock wave and reactiveinert gas interface as the reference point. In all cases, (β -α) dependent on the interface angle α. More specifically, (β -α) has a maximum around α = 15 • , and decreases rapidly to small relative angles beyond α > 45 • . Evidently, this behavior was expected, as (β -α) → 0 when α → 90 • . This must be satisfied, since this extreme case would represent a contact surface parallel to the detonation front. Any resulting transmitted shock wave, contact surface, or mixing zone would also lie parallel to the original detonation front. The θ TMZ was found to have a minor dependence on α, and to be θ TMZ > 5 • over the range of α considered (0

• ≤ α ≤ 45 • ).
It is important to note, however, that while the measured β from simulation and experiment match within some margin of error, a clear discrepancy arises between the measured θ gap and θ TMZ , for α = 45 • . These angles were measured to be on the order of θ gap ∼ 2 • and θ TMZ ∼ 5 • from simulation; Lieberman reported angles of θ TMZ = θ gap = 7 • . A perfect match is not expected due to typical challenges associated with measurement of both simulation and experiment. First of all, the experimental measurements of θ TMZ and θ gap were taken directly from a single instantaneous schlieren image. Locating the start and end points of the shock and TMZ must be done subjectively to a flow field that is unsteady and turbulent. As a result, difficulties arise in choosing satisfactory limits for θ TMZ and θ gap . Lieberman reported an experimental measurement error of ∼ ±1 • , thus we can expect θ TMZ to be as small as ∼ 6 • . Also, additional challenges ensuring the correct orientation of the detonation wave relative to the contact surface contributed to the largest source of error in the experiment. It is also likely that the presence of the diaphragm in the experiment may have influenced the growth rate of the TMZ. As a result, each experiment is likely to exhibit some variability in θ TMZ . For example, at a similar condition, where α = 44 • , θ TMZ = 4.5 ± 1 • , [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] well in line with the numerically predicted value. Further discrepancies for θ gap were found from the Lieberman experiments, [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] where 3.5 ± 1 • ≤ θ gap ≤ 7 ± 1 • for the range 41 • ≤ α ≤ 45 • . In our measurements for θ TMZ and θ gap , we report the error of the slope coefficient (measured angles) from the linear regression procedure to be on the order of ∼ ±0.5 • or less, as reported in table I, with the exception of α = 60 • where significantly less data were available for the averaging procedure.

To address this significant discrepancy in θ gap , we repeated the simulation for α = 45 • , but with γ = 1.3. In this case, a slightly larger transmitted shock angle was observed, with β = 71.6 ± 0.1 • . This also lead to a smaller θ TMZ , which was measured to be 2.7 ± 0.4 • . θ gap , on the other hand, was found to compare better to the experiments, with θ gap = 4.1 ± 0.4 • . Although γ = 1.3 better represents the inert and combusted product gases, this choice also lead to a slightly smaller mean cell size of λ sim = 2.5 mm instead of 2.8 mm.

Despite the difference in θ TMZ measured between simulation and experiment, what we find of particular interest is how the size of the TMZ at lower angles, for α ≤ 45 • , appears to have very little dependence on α itself, and more so on the

• 5.2 ± 0.3 • 2.0 ± 0.2 • 45 • 1.126 68.8 ± 0.1 • 23.8 ± 0.1 • 5.6 ± 0.4 • 1.1 ± 0.3 • 45 • 1.3 71.6 ± 0.1 • 25.6 ± 0.1 • 2.7 ± 0.4 • 4.1 ± 0.4 • 60 • 1.126 77.5 ± 0.5 • 17.5 ± 0.5 • 3.6 ± 1.2 • 0.7 ± 1.2 • 41 • (Ref. 26 ) 67 ± 1 • 26 ± 1 • 6.5 ± 1 • 3.5 ± 1 • 44 • (Ref. 26 ) 69 ± 1 • 25 ± 1 • 4.5 ± 1 • 6 ± 1 • 45 • (Ref. 26 ) 70 ± 1 • 25 ± 1 • 7 ± 1 • 7 ± 1 •
choice of γ (for which increasing values lead to a larger mean cell size and a smaller TMZ). In fact, analysis performed by Lieberman, [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] taking into account realistic thermo-chemistry, indicated that as α → 0 • , the convective Mach numbers of the shear layer, and shear layer growth in the TMZ resulting from Kelvin-Helmholtz instability, should both tend to zero. Experiments were not conducted to verify this, but nevertheless, these contradictory findings from our simulations warrant further analysis to determine the source of turbulent mixing in the observed TMZ, and to what extent velocity shear plays a role.

IV. DISCUSSION

A. Gas-dynamic analysis

To look more closely at the source of turbulent mixing in the TMZ, it is first necessary to reproduce the gas-dynamic analysis of Lieberman. [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] In this analysis, the steady solution of possible wave geometries and thermodynamic states can be determined in the frame of reference of the 'node', where the detonation front (D) meets the reactive-inert gas interface (CS), as shown in Fig. 8. Owing to the small difference in refractive index between the two quiescent gases, the solution system, for the range of α under investigation, contains a transmitted shock (T), a reflected expansion wave (RE), and the deflected contact surface (CS'). The solution is composed of 5 possible thermodynamic states. These states are the quiescent reactive gas (1), the CJ-combustion products behind the detonation wave (2), the expanded combustion products (3), the shocked inert gas (4), and the quiescent inert gas (5). The solution procedure is the following: the CJ-state (2) is first obtained across the detonation wave for a given D CJ . Then, the Prandtl-Meyer solution [START_REF] Thompson | Compressible-Fluid Dynamics[END_REF] is applied across the expansion fan to determine state (3). Finally, oblique shock jump conditions are applied to obtain state (4). In this solution methodology, the transmitted shock strength and expansion fan size must be iterated until the pressure, p, and particle path trajectories, u, are matched across states (3) and (4) (u 3 = u 4 ; p 3 = p 4 ). While Lieberman considered only the solution with temperature dependent thermodynamic properties, using Can-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. tera, we also show solutions for a calorically perfect gas using a constant γ = 1.126 to be consistent with our numerical simulations (single-γ model). For comparison, we also consider a solution where three independent, constant values for γ are known in states ( 1), (2), and (5). For this 3-γ model, we take γ 1 = 1.30, γ 4 = γ 5 = 1.40, and γ 2 = γ 3 = 1.28 (determined for each state using Cantera). The jumps across the detonation wave, shock wave, and expansion fan are described in Ref. 46 for the single-γ and 3-γ models. Additional details of the solution procedure can be found in the work of Lieberman. [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] As part of the gas-dynamic solution obtained, the transmitted shock angle (β ) for all three-methods (temperature dependent thermodynamics, and calorically perfect gases with single-and 3-γ) as a function of α are shown in Fig. 9. The transmitted shock angles measured experimentally, [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] and those predicted numerically are also included in the figure. For all values of α, both calorically perfect gas models yielded a slightly lower transmitted shock angle, β , compared to the temperature dependent thermodynamics solution. This result was consistent with the fact that the experimental transmitted shock angle, β , at α = 45 • , was also slightly larger than that obtained numerically (within 2.2 • ). It is thus sensible to attribute the differences between β exp and β sim to the calorically perfect gas assumption adopted. Notably, both calorically perfect gas models resulted in the same β values, despite the differences observed in densities and velocity magnitudes in states (3) and (4). Additional sources of discrepancies between experiments and models are later discussed in Section IV C, in the context of the TMZ growth rates.

PLEASE CITE THIS ARTICLE AS

The transmitted shock angles, β , measured from the CLEM-LES followed the same trend as the solutions obtained from the gas-dynamic analysis. However, in all cases, the β values were systematically under-predicted by a few degrees (1 to 3 • ), see Fig. transmitted shock with the top boundary making it relax at a lower angle, whereas the latter results from the wall boundary that was specified behind the detonation, which is known to develop a self-similar Taylor expansion-wave solution behind the reactive wave front. [START_REF] Fickett | Detonation Theory and Experiment[END_REF] While such expansion waves do not influence the detonation velocity, they are able to interfere with the gas-dynamic solution, thus leading to a reduction in velocities and pressure in states ( 3) and ( 4), and thus a weaker transmitted shock is required to satisfy the matching condition. It should also be pointed out that β exp at α = 45 • was also less than that obtained using temperature dependent thermodynamics. We note that the experiment likely had many other factors which may have contributed to this discrepancy. Such as the presence of the film used to separate the gases, as well as the usual velocity deficit observed in real detonations propagating in narrow channels. In the experiments of Lieberman, 26 a velocity deficit up to 5% below the CJ condition was reported.

B. Inert simulations in frame of reference of node

In addition to the full lab-frame CLEM-LES simulations conducted in this work, to isolate the contributions of Kelvin-Helmholtz instability to the TMZ growth, inert simulations in the frame of reference of the node were also carried out. By transforming the frame of reference, we eliminate any influence coming from the left boundary (wall) thereby expecting a better prediction to the gas-dynamic solution for the dynamically evolved wave angles obtained, and the CJ solution may be prescribed as an initial/boundary condition, avoiding the requirement for integrating reaction source terms. This has the added benefit of removing the influence of cellular instability on the solution.

Here, domain sizes of 1000∆ 1/2 wide by 1000∆ 1/2 high were considered, with the node centered at (x, y) = (0, 0), and the same resolution as in the lab-frame simulations was used. The domains for each simulation (α = 0, 30, 45, 60 • ) were initialized, in the frame of reference of the node, with the CJ-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073 solution for a calorically perfect single component gas 47 everywhere, i.e.,

ρ(x, y) = ρ CJ = ρ 0 D CJ D CJ -u CJ , ( 13
) p(x, y) = p CJ = p 0 + ρ 0 D 2 CJ γ + 1 , ( 14
) u(x, y) = u CJ -D CJ = p CJ -p 0 ρ 0 D CJ -D CJ , ( 15
) v(x, y) = D CJ tan α, (16)
where ρ 0 = 1, p 0 = 1/γ, and D CJ = 8. The right boundary, at x = 0, was prescribed to have the following conditions:

ρ(0, y) = ρ 0 if y > 0 ρ CJ otherwise, ( 17
)
p(0, y) = p 0 if y > 0 p CJ otherwise, (18)
u(0, y) = -D if y > 0 u CJ -D CJ otherwise, (19)
v(0, y) = D CJ tan α for all y. ( 20
)
The remaining boundaries for all variables were set to Neumann type (zero gradient). The resulting unsteady flow was allowed to evolve naturally up to t = 500 (150 μs), sufficiently long to reach a steady structure, before averaging statistics were obtained. Time-averaged density gradient fields for three different angles of α are presented in Fig. 10. Upon measuring various state properties and wave angles, we found the inert simulation solutions matched closely those from the single-γ gasdynamic solution (see Fig. 9). In addition to this agreement with the gas-dynamic model, very little shearing was observed for all values of α. Moreover, this result was insensitive to changes in numerical resolution. As shown in Section A, our LES methodology is able to reproduce the Kelvin-Helmholtz instability observed in supersonic shear layer experiments of Brown and Roshko. [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] Remarkably, by removing cellular instabilities at the detonation front, the resulting TMZ growth rates are much smaller, nearly negligible.

C. Growth rate of Turbulent Mixing Zone

To investigate further the growth rate of the TMZ, we carried out the same analysis done by Lieberman and Shepherd, [START_REF] Lieberman | Detonation interaction with an interface[END_REF] but we carefully choose the frame of reference in which shear growth rates occur. In Lieberman's analysis, the following transformation to laboratory coordinates for the two stream velocities on each side of the TMZ was applied, from which the shear layer growth rates were calculated:

U 1,x ′ = | w 3 | -| w 1 | cos (δ -π -α) (21)
and

U 2,x ′ = | w 4 | -| w 1 | cos (δ -π -α), (22)
where x ′ direction is positive away from the node, as indicated earlier in Fig. 8b. Here, w 1 , w 3 , and w 4 refer to the flow velocities in states ( 1), (3), and (4), in the frame of reference attached to the node. The angles α and δ are the angles of the interface and the deflected contact surface, respectively, also indicated in Fig. 8. Upon transformation to laboratory coordinates, the flow velocities U 1 and U 2 have equal components along the y ′ coordinate:

U 1,y ′ = U 2,y ′ = | w 1 | sin (δ -π -α). (23)
Upon solving Eqs. ( 21) and ( 22) for the entire range of α, we find that U 1,x ′ and U 2,x ′ are always negative. Although the relative velocities w 3 , and w 4 were aligned with CS', but in a direction travelling away from the node, the absolute lab-frame velocities have trajectories in the +x and +y directions (also in the -x ′ and +y ′ direction). This results in the TMZ not only shearing, but also advected forward and upwards when viewed from a fixed position. Also, while shear growth occurs in the -x ′ direction, the node travels at supersonic velocities in the +x direction. As a result, the laboratory-frame flow field is always locally unsteady, making the choice of U 1,x ′ and U 2,x ′ to measure shear layer growth inconsistent with past work involving splitter plates. [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF][START_REF] Konrad | An Experimental Investigation of Mixing in Two-Dimensional Turbulent Shear Flows with Applications to Diffusion-Limited Chemical Reactions[END_REF][START_REF] Dimotakis | The mixing transition in turbulent flows[END_REF] Instead, since conventional shear layers have always been measured and analysed as a spatially growing instability in a quasi-steady state, w 3 and w 4 should be used instead to determine the shear layer growth contribution to θ T MZ , in the node-frame of reference.

In this frame of reference, the flow becomes quasi-steady, as if a splitter plate was attached to the node itself, and more consistent to analyze with previously developed shear layer growth correlations. Therefore, to estimate the K-H contribution in a quasi-steady frame of reference, we select our shear flow velocity streams to have the following velocities:

U 1 = w 4 ( 24
)
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073

and

U 2 = w 3 . ( 25
)
By convention, we choose these such that U 1 > U 2 . [START_REF] Smits | Turbulent Shear Layers in Supersonic Flow[END_REF] Then, the convective velocity, U c , is defined as the velocity of the centre of shearing between the two streams. For low convective Mach number shear layer growth evolution, U c is approximated as

U c ≈ 1 + r √ s 1 + √ s ( 26
)
where r = U 2 /U 1 and s = ρ 2 /ρ 1 such that states 1 and 2 for the shear layer analysis correspond to the expanded burned products and shocked inert gas streams, respectively. The respective stream Mach number are then defined as

M c1 = U 1 -U c c 1 and M c2 = U c -U 2 c 2 ( 27
)
where c 1 and c 2 are the speeds of sound in the expanded burned products and shocked inert gas streams. For the special case where γ is equal and constant in both streams,

M c = M c1 = M c2 = U 1 -U 2 c 1 + c 2 . ( 28
)
Then, following Ref. 26, profiles of s, r, and M c as functions of α, obtained using the three different models described earlier, are shown in Fig. 11. We note some discrepancies between the calorically perfect gas models and the temperature dependent thermodynamics model: (i) The single-γ model always produces some deviation from the real gas solution for s, r, and M c . We thus expect the largest error in estimating the shear layer growth with this model. Note, however, that the inert simulations, in the frame of reference of the node, capture well the values obtained from the gas-dynamic analysis (also shown in the figure). (ii) The 3-γ model gave better predictions of s, r, and M c , especially for the density ratio, s. Deviations in the solution for the velocity ratio r were evident for all models, and also lead to deviations in the solution for M c . The significant difference in values obtained for r and s compared to Lieberman [START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] are due to our selected frame of reference, however the variation of M c compares well to Lieberman's model results.

Then, using the solutions for s, r, and M c , the shear layer growth rate, dδ /dx ′ , is readily obtained using the correlations of Dimotakis, 51

dδ dx ′ Mc→0 ≈C δ (1 -r)(1 + √ s) 2(1 + r √ s) × 1 - (1 - √ s)/(1 + √ s) 1 + 2.9(1 + r)/(1 -r) ( 29
)
where C δ lies in the range 0.25 < C δ < 0.45. [START_REF] Slessor | Turbulent shear-layer mixing: growth-rate compressibility scaling[END_REF] To account for compressibility effects, the correction of Papamoschou and Roshko [START_REF] Papamoschou | The compressible turbulent shear layer: an experimental study[END_REF] is applied such that with f ∞ = 0.2. Finally, the growth angle of the TMZ, due to flow shearing, may be determined from

dδ dx ′ = (1 -f ∞ )e -3M 2 c1 + f ∞ dδ dx ′ Mc→0 ( 30
θ = tan -1 dδ dx ′ (31)
which represents the angle formed as the shear layer size, δ , grows as it is advected from a reference point (i.e., a splitter plate, or the node in our case) with velocity U c . θ was determined for all three gas-dynamic models and presented in . This is much less than the Liberman model predictions, experiments, and the shear growth rates measured from the CLEM-LES simulations, all of which are compared in Fig. 12. Again, despite Lieberman's agreement, there is some inconsistency in using the laboratory-frame of reference to measure the shear-growth rates of the TMZ due to the fact that the flow field is not quasi-steady. For the Lieberman model, and our temperature dependent model both exhibit vanishing θ as α → 0. This is contrary to the CLEM-LES observations, which suggest there may be other significant contributing factors to the TMZ growth. As will be discussed in the next section, we suggest that pressure pulses associated with the cellular nature of the detonation front significantly influence the TMZ growth. Upon measuring θ in the inert node-frame simulations, using the same procedure of the lab-frame CLEM-LES simulations, good agreement with our predicted shear growth angle of θ was found. In fact, θ obtained from the inert node-frame simulations was also always less than 1 • . Despite the small difference in θ predicted from the gas-dynamic models and the inert node-frame simulation, it is reasonable to speculate that shearing due to K-H is minimal, and only contributes to a fraction of the actual TMZ growth as observed in experiments and lab-frame CLEM-LES.

D. Contribution of pressure pulses to turbulent mixing

Due to the rather small contribution of K-H instability on the observed TMZ growth, further investigation of the source of turbulent mixing in the TMZ is necessary. We investigate the influence of instabilities originating from the detonation front itself on the TMZ evolution. For this, we once again use inert node-frame simulations, as described previously in Section IV B, where irregular cellular instabilities have been removed. In this case, however, we add controlled sinusoidal perturbations to the inflow, such that boundary condition for pressure in Eq. ( 14) is replaced with

p(0, y) = p 0 if y > 0 p CJ + 15 sin (2πt/2.5) otherwise. ( 32
)
The amplitude of the perturbation for p (value of 15) was chosen to be representative of pressure pulses associated with changes in the detonation velocity up to ±25% the CJ value, which may be expected owing to the true unsteady nature of the detonation wave front. The period of the perturbation was chosen, roughly, to reflect the time it takes for a detonation wave to pass over the characteristic cell size. In Fig. 13 an instantaneous density field obtained for α = 45 • resulting from the sinusoidal forcing is shown at an instant in time once a quasi-steady state was achieved in the evolution of the TMZ. For comparison, time-averaged density fields from the inert simulation with and without pressure forcing at the boundary are also shown. Pressure fluctuations originating from the inflow have a much larger contribution to the TMZ evolution and growth than velocity shearing alone. Upon time-averaging the perturbed flow field through 300 instances in time, and measuring the growth rate using the procedure described in Section III, θ TMZ = 1.2±0.3 • for α = 45 • . Although this TMZ growth rate is still less than the lab-frame simulation at the same value of α, θ TMZ = 5.6 ± 0.4 • , both angles are much larger than the value obtained for the inert simulation without forcing (θ TMZ = 0.28 ± 0.05 • ). To test the effects of the pressure pulses, we repeated the simulation by independently doubling the pressure pulse amplitude and period. We found that by doubling the amplitude of pressure pulses, θ TMZ = 1.6 ± 0.3 • . Then, in another simulation where the period was halved, and thus more frequent, we found that θ TMZ = 2.0 ± 0.1 • . For both a doubled amplitude and frequency, θ TMZ = 2.2 ± 0.3 • . Based on this result, it appears that the observed TMZ growth is very sensitive to pressure pulses originating from the unsteady cellular detonation front, and that such pulses have a significant contribution to the observed TMZ growth rate. More specifically, the latter pulses, originating from a normal direction to the detonation wave and likely enhanced by baroclinic instabilities, account for 20 to 39% of the TMZ growth, while K-H instability on its own accounts for only 5%. The remaining contribution may originate from multi-dimensional transverse waves and irregularity associated with the cellular structure of a real detonation front, plausibly enhanced by Richtmyer-Meshkov instabilities between the transmitted oblique shock and the TMZ. From these figures, we find that these normal pressure pulses and transverse waves arguably have a much higher impact on the TMZ growth than K-H instability alone.

E. Three-dimensional effects and the critical Reynolds number

While our two-dimensional computations attempted to clarify the role of Kelvin-Helmholtz (K-H) shearing, such phenomena is actually three-dimensional in nature. To fully capture vortex-stretching in all directions, and its influence on This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. the actual TMZ growth, three-dimensional simulations are required. Furthermore, the cellular structure of the detonation wave, and the resulting transverse pressure perturbations are also three-dimensional and may influence the TMZ growth.

PLEASE CITE THIS ARTICLE AS

In Section A we showed that three-dimensional effects could increase the K-H contribution to the shear layer growth by ∼10% compared to two-dimensional simulations. However, at the scales considered here, the K-H contribution was found to be small, i.e. θ TMZ < 1 • when no external forcing of pressure perturbations is present. To investigate this further, the critical Reynolds number, Re crit as defined in Refs. 48 and 49 is an important parameter to quantify, as Konrad previously found that planar gas-phase shear layers are only able to sustain three-dimensional fluctuations when Re crit 10, 000 -20, 000 thus marking a transition to a threedimensional regime as a result of increased molecular mixing.

Re crit = ∆Uδ crit ν , (33)
where ∆U = U 1 -U 2 , the difference in velocity between the fast and slow streams, respectively, in the steady-state frame of reference. From this, the physical distance before which three-dimensional shearing is not expected to play a role can

1×10 1 1×10 2 1×10 3 1×10 4 1×10 5 1×10 6 1×10 7 0.01 0.1 1 1×10 0 1×10 1 1×10 2 1×10 3 1×10 4 1×10 5 1×10 6 x' crit (∆ 1/2 ) [mm] ∆U/(U 1 +U 2 )
Simulation, s=7 Simulation, s=0.130

Curve fit, s=7 Curve fit, s=1/7 FIG. 14. Distance along x ′ to Re crit , before three-dimensional effects from K-H shearing influence the flow. The density ratio between the slow and fast streams is given by s = ρ 2 /ρ 1 .

be found:

x ′ crit ≈ dδ dx ′ -1 δ crit . (34)
From Brown and Roshko, [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] correlations may be used to estimate the shear layer growth depending on the density ratio:

dδ dx ′ = C U 1 -U 2 U 1 +U 2 , (35)
with C = 0.51 for ρ 2 /ρ 1 = 7, or C = 0.28 for ρ 2 /ρ 1 = 1/7. Figure 14 shows the dependence of x ′ crit on the ratio s = ρ 2 /ρ 1 and the velocity ratio ∆U/(U 1 +U 2 ). In this figure, we show how the results from our shear layer simulations of Section A compare to the correlation of equation (34), when s = ρ 2 /ρ 1 = 7, which demonstrates very good agreement. For the conditions of the situation examined in this study, i.e. Fig. 11, we note that ρ 2 /ρ 1 is actually closer to 1/7 since the fast stream, in the steady-state frame of reference, corresponds to the heavier shocked inert gas.

To examine in more detail, we chose to analyse the situation most prone to K-H shearing, when M c is a maximum. From analysis conducted in Fig. 11, this occurs at α = 25 • . For this case, we repeated the inert simulation described in Section IV B, in the frame of reference of the node, but on a much larger domain. In this new simulation, shown in Fig. 15, the domain size considered was 8000∆ 1/2 (800 mm) by 11,776∆ 1/2 (1177.6 mm), nearly ten times larger than the scales observed in the Lieberman and Shepherd experiments. [START_REF] Lieberman | Detonation interaction with an interface[END_REF] The finest grid resolution was kept unchanged, to maintain consistency with our earlier results. This new simulation shows that visible K-H structures do not appear until the shear layer is about 1 m away from the node. Upon time-averaging the flow fields, the size of shear layer (δ ) vs. distance from the node, x ′ , was measured (see Fig. 16). Also shown in this figure is the Reynolds number vs. x ′ , which permits us to determine x crit when Re crit = 10, 000 from Eq. 33. For this particular case, a critical Reynolds number This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. is not achieved until the flow is 1.11m away from the node, and is shown in Fig. 14, which compares well to the correlation for s = ρ 2 /ρ 1 = 1/7. Moreover, the shear layer growth rate at this location was dδ /dx ′ = 0.019, which corresponds to θ TMZ = 1.09 • . This value is significantly greater than those predicted earlier, on much smaller domains, as critical Reynolds numbers were not achieved. With spreading rates up to dδ /dx ′ = 0.019, it is not possible to achieve x ′ crit due to K-H shearing alone within 100 mm from the node suggesting that three-dimensional shearing effects are likely not significant on the scales observed in the past experiments, [START_REF] Lieberman | Detonation interaction with an interface[END_REF] and our corresponding CLEM-LES investigations. In the vicinity of the node, i.e. within 100 mm, we argue that pressure pulses and transverse waves are the most likely contributors to the observed TMZ growth.

PLEASE CITE THIS ARTICLE AS

Finally, in spite of the fact that our statistics reveal selfsimilarity of the flow, the proximity of the shock wave to the TMZ could possibly trigger the development of instabilities in the TMZ sooner than in canonical shear flow experiments. Physically, one would expect that a portion of the acoustic disturbances / pressure waves originating from the detonation-inert interface interaction would reflect from the oblique shock wave and interact in non-trivial ways with the TMZ. This could potentially lead to the destabilization of the shear layer itself. The extent of the latter interaction, however, remains to be characterized in detail. The influence of threedimensional cellular structure effects, and also the effects of realistic and variable-γ values will also be investigated in future work.

V. CONCLUSIONS

Detonation propagation into sharp interfaces of reactive and inert gas was investigated using the CLEM-LES framework, and validated with previous experiments. [START_REF] Lieberman | Detonation interaction with an interface[END_REF][START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] The full-scale simulations reproduced the experimental flow features (quali- , which is still less than actual TMZ growth rate. tatively and quantitatively), which consisted of a cellular detonation front, a transmitted shock wave-TMZ complex, and a Mach stem. In addition, measurement of the resulting TMZ growth rates from the CLEM-LES showed good agreement with experimentally reported values (within ∼ 2 • ). Upon recasting the simulations in the frame of reference of the node, and by removing instabilities at the detonation front, shear growth rates were found to be insignificant when perceived from the node. This was further substantiated by performing a detailed gas-dynamic and shear layer growth analysis, in both, the laboratory-and node-frame of reference. Most notably, perturbing the detonation front pressure, in the inert node-frame simulations, showed that the observed shear layer growth is heavily influenced by baroclinic instabilities originating at the detonation front. However, in order to account for the full growth rate of the TMZ, as observed in experiments, transverse waves associated with the detonation front cellular structure are likely the largest contributor through additional Richtmyer-Meshkov instabilities. In future work, to better capture all of the gasdynamic features present, the combustion reaction mechanism description of the CLEM-LES methodology should be extended to include multiple steps and more realistic properties such as temperature-dependent heat capacities. Three-dimensional simulations should also be sought.

We simulated several inert shear layers in both two and three-dimensions, for different stream velocity ratios, and at various resolutions to ensure that our LES strategy adequately captures the correct shear layer growth for conditions expected in our TMZ. Since there are no chemical reactions, closure to ω is not required via the LEM subgrid. In this case, the development of the inert and turbulent shear layers is computed entirely on the large-scale LES, Eqs. ( 1) to (4). The shear layers simulated have initial conditions relevant to the conditions encountered in this study, i.e. a fast moving high temperature hot stream mixing with a slower moving low temperature cold stream, and both streams having the same molecular weight. The density ratio for all shear layer simulations was kept at ρ 2 /ρ 1 = 7, diffusion coefficients and specific heat ratio were kept the same for both streams (ν = 0.00786, Pr = 0.77, Pr t = 1, γ = 1.126 ); the initial pressure was chosen to be close to the pressure expected in the TMZ of the detonation interaction with a sharp interface at α = 45 • ( p 0 = 280 kPa). The velocity of the low temperature, high density stream (ρ 2 = 1.8kg/m 3 ) was kept constant at U 2 = 0.21 (68.4 m/s), while the velocity of the hot low density stream, U 1 , was varied. The domain size simulated was x = 1536∆ 1/2 (153.6 mm) in length by y = 1024∆ 1/2 (102.4 mm) in height, with the splitter plate located at y = 512∆ 1/2 (51.2 mm). This domain size was sufficient to capture the detailed large-scale shearing, without influencing the solution. For three-dimensional simulations, a width of z = 36∆ 1/2 (3.6 mm) spanning the splitter plate with periodic boundary conditions was specified. The spatial growth of the shear layer size, δ (x), as a function of distance x downstream from the initial splitter plate, was determined by averaging temperature fields from 750 instances in time and recording the locations where the temperatures changed from their initial value by 4%. This particular threshold was used for consistency with our measurements of shear growth in Section III B. Each simulation was run for 7500 time units (2.29 ms) prior to averaging to allow the shear layer to develop to a quasi-steady state; the averaging was carried out over the same amount of time (2.29 ms). This was sufficient to give statistically resolved measurements in all cases.

For verification, Fig. 17 shows the effect of grid resolution on the shear layer growth rate (dδ /dx) for U 2 /U 1 = 0.257, for both two-and three-dimensional simulations. For the two-dimensional simulations, the shear layer growth rate converges at a resolution of at least 2 grids per ∆ 1/2 (0.05 mm). While we did not conduct three-dimensional simulations beyond this resolution due to resource limitations, it seems that three-dimensional effects lead to increased shear growth rates; the difference in shear growth rates between the twodimensional and three-dimensional cases are less than 10%. For validation, Fig. 18 shows an example of a resolved instantaneous density gradient flow obtained from a twodimensional numerical simulation, compared directly to an experimental schlieren image, to scale, for a case where U 2 /U 1 = 0.143 and ρ 2 /ρ 1 = 7. The resolution in this case was 4 grids per ∆ 1/2 (0.025 mm). Also shown in this figure is the location of Re δ ,crit = 10, 000. It is above this local Reynolds number that the flow is able to sustain threedimensional fluctuations, [START_REF] Konrad | An Experimental Investigation of Mixing in Two-Dimensional Turbulent Shear Flows with Applications to Diffusion-Limited Chemical Reactions[END_REF][START_REF] Dimotakis | The mixing transition in turbulent flows[END_REF] and therefore influence mixing in the third dimension. Konrad [START_REF] Konrad | An Experimental Investigation of Mixing in Two-Dimensional Turbulent Shear Flows with Applications to Diffusion-Limited Chemical Reactions[END_REF] argued that even though three-dimensional effects are present beyond this Reynolds number, the shear layer continues to grow at the same linear rate as below this threshold, and that the large-scale flow structures essentially remain two-dimensional. The twodimensional simulation in Fig. 18 captures the qualitative behaviour when compared to the experiments of Brown and Roshko. [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] Finally, Fig. 19 quantitatively compares the simulation shear layer growth rates to experiments for different velocity ratios. The resolution for the two-dimensional simulations was 4 grids per ∆ 1/2 (0.025 mm) whereas the three-dimensional simulation resolution was 2 grids per ∆ 1/2 (0.05 mm). Results show that our LES strategy correctly captures the experimental shear layer growth rates expected for a wide range of velocity ratios.

FIG. 1 .

 1 FIG. 1. Initial fields for sharp interfaces for (a) α = 15 • to 60 • , and (b) α = 0 • . The dark regions represent the reactive fuel-oxidizer mixture ( Ỹ = 1), while the white regions represent the inert gas (or detonation products -Ỹ = 0). A ZND detonation profile is initialized at x = 0, with a Dirichlet type wall boundary condition imposed at the left (at x = -100∆ 1/2 ), and Neumann type boundaries elsewhere.

FIG. 2 .

 2 FIG.2. Unsteady one-dimensional laminar flame speed, obtained for different resolutions using the one-step model shown in Eq. (12), for an initial state corresponding to the post-shock conditions of M D = 0.5M CJ . The flame speed evolutions are found to converge to the solution obtained using Cantera[START_REF] Goodwin | Cantera: An objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes[END_REF] and the Wang et al. mechanism.[START_REF] Wang | A comprehensive mechanism of C2Hx and C3Hx fuel combustion[END_REF]

FIG. 3 .FIG. 4 .

 34 FIG. 3. Portion of numerical soot foils fora) ∆ = ∆ 1/2 /4, N = 16 ( ∆e f f = ∆ 1/2 /64), b) ∆ = ∆ 1/2 /8, N = 16 ( ∆e f f = ∆ 1/2 /128), c) ∆ = ∆ 1/2 /8, N = 32 ( ∆e f f = ∆ 1/2 /256), and d) ∆ = ∆ 1/2 /16, N = 16 ( ∆e f f = ∆ 1/2 /256). The cell sizes found from autocorrelation were a) λ = 5.1 mm, b) λ = 2.8 mm, c) λ = 3.0 mm, and d) λ = 2.4 mm.

A

  . Overall flow structure for α = 45 • and α = 0 • A snapshot of the resulting numerical density gradient field and an experimental schlieren image from Ref. 26 for the corresponding experiment are shown in Fig.

FIG. 5 .

 5 FIG. 5. Schlieren (density gradient) images of detonation propagation through sharp interfaces at (a) α = 45 • and (b) α = 0 • . Chemical reaction rate locations of ω > 5.9 × 10 5 kg/m 3 /s are superimposed in red. Note: The experimental viewing window was 150 mm diameter, while the height of the simulation domain was 100 mm. The experimental image was reproduced from Lieberman and Shepherd [Phys. Fluids 19, 096101 (2007)],[START_REF] Lieberman | Detonation interaction with an interface[END_REF] with the permission of AIP Publishing.

FIG. 6 .

 6 FIG. 6. Time-averaged numerical schlieren images for (a) α = 0 • , (b) α = 15 • , (c) α = 30 • , (d) α = 45 • , and (e) α = 60 • .

1 FIG. 7 .

 17 FIG.7. Averaged density (top) and temperature (bottom) profiles for α = 0 • at various positions along x. Other interface angles (α) were also found to produce self-similar profiles.

1 w 5 =w 1 w 3 w 4 FIG. 8 .

 1348 FIG. 8. Detonation refraction configuration showing states (1)-(5): (a) particle paths P1 and P2 originating from states 1 and 5; (b) wave angles. D, T, RE, CS, and CS' refer to the detonation front, transmitted shock, rarefaction fan, contact surface, and deflected contact surface (shear layer), respectively. Diagram adapted from Ref. 26.

9 .

 9 FIG. 9. Transmitted shock angles, β , as a function of interface angle α.

FIG. 10 .

 10 FIG. 10. Density gradient fields for inert simulations in the frame of reference of node (where the detonation meets the reactive-inert gas interface), for (a) α = 0 • , (b) α = 30 • , and (c) α = 60 • .

)

  FIG.11. Density ratio, s, velocity ratio r, and convective Mach number, M c as a function of interface angle α.

Fig. 12 .

 12 Fig.12. In this figure, the error bars represent the range of angles obtained for 0.25 < C δ < 0.45. For all three models and α values, θ < 1 • . This is much less than the Liberman model predictions, experiments, and the shear growth rates measured from the CLEM-LES simulations, all of which are compared in Fig.12. Again, despite Lieberman's agreement, there is some inconsistency in using the laboratory-frame of reference to measure the shear-growth rates of the TMZ due to the fact that the flow field is not quasi-steady. For the Lieberman model, and our temperature dependent model both exhibit vanishing θ as α → 0. This is contrary to the CLEM-LES observations, which suggest there may be other significant contributing factors to the TMZ growth. As will be discussed in the next section, we suggest that pressure pulses associated with the cellular nature of the detonation front significantly influence the TMZ growth.Upon measuring θ in the inert node-frame simulations, using the same procedure of the lab-frame CLEM-LES simulations, good agreement with our predicted shear growth angle of θ was found. In fact, θ obtained from the inert node-frame simulations was also always less than 1 • . Despite the small difference in θ predicted from the gas-dynamic models and the inert node-frame simulation, it is reasonable to speculate that shearing due to K-H is minimal, and only contributes to a fraction of the actual TMZ growth as observed in experiments and lab-frame CLEM-LES.

FIG. 13 .

 13 FIG. 13. Time-averaged density fields of inert simulations in the node frame of reference without (a) and with (b) pressure perturbations at the inlet. Frame (c) shows an instantaneous density field for the case with forced pressure perturbations.

1 •

 1 FIG. 15. (a) Time-averaged, and (b) instantaneous density gradient fields of large-scale inert simulations in the node frame of reference (no forcing of pressure perturbations).

  FIG.16. Shear layer growth for α = 25 • as a function of distance along the shear layer, x ′ . The TMZ growth rate eventually corresponds to θ TMZ = 1.09 • , which is still less than actual TMZ growth rate.

FIG. 17 .

 17 FIG. 17.Influence of LES grid resolution ( ∆) and multidimensionality on shear layer growth rate (dδ /dx) for U 2 /U 1 = 0.257 and ρ 2 /ρ 1= 7.

FIG. 18 . 44 Accepted

 1844 FIG.18. Numerical schlieren image of a shear layer compared to experiments of Brown and Roshko,[START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] to scale. Both the simulation and experiment have U 2 /U 1 = 0.143 and ρ 2 /ρ 1 = 7. Also shown is the numerical location along x where Re δ ,crit = 10, 000. Note, dimensions are given in ∆ 1/2 . The experimental image was reproduced from Brown and Roshko [J. Fluid Mech., 64, 775-816 (1974)],[START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] with the permission of Cambridge University Press.

2 1 w 5 =w 1 w 3 w 4 1 •

 21341 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5(γ=1.126) Perfect gas model (multiple γ) Realistic thermo-chemical properties CLEM-LES Inert simulation (node-frame) Liebermann (2007) y a) α=0°: b) α =30°: c) α =60°: This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5(γ=1.126) Perfect gas model (multiple γ) Realistic thermo-chemical properties Inert simulation (node-frame) Lieberman ((γ=1.126) Perfect gas model (multiple γ) Realistic thermo-chemical properties Inert simulation (node-frame) Lieberman (2007) model This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5Perfect gas model (γ=1.126) (Mc 1 ) Perfect gas model (multiple γ) (Mc 2 ) Perfect gas model (multiple γ) (Mc 1 ) Realistic thermo-chemical properties (Mc 2 ) Realistic thermo-chemical properties (Mc) Inert simulation (node-frame) (Mc 1 and Mc 2 ) Lieberman (2007) model This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5(γ=1.126) Perfect gas model (multiple γ) Realistic thermo-chemical properties Inert simulation (node-frame) CLEM-LES Lieberman (2007) model Lieberman (2007) experiment 10.1063/5.0113073 400 y a) α=45 • no forcing at CJ boundary (time-averaged) This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5Simulation, s=0.130 Curve fit, s=7 Curve fit, s=1/7 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5rate: dδ/dx'=0.019 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5LES grids per ∆ 1/2 (grids per 0.1mm) 2D 3D

TABLE I .

 I Shock, TMZ, and gap angles for detonation interaction with sharp interfaces, obtained numerically and compared to experiments.[START_REF] Lieberman | Detonation interaction with sharp and diffuse interfaces[END_REF] .33 ± 0.03 • 26.33 ± 0.03 • 5.3 ± 0.1 • 2.0 ± 0.1 • 15 • 1.126 43.44 ± 0.06 • 28.44 ± 0.06 • 5.4 ± 0.5 • 1.9 ± 0.2 • 30 • 1.126 57.82 ± 0.07 • 27.82 ± 0.07

	α 0 •	γ 1.126 26	β	(β -α)	θ TMZ	θ gap

Accepted to Phys. Fluids 10.1063/5.0113073

  Origins of instabilities in turbulent mixing layers behind detonation propagation into reactive-inert gas interfaces 11

	100
	10
	1
	θ (˚)
	0.1
				Perfect gas model (γ=1.126)
				Perfect gas model (multiple γ)
	0.01		Realistic thermo-chemical properties Inert simulation (node-frame)
						CLEM-LES
					Lieberman (2007) model
				Lieberman (2007) experiment
	0.001
	0	10	20	30	40	50	60
				α (˚)
	FIG. 12. TMZ growth angles (θ ) as a function of interface angle α
			.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073

ACKNOWLEDGMENTS

This research was enabled in part by high performance computing resources provided the Core Facility for Advanced Research Computing at Case Western Reserve University. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0113073

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.


C Raimbault

P Laure

G Francois

S A E Boyer

F Choquart

J F Agassant

Injection of polyurethane foaming on a porous medium

Keywords: mousse polyuréthane, feutre, injection réactive, simulation polyurethane foam, felt, reactive injection, simulation 1

. Il s'agit d'une plaque rectangulaire injectée par son centre. Le matériau poreux considéré est un feutre constitué de fibres de polyester (PET). Par ce type de calcul, on peut donner une évaluation de la profondeur de pénétration dans le milieu fibreux.

Introduction

Polyurethane foams are classically used due to their variably adjustable densities as well as their ability to produce either rigid or soft complex forms. Soft-elastic polyurethane foams are rather used in the upholstery area for car seats or mattresses, whereas hard foams have mainly been used as thermal insulation. The epoxy resins or thermoset polymers are used to produce composite. For example, they can contain random continuous glass fiber mats [2]. In other applications, a textile structure is added in the mold, and it acts as a porous medium with both large and small pores. Therefore for producing lightweight composites, polyurethane foam can be advantageously reinforced with fabrics or felts. The standard procedure is to inject the reacting polyurethane mixture through a special designed runner system in a closed mold containing the porous medium [1]. In other applications, the expanding foam must impregnate free space between fibers during this RIM-process. In the same way, the creation of sandwich design [3] in which the polyurethane foam allows increasing lightweight. In this context, the adhesion between foam and the surrounding materials is an important issue [4]. Our study devoted to automotive industries consists in associating a polyester felt to a flexible polyurethane foam. The objective is to improve both mechanical and acoustic properties of the foam. In this context, the evaluation of impregnation depth is important as it gives information on the adhesion between the foam part and the felt components. There are already numerous articles which are focused on modelling the injection-expansion of polyurethane foam in a mold [5,14,15]. However, the modelling and simulation of the combined foam expansion and impregnation process is still a challenging task.

Modelling the flow motion in porous media

The modelling without porous media has been already presented in a previous paper [5] and we present here the modifications to handle this porous medium. To model the filling of a mold cavity partially filled by a felt, the momentum, energy and reaction kinetics are required. The porous medium can be characterized by its porosity 𝜑 and permeability 𝐾 . The stress tensor is assumed to follow a Newtonian law and it has this expression as the foam is compressible:

𝝈 = -𝑝 𝐈 + 2 𝜺̇- 1 3 tr𝜺̇ 𝐈 (1)
𝜺̇ is the rate of strain tensor and tr𝜺̇ its trace, I is the identity tensor, 𝑝 is the pressure and 𝜂 the viscosity which is a complex function of temperature, gas content and curing reaction.

This simplest approach to deal with porous media is to mix the Stokes equation and the Darcy's law. That gives the Navier-Stokes-Brinkman macroscopic model which has been successfully used in simulations of flow in porous media. If the inertial term is neglected, the resulting equation will be valid in the whole computational domain:

∇. 𝝈 +  𝐾 𝒖 + 𝜌 𝒈 = 0 (2)
where the viscosity  and the permeability 𝐾 are space dependent and are defined on the whole computational domain. A large value of 𝐾 is chosen outside the porous medium to have the ratio  𝐾 ⁄ less than 10 . Note that the velocity in equation ( 2) is the average velocity through the porous medium whereas the fluid velocity inside the porous medium is equal to 𝒖 𝜑 ⁄ . Therefore, the mass balance inside the porous medium is modified as follows:

𝛁. 𝒖 = 𝜑 1 - 𝑑 𝑑𝑡
where  is the gas content defined as 𝑉 𝑉 ⁄ .

(

) 3
The equation describing the average temperature 𝑇 in the porous medium is written as:

𝜌𝑐 𝜕𝑇 𝜕𝑡 + 𝜑 𝜌 𝑐 𝒖 . ∇𝑇 = ∇ • 𝜆 ∇𝑇 + 𝜑  𝜀̅ ̇ + 𝑑𝑇 𝑑𝑡 (4)
where the heat generated by reaction kinetics is given by the term 𝑇 . In this equation, the temperature represents the thermal equilibrium between the fluid and the porous medium. Therefore, the mean value 𝜌𝑐 can be deduced from the density 𝜌 and the specific heat 𝑐 of each component by a linear mixing law [START_REF] Bruschke | [END_REF] as described by this formula:

𝜌𝑐 = 𝜑 𝜌 𝑐 + (1 -𝜑 )𝜌 𝑐 (5)
where the solid and fluid parts are labelled by the subscript p and f respectively. The equivalent thermal conductivity 𝜆 is more complicated to express [START_REF] Fournier | Optimisation et modélisation du procédé R[END_REF]. For foaming flows the Peclet number is usually very low, hence we assume that the heat dispersion effect due to the pores is insignificant. The effective estimation of this parameter can be made by performing numerical computations on a REV. The equivalent conductivity depends on the main flow direction with respect to felt orientation: if the fluid goes through the porous medium, a harmonic mixing law can be used whereas it is a linear law if the flow is parallel to the porous medium. However, if the flow is at rest, it is shown [START_REF] Fournier | Optimisation et modélisation du procédé R[END_REF] that the linear mixing gives suitable approximation. In our case the thermal regulation spends more time than foaming, we choose a linear law comparable to relation (5) for the thermal conductivity.

Another novelty corresponds to a new treatment of the pressure. Indeed, the bubble size will exceed the maximal size allowed by the porous medium. Therefore, the bubbles are cut through and this leads in the porous medium to lower gas volume fraction, higher polyurethane volume fraction and eventually higher final mass. This modelling allows a relaxation of the pressure as soon as the gas created during foaming exceeds a critical value [START_REF] Raimbault | Moussage de polyuréthane en injection sur un milieu poreux[END_REF] which is written as follow:

1 𝑝 𝑑𝑝 𝑑𝑡 = -(𝜀̅ ̇ , ) (6)
The relaxation term  also depends on a critical strain dependent term 𝜀 ̇. That reflects the fact that the gas bubbles open more easily under the effect of shear in addition of the size of void in the porous medium. This term can be adjusted by comparison with experimental measurement.

The equations are solved in the whole computational domain 𝛺 with a finite element method. More precisely, all subdomains are embedded in a unique Eulerian mesh using the immersed volume method [START_REF] Coupez | Valette In Arbitrary Lagrangian-Eulerian and Fluid-Structure Interaction[END_REF]. A mixed finite element formulation in velocity and pressure is used for Brinkman equations [START_REF] Pacquaut | Couplage Stokes/Darcy dans un cadre Level-set en grandes déformations pour la simulation des procédés d'élaboration par infusion de résine[END_REF][START_REF] Pacquaut | [END_REF]. Two level-set functions [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] are used to get implicit representation of both foam/air interface and felt domain. The different sub-domains are depicted in Figure 1.

A) Before impregnation B) After impregnation

Fig. 1: sketch of Multiphysics approach: notation and definition of different domains (A) before and (B) after percolation in the porous medium.

In this way, one gets a proper definition of the viscosity , the permeability 𝐾, the density 𝜌, the specific heat c and the thermal conductivity  as space dependent functions. A global property 𝐹 is defined using two Heaviside functions 𝐻 and 𝐻 deduced from these level set functions [START_REF] Nakhoul | [END_REF]:

𝐹 = 𝐹 𝐻 + 𝐹 1 -𝐻 1 -𝐻 + 𝜑 𝐹 𝐻 + 𝐹 1 -𝐻 + 1 -𝜑 𝐹 𝐻 (7)
where the subscripts f, p and a are respectively associated to foam, porous and air phases. Moreover, as explained in [START_REF] Nakhoul | [END_REF], the level-set associated to the fluid/air interface is transported by the velocity solution of the equations (1-3) describing the foam growth.

Results

The characteristics of injected polyurethane come from our previous article [5]. As showed in Figure 2, the mold consists of a square metallic frame (size: 900 mm, thickness whithout the felt: 16 mm) closed by two horizontal plates thermally regulated at 55 °C. The fibrous layer is positioned between the frame and the lower plate. Except for the mass of injected mixture which depends on the absence or presence of the fibrous medium, the injection parameters are identical: the injection rate is equal to 300 g.s -1 , the injection temperature is equal to 30°C and the temperature of the upper and lower walls of the mold is equal to 55°C. The variable parameters are the mass of injected mixture, the thickness of the fibrous support as well as its density, its porosity and its permeability.

Fig. 2: Vertical cut-plane of mold geometry

To study the mold filling, the case without fibers and an injected mass of 720 g is compared to the case with a low density and thick fibrous layer (5 mm). In this later case, the injected mass is increased to 740 𝑔 to fill up the cavity. The characteristic parameters of the felt are 𝜌 = 153.7 kg.m -3 , 𝜑 = 0.902, 𝐾 = 3.7.10 -10 m 2 . The mold filling dynamics are computed and shown in Figure 3. The main cavity is filled after 33 s in the presence of the porous layer and after 26 s otherwise. This shows that the presence of the porous layer slows down the filling because the foaming kinetics have not changed. Furthermore, the porous layer was impregnated with a volume of 285 cm 3 of mixture (curve in small dotted lines), which corresponds to 28% of its volume. In addition, another computation is made with a felt having thinner thickness (3 mm) and higher density (𝜌 = 407.7 kg.m -3 , 𝜑 = 0.721, 𝐾= 3.7 10 -11 m 2 ). This new computation shows a lower decrease in density for a high-density medium [START_REF] Raimbault | Moussage de polyuréthane en injection sur un milieu poreux[END_REF]. In the same way, the impregnation depth is less important (Figure 6).

Conclusions

The modelling of polyurethane foaming on a porous is studied by a multiphase approach. This modelling considers a Navier-Stokes-Brinkman flow for the foam in the porous medium. A pressure relaxation model is introduced, and it is possible to account for the decrease in pressure linked to the evacuation of the gas after the opening of the cells due to the porosity of the porous medium. Based on the foaming experiment in a fibrous material [START_REF] Raimbault | Moussage de polyuréthane en injection sur un milieu poreux[END_REF], a numerical estimate of the value of this coefficient is proposed. A sensitivity study also showed the strong influence of the permeability of the fibrous medium and the viscosity of the foam on the impregnation distance [START_REF] Raimbault | Moussage de polyuréthane en injection sur un milieu poreux[END_REF]. The model equations were implemented in REM3D® and the computations on an industrial geometry are presented.

Simulation results showed that the presence of the porous layer increases the time to fill the free cavity. Also, the foam has a lower density, and the size of pores is larger. The temperature field is in the mold The temperature in the foam is indeed higher near the porous medium, and it evolves linearly in the non-impregnated porous medium up to the temperature boundary condition imposed on the lower wall. The impregnation is compared between low density and high-density porous media. The impregnation is more important for low density fibrous media in the first third of the plate on the side of the injection point. This zone may correspond to the initial spreading distance of the injected mixture before the starting time of foaming, which can vary depending on the initial viscosity and the cream time of the foam.

Fig. 3 :

 3 Fig. 3: Evolution of the volume of foam without and with the porous layer Figures 4 compare the temperature fields in the mold thickness at the end of the foaming during the polymerization (t = 100 s) in configurations with and without the porous layer. They point out an asymmetry of temperature in the foam thickness.

Fig. 4 :Figure 5

 45 Fig. 4: Temperature field in the mold thickness at t=100 s in the cut vertical plane without and with the porous layer

Fig. 5 :

 5 Fig. 5: Density in the mold thickness at t=100s without and with the porous layer

Fig. 6 :

 6 Fig. 6: Impregnation depth along the radial component from injection point for felts having low and high densities.