Numerical Investigation of Remote Ignition in Shock Tubes

Highly resolved two- and three-dimensional computational fluid dynamics (CFD) simulations are presented for shock-tube experiments containing hydrogen/oxygen (H 2 /O 2 ) mixtures, to investigate mechanisms leading to remote ignition. The results of the reactive cases are compared against experimental results from Meyer and Oppenheim (Proc Com-bust and stfla me.2013.03.026). The results of the non-reactive case are compared against shock tube experiments, recently carried out in Duisburg and Texas. The computational domain covers the end-wall region of the shock tube and applies high order numerics featuring an all-speed approximate Riemann scheme, combined with a 5th order interpolation scheme. Direct chemistry is employed using detailed reaction mechanisms with 11 species and up to 40 reactions, on a grid with up to 2.2 billion cells. Additional two-dimensional simulations are performed for non-reactive conditions to validate the treatment of boundary-layer effects at the inlet of the computational domain. The computational domain covers a region at the end part of the shock tube. The ignition process is analyzed by fields of localized, expected ignition times. Instantaneous fields of temperature, pressure, entropy, and dissipation rate are presented to explain the flow dynamics, specifically in the case of a bifurcated reflected shock. In all cases regions with locally increased temperatures were observed, reducing the local ignition-delay time in areas away from the end wall significantly, thus compensating for the late compression by the reflected shock and therefore leading for first ignition at a remote location, i.e., away from the end wall where the ignition would occur under ideal conditions. In cases without a bifurcated reflected shock, the temperature increase results from shock attenuation. In cases with a bifurcated reflected shock, the formation of a second normal shock and shear near the slip line is found to be crucial for the remote ignition to take place. Overall, the two- and three-dimensional simulations were found to qualitatively explain the occurrence of remote ignition and to be quantitatively correct, implying that they include the correct physics.


