Abstract

Numerical simulation can be key to the understanding of the multidimensional nature of transient detonation waves. However, the accurate approximation of realistic detonations is demanding as a wide range of scales needs to be resolved. This paper describes a successful solution strategy that utilizes logically rectangular dynamically adaptive meshes. The hydrodynamic transport scheme and the treatment of the nonequilibrium reaction terms are sketched. A ghost fluid approach is integrated into the method to allow for embedded geometrically complex boundaries. Large-scale parallel simulations of unstable detonation structures of Chapman-Jouguet detonations in low-pressure hydrogen-oxygen-argon mixtures demonstrate the efficiency of the described techniques in practice. In particular, computations of regular cellular structures in two and three space dimensions and their development under transient conditions, that is, under diffraction and for propagation through bends are presented. Some of the observed patterns are classified by shock polar analysis, and a diagram of the transition boundaries between possible Mach reflection structures is constructed.

1. Introduction

Reacting flows have been a topic of on-going research since more than one hundred years. The interaction between hydrodynamic flow and chemical kinetics can be extremely complex and even today many phenomena are not very well understood. One of these phenomena is the propagation of detonation waves in gaseous media. While detonations propagate at supersonic velocities between 1000 and 2000 m/s, they inhibit nonnegligible instationary substructures in the millimeter range. Experimental observations can provide only limited insight and it is therefore not surprising that the understanding of the multidimensionality has improved little since the first systematic investigations [1, 2]. An alternative to laboratory experiments is numerical simulations of the governing thermo- and hydrodynamic equations. But the additional source terms modeling detailed nonequilibrium chemistry are often stiff and introduce new and extremely small scales into the flow field. Their accurate numerical representation requires finite volume meshes with extraordinarily high local resolution.

In this paper, we summarize a significant body of multiyear work [38] devoted to simulating multidimensional detonations with detailed and highly stiff chemical kinetics on parallel machines with distributed memory, especially on clusters of standard personal computers. The objective of all presented computations is to resolve the transverse-wave structures, that form the characteristic cellular structures intrinsic to detonation propagation in two and three space dimensions, accurately and without computational ambiguity. We sketch the design of AMROC (Adaptive Mesh Refinement in Object-oriented C++), a freely available dimension-independent mesh adaptation framework for time-explicit Cartesian finite volume methods on distributed memory machines [9, 10], and discuss briefly the locality-preserving rigorous domain decomposition technique it employs [11]. The framework provides a generic implementation of the adaptive mesh refinement (AMR) algorithm after Berger and Colella [12] designed especially for the solution of hyperbolic fluid flow problems on logically rectangular grids. The ghost fluid approach is integrated into the refinement algorithm to allow for complex moving boundaries represented implicitly by additional scalar level set variables [13, 14]. Briefly, we describe the employed numerical methods and the treatment of the non-equilibrium reaction terms. In Section 6, several two- and three-dimensional simulations of cellular detonation structures in low-pressure hydrogen-oxygen-argon mixtures simulated in purely Cartesian geometry and with embedded static boundaries are presented. The paper concludes by demonstrating that the high-resolution adaptive computations allow a detailed quantitative analysis of the observed transverse-wave patterns. It is found that nonreactive shock wave reflection theory [15], extended to thermally perfect gas mixtures, is applicable to classify the observed transient triple point structures unambiguously as transitional and double Mach reflection patterns.

2. Detonation Theory

A detonation is a shock-induced combustion wave that internally consists of a discontinuous hydrodynamic shock wave followed by a smooth region of decaying combustion. The adiabatic compression due to the passage of the shock rises the temperature of the combustible mixture above the ignition limit. The reaction results in an energy release driving the shock wave forward. In a self-sustaining detonation, shock and reaction zone propagate essentially with an identical speed 𝑑𝐶𝐽 that is approximated to good accuracy by the classical Chapman-Jouguet (CJ) theory, see for example, [16]. But up to now, no theory exists that describes the internal flow structure satisfactory. The Zel'dovich-von Neumann-During (ZND) theory is widely believed to reproduce the one-dimensional detonation structure correctly (see e.g., Figure 1), but already early experiments [1] uncovered that the reduction to one space dimension is not even justified in long combustion devices. It was found that detonation waves usually exhibit nonnegligible instationary multidimensional substructures and do not remain planar. The multidimensional instability manifests itself in instationary shock waves propagating perpendicular to the detonation front. A complex flow pattern is formed around each triple point, where the detonation front is intersected by a transverse shock. Pressure and temperature are increased remarkably in a triple point and the chemical reaction is enhanced drastically giving rise to an enormous local energy release. Hence, the accurate representation of triple points is essential for safety analysis, but also in technical applications, for example, in the pulse detonation engine. Some particular mixtures, for example, low-pressure hydrogen-oxygen with high argon dilution, are known to produce very regular triple point movements. The triple point trajectories form regular “fish-scale” patterns, so-called detonation cells, with a characteristic length 𝐿 and width 𝜆 (compare Figure 2).

Figures 2 and 3 display the hydrodynamic flow pattern of a detonation with regular cellular structure as it is known since the early 1970s, see for example, [2] or [17]. Figure 3 shows the periodic wave configuration around a triple point in detail. It consists of a Mach reflection, a flow pattern well known from nonreactive supersonic hydrodynamics (see e.g., [15, 18] for details). The undisturbed detonation front is called the incident shock, while the transverse wave takes the role of the reflected shock. The triple point is driven forward by a strong shock wave, called Mach stem. Mach stem and reflected shock enclose the slip line, the contact discontinuity.

The Mach stem is always much stronger than the incident shock, which results in a considerable reduction of the induction length 𝑙ig, the distance between leading shock and measurable reaction. The shock front inside the detonation cell travels as two Mach stems from point A to the line BC. In the points B and C, the triple point configuration is inverted nearly instantaneously and the front in the cell becomes the incident shock. Along the symmetry line AD, the change is smooth and the shock strength decreases continuously. In D, the two triple points merge exactly in a single point. The incident shock vanishes completely and the slip line, which was necessary for a stable triple point configuration between Mach stem and incident shock, is torn off and remains behind. Two new triple points with two new slip lines develop immediately after D.

