Shock Mach number influence on reaction wave types and mixing in reactive shock–bubble interaction Academic research paper on "Chemical engineering"

Contents lists available at ScienceDirect

Combustion and Flame

journal homepage: www.elsevier.com/locate/combustflame

Shock Mach number influence on reaction wave types and mixing in reactive shock-bubble interaction

CrossMaik

Felix Diegelmanna'*, Stefan Hickela,b, Nikolaus A. Adams3

a Institute of Aerodynamics and Fluid Mechanics, Technische Universität München, Garching 85748, Germany b Faculty of Aerospace Engineering, TU Delft, 2629 HS Delft, Netherlands

A R T I C L E I N F 0

Article history: Received 4 March 2016 Revised 3 August 2016 Accepted 13 September 2016

Keywords: Shock wave

Richtmyer-Meshkov instability Shock-bubble interaction Detonation Deflagration

A B S T R A C T

We present numerical simulations for a reactive shock-bubble interaction with detailed chemistry. The convex shape of the bubble leads to shock focusing, which generates spots of high pressure and temperature. Pressure and temperature levels are sufficient to ignite the stoichiometric H2 -O2 gas mixture. Shock Mach numbers between Ma = 2 . 13 and Ma = 2 . 90 induce different reaction wave types (deflagration and detonation). Depending on the shock Mach number low-pressure reactions or high-pressure chemistry are prevalent. A deflagration wave is observed for the lowest shock Mach number. Shock Mach numbers of Ma = 2 .30 or higher ignite the gas mixture after a short induction time, followed by a detonation wave. An intermediate shock strength of Ma = 2 . 19 induces deflagration that transitions into a detonation wave. Richtmyer-Meshkov and Kelvin-Helmholtz instability evolutions exhibit a high sensitivity to the reaction wave type, which in turn has distinct effects on the spatial and temporal evolution of the gas bubble. We observe a significant reduction in mixing for both reaction wave types, wherein detonation shows the strongest effect. Furthermore, we observe a very good agreement with experimental observations.

© 2016 The Authors. Published by Elsevier Inc. on behalf of The Combustion Institute.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/40/).

1. Introduction

The interaction between high-speed reactive flows and shock waves is a generic situation present in many combustion systems. Controlled application can promote mixing; uncontrolled interactions, however, can lead to undesirable heat release and thermomechanical loads. Especially in supersonic combustion, where the rapid and efficient mixing of fuel and oxidizer is crucial, as the residence time of the fuel-oxidizer mixture in the combustion chamber is only a few milliseconds [1], mixing can be enhanced sufficiently by shock-induced instabilities. The selected generic configuration of reacting shock-bubble interaction (RSBI) is representative for a large range of hydrodynamic instabilities and different reaction wave types occurring in application, and allows us to study the interaction between different effects in detail.

1.1. Hydrodynamic instabilities

Two hydrodynamic instabilities dominate in a RSBI: the Richtmyer-Meshkov instability (RMI) and the Kelvin-Helmholtz instability (KHI). RMI can enhance mixing in high-speed reactive

• Corresponding author. E-mail addresses: felix.diegelmann@aer.mw.tum.de, felix.diegelmann@gmail.com (F. Diegelmann).

flows, promote turbulent mixing and thus increase the burning efficiency of supersonic combustion engines [2]. The instability occurs at the interface between two fluids of different densities. Theoretically stated in 1960 by Richtmyer [3] and experimentally verified by Meshkov [4] in 1969, RMI can be considered as the impulsive limit of the Rayleigh-Taylor instability [5,6]. The misalignment of pressure gradient, Vp , associated with a shock wave and density gradient, V p, at the material interface causes baroclinic vorticity production at the interface. For comprehensive reviews the reader is referred to Brouillette [7] and Zabusky [8]. RMI occurs on a wide range of highly reactive environments from extremely large scales in astrophysics [9], to intermediate scales in combustion [1,10] and down to very small scales in inertial confinement fusion [11].

RMI induces velocity shear and small perturbations at the interface of the bubble, which are necessary preconditions for KHI [12]. The perturbations are amplified, eventually generating vortices at the interface accompanied by the appearance of smaller scales [7]. KHI drives the breakup of large-scale structures [13] and forces mixing [14]. Both effects are the main hydrodynamic drivers in RSBI.

1.2. Shock-induced ignition and reaction waves

Independently of the scale, RMI is accompanied by a second phenomenon in reactive gas mixtures: the shock-induced variation

http://dx.doi.org/10.1016/jxombustflame.2016.09.014

0010-2180/© 2016 The Authors. Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

of thermodynamic properties, which can lead to ignition, followed by a reaction wave. Two reaction wave types can be distinguished: deflagration and detonation. Deflagration is a subsonic diffusion-driven reaction wave that propagates through the gas mixture due to direct transfer of chemical energy from burning to unburned gas [15]. Detonation is driven by a fast chemical reaction and the associated large heat release within the reaction wave. A shock wave immediately precedes the detonation wave and preheats the gas mixture by compression [15]. The detonation wave propagates up to 108 times faster than the deflagration wave [16]. Due to the large differences in the characteristic reaction time scales, the reaction wave type has a crucial influence on the flow evolution.

Under certain circumstances a deflagration wave can transform into a detonation wave. Deflagration-to-detonation transition (DDT) is one of the most interesting unresolved problems in combustion theory. Generally, a self-propagating deflagration wave is unstable and tends to accelerate. Under specific conditions the continuous acceleration can suddenly transition into a detonation wave [17]. Liberman et al. [18] proposed a mechanism mainly driven by flame acceleration divided into three stages. The reaction front accelerates and produces shock waves far ahead of the flame. Thereafter, the acceleration decreases, shocks are formed on the flame surface and pockets of compressed and heated unburnt gas emerge (preheat zone). In the final stage the transition to detonation happens: the flame propagates into the preheat zone and produces a large amplitude pressure pulse. Increasing pressure enhances reaction rates and the feedback between the pressure peak and the reaction leads to a growth of the pressure peak, which steepens into a strong shock that, coupled with the reaction zone, finally forms an overdriven detonation wave.

Furthermore, the flame front can propagate into regions of gas that already have been compressed and preheated by preceding shock waves such as in shock-bubble interactions (SBI). The reaction rates and the heat release are enhanced in these regions, which in turn increases the pressure pulse and accelerates the transition to detonation. In general, DDT can occur in two regions: it develops from the preheated, compressed gas mixture between the leading shock wave and the flame or it arises from within the flame [19]. The latter transition process is relevant for the presented study as RSBI contains regions of irregular compression by the initial shock wave.

1.3. Reacting shock-bubble interaction

The impact of a shock wave on a reactive gas bubble allows to investigate the interaction between shock-induced hydrodynamic instabilities and ignition. The shock wave triggers RMI and the pressure and temperature increase leads to the formation of radicals, which accumulate until the gas mixture ignites. RMI, due to the misalignment of the pressure and density gradient at the bubble interface, causes the bubble to evolve into a vortex ring. Provided that the initial kinetic-energy input is sufficient, the flow develops a turbulent mixing zone through non-linear interactions of the material interface perturbations [7,8]. Upon contact, the incident shock wave is partially reflected and partially transmitted. For an Atwood number A = (p1 - p2) /(pi + p2) < 0 (the bubble gas is lighter than the ambient gas), the transmitted shock wave propagates faster than the incident shock wave. A > 0 shows the converse effect, the transmitted shock wave travels slower than the incident shock wave outside of the bubble. The transmitted shock wave focuses at the downstream pole of the bubble and collapses into a single point (shock-focusing point).

Classical inert SBI was the subject of several studies over the last decades. Haas and Sturtevant [20] investigated the interaction of shock waves propagating in air with a gas bubble filled with either helium or R-22. Their experimental results contributed to a

better understanding of the temporal bubble evolution under shock acceleration and established a new class of canonical flow configurations. These experimental findings were completed by the investigations of Quirk and Karni [21], providing detailed numerical results of shock-bubble interaction problems. They reproduced the transition from regular to irregular refraction, shock wave focusing and the formation of a jet towards the center of the bubble. For a detailed review of SBI see Ranjan et al. [22].

A new level of complexity can be added to the setup of SBI by replacing the inert gas with a reactive gas mixture. A strong shock wave can ignite the reactive gas mixture directly at the interface, whereas the additional increase of pressure and temperature in the shock-focusing point is required for ignition at lower shock Mach numbers. Two types have to be differentiated: non-premixed and premixed gas mixtures. Reacting SBI of non-premixed gas mixture was studied by Billet et al. [23]. In their setup a H2 gas bubble surrounded by air is shocked to study the influence of the volume viscosity on the bubble evolution and vorticity production. Attal et al. [24] verified the results of Billet et al. [23] and furthermore observed the formation of a double diffusion flame in the bridge region of the shocked bubble. Attal and Ramaprabhu [25] studied single-mode reacting RM in a non-premixed setup at different interface thicknesses. They observed shock-induced ignition and mixing enhancement by reshocking the propagating flame. Furthermore shock-flame interaction increases the surface area of the flame and the energy release and therefore the burning rate [26]. Massa and Jha [27] showed that small scales are damped by the shock wave and that the growth of RMI and KHI are reduced.