Shock Tubes
Shock tubes have been an important tool for many years to investigate fast reaction kinetics and to support the development of reaction mechanisms. A membrane initially separates a pressurized inert gas in the driver section (denoted as region 4) at an elevated pressure from the test mixture in the driven section (denoted as region 1). The experiment starts in the instant the membrane bursts, after which a shock (incident shock) quickly evolves due to the pressure difference. A contact discontinuity forms at the place where the driven and driver gases touch, creating a zone of compressed driven gas between contact discontinuity and incident shock (referred to as region 2) and a zone of expanded driver gas between contact discontinuity and expansion fan (region 3). Eventually, the incident shock wave will reach the end wall of the shock tube having accelerated all driven gas towards the end wall. The shock is then reflected (reflected shock), running back towards the driver gas, and bringing the driven gas to rest, thus increasing pressure and temperature (designated as region 5). The initial conditions are usually chosen such that the temperature T 5 behind the reflected shock exceeds 800 K, initiating chemical reactions, so that an ignition-delay time ig can be measured (Gaydon and Hurle 1963). It must, however, be ensured that the ignition-delay time is sufficiently short compared to the possible test time that is limited by the arrival of the contact discontinuity. At high temperatures and accordingly short ignitiondelay times ( ig < 1 ms), a nearly homogeneous thermodynamic state establishes (until the mixture ignites) behind the reflected shock and the problem can be modeled as inviscid and adiabatic (Bhaskaran and Roth 2002). These well defined initial boundary conditions allow to compare the measurements to those of a perfectly mixed reactor at constant volume.
At longer test times, several phenomena, most of them directly related to the formation of a boundary layer initiated by the motion of the incident shock, can severely impact the results and lead to huge deviations of measured ignition-delay time and the expected ignition-delay time at constant volume. The effect of the boundary layer is well known: It decelerates the (near-wall) flow behind the incoming shock and thus "removes" mass from the core flow and therefore affects the state in the core flow outside of the boundary layer. The change of state in the core-flow has been modeled by Mirels (1957), Mirels and Braun (1957) and Mirels and Mullen (1964) using perturbation theory and has been successfully applied to compute the attenuation of the incident shock or the change of state at a given location in time. According to perturbation theory, the pressure can be approximated by the superposition of the ideal pressure from an adiabatic, inviscid process and weak pressure perturbation waves. This leads to a spatially and temporally changing distribution of state quantities behind the incident shock (non-uniformities) and also affects the change of state behind the reflected shock, as the pressure variations are amplified across the reflected shock (Rudinger 1961). Assuming a purely laminar or purely turbulent boundary layer, Mirels found that the state quantities are constantly increasing between the incident shock and the contact discontinuity. While pressure and temperature fall short compared to the ideal values, the particle velocity behind the shock front decreases due to the reduced shock strength, but accelerates near the contact discontinuity. Hence, the distance between the incident shock wave and the contact discontinuity decreases and the maximum test time is reduced. In some cases, e.g., at very low pressures, the contact discontinuity can even reach the speed of the incident shock (Mirels 1966). Typically, the boundary-layer induced variations of state lead to a slow, continuous rise of pressure ( p 5 ∕ t ) at the end wall.
Another boundary-layer effect is due to the interaction with the reflected shock under conditions that promote the formation of a shock-bifurcation structure. Mark (1958) developed a simple model to predict the occurrence of the shock bifurcation and to describe the geometry and size of this structure. A brief summary is given below.
For his analysis, Mark examined the Mach number of the fluid in the boundary layer M bl and in the main flow M 2 in reflected-shock coordinates. Mark made the assumptions that the fluid in the boundary layer is in thermal equilibrium with the wall and has no velocity relative to laboratory coordinates. He found that the boundary layer Mach number M bl (as one would expect due to the velocity deficit) is below the Mach number of the main flow M 2 at moderate Mach numbers M 1 of the incident shock, as presented in Fig. 1. However, if the Mach number M 1 further increases, a critical point M * 1 is reached after which the boundary-layer Mach number M bl exceeds the Mach number M 2 of the main flow. This phenomenon is due to the cooling effect of the wall. It is also worth noting that the ratio of specific heats has a very large influence. The point at which the Mach number of the boundary layer exceeds that of the main stream is reached earlier for large values of (e.g., M * 1 = 7.4 for = 1.4, M * 1 = 3.8 for = 1.67 ). In addition, the stagnation pressure of the boundary layer p bl,st and the static pressure of the main flow p 5 behind the reflected shock were investigated and compared. At an incident Mach number M 1 close to 1, the stagnation pressure p bl,st is higher than the static pressure p 5 of the main flow and the boundary layer fluid can pass the reflected shock into the end-wall region. At higher incident shock Mach numbers, a lower crossover point M 1,L is reached, where the static pressure of the main flow p 5 exceeds the stagnation pressure p bl,st instead. The interaction of shock and boundary layer is expected to be fundamentally different at Mach numbers M 1 larger than the lower cross-over point M 1,L . If the stagnation pressure in the boundary layer is below the static pressure in the main flow, the boundary layer fluid cannot match the static pressure Fig. 1 Sketch of a bifurcation structure, showing the speed of the reflected shock V RS in laboratory coordinates, the particle speed u 2 in region 2 in laboratory coordinates, the speed of sound a 2 in region 2, and the speed of sound a bl in the boundary layer. M 2 and M bl are Mach numbers in reflected shock coordinates. Solid lines indicate locations with high density gradients. Dashed lines show particle paths. The sketch shows the bifurcation structure under flow conditions such that the Mach number of the boundary layer is less than 1 and thus no shock occurs in the boundary layer even after stagnating. Instead, as it works against the static pressure of the main flow, it is decelerated until stagnation and reversely accelerated in the direction of the reflected shock, causing the fluid to accumulate. This phenomenon is referred to as bifurcation of the reflected shock. At even higher Mach numbers of the incident shock, the Mach number of the boundary layer surpasses that of the main flow at the Mach number M * 1 . As a result, an upper cross-over point M 1,U exists, after which the boundary-layer fluid can pass the reflected shock again and no bifurcation emerges. It is important to note that the two crossover points are closer together in the case of large ratios of specific heats, while the ratio p bl,st ∕p 5 is higher at the same time (e.g., M 1,L = 1.33, M 1,U = 6.45 and p bl,st ∕p 5 ≥ 0.5 for = 1.4, M 1,L = 1.57, M 1,U = 2.8 and p bl,st ∕p 5 ≥ 0.9 for = 1.67 ). Furthermore, the assumptions that are made to calculate the limits of this phenomenon are very conservative. In fact, shock-tube experiments containing a gas or a mixture with large values of (e.g., Argon) are very unlikely to suffer from reflected shock bifurcation (Davidson and Hanson 2004).
A bifurcation structure is characterized by a triple point, where the oblique shock, the normal reflected shock and the tail shock meet. A slip line emerges from the triple point, which separates fluid from the oblique shock and the normal reflected shock. While the fluid behind the shock and along the slip is mechanically in balance, the difference in entropy provokes instabilities and triggers the formation of vortices. The flow field behind a bifurcated reflected shock is always inhomogeneous and the variations in temperature can provoke the ignition from small ignition kernels (mild-ignition) (Lipkowicz et al. 2019).

Remote Ignition
Under ideal conditions and corresponding homogeneous fields of pressure and temperature behind the reflected shock, the reactive mixture will always ignite at the end wall of the shock tube, as the ignition process depends exclusively on the time that passed after the compression by the reflected shock wave. This type of ignition is commonly referred to as strong ignition. Under real conditions, however, ignition processes were frequently observed that start from small ignition kernels at various positions in the test section (Fay 1953;Berets et al. 1950;Steinberg and Kaskan 1955) and often transition into a detonation. These ignition kernels can be located near the end wall (e.g., along the slip line in case of a bifurcated shock, which will be denoted as mild ignition) or further away, which we will refer to as remote ignition. For mild ignition to occur, the ignition-delay time in general must exceed a certain, yet unknown, limit so that flow-induced inhomogeneities can evolve in the flow field, which significantly reduces the local ignition-delay time. One criterion by Meyer and Oppenheim (1971) states that the change of ignition-delay time with respect to the change of temperature ( ig ∕ T) p must be below a critical value (e.g., ig ∕ T = −2μs /K for stoichiometric hydrogen/oxygen mixtures) such that mild ignition occurs. Many numerical studies were published studying the mild ignition phenomenon in shock tube-simulations in 2D (Ihme et al. 2013;Grogan and Ihme 2015;Oran and Gamezo 2007) and 3D (Khokhlov et al. 2015). The studies focussed on the impact of the incident shock Mach number (Oran and Gamezo 2007), the evolution of ignition kernels due to velocity fluctuations (Ihme et al. 2013), the effect of wall treatment on mild ignition (Grogan and Ihme 2015), and the development of hot spots in 3D (Khokhlov et al. 2015). However, all these studies featured cases with bifurcated reflected shocks, while remote ignition was also observed in the absence of bifurcated reflected shocks (Hanson et al. 2013). The present paper aims to shed light on the physics of remote ignition in cases with and without bifurcated reflected shocks.