3. Governing Equations

The appropriate model for detonation propagation in premixed gases with realistic chemistry are the inviscid Euler equations for multiple thermally perfect species with reactive source terms [16, 19]. These equations form a system of inhomogeneous hyperbolic conservation laws that reads𝜕𝑡𝜌𝑖𝜌+𝑖𝑢=𝑊𝑖̇𝜔𝑖,𝜕𝑡𝜕𝜌𝑢+𝜌𝑢𝑢+𝑝=0,𝑡(𝜌𝐸)+(𝜌𝐸+𝑝)𝑢=0,(1) with 𝑖=1,,𝐾. Herein, 𝜌𝑖 denotes the partial density of the 𝑖th species and 𝜌=𝐾𝑖=1𝜌𝑖 is the total density. The ratios 𝑌𝑖=𝜌𝑖/𝜌 are called mass fractions. We denote the velocity vector by 𝑢 and 𝐸 is the specific total energy. We assume that all species are ideal gases in thermal equilibrium, and the hydrostatic pressure 𝑝 is given as the sum of the partial pressures 𝑝𝑖=𝑇𝜌𝑖/𝑊𝑖 with denoting the universal gas constant and 𝑊𝑖 the molecular weight, respectively. The evaluation of the last equation requires the previous calculation of the temperature 𝑇. As detailed chemical kinetics typically require species with temperature-dependent material properties, each evaluation of 𝑇 involves the approximate solution of an implicit equation by Newton iteration [3].

The chemical production rate for each species is derived from a reaction mechanism of 𝐽 chemical reactions aṡ𝜔𝑖=𝐽𝑗=1𝜈𝑟𝑗𝑖𝜈𝑓𝑗𝑖𝑘𝑓𝑗𝐾𝑙=1𝜌𝑙𝑊𝑙𝜈𝑓𝑗𝑙𝑘𝑟𝑗𝐾𝑙=1𝜌𝑙𝑊𝑙𝜈𝑟𝑗𝑙,(2) with 𝜈𝑓/𝑟𝑗𝑖 denoting the forward and backward stoichiometric coefficients of the 𝑖th species in the 𝑗th reaction. The rate expressions 𝑘𝑗𝑓/𝑟(𝑇) are calculated by an Arrhenius law, see for example, [16].

4. Numerical Methods

We use the time-operator splitting approach or method of fractional steps [20] to decouple hydrodynamic transport and chemical reaction numerically. This technique is most frequently used for time-dependent reactive flow computations. The homogeneous Euler equations and the usually stiff system of ordinary differential equations𝜕𝑡𝜌𝑖=𝑊𝑖̇𝜔𝑖𝜌1,,𝜌𝐾,𝑇,𝑖=1,,𝐾,(3) are integrated successively with the data from the preceding step as initial condition. The advantage of this approach is that a globally coupled implicit problem is avoided and a time-implicit discretization, which accounts for the stiffness of the reaction terms, needs to be applied only local in each finite volume cell. We use a semi-implicit Rosenbrock-Wanner method [21] to integrate (3) within each cell. Temperature-dependent material properties are derived from lookup tables that are constructed during start-up of the computational code. The expensive reaction rate expressions (2) are evaluated by a mechanism-specific Fortran-77 function, which is produced by a source code generator on top of the Chemkin-II library [22] in advance. The code generator implements the reaction rate formulas without any loops and inserts constants such as 𝜈𝑓/𝑟𝑗𝑖 directly into the code.

Since detonations involve supersonic shock waves, we use a finite volume discretization that achieves a proper upwinding in all characteristic fields (cf. [4, 6, 7]). The scheme utilizes a quasi-one-dimensional approximate Riemann solver of Roe-type [23] and is extended to multiple space dimensions via the method of fractional steps, see for example, [24]. To circumvent the intrinsic problem of unphysical total densities and internal energies near vacuum due to the Roe linearization, cf. [25], the scheme has the possibility to switch to the simple, but extremely robust Harten-Lax-Van Leer (HLL) Riemann solver. Negative mass fraction values are avoided by a numerical flux modification [26]. Finally, the occurrence of the disastrous carbuncle phenomena, a multidimensional numerical crossflow instability that destroys every simulation of strong grid-aligned shocks or detonation waves completely and has been documented first by Quirk [27], is prevented by introducing a small amount of additional numerical viscosity in a multidimensional way, cf. [28]. A detailed derivation of the entire Roe-HLL scheme including all necessary modifications can be found in [3]. This hybrid Riemann solver is extended to a second-order accurate method with the MUSCL-Hancock variable extrapolation technique by Van Leer (see [24] for a detailed derivation).

4.1. Embedded Complex Boundaries

Higher-order shock-capturing finite volume schemes are most efficient on rectangular Cartesian grids. In order to consider geometrically complex moving boundaries within the Cartesian upwind scheme outlined above, we use some of the finite volume cells as ghost cells to enforce immersed boundary conditions, cf. [13] or [29]. Their values are set immediately before the original numerical update to model moving embedded walls and can be understood as part of the numerical scheme. The boundary geometry is mapped onto the Cartesian mesh by employing a scalar level set function 𝜑 that stores the signed distance to the boundary surface and allows the efficient evaluation of the boundary outer normal in every mesh point as 𝑛=𝜑/|𝜑| [30]. A cell is considered to be a valid fluid cell within the interior if the distance in the cell midpoint is positive and is treated as exterior otherwise. The numerical stencil by itself is not modified, which causes a slight diffusion of the boundary location throughout the method and results in an overall nonconservative scheme. We alleviate such errors and the unavoidable staircase approximation of the boundary with this approach effectively by using the dynamic mesh adaptation technique described in Section 5 to also refine the Cartesian mesh appropriately along the boundary. Other authors have also presented cut-cell techniques that utilize the correct boundary flux [31, 32]; however, the proposed numerical circumventions of the severe time step restriction in time-explicit schemes, that can result from small cells created by the boundary intersection, are inherently complicate and most approaches have not been extended successfully to three space dimensions yet.