In 2012, Haehn et al. [28] investigated the interaction of a shock wave with a premixed gas bubble, filled with a stoichio-metric gas mixture of hydrogen (H2 ) and oxygen (O2 ), diluted by xenon (Xe). Besides triggering hydrodynamic instabilities, such as RMI, the shock wave also increases the temperature and pressure, which in turn induces faster chemical reaction rates up to the ignition of the gas mixture. Maximum pressures and temperatures are reached when the shock passes the bubble. Subsequently, the gas mixture relaxes and the two main parameters controlling the reaction rate, temperature and pressure, decrease. The experiments of Haehn et al. [28] covered both ignition types deflagration and detonation, by varying the shock wave Mach numbers between Ma = 1. 34 and Ma = 2 . 83.

A weak shock wave with Ma = 1. 34 does not ignite the gas mixture within the experimental timeframe. Compression is not sufficient to start a self-sustaining chemical reaction. An increase of the shock strength results in an ignition followed by a deflagration wave. The reaction wave type changes for higher shock Mach numbers; Haehn et al. [28] observed a detonation wave for Ma = 2 . 83 , even before the shock wave has reached the shock focusing point. Damkohler numbers between 0.25 (Ma = 1 . 65) and 8.00 (Ma = 2 . 83) were determined. Haehn et al. [28] conclude that heat conduction plays an important role at lower shock Mach numbers, and that the Zeldovich mechanism becomes important at higher shock Mach numbers. Their conclusion is consistent with the two limiting cases of shock-induced combustion, the strong and the weak ignition [19]. Strong ignition results in a detonation essentially initiated directly by the shock wave. Weak ignition is characterized by the appearance of small flames that can undergo transition into detonation waves. Several chemiluminescence exposures are provided by Haehn et al. [28] to depict the qualitative evolution of the bubble and reaction processes. Furthermore, quantitative data for the temporal evolution of the transverse diameter of the bubble as well as for the vortex ring diameter are presented. However, the complex experimental setup implies uncertainties. Haehn et al. [28] estimate the uncertainty of the Damkohler number at the highest shock Mach number (Ma = 2 . 83) of up to 50% (Da = 8 ± 4). At the lowest shock Mach number (Ma = 1.34), 30%

of all measurements showed no ignition within the given experimental time frame. Hence, numerical studies of RSBI can provide more certainty and complementary insight into RSBI phenomena that cannot be achieved by purely experimental work.

1.4. Scope of the present work

The present numerical study complements the work of Haehn et al. [28] and continues the first numerical approach to RSBI (Diegelmann et al. [29]). The main emphasis is placed on the general temporal and spatial evolution of RSBI, the comparison with SBI, and the dependence of the bubble evolution on the reaction wave type. In our study, the shock Mach number is varied between Ma = 2 . 13 and Ma = 2 . 90 at a constant initial pressure and temperature. Besides the limiting cases of deflagration and detonation we study two special phenomena in detail, which have not been discussed before: DDT at Ma = 2 . 19 and a double detonation at Ma = 2 . 50. Haehn et al. [28] observed an effect, which they assume is either a double detonation or a reflection of measurement signals, but the experimental measurement technique did not allow a clear identification. Our present numerical study confirms the observed physical effect and gives a deeper insight into the gas composition of the two ignition spots during the induction time. Intentionally, we focus on two-dimensional configurations as they facilitate particular analysis and phenomenological investigation. Moreover, in [29] it was shown that early stages of RSBI can be well reproduced by a two-dimensional approximation.

The chemical reaction rates of most gas mixtures increase with pressure. H2 -O2 reactions, however, show a different behavior [30]. Some intermediate reaction rates are proportional to the square of the pressure, others are linearly proportional [31]. Hence, the variation of the shock Mach number, or more precisely the post-shock pressure, affects the chemical reaction process and determines the occurrence of either detonation or deflagration.

We structure the paper as follows: Section 2 outlines the governing equations, including molecular transport properties for mul-ticomponent flows and chemical reaction kinetics. Section 3 describes the computational domain and the initial conditions of our setup. General results are discussed in Section 4. The spatial and temporal evolution of the RSBI are presented. The effect of different types of reaction waves on bubble deformation are compared with each other and with their non-reacting counterparts. The chemical reaction process during shock passage until ignition is analyzed in detail. A consistent definition of the dimension-less Damkohler number is used to evaluate whether hydrodynamic or chemical reaction time scales dominate the flow field. Integral quantities, such as enstrophy or the molar mixing fraction, are estimated to assess the effect of the reaction waves on mixing of the bubble gas. In Section 5, we discuss two special cases of RSBI: First the transition of a deflagration into a detonation wave and second a simulation with a simultaneous detonation at two spots. Section 6 presents a comparison to experimental results and a critical discussion. Section 7 summarizes the key findings.

2. Numerical model

+ V • [ (E + p)u] -V • ( t • u - qc - qa) - ocT = 0

dpY d t

+ V • (puYj ) + V • Ji - oc = 0

where p is the mixture density and u the velocity vector. The identity matrix is given by S, total energy by E and pressure by p. Y are the mass fraction of species i = 1) 2 i ..., N, with N being the total number of species. The heat release 0)T and species formation and destruction in terms of individual mass rates ocj represent the chemical reaction kinetics.

2.2. Caloric and transport properties

The viscous stress tensor t for a Newtonian fluid is given by

t = 2f

1 (Vu + (Vu)T) - 3S(V • u)

with f as the mixture viscosity

EN= i fYj /M) c 2

iN= 1 Yi /Mi 1 / 2

Mi is defined as the molecular mass of each species i. The calculation of the viscosity of each species ^ is based on the Chapman-Enskog viscosity model

f = 2 i6693 • 10-

; jW fifii of

where T is the temperature and a, the collision diameter. The collision integral Q,^., [32] is defined as

V^i = A(T*)B + C exp (DT*) + E exp (FT*). (8)

with A = 1. 16145. B = -0. 14874. C = 0.52487 . D = -0. 7732 . E = 2 . 16178 . F = -2 .43787 and T* = T / (e/k) , using the Lennard-Jones energy parameter (e /k) for species i . According to the Fourier law, we define the heat conduction as

qc = -KV T,

with k as the mixture heat conductivity, which is calculated from [33]

Eil 1 kY. /M. . 2

eN= i Y /M) c 2

ki is the thermal conductivity of species i. The interspecies diffu-sional heat flux qa [34] is given by

qa = X> Jc, (11)

with hi as the individual species enthalpy. The species diffusion Ji is modeled as

Jc = Dj V Y - Y £ D jVYj. (12)

Dj describes the effective binary diffusion coefficient of species i

2.1. Navier-Stokes equations

We solve the full set of compressible reacting multicomponent Navier-Stokes equations in conservative form

dp + V • (p u) = 0

+ V • (p uu + pS - t ) = 0