Code
The simulations are carried out with the in-house code PsiPhi that has been developed at Imperial College in London and at the university of Duisburg-Essen by Kempf and coworkers (Lipkowicz et al. 2019;Inanc and Kempf 2019;Cifuentes et al. 2019;Rieth et al. 2016). PsiPhi solves the fully compressible set of Favre-filtered conservation equations for mass, momentum, total absolute internal energy and partial densities to simulate reactive flow problems: The equations feature the Favre-filtered velocity components ũ i in the ith direction, the Favre-filtered total absolute internal energy Ẽ , the viscous stress tensor ̄j i , the heat-flux density q i due to heat conduction and due to enthalpy fluxes caused by mass diffusion, the Favre-filtered mass fraction Ỹ k of the kth species, the source term ̄̇k of species k, the mixture-averaged diffusion ̄j i,k of species k, the sensible enthalpy h s,k of species k, and the correction velocity V i,c to achieve consistency between partial densities and the transported density.
Applying a LES (large-eddy simulation) filter operation to the Navier-Stokes equations leaves unclosed terms that need to be modeled. In this work, sub-filter fluxes are modeled with the eddy-viscosity approach for momentum and the eddy-diffusivity approach for scalars with a turbulent Schmidt-and turbulent Prandtl number of Sc t = Pr t = 0.7 , where the sub-grid viscosity is computed using the -model proposed by Nicoud et al. (2011). (The sigma model has been tested extensively against the static and dynamic Smagorinsky model, using the in-house code PsiPhi Rieth et al. 2014).
No modeling is required with regards to the filtered chemical source term, since gradients of scalars are small in the ignition regions. (This assumption is generally valid until combustion waves lead to strong spatial gradients of the scalars. However, this work focuses on the period up to the ignition only.) The finite-volume method (FVM) is utilized to discretize the equations on an equidistant, cartesian grid, where no local refinement or coarsening is applied to ensure a high level of consistency, even in the region behind the reflected shock, where transport and mixing must be resolved to predict weak ignition. PsiPhi uses a distributed memory domain decomposition approach, utilizing the message passing interface (MPI) and a non-blocking implementation for simultaneous computations and exchange of data, yielding high parallel efficiency. Diffusive fluxes are discretized using a 2nd order accurate central-difference scheme. The solution is advanced in time, using a low-storage explicit Runge-Kutta scheme (Williamson 1980) of 3rd order.
A wide range of velocities is present in simulations of shock-tube experiments. Hence, the all-speed approximate Riemann solver HR-Slau2 (Kitamura and Hashimoto 2016) developed by Kitamura, is used for the computation of convective fluxes, which reduces the contribution of the numerical dissipation term regarding the computation of the interface pressure in the low Mach-number limit. The states left and right to a cell interface are determined by a 5th order accurate monotonicity-preserving reconstruction scheme (MP5) by Suresh and Huynh (1997) that either reconstructs the local, one-dimensional characteristic variables or the set of primitive variables. The five-point stencil of the reconstruction scheme is also used to distinguish discontinuities from extrema such that the accuracy reduces to first order only next to discontinuities. The reconstruction scheme of Suresh and Huynh has been tested against the classical weighted essentially non-oscillatory (WENO5) scheme by Scandaliato and Liou (2012), where it proved to be more efficient and more accurate. Both of the reconstruction schemes have also been tested in a publication by Zhao et al. (2019) addressing the wave propagation errors of smooth waves, when interacting with discontinuities, where the MP5 scheme turned out to be a good compromise regarding wave propagation errors and was slightly more efficient than other schemes with a formal accuracy of 5th order. Further publications have dealt with the properties and defects of shock-capturing schemes. The reader is referred to the works of Pirozzoli (2006), Larsson (2010), Quirk (1997) and LeVeque (1998).
One of the disadvantages of using a high-order scheme and characteristic variables is the occurrence of strong oscillations under special circumstances, for example when two discontinuities interact (Harten et al. 1987). This can manifest in negative values of density and pressure. In order to avoid these unphysical solutions, recursive-order reduction (ROR) is partially applied (Gerolymos et al. 2009). A completely different problem concerns the numerical dissipation of a Riemann solver at low speeds, for example in a boundary layer. The jump of normal velocity components across a cell interface is the main contribution to numerical dissipation in the context of high-order, shock-capturing schemes, according to Thornber et al. (2008a). He proposed a simple fix (Thornber et al. 2008b) that relaxes the normal velocity components towards the arithmetic average as the low Mach-number limit is reached.
Thermochemical and transport properties of individual species and reaction-rate constants are first determined for each species with the aid of Cantera Goodwin (2009) and tabulated as a function of temperature to reduce computational effort. Reaction rates are computed during runtime, using the tabulated reaction-rate constants and the effective local concentrations. Furthermore, reaction-rate constants of fall-off reactions that depend on pressure, are tabulated not only as a function of temperature, but also as a function of effective concentration. The mixture-averaged molecular viscosity, heat conductivity, and molecular diffusion are determined by models of Wilke (1950), Peters and Warnatz (1982) and Kee et al. (2005) respectively. Direct chemistry is used, where the system of ordinary differential equations, is implicitly solved by CVODE Cohen et al. (1996) within a Strang (1968) operator-splitting framework. The reaction model FFCM-1 by Smith et al. (foundational fuel chemistry model, 29 reactions/11 species) (Smith et al. 2016) and alternatively the model by O'Conaire (40 reactions/11 species) (Oconaire et al. 2004) are used to simulate auto-ignition in hydrogen-oxygen mixtures.