For the inviscid Euler equations (1), the boundary condition at a rigid wall moving with velocity 𝑤 is 𝑢𝑛=𝑤𝑛. Enforcing the latter with ghost cells, in which the discrete values are located in the cell centers, requires the mirroring of the primitive values 𝜌𝑖, 𝑢, 𝑝 across the embedded boundary. The normal velocity in the ghost cells is set to (2𝑤𝑛𝑢𝑛)𝑛, while the mirrored tangential velocity remains unmodified. The construction of the velocity vector within the ghost cells therefore reads2𝑡𝑡𝑢=𝑤𝑛𝑢𝑛𝑛+𝑢=2𝑤𝑢𝑛𝑛+𝑢(4) with 𝑡 denoting the boundary tangential. This construction of discrete values in ghost cells (indicated by gray) for (1) is depicted in one space dimension in Figure 4.

The utilization of mirrored ghost cell values in a ghost cell center 𝑥 requires the calculation of spatially interpolated values in the point̃𝑥=𝑥+2𝜑𝑛(5) from neighboring interior cells. For instance in two space dimensions, we employ a bilinear interpolation between usually four adjacent cell values, but directly near the boundary the number of interpolants needs to be decreased, cf. Figure 5. It has to be underlined that an extrapolation in such situations is inappropriate for hyperbolic problems with discontinuities like detonation waves that necessarily require the monotonicity preservation of the numerical solution. Figure 5 highlights the reduction of the interpolation stencil for some exemplary cases close to the embedded boundary. The interpolation locations according to (5) are indicated by the origins of the red arrows. See [7] for further details.

4.2. Meshes for Detonation Simulation

Numerical simulations of detonation waves require computational meshes that are able to represent the strong local flow changes due to the reaction correctly. In particular, the shock of a detonation wave with detailed kinetics can be very sensitive to changes of the reaction behind, and if the mesh is too coarse to resolve all reaction details correctly, the Riemann Problem at the detonation front is changed leading to a wrong speed of propagation. We make a simple discretization test in order to illustrate how fine computational meshes for accurate detonation simulations in fact have to be. Figure 1 displays the density and temperature profile for the frequently studied H2 : O2 : Ar Chapman-Jouguet detonation with molar ratios of 2 : 1 : 7 in the unburned mixture with initial temperature 𝑇0=298 K and pressure 𝑝0=6.67 kPa. Figure 6 shows the distributions of 𝑌H2O and 𝑌H2O2 according to the ZND detonation model for this detonation discretized with different grids. Apparently, a resolution of 4 finite volumes per induction length (4 Pts/𝑙ig with 𝑙ig1.404 mm) would not be sufficient to capture the maximum of the intermediate product H2O2 correctly. This requires at least 5 to 6 Pts/𝑙ig; however, in triple points significantly finer resolutions can be expected. As discretizations of realistic combustors with such fine uniform meshes typically would require up to 109 points in the two- and up to 1012 points in the three-dimensional case, the application of a dynamically adaptive mesh refinement technique is imperative.

5. Adaptive Mesh Refinement