Dc = (1 -xc > j j ■ (13)

with Xj as the mole fraction of species ic Eq. (13) ensures that the interspecies diffusion fluxes balance to zero. The constitutive empirical law is used to compute the mass diffusion coefficient of a binary mixture [33]

0. 0266 T 3/ 2

Djj = fi

D,ij pc/MCjai2 '

2 , a + a j

M* j = -,-T- and aij =

ij( = J_ _i_ _L

Mi + Mj

The collision integral for diffusion aD, ¡j is given by

a,Di = A * (Tj * )B*+ C* exp (D* T( )

, ^ "»J- V 'ij / ' " V ■ 'ij .

The parameters are defined as A* = 1.06036. B* = -0 . 1561, C* = 0. 19300, D* = -0.47635 , E * = 1.03587 , F * = -1.52996 , G* = 1.76474 , H* = -3 . 89411, and T* = T/T,*ij . T*ij have been obtained from the Lennard-Jones energy parameters for species i and j as

Te„ =

(¡HU-

2.3. Equation of state

The equation of state for an ideal gas is used to close the equations

p(E, Y) = (Y - 1)E.

Y represents the ratio of specific heats of the mixture

cp - R

cp = cp,i ■ i=1

The specific gas constant of the mixture is defined by R = R/M, with R as the universal gas constant. M is the molar mass of the mixture

= J2XjMj ■

Cp i represents the specific heat coefficient

c nA —

Yi - 1

Ri ■

Rj is defined as Rj = R/Mj. The temperature is computed from P

2.4. Chemical reaction kinetics

The accurate calculation of chemical reaction kinetics is most important for the precise prediction of combustion effects, such as DDT. The review paper of Oran et al. [19] summarizes several studies about DDT, mainly operating with one-step chemical kinetics. DDT through the Zeldovich gradient mechanism was observed, arising due to the gradient of induction time within the hot spots in front of the flame, where temperature varies in the range of 60 0 to 80 0 K. A precise computation of the induction time and the corresponding heat release is therefore essential for an accurate description of DDT [18]. However, it was shown that the induction time of detailed mechanisms is larger than for one-step mechanisms [35] and also larger than the time between flame initiation and transition to detonation, which renders numerical results obtained with simple mechanisms questionable. Furthermore, important quantities of combustion such as detonation initiation and induction time in chain-branching kinetics are not correctly reproduced by one-step mechanisms [36]. Studies of Ivanov et al.

[36] reveal significant differences between the temperature gradient that leads to detonation with one-step and detailed mechanisms. For the detailed mechanism a much smaller temperature gradient is sufficient to ignite detonation, which is in accordance with the behavior of real combustible mixtures [18].

Chemical reaction kinetics are expressed by the heat release coT and species formation and destruction in terms of individual mass rates ¿Oj. The specific heat release coT is defined as

cot = -J2 Ahf,ivi,

with h0 i as the heat of formation of each species i. Mass rates ¿¿j for each species are estimated by

co* = Wj Vir Fr kfr nix* ] ^r - kbr niX* ] v r=1 \ j=1 j=1

with Nr as the number of reactions, Wj the molecular weight, rr the third body efficiency of reaction r, Xj the molar concentration, and v'jr and v" the molar stoichiometric coefficients of the reactant and the product of reaction r. vir is the net stoichiometric coefficient

Vir = vr - v'ir. (26)

The Arrhenius law is used to calculate the forward and backward reaction rates f and kbr. The forward reaction rates are defined as

k fr = A fr Tßf exp (

where Af r is the pre-exponential factor, Ef r is the activation energy, f)fr is the temperature exponent for each reaction r [37]. The backward reaction rates are calculated by using the equilibrium constants Kcr

K(r is given by

. = (P- r

/AS, exp (

with p° as a pressure of 1 atm, vr as the net change in the number of species in the reaction, AS.i as the net change in entropy and AH,i as the net change in enthalpy.

Furthermore, pressure dependent and duplicated reactions are considered; for this purpose Eq. (27) is modified. Pressure dependence is taken into account by calculating two forward reaction rates k fr. and k frm for the high-pressure and for the low-pressure limit, respectively. A blending function composed of these high- and low-pressure Arrhenius rate parameters is applied for a smooth pressure dependence. For more details on the so-called fall-off reactions, the reader is referred to Troe [38]. Duplicated reactions are considered by extending Eq. (27) to

kfr = J2Afr,Tß fr ex- (f

i=1 RT

The mechanism, which provides the parameters for the Arrhenius law, is essential for the accuracy of the numerical investigation and has to be chosen carefully. As shown by the authors in a previous publication [29], available mechanisms show large discrepancies in ignition delay time and pressure sensitivity. A certain number of intermediate reactions, third body efficiencies, duplicated and pressure dependent reactions are necessary for the accurate prediction of the reaction kinetics within the wide range of pressures and temperatures considered in this study.

We have chosen the Ó Conaire [39] reaction mechanism for the reaction rate parameters of the Arrhenius law. The mechanism is valid for a wide range of pressures (0.05-87 atm) and temperatures (298-2700 K). 8 + N species (two reactants: H2 , O2 ; 5 chain-carrying intermediates: hydrogen radical (H), oxygen radical (O), hydroxyl radical (OH), hydroperoxyl radical ( HO2 ), hydrogen peroxid (H2 O2 ); the product: hydrogen oxide (H2 O); N inert gases) and 19 intermediate reactions are considered, including duplicated and pressure dependent reactions as well as third-body efficiencies. Third-bodies absorb energy during the two-body recombination reaction and stabilize the final combination. The available modes for energy storage control the energy absorption. The third-body efficiencies of Xe are set identical to argon (Ar), which are provided by Ó Conaire [39]. As the available modes of Ar and Xe are identical, the third-body efficiencies can be assumed to be comparable. Also, the steric factor for monoatomic gases, which accounts for the geometry influence on the collision between molecules, is similar [40].

The mechanism of Ó Conaire [39] has been used widely in the recent years [41,42]. As part of a preceding validation campaign [29], the applied reaction mechanism has been compared to simpler reaction mechanisms. Accurate ignition delay times, crucial for the spatial evolution of the bubble and mixing, can only be achieved with a detailed description of chemistry by a sufficiently complex reaction mechanism.

2.5. Numerical method

The 2nd-order accurate Strang time splitting scheme [43] is used to solve the system of equations (Eq. (1-4)). The Strang splitting scheme separates the stiff terms, containing the chemical reaction kinetics (ó)T and ó)i ), from the Navier-Stokes equations. This results in a system of partial differential equations (PDE) and in a stiff system of ordinary differential equations (ODE). We use a finite-volume discretization scheme that applies a flux projection onto local characteristics for the hyperbolic part for the PDE system. The Roe matrix required for the projection is calculated for the full multi-species system [44,45]. The numerical fluxes at the cell faces are reconstructed from cell averages by the adaptive central-upwind 6th-order weighted essentially non-oscillatory (WENO-CU6) scheme [46]. The scheme uses a non-dissipative 6th-order central stencil in smooth flow regions and a non-linear convex combination of 3rd-order stencils in regions with steep gradients. Time integration is performed by the 3rd-order strongly stable Runge-Kutta scheme, developed by Gottlieb and Shu [47]. Our numerical model has been tested and validated for shock induced turbulent multi-species mixing problems at finite Reynolds numbers [13,48-50] and for shock-bubble interactions including chemistry [29]. The stiff ODE, governing the specific heat release and mass rates for each species, is separately solved by a variable-coefficient ODE solver using 5th-order backward differentiation formulas [51].

3. Computational setup

We study RSBI within a two-dimensional rectangular domain with a symmetry plane at the center axis of the bubble, see Fig. 1. Inflow boundary conditions are imposed at the left domain boundary and outflow boundary conditions at the right and upper domain boundaries. The domain size is set to 32.5 r x 10.5 r, with r as the initial bubble radius. The distance between the bubble and domain boundaries are chosen sufficiently large to avoid artifacts due to shock reflections. The Cartesian grid in the region of interest is refined by a factor of 25 compared to the coarse outer grid to reduce computational costs. A detailed grid study can be found in our previous paper. We demonstrated grid convergence by comparing four different grid resolutions, with cell sizes of Axy = 234, 117, 59 and 29 ]am in the high resolution part. The simulations are performed at a CFL-number of 0.3.

The gas bubble is filled with H2 , O2 and Xe in a stoichiometric composition of 2/1/3.67 mass fractions and surrounded by pure Nitrogen (N2 ). The bubble diameter is set to D = 2r = 0 , 04 m. The heavy inert gas Xe is used to increase the density of the bubble, leading to an Atwood number of A = 0 ,476, The gas composition of our domain and the bubble diameter are identical to the experimental setup of Haehn et al. [28]. A sharp and fully resolved interface between the bubble gas and its surrounding is defined in terms of the molar fraction of N2

tanh ((sjx2 + y2 - r)Ç) + 1

with r as the radius of the bubble and £ as parameter for controlling steepness, which is set to £ = 20,0 00. The molar fraction (X = 1 - XN;, ) inside the bubble is distributed among the three gases, ensuring the stoichiometric mixture with a relative composition of 2/1/3.67(H2 / O2 /Xe).

The shock wave is initialized on the left side of the bubble. The pre-shock state is defined by T0 = 350 K and p0 = 0 . 50 atm. The shock Mach number is varied between Ma = 2 . 13 and Ma = 2 . 90 . The post-shock thermodynamics state is given by standard Rankine-Hugoniot conditions

Pn , = Pn2

(Yn. + 1)Ma2

2 + (Yn 2 - 1)Ma2

UN2 = Ma Cn2 ^ - PN2)'

p,N2 = P0 (t + 2 ynYTT (Ma2 -1)2

where cN;. = ^/yN2 P0 /Pn. . Variables indicating post-shock conditions are marked with a prime.

Note that the initial parameters of our setup slightly deviate from the experimental pressure and temperature. To avoid that

Fig. 1. Computational domain of the RSBI, r = 0 . 02 m.

the pressure peak of the detonation front is outside of the validation range of currently available reaction mechanisms for detailed H2 -O2 reaction kinetics, we slightly decrease the initial pressure and increase the initial temperature as compared to the experiment to achieve a similar reaction behavior. We believe that it is important for further numerical investigations to operate inside the validated range of the reaction mechanism.

4. Results and discussion

The present numerical investigation of RSBI covers different reaction wave types triggered by shock Mach numbers between Ma = 2 . 13 and Ma = 2 . 90. Deflagration is induced by the lowest shock Mach number of Ma1 = 2 . 13. Increasing shock strength leads to three different types of supersonic reaction waves: Ma2 = 2 . 19 induces a deflagration, which transitions to a detonation. Ma3 = 2 .30 and Ma5 = 2 . 90 immediately cause detonations behind the shock wave at the downstream or upstream pole of the bubble, respectively. A shock Mach number of Ma4 = 2 . 50 leads to a nearly simultaneous double detonation in two bubble regions. Temporal and spatial bubble evolution, enstrophy production and mixing are strongly affected by the reaction wave type, which we discuss comprehensively in the following sections. The simulation with a shock Mach number of Ma2 = 2 . 19 is excluded from the discussion, as the global bubble dynamics are nearly identical to

RSBI at Ma3 = 2. 30. The transition process will be discussed in Section 5.1 and the double detonation at Ma4 = 2 .50 will be discussed in Section 5.2.

4.1. Global bubble dynamics

The qualitative influence of the chemical reaction kinetics on the temporal evolution of SBI is shown in Fig. 2. The contour plots of the temperature inside the bubble show the compression and propagation of the reaction front. The upper part of each sequence shows the reacting simulation, the lower part provides results for the non-reacting simulation at the same shock Mach number. For clarity we first compare simulations with shock Mach numbers Ma1 = 2 . 13 . Ma3 = 2 . 30 and Ma5 = 2 . 90.

We refer to the lower part of the first sequence of contour plots for a general description of the characteristic stages of bubble evolution. At t = 50 js and t = 86 js, the shock wave propagates through the bubble and compresses the gas mixture inside. Before the shock wave has passed, first instabilities on the interface start to arise (t = 86 jas). At t = 120 jas, the roll-up of the bubble has started, primary vortices form and secondary instabilities grow due to shear at the material interface. Finally at t = 500 js, the bubble gas shows a high degree of mixing with the surrounding gas N2 . The two main vortices are fully developed and connected over the bridge region at the upstream pole of the bubble. Secondary vortex

26 (is

50 (is

120 us

500 |js

Fig. 2. Temperature contour plots of SBI: upper parts show reacting SBI, lower parts the inert counterparts. Ma] = 2, 13 induces a deflagration wave, Ma3 = 2 , 30 and Ma5 = 2 j 90 a detonation wave.

structures are clearly visible at the outer interface. The general evolution of inert SBI and the characteristic stages are similar. The different shock strengths leads to a range of propagation velocities of the shock waves, which in turn shift the overall bubble evolution in time. Furthermore, a higher shock Mach number causes finer structures in the long-term development of the SBI.

Similarity of evolution at different shock Mach numbers cannot be observed when chemical reaction kinetics are taken into account. As the reactions are strongly pressure sensitive, the shock Mach number affects the induction time and the subsequent reaction process. The first and weakest shock wave (Ma1 = 2 . 13) leads to deflagration. The gas mixture is ignited shortly before t = 120 jas, as shown in the upper row of Fig. 2. The reaction front propagates slowly through the bubble gas; even at t = 500 jas the reaction front has not yet reached the upstream pole. Thus bubble structures of reacting and non-reacting SBI are still similar, especially the outer interfaces with evolving Kelvin-Helmholtz instabilities show the same characteristics. Merely, the overall bubble expansion increases due to isobaric heating over the reaction front.

When we increase the shock Mach number to Ma3 = 2 .30 . the bubble dynamic changes distinctly. The reactive gas mixture ignites earlier, followed by a detonation wave, depicted in the second row of Fig. 2. The supersonic reaction wave propagates within approximately At = 10 jas through the bubble. Strong heat release and density decrease result in a rapid and significant bubble expansion. When the detonation wave reaches the interface of the bubble, vorticity is generated, with the opposite sign compared to the vorticity induced by the initial shock wave. As a consequence, the growth of the secondary instabilities is decelerated, which has a significant effect on mixing. The contour plots at t = 120 jas and t = 500 jas show the bubble after the reaction wave has propagated through the reactive gas mixture. Comparison with the inert counterpart reveals different growth rates and characteristics. Furthermore, the detonation wave amplifies the N2 -jet at the symmetry plane at the downstream pole of the bubble. Hence the bridge region, connecting the two main vortices, vanishes completely.

A further increase of the shock Mach number to Ma5 = 2 . 90 shortens the induction significantly. The strong shock ignites the gas mixture directly at the upstream pole of the bubble, followed by a detonation wave. The detonation wave merges with the initial shock wave inside the bubble and subsequently propagates through the unshocked gas mixture. Comparison with the inert counterpart shows that the detonation wave propagates more than twice as fast as the initial shock wave, which significantly influences the spatial bubble evolution. Similarly to the simulation with Ma3 = 2 .30 . some of the secondary instabilities are suppressed; however, the N2 -jet is less amplified and the bridge region is preserved even in the long-term evolution. The bubble is penetrated by a single detonation wave, whereas the lower shock Mach number of Ma3 = 2 .30 induces a detonation wave that propagates in upstream direction through the pre-shock gas mixture. The interface differs from the other simulations, the KHI are not aligned along the outer interface. Furthermore, the timestep at t = 120 js shows secondary RMI arising from the KHI. Different evolution of primary and secondary instabilities, the bridge region, and the spatial bubble expansion have a significant effect on integral quantities, such as mixing and enstrophy production.

The vortex Reynolds number fter = r0 /v, defined by Glezer [52], where r0 is the initial deposition of vorticity and v the kinematic viscosity of the interface, amounts to 1. 4 x 105 for Ma = 2 . 13 and increases with higher shock Mach numbers up to 2 .4 x 105 (Ma = 2 . 90). The critical Reynolds number for transition from laminar to turbulent flow is 104 to 2 x 104 (Dimotakis [53]). All configurations exceed this mixing-transition Reynolds number by at least one order of magnitude.

Time [s]

Fig. 3. Normalized transverse bubble diameter for different shock Mach numbers. — : reaction; - : no reaction; Q : Ma. = 2 . 13 , A : Ma3 = 2.30 , □ : Ma5 = 2. 90.

4.2. Transverse bubble diameter

The transverse bubble diameter Ay = Ay/D0 normalized by the initial bubble diameter D0 is used to measure the impact of chemical reaction processes on the large-scale evolution of the bubble. The bubble diameter Ay is measured based on a threshold value of the xenon mass fraction of YXe = 0 . 01.

Figure 3 shows the temporal evolution of Ay for the inert and the reacting simulations at three different shock Mach numbers. For the inert simulations (dashed-dotted lines) we observe a nearly linear increase in the bubble diameter. Some variation can be found in the long-term evolution: the weaker the shock strength the smaller the streamwise expansion, which leads to a slightly larger transverse bubble diameter. At the highest shock Mach number of Ma5 = 2 . 90 , the evolution levels out much earlier than for the other shock Mach numbers. Note that the roll-up of the primary vortices leads to a wave-like temporal growth of Ay . In the long-term evolution the main vortices are fully developed and the bubble gas rotates around the vortex cores, which results in a flattening of the transverse bubble diameter evolution.

The deflagration wave triggered at a shock Mach number of Ma1 = 2 .13 affects the normalized transverse bubble diameter only with respect to the long-term evolution. The propagation velocity of the reaction front is low compared to the evolution of the hydrodynamic instabilities. The density increase over the reaction front accompanied by a spatial expansion leads to a slight divergence from the inert counterpart after t = 300 jas. Simulations with shock Mach numbers of Ma3 = 2. 30 and Ma5 = 2. 90 exhibit supersonic detonation waves, which have a significantly stronger effect on the normalized transverse bubble diameter. The rapid reaction leads to an instantaneous expansion of the reacted gas and a sudden increase of Ay up to 175% of the initial bubble diameter. The detonation wave at Ma3 = 2 .30 propagates in upstream and orthogonal direction of the flow field. Thus, an overshoot in the transverse bubble diameter is visible at t = 180 js, which decreases over time to a slightly lower value. The higher shock Mach number of Ma5 = 2. 90 shows an earlier increase of the transverse bubble diameter and levels out at a higher value compared to the simulation with a shock Mach number of Ma 3 = 2. 30. The higher compression leads to a higher temperature and therefore to a larger spatial expansion of the bubble gas.

4.3. Identification of the reaction wave type

To outline the different features of deflagration and detonation waves we analyze the evolution of characteristic parameters across

■10"

| 106 h

-2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1

x-dircctioo [m] -10-2

(a) Mai = 2.13 (deflagration)

■10"

°'5 Jjj û

-2.7 -2.6 -2.5 -2.4 -2.3

x-direction [m] -10-2

(b) Ma3 = 2.30 (detonation)

Fig. 4. Pressure and gas composition across the fully developed reaction wave front. — : pressure; - : YH.O ; — : YHOl .

the reaction front. Figure 4 shows the pressure and the mass fractions of the radical HO2 and the product gas H2 O across the reaction waves for Ma1 = 2 , 13 (deflagration) and Ma3 = 2 , 30 (detonation) at a time and position, when the reaction wave is fully developed and the initial shock wave has passed the bubble. (t = 90 jas for the detonation wave and t = 300 jas for the deflagration wave.)

Figure 4 (a) shows data for the deflagration reaction wave. The reaction front propagates from right to left and its location coincides with the peaks of the HO2 mass fraction. Accompanied by the peak a rapid increase of the product gaseous H2 O is apparent. The pressure is not affected by the chemical reaction and remains constant across the reaction front.

A different evolution is observed for the supersonic reaction wave at Ma3 = 2 ,30 , see Fig. 4(b). In addition to the increase of the product gaseous H2 O and the intermediate species HO2 across the reaction front, also the pressure exhibits a pronounced peak, which is caused by the shock wave preceding the detonation wave. The pressure decreases behind the shock wave but levels out at a larger value compared to the deflagration wave. Moreover, the amplitudes of the HO2 mass fraction peaks in the reaction zone differ. The detonation wave shows a significantly higher amount of HO2 and a breakdown across the reaction zone indicating that the third explosion limit is crossed and high-pressure reactions dominate [54]. Both reaction waves result in the same amount of the product gas H2 O.

Figure 5 shows the mass fraction of H2 O2 for the two different reaction wave types. The spatial coordinate £ denotes the distance from the reaction front, which propagates from right to left. The peak of YH2 Ot indicates a reaction above the third explosion limit.

HO2 collides with H2 forming either H2 O2 or directly H2 O. [54]. The strong reduction of Y H2 O2 on the right side of the plot identifies the bubble interface and consequently the boundary of the reaction zone.

4.4. Damkohler number

The flow field of RSBI is affected by hydrodynamic effects and chemical reaction kinetics. The Damkohler number, defined as the ratio of the hydrodynamic and chemical reaction time scales,

Da = ^ , (35)

indicates which effect dominates. Da > 1 characterizes a flow field mainly driven by chemical reactions, Da < 1 implies a domination of the hydrodynamic effects. The two time scales are defined as follows:

Th = M •

Tr = Tign +

The characteristic hydrodynamic time scale th is defined by the total vorticity W, averaged from the first contact of the shock wave with the bubble until the reaction wave has propagated through the bubble. The chemical reaction time scale tr consists of two time intervals: tign is the period from the first contact of the shock with the bubble until ignition, and D0 /(2VRW) is the time the reaction wave needs to propagate through half of the initial bubble shape with D0 as the initial bubble diameter. The propagation velocity of the deflagration wave is sensitive to temperature and pressure, considered by a power law expression as introduced by Dehoe [55]

V»=a ) ( P0)

Fig. 5. H2 O2 mass fraction across the fully developed reaction front. (deflagration) and — Ma, = 2,30 (detonation).

: Ma. = 2. 13

with SLt denoting the laminar burning velocity at reference conditions (T0 and p0), available in recent literature [55]. Temperature T and pressure p are taken from the hot spot shortly before ignition. The parameters /31 and f) 2 are 1.54 and 0.43, respectively, as we are dealing with a stoichiometric mixture [56]. We also computed the propagation velocity directly from the simulation and found good agreement to the literature. As the varying temperature and pressure distributions inside the bubble lead to a range of different propagation velocities, we decided to use the velocity calculated from Eq. (38). The obtained velocity proved to be a reasonable estimated value for the propagation velocity. For detonation the propagation velocity of the reaction wave is more stable and therefore determined directly from the simulation. Table 1 provides

Table 1

Damkohler numbers and characteristic time scales for different shock Mach numbers.

Ma [-] t h [f ] Tr [S ] Da [ -]

2.13 4. 926 > 10 -4 1. 739 > 10 -3 0.283

2.19 1. 227 > 10 -4 9. 291 > 10 -5 1.321

2.30 1. 094 > 10 -4 8. 278 > 10 -5 1.322

2.50 1 . 286 > 10 -4 7. 499 > 10 -5 1.715

2.90 3. 086 > 10 -4 4. 499 > 10 -5 6.859

the Damköhler numbers for the different initial shock Mach numbers, including the ones for Ma2 = 2 . 19 and Ma4 = 2 . 50 , which are discussed in Section 5.

The deflagration wave induced at a low shock Mach number of Ma1 = 2 . 13 leads to a Damköhler number of Da = 0.283. The hy-drodynamic time scale dominates the flow field, bubble evolution and the growth of secondary instabilities are mainly driven by hy-drodynamic effects. As shown in Fig. 2, the influence of the reaction front is minor.

At higher shock Mach numbers, chemical reactions start to play a crucial role for the overall bubble dynamics. The fast propagation velocity of the detonation wave shortens the chemical reaction time scale tr and leads to a Damköhler number of Da = 1.322 for a shock Mach number of Ma3 = 2 .30 . The primary vortex region and the growth of secondary instabilities are directly affected by the detonation wave. At the highest shock Mach number of Ma5 = 2 . 90 , evolution is entirely dominated by the detonation wave. After it merges with the initial shock wave, it determines the spatial evolution of the gas bubble. RMI, the primary vortices and the secondary instabilities emerge under the influence of the reaction wave. The ignition at the upstream pole of the bubble shortens the reaction time scale. The single reaction wave reduces the vorticity production compared to the detonation wave, which originates at the downstream pole of the bubble and therefore increases the hy-drodynamic time scale. The reaction wave dominates the flow field, which finds expression in a significant increase of the Damköhler number up to Da = 6 .859 . At a lower shock Mach number of Ma3 = 2 .30 . the early RMI evolves in the unburnt gas mixture of the bubble, which leads to a lower Damköhler number compared to the simulation at Ma5 = 2. 90.

4.5. Enstrophy generation We use the enstrophy e = j m 2 dxdy

to determine the influence of the different reaction waves on the vorticity production. Figure 6 outlines the enstrophy over time for reacting and inert SBI. The enstrophy is zero until the shock wave reaches the upstream pole of the bubble. Baroclinic vorticity production leads to an increase during the shock wave passage. A first local maximum in enstrophy is reached after the shock has passed half of the bubble, an effect that can be observed for all simulations. Thereafter, a slight decay is visible, followed by another increase due to shock focusing and shock reflections at the interface. The enstrophy gradually decays after the passage of the shock. The same pattern is observed for all inert simulations, independently of the shock Mach number; only overall enstrophy levels differ as stronger shock waves generate more enstrophy.

The deflagration wave induced by a shock Mach number of Ma1 = 2 . 13 has no noticeable influence on the enstrophy, the variation between the reacting and inert simulations is negligible. The detonation waves induced by a shock Mach number of Ma3 = 2 . 30 produce significant amounts of additional vorticity, which leads to

o ° o

Fig. 6. Enstrophy. — : reaction; 2.30 . □ : Ma5 = 2. 90.

: no reaction; Q : Ma] = 2. 13 , A : Ma3 =

a distinct enstrophy peak, see Fig. 6. The elevated enstrophy levels persist for about 50 |as. The highest shock Mach number of Ma5 = 2 . 90 shows a different behavior. As the ignition spot is located at the upstream pole of the bubble and as the mixture ignites immediately after the first contact of the shock wave, enstro-phy production is dominated by the detonation wave. Therefore, we have two enstrophy peaks of similar magnitude, one during the shock wave passage of the upstream part of the bubble and one during the passage of the downstream part. Thereafter, similar to the simulation at Ma3 = 2 .30 . the enstrophy decays faster compared to their inert counterparts. The inert SBI, shown in the contour plots in Fig. 2, is characterized by several small vortices at the outer interface, even in the long-term evolution. The reacting SBI shows a much smoother flow field with fewer vortices. The detonation waves of reacting SBIs decelerate the growth of secondary instabilities and reduces the appearance of smaller vortices as it induces vorticity with opposite sign compared to the vortic-ity produced by the initial shock wave. Furthermore the increased diffusion across the reaction front damps the growth rate of secondary instabilities. The shock wave Mach number of Ma5 = 2 . 90 reveals an additional effect: enstrophy production during the first half of the shock wave passage for the reacting SBI is higher than for the inert SBI. Parts of the detonation wave are reflected, when it merges with the initial shock wave. The reflected wave produces additional vorticity at the internal interface inside the bubble. During the second half of the shock wave passage, the inert SBI shows a larger enstrophy production. The density gradient across the shock wave is higher than the gradient across the detonation wave, leading to a higher enstrophy production. The vorticity of the reflected shock wave in a reacting SBI, increasing the enstro-phy during the first part of the shock wave passage, has already decayed at this stage of SBI.

4.6. Mixing in RSBI

The shock-bubble interaction provides a complex flow field, where RMI and KHI induce local spots of high mixing rates. Tomkins et al. [57] identified three main regions of mixing: the main vortices, the outer interface including KHI and the bridge region, which connects the two main vortices. The latter contributes up to 40% to the mixing. To estimate the impact of the reaction waves on the mixing, we use the molecular mixing fraction (MMF), defined by Danckwerts [58] as

&(t ) =

fZ Xn.Xxe)dx fZ (Xn2 )(Xxe )dx