Simulation Setup and Experimental Facilities
Simulations are conducted with non-reactive mixtures (NR) and reactive hydrogen-oxygen mixtures (R) at low pressure. The results are compared to shock-tube experiments from Berkeley (B) (Meyer and Oppenheim 1971), Duisburg (D), Stanford (S) (Hanson Table 1 Overview of two-, and three-dimensional simulations performed in the scope of this work D and W depict the geometry, where D is the diameter (or the height regarding ducts) and W is the width of the numerical domain in k-direction. The numerical grid resolution is denoted by Δ with the number of cells n I , n J , and n K in the corresponding directions The numerical domain of the main simulations covers the end part of the shock tube (13.5-132 cm) to allow a higher numerical resolution. In order to have realistic profiles of scalar quantities and velocities behind the incident shock wave in terms of a suitable initial solution, precursor simulations from the time of membrane rupture are performed. After the Mach number of the incident shock has reached the target value, a part of the solution behind the incident shock is stored and applied as an initial condition for the following main run. An isothermal no-slip boundary condition is used at all boundaries, except for the inlet of the numerical domain, located on the "left" of the numerical domain.
Simulations that cover only the end part of the shock tube require that the evolution of the state variables and the velocity at the inlet must be modeled. This concerns both the variation within the boundary layer and its growth, but also the variation of the state variables in the core flow. The perturbation theory by Mirels (1957) and Mirels and Braun (1957) is used to compute the variation of quantities outside the boundary at the Table 2 Overview of two-, and three-dimensional simulations performed in the scope of this work The shock Mach number M refers to the value just before the reflection. p 1 is the pressure of the initial quiescent gas in the driven section with the pressure p 5 and T 5 behind the reflected shock. The ideal ignitiondelay time ig,0 is the result of zero-dimensional reactors at constant volume/energy and ig is the result from the two-, and three-dimensional simulations. To compute the source terms of individual species caused by chemical conversion, either the O'Conaire mechanism or the foundational fuel chemistry model are used  Mirels (1955) by solving the Blasius differential equations in a shock-fixed frame, while the solution for the turbulent boundary layer is computed according to the equations by Petersen and Hanson (2001). The resulting profiles of the axial flow field and the temperature field are tabulated as function of time and are applied in terms of a Dirichlet boundary condition.

Non-reactive Cases
Non-reactive cases are simulated to compare the temporal pressure variation at the end wall to experimental measurements, thus validating the code and the modeling of boundarylayer effects. Figure 2 shows numerical schlieren visualizations after the reflection of the shock with argon ( Fig. 2a) and with nitrogen ( Fig. 2b) as test gases, by evaluating the absolute gradient of density. When nitrogen is used, a pronounced bifurcation of the reflected shock is present. The shear layer between reversed fluid and fluid that passes the oblique shock, produces turbulent kinetic energy, while vortices form along the slip line. As a result, the fields of state and velocity are highly inhomogeneous and the conditions are not ideal for shocktube experiments. Argon, on the other hand, typically suppresses bifurcation. The reflected shock is then curved due to a higher propagation speed of the reflected shock within the turbulent boundary layer. The absence of the bifurcated shock-induced vortices and shear layers leads to much smoother distributions of state quantities, which is a prerequisite for meaningful results from shock-tube experiments. However, the variation of state variables along the center line, introduced by the development of the boundary layer, still affects the state behind the reflected shock in space and time, especially because the variations are amplified by the reflected shock (Rudinger 1961). Fig. 2 Instantaneous numerical schlieren visualizations from two-dimensional simulations NRD1 (a) and a two-dimensional simulation with nitrogen as driven gas (b), to illustrate the effect of shock bifurcation on the region behind the reflected shock. The reflected shock is moving to the "left", away from the end wall on the "right". The schlieren visualizations are generated by computing the absolute gradient of density, followed by division with the maximum value and subsequent application of the decadic logarithm. The end wall of the shock tube is located on the right ( x = 0 mm) Figure 3 compares the evolution of pressure at the end wall for all non-reactive simulations against the respective measurements. The cases cover initial pressures from 44 to 120 mbar with pressures behind the reflected shock ranging from 1600 to 3200 mbar at Mach numbers of the incident shock from 2.3 to 2.8.

(b) (a)
Very good agreement is achieved in most of the cases, especially in terms of shock tubes with a large diameter (Fig. 3a, b) or at high pressure (Fig. 3c). At low pressure or small shock-tube diameters (e.g., Fig. 3d), deviations appear after 1 ms. Deviations however are expected, as the perturbation theory relies on the assumption that the thickness of the boundary layer is negligible compared to the height/diameter of the shock tube. At low pressure and/or small shock-tube diameters, these assumptions are . 3 Pressure histories from the non-reactive cases (dark black) and from the corresponding experiments (light orange), as well as temperature histories from simulation (dark purple) and from experiments (light yellow), evaluated at the end wall of the shock tube easily violated. Nevertheless, the results are initially consistent with those of the experiments, which would not have been the case with a primitive inflow condition neglecting the evolution of the boundary layer. The remaining deviations can be partially attributed to shock-contact surface interaction or the arrival of the expansion wave, effects that are not considered in the simulations. Strikingly, simulations and experiments (Fig. 3a, b, e-g) both show an initial decrease of pressure followed by a linear increase, which is in contrast to the usual expectation of a purely linear increase of pressure. This assumption of a linear increase of pressure is also incorporated in many low-order models for reactors that are used for the validation of reaction mechanisms, thus neglecting the observed behaviour could lead to large errors. It is obvious to attribute the unexpected drop of pressure to transition effects of the boundary layer. Therefore, the observed characteristic evolution of pressure is only expected at low pressures behind the incident shock, when the laminar boundary layer can not be neglected. Figure 4 presents stacked center-line profiles of both pressure and temperature for simulation NRD1, which is in excellent agreement with the experiment. The normalized pressure increase at the end of the shock tube is linear at a value of p 5 ∕ t∕p 5 ≈ 3.6 %/ms and presents the only simulation with a linear pressure evolution in this study. In contrast Fig. 4 Stacked center-line profiles of simulation NRD1 of temperature (top) and pressure (bottom), illustrated as a function of end-wall distance x and time after shock reflection t. The end wall of the shock tube is located on the right ( x = 0 cm) to the results shown in Fig. 3, the surfaces illustrate the entirety of changes in time and space, where the profile for x = 0 cm (end wall) in the lower panel refers to the solution in Fig. 3c. According to the results, the strength of the reflected shock increases, while travelling upstream (away from the end wall), which is reflected in higher pressures and temperatures behind the shock. However, the evolution of temperature and pressure at a fixed location is very different directly at the end wall, compared to locations further away. This is well illustrated by the fact that for a fixed time of t = 2.5 ms, the pressure decreases with the distance from the end wall whereby the temperature increases. While pressure and temperature are connected by isentropic relations at the end wall, this is clearly not the case further away from the end wall, a circumstance which is to be led back among other things to the variation of entropy by shock attenuation. The temperature increases continuously with distance and time such that the temperature maximum of the presented data is reached for x = 40 cm and t = 2.5 ms and is significantly larger ( ≈ 50 K) than the temperature at the end wall at the same time. Such a temperature distribution could lead to a remote ignition, if a reactive mixture were used instead. Figure 5 also presents stacked center-line profiles, but for simulation NRD3, a simulation that matches the experimental results although the initial pressure is considered very Fig. 5 Stacked center-line profiles of simulation NRD3 of temperature (top) and pressure (bottom), illustrated as a function of end-wall distance x and time after shock reflection t. The end wall of the shock tube is located on the right ( x = 0 cm) low. The differences of pressure and temperature surfaces from the previous case (Fig. 4) are apparent. In this case, a characteristic valley forms both in pressure and in temperature, independent of the end-wall distance. While the evolution at a fixed location is qualitatively similar, the strength of the changes (gradients) decreases with the wall distance. As in the previous case, pressure decreases with wall distance at a simulation time of t = 2.5 ms, whereas the temperature increases, while the distribution of both the quantities along the center line is much more homogeneous compared to the previous case. The pressure and temperature distributions at these low pressure levels are over all very complex and the values vary strongly in time. At this point we want to emphasize that it will be very important to quantify such effects in low-pressure experiments to be able to interpret the measurement results.
Slice-integrated profiles of the mass flux per unit depth, are presented in Fig. 6 for the cases NRD1 and NRD3 and for different times. Without viscous effects and heat losses, a reflected shock of constant strength forms, such that the fluid behind the reflected shock  NRD3 (b). Colors indicate the time that has passed since the reflection of the shock wave. The end wall of the shock tube is located on the right ( x = 0 mm) wave is instantly at rest. The local distribution of the state variables of the investigated cases in contrast leads to a change in shock strength and, for example, to a change in the momentum of the fluid behind the reflected shock, as presented in the panels (a) and (b) of Fig. 6. In the case of NRD1, this means that the fluid behind the reflected shock wave still has a residual momentum. It is eventually brought to a rest, accompanied by an increase in pressure and temperature. In the case of NRD2, however, the fluid behind the reflected shock has a negative momentum, hence moving in the direction of the reflected shock. The gas behind the reflected shock thus expands, decreasing the temperature. A onedimensional inviscid simulation of case NRD1, presented in supplementary material, confirms that the fluctuations visible in panel (a) do not stem from the applied algorithms, but instead are linked to the transition from the laminar to the turbulent boundary layer and occur first at the inlet, where artificial turbulence is created in the boundary layer. The much higher pressure in case NRD1, in contrast to that of case NRD3, causes a very early transition, which is why no fluctuations are visible in panel (b), since the boundary layer has not yet turned over at this point.

Remote Ignition Simulated in 2D
The observed agreement of experiments and simulations, both qualitative and quantitative, indicates that the most important phenomena including boundary-layer effects have been simulated successfully. Therefore the code can be used to also examine remote ignition events. Figure 7 presents temperature fields at different instances of simulation RS1. Since the mixture ignites simultaneously in a region near the end wall, this ignition can be (a) Fig. 7 Instantaneous temperature fields of simulation RS1 at different times after the reflection of the shock wave, illustrating strong ignition and the subsequent formation of a strong wave. The end wall of the shock tube is located on the right ( x = 0 mm) classified as a strong ignition. The pressure increase resulting from the combustion is particularly strong due to the closed end of the shock tubes, and a strong "left"-running wave is formed.
As highlighted in Table 2, the Mach number in case RS2 and case RS3 is slightly lower than that in case RS1, but this small difference is sufficient for the ignition to take a  Oconaire et al. (2004)] at different times after the reflection of the shock, illustrating mild ignition remotely from the end wall. The end wall of the shock tube is located on the right ( x = 0 mm) different course of events, as shown in Figs. 8 and 9. Both simulations (RS2, RS3) use the same initial-and boundary conditions and differ only in the mechanism used for solving the chemistry (FFCM-1 in terms of simulation RS2 and O'Conaire in terms of RS3).
According to Fig. 8a, the ignition starts from small ignition kernels, located at a distance of approximately 500 mm away from the end-wall. More ignition kernels appear at t = 3.04 ms between end-wall distances of 400 and 600 mm before the whole mixture ignites remotely at t = 3.08 ms. Also, an increase of temperature is observed at the end wall, which suggests that these conditions mark the transition from a strong to a remote ignition. This means that the reduction of the ignition-delay time caused by fluid dynamics away from the end wall just compensates for the delayed compression by the reflected shock away from the end wall.
If the reaction mechanism of O'Conaire is used instead, the ignition event slightly deviates from the previous result, as can be seen in Fig. 9, and attributed to uncertainties of the reaction mechanisms at low temperatures. This time, ignition kernels are already visible after 2.93 ms and are located even further away from the end wall at a distance of 600 mm, while the mixture is consumed more rapidly. However, the largest deviation in comparison to simulation RS2 concerns the region near the end wall, where no significant temperature increase is observed. Hence, the competition of the characteristic time scales favors the remote ignition event. According to Table 2, boundary-layer effects in each of the simulations (RS1, RS2, RS3) greatly reduce the ignition-delay time ig compared to the ideal ignition-delay time ig,0 as obtained from low-order simulations. The reduction is particularly pronounced in cases RS2 and RS3, where remote ignition occurs.
The local heat-release rate ̇H R plays an obvious and important role in the ignition process as it results from chemical conversion and accelerates it at the same time. Figures 10  and 11 present scatter plots of local heat-release rate over temperature, colored with the end-wall distance. Figure 10 shows the result for simulation RS1 and thus for the case of a strong ignition. As expected, the heat-release rates are initially highest at the end wall. Nothing changes subsequently in this overall picture despite higher temperatures at greater Fig. 10 Scatter plot of local heat-release rate ̇H R over temperature T and colored with the respective endwall distance x for simulation RS1 distances from the wall. A completely different picture emerges for Simulation RS2. Figure 11b presents the distribution of local heat-release rate of simulation RS2 at the same time as Fig. 10a for simulation RS1. The heat-release rates at the end wall in this case are an order of magnitude below those observed in simulation RS1. At a time of 2 ms after the reflection of the shock, this gap has further increased, and the heat-release rates at the end wall of simulation RS1 now exceed those of simulation RS2 by two orders of magnitude. In contrast to simulation RS1, simulation RS2 shows more concentrated distributions of heatrelease rates prior to ignition, again emphasizing that flow-induced temperature inhomogeneities upstream made up for the delayed compression.
In order to illustrate how these temperature variations affect the localization of ignition, we determine the expected time of ignition after shock-wave reflection (i.e., a "local" ignition delay time) in a post-processing step. For each numerical cell of the computational domain, the instantaneous thermochemical state is utilized to estimate the related ignition-delay time based on the assumption of isochoric 0D reactors. The results are therefore decoupled from convection and diffusion. The first panel of Fig. 12 presents the field of expected time of ignition, 0.75 ms after the reflection of the shock for simulation RS1. A general spatial gradient of the expected time of ignition is recognizable in axial direction, favouring ignition near the end wall, whereby the expected ignition times vary strongly, specifically near the walls of the shock tube. As the process progresses, these local gradients disappear at the end wall, so that the ignition is globally initiated instead of an ignition from smaller kernels. Results of the same type of post-processing, but for simulation RS2, are shown in Fig. 12d-f. In contrast to the results from simulation RS1, it is not possible to predict where the ignition will take place based on the result 0.34 ms after the reflection of the shock. The field of expected ignition times remains heterogeneous until the point of ignition, with deviations between the shortest and highest expected ignition time of about 0.3 ms. The small kernels with the shortest expected ignition times dictate the subsequent ignition process, starting at a distance of 500 mm. In order to facilitate the interpretation and to quantify the flow-induced reduction of expected ignition time as a function of axial location, the fields of expected ignition time are averaged in the vertical direction ( < ⋅ > ), followed by a filter operation ( ⋅ ) to eliminate fluctuations. Figure 13 presents the results for simulation RS1 on the left and RS2 on the right. Solid lines in the upper panels show the averaged and filtered profiles of expected ignition time at time t n . The dashed lines show projected profiles t * n resulting from shifting the previous profile t n−1 by the difference in simulation time between the samples, i.e., by Δt = t n − t n−1 . In this discussion, it is worth pointing out that the expected time of ignition increases monotonically with x for RS1, whereas it is almost constant for RS2, making the location of ignition a lot more sensitive to small perturbations. A flow-induced reduction  Fig. 12 Expected time of ignition, using instantaneous values from simulation RS1 (a-c) and from simulation RS2 (d-f) as initial condition for isochoric 0D-reactors. The end wall of the shock tube is located on the right ( x = 0 mm) of ignition time Δ fi is introduced by subtracting the computed profile of expected ignition time from the projected profile. This variable approximates the reduction of ignition time due to fluid dynamics within a given time interval. It is interesting to note that fluid dynamics reduce the ignition time independent from the axial location. Nevertheless, the trend can be observed that the expected ignition time is reduced more at greater distance to the end wall. In case of simulation RS2, the initially flat profile of expected ignition time thus gets altered by the fluid dynamics such that remote ignition occurs.

Remote Ignition Simulated in 3D
While remote ignition in the previous cases is governed by effects due to the formation of a boundary layer, remote ignition can also be related to the flow field evolving behind bifurcated shocks. In contrast to the simulations in 2D, no complex inlet treatment is applied with regards to simulations RB1 and RB2, as the bifurcation-induced effects are studied exclusively. Instead, the solution resulting from the ideal shock relations is applied. Numerical schlieren visualizations, shortly after the ignition, are presented in Fig. 14.
The upper image shows the result of simulation RB1, where the mixture ignites in the core of the shock tube at a distance of approximately 70 mm away from the end wall. Reasonable agreement with the experiments by Meyer and Oppenheim (1971) has been achieved with an ignition-delay time of 281 μ s (approximated by the time, when the  (Table 2), the ignition-delay time from this simulation is reduced by a factor of four by taking into account the shock-boundary layer interaction.
The ignition locations in simulation RB2 are very different, as the ignition starts near the corners close to the slip line before a larger volume spanning the entire cross section of the shock tube ignites. Many of the schlieren photographs by Meyer and Oppenheim (1971) showed very similar ignition locations at comparable temperatures behind the reflected shock. The reasons for the remote ignition in case RB1 were discussed in detail before (Lipkowicz et al. 2019), but are described briefly again for the sake of completeness.
One of the striking features is the formation of a second normal shock after the bifurcation structure has grown significantly. The second shock is clearly visible in the upper panel of Fig. 14 Table 2 in terms of case RB1 each of which are increasing the local temperature and are thus moving at a higher speed than the waves upstream. Eventually, they will catch up with other waves and form the second shock wave. Similar shock-wave patterns were observed in two-dimensional simulations as reported by Weber et al. (1995). Figure 16 introduces instantaneous fields of pressure, temperature, and hydroperoxyl (HO 2 ) complemented by center-line plots of pressure and temperature before the ignition has taken place, to fully understand the physics. According to Fig. 15, hydroperoxyl can be used as an auto-ignition marker, since the mass fraction increases monotonically while the growth is almost perfectly exponential over a wide range. In contrast to the previously demonstrated auto-ignition marker, where the expected time of ignition is evaluated based on instantaneous fields, this marker and the use of a logarithmic scale reveals features that are not visible in the other marker field. Fig. 16 Instantaneous fields of pressure, temperature and HO 2 species mass fraction from simulation RB1, as well as center-line plots of pressure and temperature at different times from simulation RB1. The end wall of the shock tube is located on the right ( x = 0 mm) The existence of the second normal shock (travelling at a similar speed compared to the reflected shock) implies an acceleration of the fluid behind the reflected shock (referred to as the core fluid) to supersonic speeds in a coordinate frame that is fixed to the reflected shock. The flow field behind the reflected shock is determined by displacement effects of outer-core fluid and pressure variations within the bifurcation structure. The displacement effect of the outer-core fluid is caused by the oblique shock, as the fluid that passes the oblique shock gains momentum in vertical direction. Combined, these effects force the fluid in the core to follow a convergent-divergent streamline pattern which forms a Lavalnozzle shaped stream tube, which is in fact well illustrated by the interface between coreand outer fluid along the slip line in the temperature field. Initially, when the bifurcation structure is small and the cross-section areas, characterizing the Laval-nozzle like flow and which are bounded by the slip line, are correspondingly large, the variations of velocity and state are also small. However, due to the growth of the bifurcation, these variations will increase. Since the pressure drops in the core flow, as a result of the velocity gain, fluid will be also accelerated from the end wall towards the reflected shock. Once the velocity in the core-flow has nearly reached super sonic speed, these pressure waves will start to "pile up" and eventually form the second shock. This evolution is well captured by the history of center-line plots in Fig. 16. The additional production of entropy, caused by the second (b) (c) (a) Fig. 17 Expected time of ignition, using instantaneous values from simulation RB1 as initial condition for isochoric 0D-reactors. The end wall of the shock tube is located on the right ( x = 0 mm) shock, is clearly evident by a sustained offset of temperature compared to the temperature at the end wall.
Apparently, this temperature variation has huge implications regarding the distribution of ignition-delay times, as presented in Fig. 17. Comparing these results to those without the occurrence of a bifurcation, it is noticeable that the values now differ greatly from each other, namely between 0.3 and 0.9 ms in the first panel, which presents the result 99 μ s after the reflection of the shock. The smallest values are within a range of 20-40 mm from the end wall. After another 80 μ s, the minimum has shifted further to the "left", now in a distance of about 70 mm to the end wall. Due to the high sensitivity of the ignition-delay time to temperature changes, the increased temperature in this region compensates for the delay regarding the compression of the reflected shock. Further away from the end wall (around 70 mm), a constriction is visible (also present in Fig. 16), where fluid is forced from the walls towards the core. The location coincides approximately with the formation of the second normal shock. In the following, the cold fluid from the walls mixes with that from the core and prevents ignition in this region, as can be seen in the third panel of Fig. 17. According to Fig. 16, the temperature to the "left" of the constriction is even higher, while this temperature difference does not compensate for the time elapsed for the (a) (b) (c) Fig. 18 Expected time of ignition, using instantaneous values from simulation RB2 as initial condition for isochoric 0D-reactors. The end wall of the shock tube is located on the right ( x = 0 mm) shock wave to process the gas upstream. Therefore, the mixture ignites to the "right" of the mixing zone and much earlier than the ignition at the end wall would have occurred.
The flow characteristics of simulation RB2 are very similar to those of simulation RB1. However, since the second normal shock has not even formed at the time of ignition in case of simulation RB2 and since the ignition starts first from small kernels near the slip line, the gas dynamics responsible for remote ignition in this case must be very different. Figure 18 presents the expected time of ignition for simulation RB2.
Here, the conditions behind the reflected shock promote a faster ignition compared to simulation RB1. Moreover, the field in the first panel appears much more homogeneous, and the results are overall more comparable with the results of a strong ignition from simulation RS1. A wave pattern is also conspicuous, whereby the orientation of the waves suggests that their origin lies in the bifurcation structure.
In order to investigate the physics of the ignition mechanism in this particular case, Fig. 19 shows instantaneous fields of the dissipation rate of kinetic energy, as well as entropy of the mixture. Aside from the expected high dissipation rates within the bifurcation and within the boundary layer, high values are also observed along the slip line and in close proximity to the slip line. Especially after the break up of the slip line and the subsequent formation of vortices, a broad region with very high dissipation values of more than 10 7 W/kg is present. The high dissipation rates are also reflected in the field of entropy, where the values are particularly high in the zones of the following ignition kernels. Assuming that a fluid element in close proximity to the slip line is affected by shear with a dissipation of 5 × 10 7 W/kg, over a period of 100 μ s, the estimated temperature increase based on these values would be 5 K, enough to explain a slightly earlier ignition according to the temperature sensitivity of ignition-delay times.
(a) (b) Fig. 19 Instantaneous fields of dissipation of kinetic energy (a) and entropy (b) from simulation RB2. The end wall of the shock tube is located on the right ( x = 0 mm) Compared to the ideal ignition-delay time ig,0 of 142 μ s (Table 2), the ignition-delay time of 140 μ s from this simulation turns out to be almost identical, which is also supported by the strong ignition in the entire volume shortly after. However, it is conceivable that the ignition mechanism proposed here could lead to a more significant reduction in other cases.

Conclusion
Two-dimensional computational fluid dynamics simulations of shock-tube experiments with non-reactive mixtures, as well as two-and three-dimensional simulations of shocktube experiments with hydrogen-oxygen mixtures were carried out using high-order numerics and detailed chemistry.
The non-reactive cases confirmed the measurements of pressure and temperature at the end wall in a case with a linear increase of pressure and in several cases, where the pressure initially decreased. This phenomenon of a decreasing pressure, followed by the transition into a linear increase, can be clearly attributed to the transition of the boundary layer, whose effects were modeled at the inlet of the computational domain. Taking these effects into account will be very important for low-order models to provide reliable results at low pressure. A strong increase of temperature with wall distance was observed, likely resulting from shock attenuation and decoupled from the pressure.
The event of remote ignition in an argon-diluted hydrogen-oxygen mixture was reproduced using two different reaction mechanisms, while another simulation at a slightly higher Mach number showed ignition in strong ignition mode. Differences of the results (using different reaction mechanisms) qualitatively and quantitatively underline the large uncertainties of reaction mechanisms at low temperatures and the need for improvement. Clearly, the remote ignition phenomenon in these cases was caused by the variation of temperature along the center line, resulting from the attenuation of the incident shock.
The remote ignitions, observed in the three-dimensional simulations, however, were caused exclusively by the complex fluid dynamics behind the reflected shock, as boundary layer effects were not modeled at the inlet of the computational domain. If the ignitiondelay time is sufficiently long, a second normal shock might evolve before ignition takes place, leading to additional formation of entropy and regions of higher temperature, thus reducing the local ignition-delay time. In the first case investigated, this mechanism led to remote ignition, whereas remote ignition was initiated in the second case, before the second normal shock had formed. Instead, shear in close proximity to the slip line caused a temperature increase, sufficient for a premature ignition from small kernels.