In order to supply the required temporal and spatial resolution efficiently, we employ the block-structured adaptive mesh refinement (AMR) method after Berger and Colella [12], which is tailored especially for hyperbolic conservation laws on logically rectangular finite volume grids. We have implemented the AMR method in a generic, dimension-independent object-oriented framework in C++. It is called AMROC (Adaptive Mesh Refinement in Object-oriented C++) and is free of charge for scientific use [9] (Note that the most recent version V2.0 of AMROC is part of the Virtual Test Facility software [10], available at http://www.cacr.caltech.edu/asc.). An efficient parallelization strategy for distributed memory machines has been found, and the codes can be executed on all systems that provide the MPI library.

5.1. Structured AMR Method

Instead of replacing single cells by finer ones, as it is done in cell-oriented refinement techniques, structured AMR methods follow a patch-oriented approach. Cells being flagged by various error indicators (shaded in Figure 7) are clustered with a special algorithm [33] into nonoverlapping rectangular grids. Refinement grids are derived recursively from coarser ones, and a hierarchy of successively embedded levels is thereby constructed, see for example, Figure 7. All mesh widths on level 𝑙 are 𝑟𝑙-times finer than on level 𝑙1, that is, Δ𝑡𝑙=Δ𝑡𝑙1/𝑟𝑙 and Δ𝑥𝑛,𝑙=Δ𝑥𝑛,𝑙1/𝑟𝑙 with 𝑟𝑙2 for 𝑙>0 and 𝑟0=1, and a time-explicit finite volume scheme (in principle) remains stable on all levels of the hierarchy. The recursive integration order visualized in Figure 8 is an important difference to usual unstructured adaptation strategies and one of the main reasons for the high efficiency of the approach.

The numerical scheme is applied on level 𝑙 by calling a single-grid routine in a loop over all subgrids. The subgrids become computationally decoupled by employing additional ghost cells around each computational grid. Three types of different ghost cells have to be considered in the sequential case, see Figure 9. Cells outside of the root domain are used to implement physical boundary conditions. Ghost cells overlaid by a grid on level 𝑙 have a unique interior cell analogue and are set by copying the data value from the grid, where the interior cell is contained (synchronization). On the root level, no further boundary conditions need to be considered, but for 𝑙>0 also internal boundaries can occur. They are set by a conservative time-space interpolation from two previously calculated time steps of level 𝑙1.

Beside a general data tree that stores the topology of the hierarchy, the AMR method requires at most two regular arrays assigned to each subgrid. They contain the discrete vector of state for the actual and updated time step. The regularity of the data allows high performance on vector and superscalar processors and cache optimizations. Small data arrays are effectively avoided by leaving coarse level data structures untouched when higher-level grids are created. Values of cells covered by finer subgrids are overwritten by averaged fine grid values subsequently. This operation leads to a modification of the numerical stencil on the coarse mesh and requires a special flux correction in cells abutting a fine grid. The correction replaces the coarse grid flux along the fine grid boundary by a sum of fine fluxes and ensures the discrete conservation property of the hierarchical method at least for purely Cartesian problems without embedded boundaries. See [12] or [3] for details.

5.2. Parallelization

Up to now, various reliable implementations of the AMR method for single processor computers have been developed, see for instance [34] or [35]. Even the usage of parallel computers with shared memory is straight-forward, because a time-explicit scheme allows the parallel calculation of the gridwise numerical update, cf. [33]. However, the question for an efficient parallelization strategy becomes more delicate for distributed memory architectures since on such machines the costs for communication can not be neglected. Due to the technical difficulties in implementing dynamical adaptive methods in distributed memory environments, only few parallelization strategies have been considered in practice yet, cf. [36] or [37].

In the AMROC framework, we follow a rigorous domain decomposition approach and partition the AMR hierarchy from the root level on. The key idea is that all higher level domains are required to follow this “floor-plan”. A careful analysis of the AMR algorithm uncovers that the only parallel operations under this paradigm are ghost cell synchronization, redistribution of the AMR hierarchy, and the application of the previously mentioned flux correction terms. Interpolation and averaging, but in particular the calculation of the flux corrections, remain strictly local, cf. [3, 11]. Currently, we employ a generalization of Hilbert's space-filling curve [38] to derive load-balanced root level distributions at runtime. The entire AMR hierarchy is considered by projecting the accumulated work from higher levels onto the root level cells. Although rigorous domain decomposition does not lead to a perfect balance of workload on single levels, good scale-up is usually achieved for moderate CPU counts. Figure 10 shows a representative scalability test for a three-dimensional spherical shock wave problem for the numerically inexpensive Euler equations for a single polytropic gas without chemical reaction. Roe's approximate Riemann solver within the multidimensional Wave Propagation Method [39] is used as efficient single-grid scheme. The test was run on a Pentium-4-2.4 GHz dual processor cluster connected with Quadrics Interconnect. The base grid has 323 cells and two additional levels with refinement factors 2 and 4. The adaptive calculation uses ~7.0 M cells in each time step instead of 16.8 M cells in the uniform case. The calculation on 256 CPUs employs between ~1,500 and ~1,700 subgrids on each level. Displayed are the average costs for each root level time step, which involves two time steps on the middle level and eight on the highest. All components of the dynamically adaptive algorithm, especially regridding and parallel redistribution are activated to obtain realistic results. Although we utilize a single-grid update routine in Fortran 77 in a C++ framework with full compiler optimization, the fraction of the time spent in this Fortran routine is 90.5% on four and still 74.9% on 16 CPUs. Hence, Figure 10 shows a satisfying scale-up for at least up to 64 CPUs.

AMROC's hierarchical data structures are derived from the DAGH (Distributive Adaptive Grid Hierarchies) package [37] and are implemented completely in C++. A redesign of major parts of DAGH was necessary to allow the AMR algorithm as it was described in the previous section. Currently, AMROC consists of approximately 46,000 lines of code in C++ and approximately 6,000 lines for visualization and data conversion.

6. Numerical Results

Self-sustaining CJ detonations in low-pressure hydrogen-oxygen mixtures with high-argon dilution are ideal candidates for fundamental detonation structure simulations, because they are known to produce very regular detonation cell patterns [2]. In the following, we will simulate such mixtures with the objective of resolving the Mach reflection structures around triple points with high accuracy. Although the detailed hydrodynamic shock wave structure of periodically oscillating detonations in rectangular two-dimensional channels has been fairly well established by now [3, 40, 41], numerous open questions remain for non-rectangular geometries and three space dimensions.

6.1. Two-Dimensional Cellular Structure

We start by reproducing the detonation structure result of Hu et al. [41], utilizing the dynamically adaptive computing approach presented in the previous sections. As initial conditions, the analytical solution according to the one-dimensional ZND theory for the CJ detonation of Figure 1 is extended to multiple space dimensions, and transverse disturbances are initiated by placing a small rectangular unreacted pocket behind the detonation front, cf. [40] or [3]. Throughout this paper, only the hydrogen-oxygen reaction mechanism extracted from the larger hydrocarbon mechanism assembled by Westbrook [42] is used. The mechanism consists of 34 elementary reactions and considers the 9 species H, O, OH, H2, O2, H2O, HO2, H2O2, and Ar. According to the ZND theory, the induction length, the distance between leading shock and head of reaction zone in one space dimension, is 𝑙ig1.404 mm for this mechanism in the above configuration. The detonation velocity is ~1626.9 m/s.

Our shown dynamically adaptive computation has effectively four times higher resolution (44.8 Pts/𝑙ig) than earlier calculations [40] and is similarly resolved as the result by Hu et al. [41] on a uniform mesh. Unfortunately, no technical details are reported for this simulation. In our case, the calculation was run on a small Beowulf-cluster of 7 Pentium 3-850 MHz-CPUs connected with a 1 Gb-Myrinet network and required 2150 h CPU time (12.8 days wall time).

We simulate a physical time of 800 μs. After an initial period of ~200 μs, the simulation starts to exhibit very regular detonation cells with oscillation period ~32 μs. Exploiting the spatial periodicity, we concentrate on the high-resolution simulation of a single cell. The calculation is done in a frame of reference attached to the detonation and requires just the computational domain 10 cm × 3 cm. The adaptive run uses a root level grid of 200×40 cells and two refinement levels with 𝑟1,2=4, where ~3550 root level time steps (automatically chosen based on 𝐶CFL0.95) with average Δ𝑡00.225μs are taken (All computations herein used time steps automatically chosen based on the largest Courant-Friedrichs-Lewy (CFL) number 𝐶CFL encountered throughout the last root level time step. In evaluating 𝐶CFL, our practical implementation considers all hierarchical time steps, yet only adjustment of Δ𝑡0 is permitted. This choice avoids having to deal with additional steps during the recursive integration procedure, cf. Figure 8. If the maximal CFL number exceeds a given threshold, the software has the capability of repeating the entire root level time step.) A physically motivated combination of scaled gradients and heuristically estimated relative errors is applied as adaptation criteria. These criteria are used similarly to generate the refinement of all computations throughout this paper (see [3, 7] for details). Note that in order to simulate the shock wave patterns around the triple points, it suffices to concentrate the computational expense on the detonation front. The flow behind a detonation is basically supersonic, and resolution can be dropped downstream without affecting the approximation of the front structures. Two typical snapshots with the corresponding refinement close to the end of the simulation, when the oscillation is perfectly regular, are displayed in the Figures 11 and 12.

A breakdown of the computational time spent in the computationally most expensive parts for this simulation is given in Table 1. The numerical integration is split into the two main steps of the fractional step method: the fluid dynamic update and the integration of the reactive source term. A ratio between these portions of 1.1 is in good agreement with timings reported by Oran and coworkers [40] and underlines that a reasonable load-balancing is achieved without considering the number of subcycling time steps in the source term integration in the workload estimation explicitly. The third important portion in Table 1 is the time for setting the ghost cell values at subgrid boundaries, which is dominated (~90%) by parallel synchronization.

6.2. Three-Dimensional Cellular Structure

On 24 Athlon-1.4 GHz double-processor nodes (2 Gb-Myrinet interconnect), our approach allowed a sufficiently resolved computation of the three-dimensional cellular structure of a hydrogen-oxygen detonation. The maximal effective resolution of this calculation is 16.8 Pts/𝑙ig and the run required 3800 h CPU time (3.3 days wall time). Our adaptive results are in perfect agreement with the calculations by Tsuboi et al. [43] for the same configuration obtained on a uniform mesh on a superscalar vector machine. A snapshot of the regular two-dimensional solution of the preceding section is used to initialize a three-dimensional oscillation in the 𝑥2-direction and disturbed with an unreacted pocket in the orthogonal direction. We use a computational domain of the size 7 cm × 1.5 cm × 3 cm that exploits the symmetry of the initial data, but allows the development of a full detonation cell in the 𝑥3-direction. Again, we simulate a physical time of 800 μs. The AMROC computation uses a two-level refinement with 𝑟1=2 and 𝑟2=3 on a base grid of 140×12×24 cells (cf. Figure 13) and utilizes between 1.3 M and 1.5 M cells, instead of 8.7 M cells like a uniformly refined grid. 3431 root level time steps with average size Δ𝑡00.233μs are taken. Note how Table 2 clearly reflects the increased expense in solving the hydrodynamic equations in three space dimensions.

After a settling time of about 20 periods (𝑡600μs), a regular cellular oscillation with identical strength in 𝑥2- and 𝑥3-direction can be observed. In both transverse directions, the strong two-dimensional oscillations is present and forces the creation of rectangular detonation cells of 3 cm width. Like in the two-dimensional case, the oscillation period is ~32 μs. The transverse waves form triple point lines in three space dimensions. During a complete detonation cell, the four lines remain mostly parallel to the boundary and hardly disturb each other, cf. Figure 14. The characteristic triple point pattern can therefore be observed in all planes perpendicular to a triple point line. Each triple point line is driven forward by a Mach stem line into an incident shock region. A region that is a Mach stem in one direction can be either a Mach stem or an incident shock in the orthogonal direction. The same is true for the incident shock, and we consequently have three different types of rectangular shock regions at the detonation front: Mach stem-Mach stem (MM), Incident-Incident (II), and mixed type regions (MI) [44]. An MM region is bounded only by Mach stems and expands in the 𝑥2- and 𝑥3-direction. An II rectangle is formed only by incident shocks and shrinks in both directions. The MI region is of mixed type and expands in one direction, while it shrinks in the other. Figure 15 displays the different rectangular shock regions for a specific snapshot. It depends on the initial flow field whether the upper or the lower situation of Figure 14 appears in the simulation. Unlike Williams et al. [44], who presented a similar calculation for an overdriven detonation with simplified one-step reaction model, we notice no phase-shift between both transverse directions. In all our computations, only this regular three-dimensional mode or a purely two-dimensional mode with triple point lines just in 𝑥2- or 𝑥3-direction did occur. The three-dimensional mode of propagation, called “rectangular-mode-in-phase”, has also been found in experiments with hydrogen-oxygen CJ detonations [45].

In order to verify the independence of the numerical solution from the computational grid, we carry out a similar three-dimensional simulation on a grid rotated by 45 degree and utilize periodic boundary conditions in the 𝑥2- and 𝑥3-directions. The domain has the size 7 cm × 32 cm × 32 cm and the base grid the resolution 140×34×34. With a two-level refinement with factors 𝑟1=2 and 𝑟2=3 and time step parameters as before, the simulation has the same effective resolution as the previous one and is consequently about four times more expensive. The lower graph of Figure 16 shows a typical snapshot that corresponds well to the situation shown in the upper graph, which illustrates that the “diagonal-mode-in-phase”, sometimes observed in experiments, is only a different, boundary-condition-dependent, appearance of the same three-dimensional mode of oscillation. Diagonal and rectangular oscillation are simultaneously present in the periodic situation depicted in Figure 15, and it depends only on the size and orientation of the observation window which oscillation type is actually seen. The graphs of Figure 16 show schlieren plots of the mass fraction for 𝑌OH overlaid by a blue isosurface of the density, which visualizes the induction length 𝑙ig. The induction length shows significantly larger variations than in the two-dimensional case. The oscillation period is identical in two- and three space dimensions which underlines that the basic two-dimensional instability is exactly preserved in three space dimensions, but that its manifestation in the hydrodynamic flow field is different.

6.3. Structure of Diffracting Detonations

Experiments have shown that the behavior of planar CJ detonations propagating out of tubes into unconfinement is determined mainly by the width of the tube. For square tubes, the critical tube width has been found to be of the order of 10-times the cell width, that is, 10 λ [17]. For widths significantly below 10 λ, the process of shock wave diffraction causes a pressure decrease at the head of the detonation wave below the limit of detonability across the entire tube width. Hydrodynamic shock and reaction front decouple and the detonation decays to a shock-induced flame. This observation is independent of a particular mixture. While the successful transmission of the detonation is hardly disturbed for tubes widths ≫10 λ, a backward-facing wave reignites the detonation in the partially decoupled region for widths of ~10 λ and creates considerable vortices.

We are interested in the decoupling of shock and reaction and also in the reignition phenomenon. Therefore, we simulate the two-dimensional diffraction of the H2 : O2 : Ar CJ detonation for the tube widths 8 λ and 10 λ. A periodically reproduced snapshot of the regular oscillating detonation of Section 6.1 propagating into unreacted gas at rest is used as initial condition. This is a reasonable idealization for the flow situation in real detonation tubes directly before the experimental setup. The symmetry of the problem is exploited by simulating just one half.

The adaptive simulations utilize a base grid of 508×288 cells and use four levels of refinement with 𝑟1,2,3=2, 𝑟4=4. The calculations correspond to a uniform computation with ~150 M cells and have an effective resolution of 25.5 Pts/𝑙ig in the 𝑥1-direction (with respect to the initial detonation). Both runs are stopped at 𝑡end=240μs after 730 root level time steps, with average step size Δ𝑡00.329μs, when the flow situations of interest are clearly established. The complete decoupling of shock wave and flame front for 8 λ is visible in the upper graph of Figure 17, while the lower graph clearly exhibits the re-ignition wave. It is interesting to note that the re-ignition wave itself is a detonation. The triple point tracks for 10 λ (not shown) uncover that it has developed out of the transverse wave of an initial triple point.

The enormous efficiency of the refinement is visualized in the upper graph of Figure 18 for the setup with tube width 10 λ. At 𝑡end, the calculation shown in Figure 18 uses ~3.0 M cells on all levels, where ~2.4 M cells are inside one of the 2,479 grids of the highest level. The bottom graph of Figure 18 displays the complexity of the domain decomposition (indicated by color) to the 48 processors. These simulations were also run on the same parallel computer system as the three-dimensional computations and required ~4500 h CPU (3.9 days wall time).

6.4. Propagation through a Pipe Bend

In order to demonstrate the enormous potential of the entire approach even for non-Cartesian problems, we finally discuss a case study that combines highly efficient dynamic mesh adaptation with the embedded boundary method sketched in Section 4.1. A two-dimensional regularly oscillating detonation similar to Section 6.3 but with initial pressure 𝑝0=10.0 kPa and 𝜆=1.6 cm is placed into a pipe of width 5 λ. The pipe bends at different angles and has the inner radius 9.375 λ. First, we consider the bending angle 𝜑=60°. When the detonation propagates through the bend, it gets compressed and consequently overdriven near the outer wall, but a continuous shock wave diffraction occurs near the inner wall. Similar to the simulation of Section 6.3, the diffraction causes a pressure decrease below the limit of detonability that leads to a continuous decoupling of shock and reaction front. This effect is visible particularly in the second graph of Figure 21 at 𝑡=150μs. The detonation exits the bend before the decay to a flame occurs across the entire tube width. A reignition wave arises from the successfully transmitted region and reinitiates the detonation in the decoupled area, see third graph of Figure 21 and first schlieren graph of Figure 19 for 𝑡=200μs. It propagates in the direction normal to the pipe middle axis and causes a strong shock wave reflection as it hits the inner wall, compare last graph of Figure 21 and second schlieren graphs of Figure 19 for 𝑡=250μs. After the partial reignition event, the front structure converges back to the initial oscillation with five regular cells across the tube width. Figure 20 shows the simulated triple point trajectories through the course of the entire simulated time of 𝑡end=400μs and the convergence back to the equilibrium state can be inferred.

This simulation uses a base grid of 1200×992 cells, four levels of refinement with 𝑟1,2,3=2, 𝑟4=4, and has an effective resolution of 67.6 Pts/𝑙ig. Approximately 7.1 M to 3.4 M cells on all and 4.8 M to 1.8 M cells on the highest level are used instead of ~1,219 M in the uniform case. The number of subgrids on the highest level varies between ~2,600 and ~800. The calculation was run on 64 dual-processor nodes of the ASC Linux cluster at Lawrence Livermore National Laboratory and required nevertheless ~70,000 h CPU (~23 days wall time), due to the necessary high local resolution in space and time. Approximately 4000 time steps were taken on the root level, yielding an average root level time step of Δ𝑡00.1μs, but note that due to automatic time step steering based on 𝐶CFL0.6 the used root level time steps vary between ~0.085 μs and ~0.120 μs. The extraordinary high efficiency of our mesh refinement approach in capturing only the essential features near the detonation front is illustrated by the right graphs of Figure 19. The effective resolution at the detonation front is so high that the simulation succeeds in capturing even secondary triple points throughout the entire run. Figure 22 shows a zoom into a triple point as it enters the bend and all features are fully resolved.

6.5. Triple Point Analysis

The high resolution of the last simulation admits a considerable refinement of the triple point pattern introduced in Section 2. We start by analyzing the shock reflection pattern around a triple point for the perfectly regularly oscillating detonation before entering the bend, which obviously also corresponds to the flow situation of Section 6.1. In order to carry out a quantitative Mach reflection pattern analysis, it is necessary to map the velocity field of the simulation into a frame of reference attached to the triple point. However, the reliable estimation of the triple point speed 𝑣 from a single output time step is nonapparent. When evaluating 𝑣, we take advantage of the fact that the triple point is formed at the tip, where the Mach stem intersects the incident region, and that the oblique shock relations [15] between two points in regions A and B (cf. Figure 24) close to the triple point must hold true. We require only the two relations 𝜌𝐴𝑢𝐴𝜙sin𝐵=𝜌𝐵𝑢𝐵𝜙sin𝐵𝜃𝐵𝑝,(6)𝐴+𝜌𝐴𝑢2𝐴sin2𝜙𝐵=𝑝𝐵+𝜌𝐵𝑢2𝐵sin2𝜙𝐵𝜃𝐵.(7) Inserting (6) into (7) allows the elimination of 𝑢𝐵sin(𝜙𝐵𝜃𝐵), which yields 𝑢𝐴=1sin𝜙𝐵𝜌𝐵𝑝𝐵𝑝𝐴𝜌𝐴𝜌𝐵𝜌𝐴.(8) As the gas is initially at rest, the triple point velocity is 𝑣=𝑢𝑎, and 𝜙𝐵, the angle of inflow, is given as the angle between Mach stem front and the triple point trajectory, which can be measured from visualizations such as those shown in Figures 25 and 26.

Figure 23 zooms into the situation shortly before the collision of two triple points. The enlargement shows clearly that the shock wave pattern around each triple point is of double Mach reflection (DMR) (aka “strong”) type. As depicted in Figure 24, the essential flow features around triple point T are inflow (A), region B behind the Mach stem M, transverse wave region (C), and region D downstream of the incident shock I. Regions B and C are separated by the slip line S, a contact discontinuity. Regions C and D are separated by the reflected shock R. The characteristic for the DMR pattern is a high supersonic velocity in region C that leads to the formation of a further shock R′ creating a secondary triple point T′ on the transverse wave. The very weak slip line S′, emanating from T′ between regions E and F, can hardly be inferred. Note that the reflected shock R is the incident shock I′ for the secondary triple point T′.

Since the triple point is far ahead of the reaction region, changes in mixture are neglected in evaluating the Mach number 𝑀 in the triple point pattern. As it can be expected in a DMR [15], 𝑀𝐶 is clearly greater than 1. A further important quantity for triple point structures is the strength 𝑆 of the transverse wave [19] that is defined as𝑝𝑆=𝐶𝑝𝐷𝑝𝐷.(9)

For the present computation, 𝑆 is rather constant throughout one regular detonation cell and varies marginally around 0.65.

6.6. Triple Point Structures under Transient Conditions

The diffraction and compression of the detonation wave at the bend leads to changes of the transverse wave strength 𝑆. While the state 𝐶 behind the reflected shock R remains initially unchanged, the change in geometry alters the incident state 𝐷. Near the inner bend wall (upper boundary in Figure 25), the incident pressure 𝑝𝐷 drops, leading to a considerable increase in 𝑆. The resulting triple point structure is of strengthened DMR type with larger detonation cells (cf. upper graph of Figure 25) and increased Mach number 𝑀𝐶. In the compression region near the outer wall (lower boundary in Figure 25), however, 𝑝𝐷 increases drastically leading to a considerable decrease of 𝑆. As a consequence, the flow in region C decelerates and the secondary shock R′ between C and E is no longer necessary for a stable configuration. The triple point now exhibits a transitional Mach reflection (TMR) or “weak” pattern. Characteristic for the TMR structure is that the flow in region C is just barely supersonic with respect to T [15]. TMR structures are also found in the diffraction region when a sufficiently large angle 𝜑 causes partial detonation failure.

The lower graph of Figure 25 depicts the interesting case of triple point reinitiation behind the bend for 𝜑=15°. The new triple point in the center of the graph is formed as a weak structure (TMR) but strengthens into a DMR after its first collision. Figure 26 visualizes the quenching of triple points at the outer bend wall that occurs for all larger values of 𝜑. As can be inferred from the lower graph of Figure 26, failing triple points with vanishing trajectories seem to exhibit a weak structure.

6.7. Transition Criteria of Mach Reflection Phenomena

The two-dimensional oblique jump conditions [15] are satisfied across each discontinuity depicted in the sketch of Figure 24. We solve the resulting nonlinear system of equations numerically under the assumption that the composition of the H2 : O2 : Ar mixture is invariant, but note that the temperature dependency of the material properties is fully considered. It is well known (cf. [15]) that changes in Mach reflection type occur when the flow in certain regions of the triple point structure becomes supersonic in the frame of reference of the triple point T or T′. Any Mach reflection requires a supersonic flow in region B measured in the frame of reference of the primary triple point, that is, 𝑀TB>1, while for 𝑀TB<1 a regular reflection (RR) occurs. A single Mach reflection (SMR) pattern encompassing just the four states A to D is found for 𝑀TC<1; for 𝑀TC>1 a transitional or a double Mach reflection pattern occurs. In a DMR, the flow in C is also supersonic in the frame of reference of the secondary triple point T′ (𝑀TC>1), however, it is subsonic with 𝑀TC<1, in a TMR. Consequently, the secondary triple point T′ is not fully developed in a TMR and the discontinuities R′ and S′ are absent. For given inflow velocity 𝑣, the transition points 𝑀TB/C=1 are unique and can be determined rather easily by solving the corresponding set of oblique shock relations. The evaluation of the TMR/DMR transition point 𝑀TC=1 is more involved as knowledge of the velocity vector 𝑎 of T′ relative to T is required. The oblique shock relations between C and D readily yield 𝑎𝑛0; however, the tangential velocity 𝑎𝑡 of T′along R is a free parameter that needs to be specified separately.

Figure 27 visualizes the transition boundaries of the different shock wave reflection phenomena evaluated numerically for varying inflow Mach number. Since the estimation of the exact minimally possible value of 𝑎𝑡 in a DMR is a considerable problem in itself, the TMR/DMR transition line is evaluated for the fixed value 𝑎𝑡=100 m/s, which puts a reasonable upper transition boundary on the TMR region (dotted in Figure 27). Specially geared to detonations, results are displayed for the relative transverse wave strength 𝑆. We use the simulation results discussed in Sections 6.5 and 6.6 as a first verification for the found quantitative transition boundaries. Throughout an unperturbed detonation cell only the DMR structure occurs, which is confirmed by plotting (with stars) the point at the end of a cell (cf. Figure 23) and an additional one immediately after triple point collision. The open circles mark the triple point reinitiation visualized in the lower graph of Figure 25, and the transition boundaries clearly confirm the transition from TMR to DMR. The failing triple point, shown in the lower graph of Figure 26, is marked with closed triangles indicating that the SMR pattern seems unstable for a detonation substructure and occurs only immediately before the triple point disappears.

7. Conclusions

We have described an efficient solution strategy for the numerical simulation of gaseous detonations with detailed chemical reaction. Beside the application of the time-operator splitting technique and the construction of a robust high-resolution shock capturing scheme, the key to the high efficiency of the presented simulations is the application of a block-structured adaptive mesh refinement method and its parallelized implementation in our software framework AMROC. AMROC provides the required high local resolution dynamically and follows a parallelization strategy tailored especially for the emerging generation of distributed memory architectures. An embedded boundary method utilizing internal ghost cells extends the framework effectively to non-Cartesian problems.

Two- and three-dimensional detonation structure simulations for hydrogen-oxygen-argon mixtures have been discussed in detail, where all temporal and spatial scales relevant for the complex process of detonation propagation and cellular structure formation from subscale Mach reflection patterns were successfully resolved. All presented simulation results have been obtained on Linux-Beowulf-clusters of moderate size. The reliability of the mesh adaptation method has been verified for periodically oscillating detonation structure simulations in a frame of reference attached to the detonation front. While the savings from dynamic mesh adaptation are necessarily modest for these types of simulations, they have been demonstrated to be enormous for detonation structure simulations in more complex geometries that have to be carried out in an Eulerian frame of reference. Depending on the required maximal resolution, increasingly larger savings in mesh size can be achieved. The large-scale two-dimensional structure simulations in pipe bends exhibit reductions of the finest mesh size of at least a factor of 250 and up to 680. These computations have an effective resolution that captures even secondary triple points reliably throughout the whole simulation and allows the unambiguous classification of the Mach reflection structures around triple points. Under two-dimensional transient conditions, transitional (aka “weak”) and double Mach reflection patterns (aka “strong” structures) have been observed.

Nonreactive but thermally perfect shock wave reflection theory was used successfully to predict the transition boundaries between Mach reflection patterns depending on the inflow velocity. Exemplary transient triple points structures from the simulations are compared with the theoretical results and good agreement between observed and predicted Mach reflection type is found. Future work in understanding the occurrence of Mach reflection patterns in self-sustained detonation waves will concentrate on deriving a physical estimate for the smallest possible relative tangential velocity 𝑎𝑡 between primary and secondary triple point in a double Mach reflection in real gases to complete the theoretical model and on comparing the substructures of the regularly oscillating situations in two and three space dimensions. First detailed investigations [8] indicate that the fully three-dimensional case shows lower values for 𝑆 and therefore a transitional rather than a double Mach reflection as basic propagation pattern. This result also suggests that fully predictive detonation simulations in general will have to be carried out in three space dimensions. Since it was demonstrated that our computational approach permits the accurate numerical simulation of detonation structures in technically relevant two-dimensional devices and the study of idealized three-dimensional configurations on current capacity computing systems, this ultimate goal might already be attainable on available petaFLOPS supercomputers.

As an outlook on computational detonation research in general, we finally want to mention that the quantification of the influence of diffusive processes on the detailed structure of the primary slip line and its extension (denoted F and G in Figure 3) in Mach reflection patterns in detonation waves has recently received some attention. This topic is of particular interest for mixtures that lead to very irregularly oscillating cellular structures, for example, hydrocarbon fuels, cf. [46, 47]. The slip line is formed by a shear layer flow that is unstable to perturbations (which is commonly known as a Kelvin-Helmholtz instability). A fully converged solution, that predicts the physically accurate behavior of this shear layer at all subscales, cannot be obtained with an inviscid model. For the resolutions and second-order scheme used in this paper, stability and structure of the triple point slip lines are still clearly determined by the numerical diffusion introduced by the approximation; yet under increasingly finer grid resolution a point is invariably reached from which on the correct physical diffusion needs to be considered. Instead of the Euler equations (1), the fully reactive Navier-Stokes equations (cf. [16]) need to be employed, which can be accomplished technically rather easily for detonation wave simulation by supplementing the presented numerical method with a conservative approximation of the physical diffusion fluxes and incorporation of an additional stability condition. Utilizing such an approach in AMROC in combination with an upwind scheme of very high order, we have recently succeeded in performing a large-scale two-dimensional adaptive computation that resolves all physical length scales along the slip line of an isolated transitional Mach reflection pattern that is created by the reflection of a planar detonation wave at an inclined wedge [48]. Note that this pseudo-DNS (A true direct numerical simulation (DNS) would require accurate resolution of all physical length scales involved. This is presently not feasible for the thickness of the shock waves, which is of the order of three to five mean free path lengths.) result has been obtained for the simplified case of a single exothermic reaction between two calorically perfect gases and not for a realistic chemical reaction mechanism, as used throughout this paper.

Acknowledgments

This work has been sponsored by the Office of Advanced Scientific Computing Research, U.S. Department of Energy (DOE) and was performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC under Contract no. DE-AC05-00OR22725, and by the DFG high priority research program “Analysis and Numerics of Conservation Laws”, grant Ba 840/3-3. The large-scale computations of Sections 6.4 to 6.7 have been carried out while the author was at the California Institute of Technology and was supported by the ASC program of the Department of Energy under subcontract no. B341492 of DOE contract W-7405-ENG-48.