r'' ...........................

Time [s]

■10"

Fig. 7. Molecular mixing fraction. — : reaction; — : no reaction; O : Ma. = 2 . 13,

A : Ma3 = 2 . 30 . □ : Ma, = 2 . 90.

The molar mixing fraction can be interpreted as the ratio of molecular mixing to large-scale entrainment by convective motion. We plot the temporal evolution of ©(f) for reacting and inert simulations in Fig. 7.

The inert simulations show a linear growth of the molar mixing fraction, the slope increases with higher shock Mach number. An increase in shock strength leads to higher enstrophy production and faster growth of secondary instabilities, which enhance mixing. The reacting counterparts show a different behavior. In general, the mixing is reduced by the reaction waves, independent of their type. The deflagration wave induced by the lowest shock Mach number leads to a decrease of up to 30%. Mixing is affected

after approximately f = 270 |as, when the reaction wave reaches parts of the interface and the main vortices. However, the bridge region remains unaffected, as the propagation velocity of the deflagration wave is too low.

The detonation wave affects all three main mixing regions. The MMF is reduced by up to 50% for Ma3 = 2 .30 as well as for Ma5 = 2 . 90. Besides the reduction of mixing in the region of the main vortices and at the interface, which are already affected by a deflagration wave at a lower shock Mach number, the bridge region is also influenced by the detonation wave. The detonation waves decelerate the growth of secondary instabilities. Especially the bubble evolution at the highest shock Mach number, see Fig. 2, shows damping of fine structures, which explains the higher reduction of the MMF.

Figure 8 outlines the mixing progress in the long-term evolution at f = 500 |as for three different shock Mach numbers. The inert SBIs show already a high degree of mixing, whereas the reacting SBIs are characterized by large areas of unmixed bubble gas. The inert SBIs show higher mixing for increasing shock Mach number, which is in accordance to the MMF plotted in Fig. 7. The reacting SBIs follow a similar trend: For low shock Mach numbers of Ma = 2 . 13 (Fig. 8(a)) and Ma = 2 .30 (Fig. 8(b)) we observe larger areas of unmixed bubble gas, whereas the highest shock Mach number of Ma = 2 . 90 (Fig. 8(c)) shows a higher degree of mixing.

5. Special cases

Two simulations of RSBI contain hydrodynamic and chemical features that have to be discussed in detail. We observed DDT for Ma2 = 2 . 19 and a double detonation with different reaction branches for Ma4 = 2 . 50 .

Fig. 8. Mixing of N. and Xe for inert and reacting SBI at t = 500 |is for different shock Mach numbers. Upper parts show the reacting SBI, lower parts the inert counterpart. Ma2 = 2.19:

P [Pa] : I

2.0E+05

2.0E+06

4.0E+06

Fig. 9. Contour plots of RSBI at Ma2 = 2 . 19 during deflagration-to-detonation transition: upper parts shows the temperature, lower parts the pressure.

Fig. 10. Detailed contour plot of RSBI at Ma. = 2. 19 shortly after DDT: upper part shows the temperature, lower part the pressure. The arrow shows the line slice for the data of Fig. 11.

5.1. Deflagration-to-detonation transition (Ma2 = 2. 19)

Figure 9 shows contour plots for Ma 2 = 2 . 19 during the early stage of ignition. Temperature (upper part of the plot) and pressure (lower part) are parameters that illustrate the ignition and transition process. At t = 93. 0 jas the shock wave still propagates through the bubble gas, and instabilities at the interface start to grow. The gas mixture ignites near the downstream pole of the bubble at approximately t = 93 . 2 js. A subsonic deflagration wave propagates through the bubble gas until t = 94 . 9 js. At t = 94 . 9 js the transition into a detonation wave starts in the lower region of the reaction front. The remaining contour plots shows a fast growth of the supersonic detonation wave out of the deflagration front and a characteristic steep pressure peak across the detonation reaction wave.

Detailed contour plots to illustrate how the detonation wave evolves from the deflagration wave are outlined in Fig. 10. The upper part of the figure shows the temperature and the lower part the pressure. A temperature isoline (T = 200 0 K) visualizes the reaction front.

The temporal evolution of the characteristic thermodynamic properties during the transition is plotted in Fig. 11. The plots show the variations of pressure, the radical H concentration (Fig. 11(a)) and the temperature (Fig. 11(b)) for seven timesteps from t = 94 . 6 js until t = 96 . 0 js. The coordinate £ is obtained by a rotation of the original coordinate system. As a result of the trans-

Fig. 12. Reaction front marked by black temperature isocontour, shortly before transition to detonation. High-pressure region denoted by red isocontour (p > 5 bar).

formation £ coincides with the propagation direction of the reaction wave. The reaction wave propagates from the right to the left. We observed an increase of the pressure peak accompanied by a decrease of H in the reaction zone, which is characteristic for DDT. The H radical appears within the deflagration wave, following the detonation wave. The peak of H can be found behind the shock wave, after the sudden increase in temperature and pressure. Furthermore, DDT can be identified by the steepening of the temperature profile, its peak at the reaction front and the higher temperature of the product gas. These findings are consistent with observations of Ivanov et al. [36] and Liberman et al. [15,18].

The deflagration front propagates in semicircular direction through the bubble gas, hence the question arises why the transition to detonation occurs at specific areas of the reaction front. The reason can be found by the detailed analysis of the pressure distribution in front of the reaction front, where DDT occurs. Figure 12 shows the reaction front shortly before the transition to detonation. The reaction zone is indicated by two black isocontours, the red isocontour depicts the region of high pressure. Due to shock focusing of the initial shock wave at the downstream pole of the bubble, a region of high pressure ( p > 5 bar) exists. The reaction front propagates into this region and high-pressure reactions are promoted, which support the transition to detonation.

As the deflagration wave persists only for a few microseconds, the overall bubble evolution of the RSBI at Ma2 = 2 . 19 is similar to that at a shock wave Mach number of Ma3 = 2. 30. The normalized transverse bubble diameter, the enstrophy production and the molar mixing fraction are nearly identical. The Damkohler numbers

(a) — : Pressure; — : radical H concen- (b) — : Temperature; — : radical H con-

tration

centration

Fig. 11. Pressure, temperature and radical H concentration in the reaction front during the transition to detonation for seven conservative time steps with At = 0 . 2 js starting at t = 94 . 6 js.

Ma* = 2.50:

t= 74 (is 80 (is 90 (is 120 ns 500ns

Fig. 13. Contour plots of RSBI at Ma4 = 2 . 50 : upper parts shows the reacting SBI, lower parts the inert SBI.

amount to Da = 1. 321 (Ma 2 = 2 . 19) and Da = 1.322 (Ma3 = 2 . 30) indicating a chemically dominated flow field. To ensure that the observed DDT is not caused by numerical artifacts, the RSBI with a shock Mach number of Ma2 = 2 . 19 was repeated at a coarsened and refined grid resolution. We used the same grids that were already applied for the grid convergence study in our previous paper [29] with cell sizes of Axy = 234, 117 and 59 pm. All simulations show the same induction time, ignition spot and location of DDT.

5.2. Double detonation (Ma4 = 2.50)

The shock wave of Ma4 = 2 50 induces two detonation waves. One originates near the downstream pole of the bubble and one at the upstream pole. Figure 13 shows the ignition spots, the propagation and interaction of the detonation waves. The two reaction waves propagate towards each other, which leads to a rapid consumption of the reactive bubble gas. The reaction time scale is significantly shortened, leading to an increase of the Damköhler number to Da = 1. 715 . see Table 1. Specific conditions have to be satisfied to cause a double detonation: The ignition delay time in the first reaction spot at the upstream pole of the bubble has to be longer than in the second spot near the downstream pole. However, the ignition has to be nearly simultaneous in the absolute timeframe.

The reaction region at the upstream pole is characterized by a slow increase of the intermediate gas mass fractions directly behind the shock wave, beginning at t = 28 ps. After an induction time of about 47 ps the gas mixture ignites and forms a detonation wave. The second ignition spot shows a different behavior; a strong reaction is directly induced after the shock wave has passed at t = 68 ps. The induction time only amounts to t = 6 ps, which is much shorter than for the first ignition spot. The higher compression and temperature at the downstream pole of the bubble lead to a faster ignition. In the absolute timeframe both spots trigger ignition, followed by a detonation wave at approximately the same time, which leads to a double detonation.

Figure 14 outlines the temporal evolution of the intermediate species H . O . OH during the induction time in the two reaction zones. The solid lines denote the ignition spot at the upstream pole of the bubble, the dashed lines show data for the downstream pole. As the chemical reactions are highly pressure sensitive the higher temperature and pressure at the downstream pole leads to a faster production of the radicals and to a shorter ignition delay time, compared to the upstream pole.

6. Discussion

6.1. Transverse bubble diameter

Our setup is based on the experimental investigation of Haehn et al. [28]. They observed, similarly to our work, different reaction wave types by varying the shock Mach number between Ma = 1 34

10° 10"2 10-1 10"6

S lo-«

\ !

5 6 Time [s]

■10"5

Fig. 14. Maximum mass fraction of intermediate products inside the two ignition spots. O = H. O = O. A = OH. — : 1st ignition spot at the upstream pole; — : 2nd ignition spot at the downstream pole.

and Ma = 2 . 83. In the following section, we compare the results of the two studies. We are well aware that a comparison of two-dimensional simulations with the experiment of Haehn et al. [28] is only reasonable along early stages of evolution [29]. Haehn et al. [28] provide the Damköhler number, the transverse normalized bubble and the main vortex diameter for several experimental setups. We will use the transverse bubble diameter to compare the general evolution of the bubble and the influence of the reaction waves. The experiments exhibit deflagration for shock Mach numbers of Ma = 1 63 and Ma = 2 07 and detonation for a shock Mach number of Ma = 2 . 83. Haehn et al. [40] also observed detonations for shock Mach numbers smaller than Ma = 2 . 83 . however, without providing quantitative data. We compare the experimental results for Maexp = 2 . 07 and Maexp = 2 . 83 with our results obtained at Manum = 2 . 13 and Manum = 2 . 90. Figure 15 shows the evolution of the normalized transverse bubble diameter for a RSBI with either a deflagration wave (a) or a detonation wave (b). The normalized time t* follows the definition of Haehn et al. [40]

t f = T., (41)

where t is the time measured from the first contact of the shock wave with the bubble. tn is defined as D0 /W. , with W. as the incident shock wave speed.

As discussed in Section 4.2, the slowly propagating deflagration wave leads only to a slight increase of the bubble diameter. Inert and reacting SBI show a similar evolution, confirmed by the experimental as well as the numerical results in Fig. 15(a) . The numerical data show a larger normalized bubble diameter in the long-term evolution, which is attributed to two-dimensional effects.

Also the data for SBIs, which induce a detonation wave, are in very good agreement. The propagation velocity of the reaction wave, the spatial expansion of the bubble and the peak of the

Fig. 15. Normalized transverse bubble diameter. Comparison between simulations and experimental data of Haehn et al. [28]. Reaction: — and No reaction: — and o. (a) Deflagration Ma = 2. 13/ 2 .07 [28], (b) Detonation Ma = 2 .90/2.83 [28].

transverse bubble diameter are nearly identical. Figure 15(b) reveals a good match of the steep bubble diameter increase, indicating that the propagation of the reaction wave inside the gas bubble is well reproduced by the two-dimensional simulations. Similarly to the deflagration setup, the bubble expansion differs from the experimental results in the long-term evolution.

6.2. Double detonation

As mentioned before, Haehn et al. [28] also performed experiments at intermediate shock strengths. For one specific setup, their chemiluminescence signal showed two bright spots, indicating two ignition points. They provided two explanations. The first one suggests that the second bright spot may be a reflection of the first initial combustion signal from the downstream surface of the un-shocked bubble interface. The second explanation assumes that simultaneous detonations at two specific points in the compressed bubble are possible as the induction time of the second reaction spot is shortened by the higher compression of the shock focusing. A proper combination of induction times thus can lead to simultaneous ignitions. Our numerical results for Ma4 = 2 . 50 support this latter hypothesis.

6.3. Limitations and critical discussion

The comparison of our two-dimensional simulations with the experimental data of Haehn et al. [28] has some limitations with respect to the bubble evolution.

As shown in Fig. 15, the transverse expansion of the bubble gas deviates from the experimental results in the long-term evolution. Recent studies of Wang et al. [59] observe the same with respect to the transverse bubble diameter of two- and three-dimensional SBIs in their experiments. The vortex stretching term of the vorticity equation vanishes for a two-dimensional simulation. The missing of this term affects the expansion and increases the transverse expansion in the long-term evolution [60], which explains the deviation of the transverse bubble diameter for t* > 6. Hejazialhosseini et al. [60] investigated the vortex dynamics in three-dimensional inert SBI and showed that the growth rate of the vortex stretching term increases significantly in the long-term evolution and contributes to a decrease of the transverse bubble diameter, an effect missing in two-dimensional simulations.

Furthermore, the shock-focusing in three dimensions is stronger, which shortens the ignition delay time. To compensate for this effect, we simulate RSBI at a slightly higher shock Mach number to achieve the same ignition delay time. The validity of two-dimensional simulations containing shock-induced instabilities has also been shown by Peng et al. [61] in their study of vortex-accelerated secondary baroclinic vorticity deposition. Klein et al. [62] investigated the interaction between a sphere and a shock wave at high shock Mach numbers. They compared the two-dimensional results with experimental data and observed good accordance in the radial and axial width of the shocked sphere. Both studies achieved very good agreement between two-dimensional RMI and experiments, even in the long-time dynamics.

Nevertheless some effects are not resolved by our simulation such as the onset of turbulence. Three-dimensional vortex rings tend to become unstable and vortex stretching may eventually result in broad-band turbulence [63]. This production mechanism is suppressed in a two-dimensional simulation. Three-dimensional effects cannot be neglected, however they become important only in the long-term evolution. Niederhaus et al. [64] studied SBIs and showed that three-dimensional effects affect the total enstrophy only at late times. Miles et al. [65] support this assumption, as they also observed no significant differences of the early growth rates of shock-induced instabilities between three- and two-dimensional simulations. Further studies report that vortex stretching affects only the long-term evolution of the mixing rate [60]. These observations support the integrity of our results, as the chemical reaction and its interaction with the hydrodynamic instabilities occur in the early stage of SBI. The very good agreement between our numerical results and the experimental data of Haehn et al. [28] indicate that three-dimensional effects may not be very significant for the specific phenomena considered here. Hence, we are confident to provide valid and reliable results for the investigation of reaction wave characteristics and its influence on the global bubble dynamics. In particular the good agreement of the normalized transverse bubble diameter for the detonating RSBI with experimental data indicates that essential mechanisms are reproduced. Furthermore, the detection and analysis of a double detonation in the simulations supports the experimental observations of Haehn et al. [28].

7. Conclusion

We have presented results of a reacting shock-bubble interaction at different shock Mach numbers with detailed H2 -O2 chemical reaction kinetics. A gas bubble filled with a reactive stoichio-metric gas mixture of H2 , O2 and Xe is penetrated by a shock wave with Mach numbers between Ma = 2 . 13 and Ma = 2 . 90. The planar shock wave propagates through the domain and interacts with the cylindrical density inhomogeneity, inducing Richtmyer-Meshkov instabilities. The convergent shape of the bubble focuses the shock, which triggers ignition of the bubble gas. Depending on the shock Mach number, the pressure sensitive H2 -O2 gas mixture shows different induction times, ignition spots and reaction wave types, which strongly affect the spatial bubble evolution and the mixing process. A weak shock wave induces a deflagration wave, higher shock Mach numbers drive high-pressure reactions, resulting in a detonation wave.

We showed that the variation of the shock Mach number covers several reaction wave types with different impact on the mixing process. A deflagration wave has a minor influence on the global bubble evolution and leads to a flow field dominated by hydro-dynamic effects (Da ^ 0.28). The growth of secondary instabilities is partially affected, which decreases mixing by about 30%. Higher shock Mach numbers and the subsequent detonation waves lead to a chemically driven flow field (Da > 1.32). The supersonic reaction

wave propagates rapidly through the reactive gas mixture and affects all bubble regions. Secondary instabilities are suppressed, the bridge region is disturbed and mixing is reduced by up to 50%. The detonation wave induces an additional peak in the enstrophy production, followed by a faster decay.

A shock Mach number Ma = 2 . 50 reveals a particular phenomenon. The bubble gas ignites simultaneously at two spots, leading to two detonation waves that propagate towards each other. At an intermediate shock Mach number of Ma = 2 . 19 . we observe a deflagration wave that transitions into a detonation wave shortly after ignition. Comparison with experimental data shows a very good agreement in terms of spatial expansion, reaction wave types and propagation velocities.

Acknowledgment

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (http://www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, http://www.lrz. de).

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 667483).

Supplementary material

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.combustflame.2016.09.014

References

[1] J. Yang. T. Kubota. E. Zukoski. Applications of shock-induced mixing to supersonic combustion, AIAA J. 31 (1993) 854-862.

[2] F. Marble. E. Zukoski. J. Jacobs. G. Hendricks, I. Waitz. Shock enhancement and control of hypersonic mixing and combustion, AIAA 26th Joint Propulsion Conference (1990).

[3] R.D. Richtmyer. Taylor instability in shock acceleration of compressible fluids, Commun. Pure Appl. Math. 13 (1960) 297-319.

[4] E.E. Meshkov. Instability of the interface of two gases accelerated by a shock wave, Fluid Dyn. 4 (1969) 101-104.

[5] L. Rayleigh. Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density, Proc. Lond. Math. Soc. 14 (1883) 170-177.

[6] G. Taylor. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. Part 1. Waves on fluid sheets, Proc. R. Soc. Lond. A Math., Phys. Sci. 201 (1950) 192-196.

[7] M. Brouillette. The Richtmyer-Meshkov Instability, Annu. Rev. Fluid Mech. 34 (2002) 445-468.

[8] N.J. Zabusky, Vortex paradigm for accelerated inhomogeneous flows: visiomet-rics for the Rayleigh-Taylor and Richtmyer-Meshkov environments, Annu. Rev. Fluid Mech. 31 (1999) 495-536 .

[9] W.D. Arnett. The role of mixing in astrophysics, Ap. J. Suppl. 127 (20 00) 213-217.

[10] A.M. Khokhlov. E.S. Oran. G.O. Thomas. Numerical simulation of deflagration-to-detonation transition: the role of shock-flame interactions in turbulent flames, Combust. Flame 117 (1999) 323-339.

[11] R.S. Craxton. K.S. Anderson. T.R. Boehly. V.N. Goncharov. D.R. Harding. J.P. Knauer. R.L. McCrory. P.W. McKenty. D.D. Meyerhofer. J.F. Myatt. A.J. Schmitt, J.D. Sethian. R.W. Short. S. Skupsky, W. Theobald. W.L. Kruer. Direct-drive inertial confinement fusion: a review, Phys. Plasmas 22 (2015).

[12] P.G. Drazin. Introduction to hydrodynamic stability, Cambridge University Press, 2002 .

[13] V.K. Tritschler. S. Hickel. X.Y. Hu. N.A. Adams. On the Kolmogorov inertial subrange developing from Richtmyer-Meshkov instability, Phys. Fluids 25 (2013) 071701.

[14] E. Oran. J. Boris. Numerical simulation of reactive flow, Cambridge University Press, 2005 .

[15] M.A. Liberman. Introduction to physics and chemistry of combustion: explosion, flame, detonation, Springer, 2008.

[16] W. Fickett. W.C. Davis. Detonation: theory and experiment, Dover Publications, 2010 .

[17] J.H.S. Lee, The detonation phenomenon, Cambridge University Press, 2008.

[18] M.A. Liberman. M.F. Ivanov. A.D. Kiverin. M.S. Kuznetsov. A.A. Chukalovsky. T.V. Rakhimova, Deflagration-to-detonation transition in highly reactive combustible mixtures, Acta Astronaut. 67 (2010) 688-701.

[19] E. Oran, V. Gamezo, Origins of the deflagration-to-detonation transition in gas-phase combustion, Combust. Flame 148 (2007) 4-47.

[20] J.-F. Haas. B. Sturtevant, Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities, J. Fluid Mech. 181 (1987 ) 41-76.

[21] J.J. Quirk. S. Karni , On the dynamics of a shock-bubble interaction, J. Fluid Mech. 318 (1996) 129-163.

[22] D. Ranjan. J. Oakley. R. Bonazza. Shock-bubble interactions, Annu. Rev. Fluid Mech. 43 (2011) 117-140.

[23] G. Billet. V. Giovangigli. G. de Gassowski. Impact of volume viscosity on a shock-hydrogen-bubble interaction, Comb. Theory Model. 12 (2008) 221-248.

[24] N. Attal. P. Ramaprabhu. J. Hossain, V. Karkhanis. M. Uddin. J. Gord, S. Roy. Development and validation of a chemical reaction solver coupled to the flash code for combustion applications, Comput. Fluids 107 (2015) 59-76.

[25] N. Attal, P. Ramaprabhu , Numerical investigation of a single-mode chemically reacting Richtmyer-Meshkov instability, Shock Waves 25 (2015) 307-328.

[26] A.M. Khokhlov. E.S. Oran. A.Y. Chtchelkanova. J.C. Wheeler. Interaction of a shock with a sinusoidally perturbed flame, Combust. Flame 117 (1999) 99-116.

[27] L. Massa. P. Jha. Linear analysis of the Richtmyer-Meshkov instability in shock-flame interactions, Phys. Fluids 24 (2012) 056101.

[28] N. Haehn, D. Ranjan, C. Weber, J. Oakley, D. Rothamer , R. Bonazza, Reacting shock bubble interaction, Combust. Flame 159 (2012) 1339-1350.

[29] F. Diegelmann, V. Tritschler , S. Hickel, N. Adams, On the pressure dependence of ignition and mixing in two-dimensional reactive shock-bubble interaction, Combust. Flame 163 (2016) 414-426.

[30] J.G. Ogilvie, The hydrogen-oxygen second explosion limit. A physical chemistry experiment, J. Chem. Educ. 48 (1971) 342-344.

[31] B. Lewis, G. von Elbe, Combustion, flames and explosions of gases, Academic Press, Inc., 1987.

[32] P.D. Neufeld. A.R. Janzen. R.A. Aziz. Empirical equations to calculate 16 of the transport collision integrals omega(l,s)* for the Lennard-Jones (12-6) potential, J. Chem. Phys. 57 (1972) 1100-1102.

[33] B.E. Poling. J.M. Prausnitz. J.P. O'Connell. The properties of gases and liquids, McGraw-Hill, 2001.

[34] A.W. Cook, Enthalpy diffusion in multicomponent flows, Phys. Fluids 21 (2009) 055109.

[35] M. Kuznetsov, M.A. Liberman, I. Matsukov, Experimental study of the preheat zone formation and deflagration to detonation transition, Combust. Sci. Tech-nol. 182 (2010) 1628-1644.

[36] M.F. Ivanov, A.D. Kiverin, M.A. Liberman, Hydrogen-oxygen flame acceleration and transition to detonation in channels with no-slip walls for a detailed chemical reaction model, Phys. Rev. E 83 (2011) 056313.

[37] T. Poinsot, D. Veynante, Theoretical and numerical combustion, R T Edwards, 2001.

[38] J. Troe, Predictive possibilities of unimolecular rate theory, J. Phys. Chem. 83 (1979) 114-126.

[39] M. O Conaire. H.J. Curran. J.M. Simmie , W.J. Pitz. C.K. Westbrook. A comprehensive modeling study of hydrogen oxidation, Int. J. Chem. Kinet. 36 (2004) 603-622.

[40] N.S. Haehn, Experimental investigation of the reactive shock-bubble interaction, University of Wisconsin-Madison, 2012 Ph.D. thesis.

[41] P.J.M. Ferrer. R. Buttay. G. Lehnasch. A. Mura. A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers, Comput. Fluids 89 (2014) 88-110.

[42] O.P. Korobeinichev. T.A. Bol'shova. Applicability of Zel'dovich's theory of chain propagation of flames to combustion of hydrogen-oxygen mixtures, Combust., Explos., Shock Waves 45 (2009) 507-510.

[43] G. Strang. On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506-517 .

[44] P.L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43 (1981) 357-372.

[45] B. Larouturou. L. Fezoui. On the equations of multi-component perfect or real gas inviscid flow, Lecture notes in mathematics, 1402, Springer Berlin Heidelberg, 1989, pp. 69-98.

[46] X.Y. Hu, V.K. Tritschler, S. Pirozzoli, N.A. Adams, Dispersion-dissipation condition for finite difference schemes, arXiv:1204.5088 (2012).

[47] S. Gottlieb. C.-W. Shu. Total variation diminishing Runge-Kutta schemes, Math. Comput. 67 (1998) 73-85.

[48] V.K. Tritschler. A. Avdonin. S. Hickel. X.Y. Hu. N.A. Adams, Quantification of initial-data uncertainty on a shock-accelerated gas cylinder, Phys. Fluids 26 (2014a) 026101.

[49] V.K. Tritschler, B.J. Olson, S.K. Lele, S. Hickel, X.Y. Hu, N.A. Adams, On the Richt-myer-Meshkov instability evolving from a deterministic multimode planar interface, J. Fluid Mech. 755 (2014b) 429-462 .

[50] V.K. Tritschler, M. Zubel, S. Hickel, N.A. Adams, Evolution of length scales and statistics of Richtmyer-Meshkov instability from direct numerical simulations, Phys. Rev. E 90 (2014c) 063001.

[51] P.N. Brown. G.D. Byrne. A.C. Hindmarsh. VODE: a variable-coefficient ODE solver, SIAM J. Sci. Stat. Comput. 10 (1989) 1038-1051.

[52] A. Glezer. The formation of vortex rings, Phys. Fluids 31 (1988) 3532-3542.

[53] P.E. Dimotakis. The mixing transition in turbulent flows, J. Fluid Mech. 409 (20 0 0) 69-98.

[54] O. Oldenberg. H.S. Sommers. The thermal reaction between hydrogen and oxygen III. The temperature coefficient of the steady thermal reaction, J. Chem. Phys. 9 (1941) 432-438.

[55] A.E. Dahoe. Laminar burning velocities of hydrogen-air mixtures from closed vessel gas explosions, J. Loss Prev. Process Ind. 18 (2005) 152-166.

[56] T. Iijima. T. Takeno. Effects of temperature and pressure on burning velocity, Combust. Flame 65 (1986) 35-43.

[57] C. Tomkins, S. Kumar. G. Orlicz. K. Prestridge. An experimental investigation of mixing mechanisms in shock-accelerated flow, J. Fluid Mech. 611 (2008) 131-150.

[58] P.V. Danckwerts. The definition and measurement of some characteristics of mixtures, Ap. Sci. Res., Sect. A 3 (1952) 279-296 .

[59] X. Wang. T. Si. X. Luo. J. Yang, Generation of air/SF6 interface with minimum surface feature by soap film technique, 29th International Symposium on Shock Waves 2 (2015), pp. 1065-1070.

[60] B. Hejazialhosseini, D. Rossinelli, P. Koumoutsakos, Vortex dynamics in 3D shock-bubble interaction, Phys. Fluids 25 (2013) 110816.

[61] G. Peng. N.J. Zabusky. S. Zhang. Vortex-accelerated secondary baroclinic vor-ticity deposition and late-intermediate time dynamics of a two-dimensional Richtmyer-Meshkov interface, Phys. Fluids 15 (2003) 3730-3744.

[62] R. Klein, K.S. Budil, T.S. Perry , D.R. Bach, The interaction of supernova remnants with interstellar clouds: experiments on the nova laser, Astrophys. J. 583

(2003) 245-259.

[63] P.A. Davidson, Turbulence: an introduction for scientists and engineers, Oxford University Press, 2004.

[64] J.H. Niederhaus. J.A. Greenough. J.G. Oakley. D. Ranjan. M.H. Anderson. R. Bonazza, A computational parameter study for the three-dimensional shock-bubble interaction, J. Fluid Mech. 594 (2008) 85-124.

[65] A.R. Miles. B. Blue. M.J. Edwards. J.A. Greenough. J.F. Hansen, H.F. Robey. R.P. Drake, C. Kuranz, D.R. Leibrandt, Transition to turbulence and effect of initial conditions on three-dimensional compressible mixing in planar blast-wave-driven systems, Phys. Plasmas 12 (2005) 056317.