Meso-resolved simulations of shock-to-detonation transition in nitromethane with air-filled cavities

TL;DR

This study uses meso-resolved simulations to explore how cavities in nitromethane influence shock-to-detonation transition, revealing different initiation behaviors and the impact of cavity distribution on detonation sensitivity.

Contribution

The paper introduces a GPU-enabled, meso-resolved simulation approach for modeling SDT in nitromethane with explicit cavity resolution without empirical models, highlighting the effects of cavity distribution.

Findings

Random cavities accelerate SDT more than regular arrays.

Heterogeneous mixtures show shorter overtake times at lower shock pressures.

Simulation results align with experimental observations of SDT behaviors.

Abstract

Two-dimensional, meso-resolved numerical simulations are performed to investigate the complete shock-to-detonation transition (SDT) process in a mixture of liquid nitromethane (NM) and air-filled, circular cavities. The shock-induced initiation behaviors resulting from the cases with neat NM, NM with an array of regularly spaced cavities, and NM with randomly distributed cavities are examined. For the case with randomly distributed cavities, hundreds of cavities are explicitly resolved in the simulations using a diffuse-interface approach to treat two immiscible fluids and GPU-enabled parallel computing. Without invoking any empirically calibrated, phenomenological models, the reaction rate in the simulations is governed by Arrhenius kinetics. For the cases with neat NM, the resulting SDT process features a superdetonation that evolves from a thermal explosion after a delay following the passage of the incident shock wave and eventually catches up with the leading shock front. For the cases wherein mesoscale heterogeneities are explicitly considered, a gradual SDT process is captured. These two distinct initiation behaviors for neat NM and heterogeneous NM mixtures agree with experimental findings. Via examining the global reaction rate of the mixture, a unique time scale characterizing the SDT process, i.e., the overtake time, is measured for each simulation. For an input shock pressure less than approximately $9.4~\mathrm{GPa}$, the overtake time resulting from a heterogeneous mixture is shorter than that for neat NM. This sensitizing effect is more pronounced for lower input shock pressures. A random distribution of cavities is found to be more effective in enhancing the SDT process than a regular array of cavities. Statistical analysis on the meso-resolved simulation data provides more insights into the mechanism of energy release underlying the SDT process.

Tables

Table 1. Table 1: Parameters for the equation of state of air

Parameter	Value (unit)
${\gamma }_1$	$1.4$ (-)
$c_{v,1}$	$718$ ($J$ ${\mathrm{kg}}^{-1}{m}^{-3}$)

Table 2. Table 2: Parameters for the equation of state of liquid NM

Parameter	Value (unit)
${\Gamma }_{0,2}$	$1.19$ (-)
$𝒜$	$0.819$ (GPa)
$ℬ$	$1.51$ (GPa)
$𝒞$	$4.53$ (-)
$𝒟$	$1.42$ (-)
${\rho }_{0,2}$	$1134$ ($\mathrm{kg}$ $m^{-3}$)
$c_{v,2}$	$1714$ ($J$ ${\mathrm{kg}}^{-1}{m}^{-3}$)

Table 3. Table 3: Parameters for the reaction law of liquid NM

Parameter	Value (unit)
$T_a$	$11350$ (K)
$C$	$2.6\times {10}^9$ ($s^{-1}$)
$Q$	$4.46\times {10}^6$ ($J$ ${\mathrm{kg}}^{-1}$)

Table 4. Table 4: Numerical convergence tests for the case with an array of regularly spaced cavities of ϕ = 8 % italic-ϕ percent 8 \phi=8\% and d c = 100 ​ μ ​ m subscript 𝑑 c 100 𝜇 m d_{\mathrm{c}}=100~{}\mu\mathrm{m} subjected to an input shock pressure of 7.81 ​ GPa 7.81 GPa 7.81~{}\mathrm{GPa} .

$dx$	$4\mu m$	$2\mu m$	$1\mu m$	$0.5\mu m$
Overtake time	$3.77\mu s$	$3.05\mu s$	$2.88\mu s$	$2.88\mu s$

Table 5. Table 5: Numerical convergence tests for the case with randomly distributed cavities of ϕ = 5 % italic-ϕ percent 5 \phi=5\% and d c = 100 ​ μ ​ m subscript 𝑑 c 100 𝜇 m d_{\mathrm{c}}=100~{}\mu\mathrm{m} subjected to an input shock pressure of 7.81 ​ GPa 7.81 GPa 7.81~{}\mathrm{GPa} .

$dx$	$4\mu m$	$2\mu m$	$1\mu m$
Overtake time	$3.68\mu s$	$3.21\mu s$	$2.99\mu s$

Equations

$$\begin{array}{r}\frac{\partial z_1{\rho }_1}{\partial t}+\nabla \cdot \left(z_1{\rho }_1\mathbf{u}\right)=0\\ \frac{\partial z_2{\rho }_2}{\partial t}+\nabla \cdot \left(z_2{\rho }_2\mathbf{u}\right)=0\\ \frac{\partial }{\partial t}\left(\rho \mathbf{u}\right)+\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}\right)+\nabla p=0\\ \frac{\partial }{\partial t}\left(\rho E\right)+\nabla \cdot \left[\left(\rho E+p\right)\mathbf{u}\right]=-z_2{\rho }_{2}KQ\\ \frac{\partial z_1}{\partial t}+\mathbf{u}\cdot \nabla z_1=0\\ \frac{\partial }{\partial t}\left(z_2{\rho }_2\lambda \right)+\nabla \cdot \left(z_2{\rho }_2\lambda \mathbf{u}\right)=z_2{\rho }_{2}K\end{array}$$

$$p({\rho }_2,e_2)=p_{\mathrm{ref},2}({\rho }_2)+{\rho }_2{\Gamma }_{0,2}\left[e_2-e_{\mathrm{ref},2}({\rho }_2)\right]$$

$$p_{\mathrm{ref},2}({\rho }_2)=\mathcal{A}{\left(\frac{{\rho }_{0,2}}{{\rho }_2}\right)}^{-\mathcal{C}}-\mathcal{B}{\left(\frac{{\rho }_{0,2}}{{\rho }_2}\right)}^{-\mathcal{D}}$$

$$\begin{array}{rl}{e}_{\mathrm{ref},2}({\rho }_2)=& \frac{-\mathcal{A}}{{\rho }_{0,2}(1-\mathcal{C})}\left[{\left(\frac{{\rho }_{0,2}}{{\rho }_2}\right)}^{1-\mathcal{C}}-1\right]\\ +& \frac{\mathcal{B}}{{\rho }_{0,2}(1-\mathcal{D})}\left[{\left(\frac{{\rho }_{0,2}}{{\rho }_2}\right)}^{1-\mathcal{D}}-1\right]\end{array}$$

$$T_i-T_{\mathrm{ref},i}({\rho }_i)=\frac{p-p_{\mathrm{ref},i}({\rho }_i)}{{\rho }_i{\Gamma }_i({\rho }_i,T_i)c_{\mathrm{v},i}({\rho }_i,T_i)}\mathrm{for}i=1,2$$

$$K=\frac{\partial \lambda }{\partial t}=-\lambda C\exp (-T_{\mathrm{a}}/T_2)$$

$${\delta }_{\mathrm{c}}=\sqrt{\frac{\pi {d_{\mathrm{c}}}^2}{4\varphi }}$$

$$N_{\mathrm{tot}}=\frac{4\varphi LW}{\pi {d_{\mathrm{c}}}^2}$$

$$\bar{p}(x,t)=\frac{1}{W}{\int }_0^{W}p(x,y,t)\mathrm{d}y$$

$$\bar{\rho }(x,t)=\frac{1}{W}{\int }_0^W\rho (x,y,t)\mathrm{d}y$$

$$\bar{u}(x,t)=\frac{1}{W\bar{\rho }(x,t)}{\int }_0^W\rho (x,y,t)u(x,y,t)\mathrm{d}y$$

$$M_{2,\mathrm{r}}(t)={\int }_0^W{\int }_0^L{\rho }_2{z}_2\lambda \mathrm{d}x\mathrm{d}y$$

$$\left|{\bar{\lambda }}_t\right|=\left|\frac{\mathrm{d}{M}_{2,\mathrm{r}}}{\mathrm{d}t}/M_{2,\mathrm{tot}}\right|$$

$$\dot{q}=-z_2{\rho }_{2}KQ$$

$$10^n\mathrm{GW}/{\mathrm{m}}^3\le \dot{q}<10^{n+\Delta n}\mathrm{GW}/{\mathrm{m}}^3$$

$$f(\dot{q},t)=\frac{m(\dot{q},t)}{M_{2,\mathrm{tot}}\Delta n}$$

$${\rho }_{\mathrm{N}}=\frac{6\varphi }{\pi {d_{\mathrm{c}}}^3}$$

$$u_{\mathrm{s}}(t)=\frac{x_{\mathrm{s}}(t+\Delta \tau /2)-x_{\mathrm{s}}(t-\Delta \tau /2)}{\Delta \tau }$$

$$a_{\mathrm{s}}(t)=\frac{u_{\mathrm{s}}(t+\Delta \tau /2)-u_{\mathrm{s}}(t-\Delta \tau /2)}{\Delta \tau }$$

Full text

I Introduction

Condensed-phase explosives can be significantly sensitized to shock initiation by the introduction of mesoscale heterogeneities. This sensitizing effect is attributed to the formation and growth of localized high-temperature regions, or the so-called “hot spots”, as an incident shock wave interacts with the heterogeneities. Since the concept of hot-spot sensitization was first proposed by Bowden and Yoffe in the 1940s Bowden and Yoffe (1948), a substantial body of research has been carried out to elucidate the underlying mechanisms.

In contrast to most solid and slurry explosives with intrinsic heterogeneities that are highly irregular in shape, size, and spatial distribution, the size and number density ($\rho_{\mathrm{N}}$) of heterogeneities added into a liquid explosive matrix is relatively controllable. Taking advantage of this feature, the use of liquid nitromethane (NM) mixed with micron-sized solid particles or hollow spheres has been explored to experimentally probe the sensitizing effect of mesoscale heterogeneities on shock-to-detonation transition (SDT) in high explosives. Campbell et al. reported in 1961 that the build-up of a detonation in NM with trapped bubbles of gas or pieces of solid is significantly faster than that in neat NM. Campbell, Davis, and Travis (1961) Later, via in-situ measurement of particle velocity, it has been revealed that the SDT process in NM sensitized by the addition of silica beads or glass micro-balloons (GMBs) features a reaction front growing in the vicinity of the incident shock Sheffield, Engelke, and Alcon (1989), in contrast to the SDT in neat NM wherein the reaction front is initiated significantly behind the incident shock, evolves into a superdetonation, and eventually overtakes the incident shock.Chaiken (1960); Campbell, Davis, and Travis (1961) To further interpret these macroscopic measurements, one needs to gain direct insights into the mechanisms of how hot spots are formed, evolve, and collectively contribute to the energy release process. Thus, studies focused on illustrating the detailed wave dynamics arising from the interaction of an incident shock wave with mesoscale heterogeneities have been motivated.

The wave and interface morphologies evolving over the course of shock-induced collapse of cavities have been experimentally observed using high-speed photography. Dear and Field (1988); Dear, Field, and Walton (1988); Field et al. (1992); Field (1992); Bourne and Field (1992, 1999a, 1999b); Bourne and Milne (2003); Swantek, Shukla, and Austin (2010); Swantek and Austin (2010) Wave interactions resulting from the collapse of neighboring cavities have been revealed by Dear and Field and later researchers.Dear and Field (1988); Bourne and Field (1999a, a); Swantek and Austin (2010) To enable visualization, millimeter-sized or larger cavities were considered in these experimental studies while typical GMBs added into liquid NM, for example, are of diameters on the order of $10$-$100$ $\mu\mathrm{m}$. It is rather challenging, although attempted by Bourne and Milne Bourne and Milne (2003), to experimentally resolve the temperature field arising from the collapse on the spatial scale of a single cavity or average spacing between cavities. Therefore, to understand the formation and evolution of hot spots, shock-induced collapse of micron-sized cavities has been otherwise investigated via numerical simulations.

One- and two-dimensional simulations of a cavity subjected to a planar incident shock wave were performed in the 1960s. As only inert scenarios were considered by Evans et al. Evans, Harlow, and Meixner (1962), Mader, for the first time, demonstrated how the formation of hot spots triggers a rapid burnout of the reactive material surrounding a cavity or a solid pellet Mader (1963, 1965). Three-dimensional simulations of shock-induced collapse of an array of $91$ regularly staggered, micron-sized cavities in NM were performed by Mader in the 1980s.Mader (1984) The simulation results showed that hot spots of increasingly greater temperature and size are formed as the leading shock passes across layers of cavities. Recent computational effort has been dedicated to revealing finer detail of the resulting temperature field from shock–cavity interactions.Ball et al. (2000); Swantek, Shukla, and Austin (2010); Swantek and Austin (2010); Hawker and Ventikos (2012); Lauer et al. (2012); Ozlem et al. (2012); Michael and Nikiforakis (2014); Kapila et al. (2015); Betney et al. (2015); Apazidis (2016); Michael and Nikiforakis (2019a, b, c) Two- and three-dimensional simulations of the collapse of a single air-filled cavity in liquid NM have been performed by Michael and Nikiforakis.Michael and Nikiforakis (2019a, b) These authors demonstrated that hot spots are not only formed due to the penetration of a micro-jet along the centerline of the cavity, but also generated behind Mach stems occurring away from the centerline. These simulation results indicate that hot spots are of irregular shapes and undergo complex evolutions. Thus, it is rather hard to measure the size of a hot spot, and it might not be informative to characterize a hot spot by its temperature without assessing the mass of reactive material associated with it. After revealing detail on mesoscales, how can one computationally examine the mechanism of a large number of hot spots affecting the SDT process in an energetic medium? Two strategies prevail nowadays: Meso-informed and meso-resolved simulations.

In meso-informed simulations, an explosive mixture with mesoscale heterogeneities is treated as a homogeneous medium with a phenomenological (or so-called “meso-informed”) model—the Ignition and Growth (IG) model proposed by Lee and Tarver Lee and Tarver (1980) is the most commonly adopted—that governs the rate of energy release for the given local flow and thermodynamic properties. The computational studies using a variety of phenomenological reaction models have recently been reviewed in detail by Handley et al.Handley et al. (2018) Conventional strategies to develop a meso-informed reaction model consists of, first, constructing the ignition and growth rate terms based on some micro-mechanistic models for single-cavity collapse, and then calibrating the values of the governing parameters against empirical data. Novel approaches using machine-trained surrogate models for the IG rates have been recently applied by Sen et al. to simulate SDT processes in HMX (cyclotetramethylene-tetranitramine, also called octogen). Sen et al. (2018) These surrogate models are trained by ensembles of simulation data of single-cavity collapse under various conditions (e.g., incident shock strength, pulse duration, cavity size and shape, etc.). Such simulations, however, are unlikely able to adequately describe the mutual influence among a large number of randomly distributed hot spots. Hence, meso-resolved simulations might be a more pertinent approach to examine the collective effect of hot spots.

In meso-resolved simulations, a statistically significant amount of mesoscale heterogeneities is explicitly considered; instead of invoking a phenomenological model, the reaction rate is determined by a model based on detailed or reduced chemical kinetics, i.e., activated kinetics in an Arrhenius form. With advanced computing technology, to perform meso-resolved simulations of detonation phenomena in high explosives has recently become feasible. Via performing meso-resolved simulations, Rai et al. Rai and Udaykumar (2015) and Kim et al. Kim et al. (2018) demonstrated that the characteristics of mesoscale features determine the macroscopic ignition behavior of HMX and polymer-bonded explosives (PBXs), respectively. These simulations, however, have not yet been extended to time or length scales over which a complete SDT process occurs.

In the current study, two-dimensional, meso-resolved simulations of complete SDT processes in liquid NM mixed with air-filled cavities are performed. The resulting initiation behaviors from the cases with randomly distributed cavities are compared to those from the cases with a regular array of cavities and neat NM. These computations are made amenable through the implementation of a diffuse-interface approach to treat two immiscible fluids Michael and Nikiforakis (2016) and a CUDA-based (CUDA is an acronym for Compute Unified Device Architecture) parallel computing on general-purpose graphic processing units (GPGPUs). The objective of this work is to explore the capability of meso-resolved simulations in capturing the experimentally observed SDT behaviors of GMB-sensitized NM without invoking a calibrated IG reaction rate model. Via performing statistical analysis on the full information of the resulting flow field, a further elucidation of the hot-spot sensitization mechanism is attempted.

This paper is organized as follows. In Sect. II, the governing equations and various scenarios considered in the simulations are stated. Section III describes the herein used numerical methodology in detail. Simulation results and analysis are presented in Sect. IV. The findings based on the results, limitations of the current work, and future development of meso-resolved simulations to further study the SDT phenomena in NM-GMB mixtures are discussed in Sect. V and summarized in Sect. VI. Numerical and statistical convergence studies are reported in Appendices A and B, respectively.

II Problem statement

II.1 Governing equations

The dynamic system at hand is simulated based on the Euler equations that are augmented to treat a reactive, multi-phase flow with two immiscible materials (i.e., liquid NM and air). Michael and Nikiforakis (2016) The use of this mathematical framework has previously been explored by Michael and Nikiforakis to simulate single cavity collapse in liquid NM. Michael and Nikiforakis (2014, 2019a, 2019b, 2019c) The governing equations are formulated as follows,

[TABLE]

The air within cavities is considered as phase $1$ and the liquid NM is considered as phase 2, which are denoted as subscripts “$1$” and “$2$”, respectively. The volume fractions of air and NM are represented by $z_{1}$ and $z_{2}$, respectively, where $z_{1}+z_{2}=1$. The total density $\rho$ of the mixture can thus be calculated as $\rho=z_{1}\rho_{1}+z_{2}\rho_{2}$. The total specific energy is defined as $E=\frac{1}{2}u^{2}+e$, where $e$ is the specific internal energy of the mixture, i.e., $\rho e=z_{1}\rho_{1}e_{1}+z_{2}\rho_{2}e_{2}$. The reaction progress variable $\lambda$ evolves from $1$ to [math], representing the mass fraction of unreacted NM. The time rate of change of $\lambda$ is denoted as $K$. The specific energy release of NM is represented by $Q$. At material interfaces, these two immiscible materials are considered to be in mechanical equilibrium but not thermal equilibrium, i.e., $p_{1}=p_{2}=p$, $\mathbf{u}_{1}=\mathbf{u}_{2}=\mathbf{u}$, but $T_{1}$ and $T_{2}$ might not be equal.

Since the thermal diffusivity of liquid NM $\alpha$ is on the order of $0.1~{}\mathrm{\mu m}^{2}/\mathrm{\mu s}$, the characteristic time scale of heat diffusion over a cavity diameter of $100~{}\mathrm{\mu m}$, i.e., $\tau_{\mathrm{d}}={d_{\mathrm{c}}}^{2}/\alpha$, is on the order of nearly $1\times 10^{5}\mathrm{\mu s}$. The time scale of the SDT process is expected to be on the order of $\mathrm{\mu s}$ for the range of input shock pressure ($\sim 10~{}\mathrm{GPa}$) in this study. It is therefore justifiable in this model to neglect heat diffusion, which is too slow to significantly affect the hot-spot-triggered combustion. Due to the same reason, heat diffusion was also neglected by MenikoffMenikoff (2011) in modeling hot spot formation in NM mixed with solid beads. The effect of viscous heating in this system consisting of liquid- and gas-phase materials is in general less significant than the viscoplastic heating effect due to intergranular friction in polycrystalline solid explosives. Also, viscous heating is speculated to significantly increase the hot spot temperature for low incident shock pressures.Field (1992) In this study, since relatively high input shock pressures ($>7~{}\mathrm{GPa}$) are considered, material viscosity is neglected in this model in order to carry out more focused parametric studies.

II.2 Equations of state (EoS)

In this work, the unreacted NM and its products are described by the same Cochran-Chan EoS Cochran and Chan (1979); Saurel, Petitpas, and Berry (2009), which is expressed as follows,

[TABLE]

where the reference pressure $p_{\mathrm{ref,2}}$ is given by

[TABLE]

the reference energy $e_{\mathrm{ref,2}}$ is given by

[TABLE]

and $\Gamma_{0,2}$ is the Grüneisen coefficient corresponding to the initial state of NM at $\rho_{2}=\rho_{0,2}$. This EoS has been broadly used in the literature. Saurel, Petitpas, and Berry (2009); Shukla, Pantano, and Freund (2010); Genetier, Osmont, and Baudin (2014); Michael and Nikiforakis (2019a, b, c) Air inside the cavities is governed by the ideal gas law, i.e., $e_{1}=p/(\gamma_{1}-1)\rho_{1}$, where $\gamma_{1}$ is the ratio of specific heat capacities of air. The EoS of the ideal gas law can be equivalently expressed in the Mie-Grüneisen form with $p_{\mathrm{ref,1}}=0$ and $e_{\mathrm{ref,1}}=0$. The values of $\gamma_{1}$ and the constant-volume specific heat capacity $c_{\mathrm{v,1}}$ of air are given in Table 1.

The temperature of each material needs to be calculated separately according to its equation of state. A general expression for calculating the temperature of each material based on the first-order Taylor expansion from the reference curve should be as follows

[TABLE]

where the constant-volume specific heat capacity $c_{\mathrm{v},i}$ and Grüneisen coefficient $\Gamma_{i}$ are functions of both density and temperature.Winey et al. (2000) The Cochran-Chan EoS used for liquid NM is based on the reference curves of isotherms, thus, $T_{\mathrm{ref},2}(\rho_{2})=T_{\mathrm{ref},2}$.Cochran and Chan (1979) The parameters, $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$, and $\mathcal{D}$, in the functions for reference pressure and energy are calibrated using the approach proposed by Cochran and Chan Cochran and Chan (1979) to reproduce the Hugoniot relation between the experimental data of shock and particle velocities. The parameters $\mathcal{A}$ and $\mathcal{B}$ have dimensions of pressure, while $\mathcal{C}$ and $\mathcal{D}$ are dimensionless. In this study, the parameter values for this EoS (summarized in Table 2) are adopted from the literature Michael and Nikiforakis (2019a, b) imposing a temperature $T_{2}=298~{}\mathrm{K}$ for liquid NM at its initial density $\rho_{2}=\rho_{0,2}$, which is a physically accurate value for the initial temperature of NM. Thus, the Cochran-Chan isothermal reference curve is at a temperature $T_{\mathrm{ref,2}}=0~{}\mathrm{K}$.

As described in this section, two major simplifications have been made in this model: (1) The Cochran-Chan EoS calibrated for liquid NM is also used to describe the reaction products; (2) the constant-volume specific heat capacities for both materials and the Grüneisen coefficient for NM are considered to be constant values corresponding to the initial state of the materials, i.e., $c_{\mathrm{v},i}(\rho_{i},T_{i})=c_{\mathrm{v},i}$ and $\Gamma_{2}(\rho_{2},T_{2})=\Gamma_{0,2}$, respectively. Simplification (1) is made to avoid the iterative root-finding procedure to determine the volume fractions of reactant and products in the reaction zone, thus, ensuring the robustness and time-efficiency of the simulations to explicitly resolve a large number of heterogeneities. It is of importance to note that both of these simplifications impose a limitation on the quantitative accuracy of the model, which is further discussed in Sect. V.5.

II.3 Reaction rate model

The reaction rate $K$ of liquid NM is governed by single-step Arrhenius kinetics as follows,

[TABLE]

where $T_{\mathrm{a}}$ is the activation temperature and $C$ is the pre-exponential factor. The value of $Q$ was calculated by Michael and NikiforakisMichael and Nikiforakis (2019b) following the approach proposed by Arienti et al. Arienti, Morano, and Shepherd (2004): Via varying the value of $Q$, the reactive Hugoniot curve can be shifted on the $p$-$v$ plane to match as closely as possible the CJ pressure of $12.5~{}\mathrm{GPa}$ reported by Dattelbaum et al. Dattelbaum et al. (2009). The value $Q=4.46\times 10^{6}\mathrm{J}/\mathrm{kg}$ has thus been obtained. Based on the velocity interferometer data obtained for the range of shock pressure from $7.5$ to $9.5~{}\mathrm{GPa}$, HardestyHardesty (1976) obtained $C=2.6\times 10^{9}~{}\mathrm{s}^{-1}$ and $T_{\mathrm{a}}=11500~{}\mathrm{K}$ for a reaction rate in the first-order, single-step Arrhenius form. Fixing the value of $C$ as calibrated by HardestyHardesty (1976), $T_{\mathrm{a}}$ was re-calibrated by Michael and NikiforkaisMichael and Nikiforakis (2019b) using the obtained value of $Q$ and the same EoS as that used in this study to match the in-situ gauge data of particle velocity reported by Sheffield et al. Sheffield et al. (2006) for neat NM shocked at $9.1~{}\mathrm{GPa}$. Thus, the value of $T_{\mathrm{a}}=11350~{}\mathrm{K}$ has been obtained and used in this work. The values for $C$, $T_{\mathrm{a}}$, and $Q$ are summarized in Table 3.

II.4 Initial configuration of the problem

Three different scenarios are considered in this study to examine the SDT behaviors: 1) Neat NM, 2) NM mixed with an array of regularly spaced circular cavities, and 3) NM mixed with randomly distributed cavities. The schematic illustration of these three scenarios is shown in Fig. 1. For all of the three scenarios, a rightward-moving incident shock wave is implemented via initially placing a shocked region near the left end of the domain with an inflow (transmissive) boundary condition applied upon the left boundary. Within this shocked region, the initial values of density, pressure, and particle velocity satisfy the Hugoniot relation for neat NM given the initial pre-shock conditions (see Sect. 5.1 in Michael’s PhD thesisMichael (2013)). The reaction progress variable is set to be [math] in this shocked region so that this amount of initially shocked NM does not contribute to the global energy release of the mixture. The longitudinal length of the shocked region $L_{\mathrm{s}}$ (in the propagation direction of the incident shock) is chosen to be $0.5~{}\mathrm{mm}$ for all of the simulations reported in this paper. The overall length of the domain ahead of the incident shock front $L$ is chosen to be sufficiently long so that the complete SDT process can be captured for the corresponding incident shock strength.

For the scenario with neat NM, since the resulting wave dynamics is expected to be one-dimensional, the transverse width of the domain $W$ is arbitrarily chosen to be $0.3~{}\mathrm{mm}$. For the case with regularly spaced cavities, as shown in Fig. 1(b), each cavity is of the same diameter $d_{\mathrm{c}}$ and spaced from its neighboring cavities by an equal distance $\delta_{\mathrm{c}}$. Specifying the overall porosity of the mixture $\phi$ (i.e., initial volume fraction of the air-filled cavities) and cavity diameter, the spacing between each two adjacent cavities can be calculated as follows,

[TABLE]

The transverse width of the domain for the scenario of regularly distributed cavities is set equal to $\delta_{\mathrm{c}}$ and periodic boundary conditions are applied on the top and bottom boundaries of the domain so that a two-dimensional array of equally spaced cavities can be simulated.

In the third scenario as shown in Fig. 1(c), the cavities are randomly distributed in a medium of NM. Knowing the total volume of the reactive mixture, i.e., $V_{\mathrm{tot}}=L\,W$, the total number of cavities in the domain can be calculated as

[TABLE]

In each simulation, the cavities are of the same size. The initial pressure and density of the unshocked air and NM are $p=10^{5}~{}\mathrm{Pa}$, $\rho_{1}=1.2~{}\mathrm{kg}/\mathrm{m}^{3}$, and $\rho_{2}=1134~{}\mathrm{kg}/\mathrm{m}^{3}$, respectively, for all of the simulations in this study. The total mass of reactive NM in the system can thus be calculated as $M_{2,\mathrm{tot}}=\rho_{2}(1-\phi)L\,W$.

III Numerical methodology

The simulation code used to solve the two-dimensional, reactive Euler equations is based upon a uniform Cartesian grid. The MUSCL-Hancock scheme with the van Leer non-smooth slope limiter and a Harten-Lax-van Leer-contact (HLLC) approximate solver for the Riemann problem were used. Toro (2009) The Strang splitting method was adopted to treat separately the hydrodynamic process and the reactive process. This numerical scheme is thus of second-order accuracy in space and time.

The simulation code was implemented in NVIDIA’s CUDA programming language to take the advantage of Graphic-Processing-Unit- (GPU-) enabled parallel computing. Each simulation was performed with a $16$-GB NVIDIA Tesla P100 GPU. The use of GPU-accelerated computing platforms has been explored in several studies on gaseous detonations. Morgan (2013); Kiyanda et al. (2015); Mi et al. (2017, 2018) The current code is a further development of the gaseous-detonation code to simulate detonations in a multiphase energetic system. This code is of a hybrid nature, consisting of parts that are executed on both CPU (Central Processing Unit) and GPU. Under a Godunov-type numerical scheme, CUDA kernel functions were implemented to parallelize the most compute-intensive parts within the solver for the reactive Euler equations, including computation of fluxes, reactive source term integration, and the update of conserved variables. In order to minimize dataflow between the global memory of a GPU and its streaming multiprocessors, the whole procedure of flux computation was implemented in such a way using the on-chip, shared memory of a GPU. Data transfer between host memory and GPU global memory is only required for domain initialization and data output at prescribed times. The code for post-simulation analysis was implemented in a serial fashion and executed on the host CPU.

A diffuse-interface approach was used in this work where two immiscible materials are separated by a thin numerical diffusion layer. Comparing to interface-tracking methodologies wherein a sharp material interface is resolved, a diffuse-interface approach demands a significantly lower computational cost. A sufficient grid resolution is however required to ensure a sufficiently thin diffusion layer at material interfaces. The validation of using this approach to simulate a single or a small number of air-filled cavities or solid particles under shock loading in liquid NM has been demonstrated in several previous studies.Michael (2013); Michael and Nikiforakis (2016, 2019a, 2019b); Michaell2019Book In this work, since a large number of cavities are present in the simulations, the total area covered by numerical diffusion zones was greater than that in previous studies. In order to avoid any numerical issues that can probably arise from entangled diffuse interfaces from neighboring collapsed cavities, a minimum spacing distance of $2.5d_{\mathrm{c}}$ was imposed to the random distribution of cavities. Further validation tests have been performed, and the computational grid-size independence of the key results has been verified and reported in Appendix  A.

IV Results and analysis

For the cases with air-filled cavities considered in this study, the overall porosity and cavity diameter were selected to be $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$. This selected value of cavity diameter is very close to the range of mean diameter of GMBs used to sensitize gelled NM as reported in a number of experimental studies.Gois, Campos, and Mendes (1996); Gois, Campos, and Plaksin (2002); Dattelbaum et al. (2010a, b); Higgins (2014); Higgins, Loiseau, and Mi (2018); Loiseau et al. (2018) The detailed comparison between these selected parameter values to realistic conditions of experimental studies on SDT in GMB-sensitized NM will be discussed in Sect. V.5. Sample results showing the evolution of wave structure for the scenarios with neat NM, regularly spaced cavities, and randomly distributed cavities are shown in Sect. IV.1 as contour plots of NM density, i.e., $z_{2}\rho_{2}$. To illustrate the difference in SDT behaviors in neat and heterogeneous NM mixtures, in Sect. IV.2 and IV.3, the corresponding spatio-temporal profiles of pressure and particle velocity, respectively, are compared. In Sect. IV.4, the global reaction rate as a function of time is shown for each scenario. Based on the time history of the global reaction rate, a characteristic time scale of the resulting SDT process, or the so-called “detonation overtake time”, can be measured. The simulation results of overtake time as a function of incident shock pressure are presented and compared to available experimental data of overtake time for GMB-sensitized, gelled NM in Sect. IV.5. Statistical analysis on the simulation data has been performed and reported in Sect. IV.6 to further investigate the collective effect of hot spots. To obtain the results reported in this section, simulations were performed at a numerical resolution of $\mathrm{d}x=\mathrm{d}y=1~{}\mu\mathrm{m}$, ensuring $100$ computational points across the diameter of a cavity. With evidence provided in the Appendices A and B, it has been verified, respectively, that further refining the computational grid or increasing the domain size in the $y$-direction would not alter the key findings in this study. It is worthwhile to remark that, as calculated by Michael and NikiforkaisMichael and Nikiforakis (2019b) (shown in Fig. 1 of their paper) using the same model for neat liquid NM, the steady Zel’dovich-von Neumann-Döring (ZND) detonation profile has a reaction zone length of approximately $200~{}\mu\mathrm{m}$, which is twice the cavity diameter $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ selected in this paper. At the nominal grid resolution of $\mathrm{d}x=1~{}\mu\mathrm{m}$ used in this study, there are thus nearly $200$ computational cells per reaction zone length of the ZND solution.

IV.1 Wave structure

The evolution of the wave structure over the course of an SDT process in neat NM, mixture with an array of regularly spaced cavities, and a mixture with randomly distributed cavities is shown in Figs. 2, 3, and 4, respectively. In these figures, the color contours represent the spatial distribution of NM density ($\rho_{2}z_{2}$). Each rectangular subfigure is a snapshot of the entire two-dimensional computational domain at a specific time. These subfigures are arranged from bottom to top in a chronological order. The incident shock pressure for the selected sample results is $7.81~{}\mathrm{GPa}$. The length of the domain is $24~{}\mathrm{mm}$.

As shown in Fig. 2, the wave structure resulting from the case with neat NM is one-dimensional. The time interval between each two consecutive snapshots is $0.3~{}\mu\mathrm{s}$. From the beginning ($t=0$), the incident shock wave can be identified as a rightward-moving abrupt change in density. From $t=3.0~{}\mu\mathrm{s}$ onwards, a high-density region appears significantly downstream from the incident shock front, e.g., at $t=2.7~{}\mu\mathrm{s}$, the incident shock has reached approximately $x=13~{}\mathrm{mm}$ and this downstream high-density region is at $x=6~{}\mathrm{mm}$. This high-density front moves rightwards with a trailing gradient and catches up with the incident shock slightly after $t=4.2~{}\mu\mathrm{s}$.

Snapshots showing the evolution of wave structure resulting from the scenario wherein a regular array of air-filled cavities embedded within a reactive matrix of NM are provided in Fig. 3. Zoom-in views of the flow field near the leading shock front at two different times, i.e., $t=2.4~{}\mu\mathrm{s}$ and $4.2~{}\mu\mathrm{s}$, are also shown in Fig. 3. For $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$, the spacing between each two neighboring cavities is $\delta_{\mathrm{c}}=0.313~{}\mathrm{mm}$. The resulting wave structure is two-dimensional and symmetric about the centerline along which the centers of the cavities are located. From $t=0$ to $t=2.7~{}\mu\mathrm{s}$, a region with non-uniformly distributed high-density loci (with $z_{2}\rho_{2}>2000~{}\mathrm{kg}/\mathrm{m}^{3}$) develops immediately behind the incident shock and spans over a significantly long distance downstream. The zoom-in view at $t=2.4~{}\mu\mathrm{s}$ reveals more details of this high-density region coupled to the leading shock front. From $t=3~{}\mu\mathrm{s}$ onwards, high-density loci can be found only in the vicinity of the incident shock interacting with a collapsing cavity (as shown in the zoom-in view at $t=4.2~{}\mu\mathrm{s}$). The high-density wave front that develops significantly downstream from the incident shock in the case with neat NM is not observed in the case where a regular array of cavities is present.

The evolution of the wave structure for the scenario with a spatially random distribution of cavities is shown in Fig. 4. The overall wave structure is qualitatively similar to that observed in the case with regularly distributed cavities: A region with high-density loci is closely attached to the incident shock from $t=0$ to approximately $t=2.7~{}\mu\mathrm{s}$; afterward, high-density loci are only present at the leading shock front as shown in the zoom-in view at $t=4.2~{}\mu\mathrm{s}$. Due to the intrinsic random distribution of the cavities, the resulting wave structure is highly irregular compared to that for the case with a regular distribution of cavities. The feature observed in homogeneous NM—a downstream developed high-density wave front catching up with the leading shock—is not observed here. Further analysis will be performed to elucidate the different behaviors resulting from neat and cavity-laden NM mixtures.

IV.2 $x$-$t$ diagrams

To further analyze the simulation results of the SDT process in neat NM, color contours of pressure and particle velocity in the $x$-direction for the case with an incident shock of $7.81~{}\mathrm{GPa}$ are plotted in an $x$-$t$ diagram as shown in Figs. 5(a) and 6(a), respectively. The trajectory of the incident shock wave can be identified as a sharp interface separating a low-pressure (black) region on the right and a relatively high-pressure (dull red) region on the left starting from the bottom of this diagram at $t=1~{}\mu\mathrm{s}$. The trajectory of the interface between the inflow inert material and the explosive is denoted by white dashed line in Figs. 5(a). In the shocked region behind the incident shock wave, the trajectory of another shock wave can be identified in Fig. 5(a) as a sharp interface with a very high-pressure (bright yellow) region on its left originating at roughly a locus of $x=5~{}\mathrm{mm}$ and $t=2.8~{}\mu\mathrm{s}$ on the trajectory of the inflow material interface. The pressure behind (on the left of) this second shock evidently increases over time from $t=2.8~{}\mu\mathrm{s}$ to $3.4~{}\mu\mathrm{s}$. A fan-shaped region of decreasing pressure in the negative $x$-direction is trailed by this second shock wave. The trajectories of the incident shock wave and the second shock wave intersect at $t=4.22~{}\mu\mathrm{s}$, suggesting that the second shock wave overtakes the incident shock wave at this time. Thus, $t=4.22~{}\mu\mathrm{s}$ can be considered as the simulation result of the overtake time of the SDT process in neat NM for an incident shock pressure of $7.81~{}\mathrm{GPa}$. After the overtake, a single shock front trailing a decreasing pressure gradient can be found near the top-right corner of this diagram.

As shown in Figs. 3 and 4, the wave dynamics resulting from the cases with cavities are two-dimensional. In order to examine the overall SDT process in an $x$-$t$ diagram, at a specific time, a one-dimensional profile (in the $x$-direction) of averaged flow and thermodynamic properties can be obtained via performing spatial averaging over the transverse direction (i.e., the $y$-direction) of the domain on the simulation data of the two-dimensional flow field. At time $t$, two-dimensional fields of pressure $p(x,y,t)$, density $\rho(x,y,t)$, and particle velocity in the $x$-direction $u(x,y,t)$ have been obtained from simulations. The one-dimensional profiles of average pressure $\bar{p}(x,t)$ and average density $\bar{\rho}(x,t)$ can be calculated as follows, respectively,

[TABLE]

and

[TABLE]

where $W$ is the width of the domain in the $y$-direction. With $\bar{\rho}(x,t)$ obtained, the density-weighted average profile of particle velocity $\bar{u}(x,t)$ be calculated as follows,

[TABLE]

The color contours of $\bar{p}$ for the cases with regularly and randomly distributed cavities subjected to an incident shock of $7.81~{}\mathrm{GPa}$ are thus plotted in an $x$-$t$ diagram as shown in Fig. 5(b) and (c), respectively. The results of $\bar{u}$ for the cavity-laden cases are used to obtain the time histories of particle velocity along particle paths (i.e., at Lagrangian positions) in Sect. IV.3.

Comparing the plots shown in Fig. 5(b) and (c), the overall SDT behaviors resulting from both heterogeneous mixtures with regularly and randomly distributed cavities are qualitatively similar: Only one leading shock trajectory, i.e., a sharp interface separating a low-pressure region on the right and a high-pressure region on the left, can be identified. The slope of this shock trajectory gradually changes over the time span of these diagrams. A smooth, but rather evident, turning point in the shock trajectory can be found roughly from $t=2.8~{}\mu\mathrm{s}$ to $3.2~{}\mu\mathrm{s}$ for the case with regularly spaced cavities, and $t=2.4~{}\mu\mathrm{s}$ to $3~{}\mu\mathrm{s}$ for the case with randomly spaced cavities. The pressure behind this shock wave gradually changes. No second shock wave with a region of pressure significantly greater than $20~{}\mathrm{GPa}$ can be seen on the $x$-$t$ plot for the cavity-laden cases. For the case with a regular distribution of cavities, as shown in Fig. 5(b), some regularly spaced fluctuations due to shock-cavity interactions can be observed in the contour plot of average pressure.

IV.3 Time history of particle velocity

The time histories of particle velocity along a series of particle paths spanning through the domain in the $x$-direction are obtained from the $x$-$t$ data of particle velocity $u$ (or $\bar{u}$ for the cavity-laden cases). Since the overall SDT behaviors resulting from the cases with regularly and randomly distributed cavities are similar, only the results of particle velocity histories for neat NM and a random distribution of cavities are reported in this section. The contour plots of $u(x,t)$ and $\bar{u}(x,t)$ resulting from the cases with neat NM and a random distribution of cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ subjected to an incident shock pressure of $7.81~{}\mathrm{GPa}$ are shown in Fig. 6(a) and (b), respectively. The red curves shown on these $x$-$t$ diagrams indicate the particle paths that are numerically tracked using forward Euler methods.

The time histories of $u$ for the case with neat NM along a number of selected particle paths are shown in Fig 7. The location of each curve on the axis of $x_{0}$ indicates the initial position of the corresponding particle path. Subfigures (a) and (b) of Fig 7 show the histories along the particle paths that are initially in the first (from $x=0$ to $12~{}\mathrm{mm}$) and second (from $x=12~{}\mathrm{mm}$ to $24~{}\mathrm{mm}$) halves of the domain, respectively. In each time history of $u$ shown in Fig. 7(a), there is an abrupt increase in $u$ at an early time followed by a period of constant velocity. At the end of this plateau, a second peak in $u$ occurs. The value of this second peak increases, and the time gap between the first increase and the second peak decreases from one particle path to the next. As shown in Fig 7(b), along the particle path initially located at $x=20~{}\mathrm{mm}$, the second peak in $u$ catches up to the first increase and merges into one abrupt increase.

For the case with randomly distributed cavities, as shown in Fig. 8, only one abrupt increase in $\bar{u}$ can be identified in the time history of $u$ along each particle path. This abrupt increase is followed by a gradual increase and decrease in particle velocity exhibiting a hump-shaped behavior. The peak value of this hump increases and its location becomes closer to the abrupt increase from one particle path to the next. The time history eventually evolves to a shape consisting of an abrupt increase followed by a gradually decreasing tail. As revealed by the contour plots in $x$-$t$ diagrams and the time histories of particle velocity, the SDT resulting from a cavity-laden NM mixture is a gradual process. It is rather inconvenient to measure a unique time scale to characterize this smooth transition process in heterogeneous NM, which is comparable to the overtake time for an SDT process in neat NM. Further analysis is required to obtain this characteristic time scale.

IV.4 Global reaction rate

Via analyzing the overall reaction rate of the mixture, this characteristic time, or the overtake time, for an SDT process in cavity-laden heterogeneous NM mixtures can be well defined, measured, and compared to that for the cases with neat NM. At each time step, the total mass of reactive NM remaining in the entire domain can be calculated via the following integration,

[TABLE]

The consumption rate of reactive NM in the entire domain can be obtained by differentiating $M_{2,\mathrm{r}}(t)$ over time, i.e., $\mathrm{d}M_{2,\mathrm{r}}/\mathrm{d}t$. For the simulations with the same incident shock strength but different scenarios of cavity distribution, the domain length $L$ is maintained the same, but domain width $W$ is different. Therefore, in order to compare the overall reaction progress resulting from cases with different domain widths on the same basis, the corresponding $\mathrm{d}M_{2,\mathrm{r}}/\mathrm{d}t$ needs to be normalized with respect to initial total mass of reactive NM. The normalization of the global reaction rate $\left|\bar{\lambda}_{t}\right|$ can be performed as follows,

[TABLE]

Since the consumption rate of reactive NM is negative, the absolute value of $\bar{\lambda}_{t}$ is taken so that the corresponding results can be plotted as positive values.

The normalized global reaction rate as a function of time for the case of neat NM with an incident shock pressure of $7.81~{}\mathrm{GPa}$ is plotted as the blue curve in Fig. 9(a). As shown on this figure, for the case of neat NM, the reaction rate is of an insignificant value before $t=2.4~{}\mu\mathrm{s}$. A rapid increase in $\left|\bar{\lambda}_{t}\right|$ occurs at $t=2.6$-$2.8~{}\mu\mathrm{s}$. After $t=3.2~{}\mu\mathrm{s}$, this rapid increase slows down until approximately $t=4.2~{}\mu\mathrm{s}$ where the curve of $\left|\bar{\lambda}_{t}\right|$ abruptly increases to its peak value. After reaching its peak, $\left|\bar{\lambda}_{t}\right|$ undergoes an abrupt decrease. From $t=4.3~{}\mu\mathrm{s}$ onwards, the reaction rate keeps decreasing but at much slower pace until the end of the simulation at $t=4.8~{}\mu\mathrm{s}$. The time at which the curve of $\left|\bar{\lambda}_{t}\right|$ reaches its peak is measured as $t=4.22~{}\mu\mathrm{s}$ from Fig. 9, agreeing with the measurement of the time when the second shock wave overtakes the incident shock front from the $x$-$t$ diagram plotted in Fig. 5(a). Thus, the time that has elapsed since the incident shock enters the reactive NM (at $t=0$) until the global reaction rate reaching its peak value can be considered as the overtake time for the resulting SDT process.

The results of $\left|\bar{\lambda}_{t}\right|$ as a function of time for the cases with a regular array of cavities and randomly distributed cavities subjected to an incident shock pressure of $7.81~{}\mathrm{GPa}$ are plotted as the green and red curves in Fig. 9(a), respectively. For both of these cavity-laden scenarios, the resulting reaction rates start to significantly increase very shortly after the incident shock enters the reactive NM mixture. As shown in Fig. 9(a), the results of $\left|\bar{\lambda}_{t}\right|$ for the cases with regularly and randomly distributed cavities remain approximately the same until $t=1.2~{}\mu\mathrm{s}$, and separate afterward: The global reaction rate for the case with a random distribution increases more rapidly than that for a regular array of cavities. The curve of $\left|\bar{\lambda}_{t}\right|$ for a random distribution reaches its peak at $t=2.67~{}\mu\mathrm{s}$ while that for the case of a regular array peaks at $t=2.88~{}\mu\mathrm{s}$. These times are thus measured as the overtake times for the SDT processes occurring in cavity-laden NM mixtures. After attaining the peaks, both curves gradually decrease to a quasi-plateau from roughly $t=3.3~{}\mu\mathrm{s}$ onwards.

Fluctuations can be observed on the curve of $\left|\bar{\lambda}_{t}\right|$ for the case with a regular array of cavities (i.e., green curve in Fig. 9(a)) throughout the entire transition process captured by the simulation. These fluctuations are of an equal frequency, i.e., local peak-to-valley time interval. The amplitude of these fluctuations is greater after the global reaction rate decreases to the quasi-plateau value. The curve of $\left|\bar{\lambda}_{t}\right|$ for the case with randomly distributed cavities (i.e., red curve in Fig. 9(a)) does not exhibit any significant fluctuations until it decreases to its quasi-plateau from $t=3.3~{}\mu\mathrm{s}$ onwards.

For a significantly greater input shock pressure, $11.51~{}\mathrm{GPa}$, the results of $\left|\bar{\lambda}_{t}\right|$ for the case of neat NM and cavity-laden heterogeneous NM mixtures are plotted in Fig. 9(b). The curves corresponding to neat NM and heterogeneous mixtures are nearly the same up to $t=0.2~{}\mu\mathrm{s}$. The resulting $\left|\bar{\lambda}_{t}\right|$ from the neat NM case increases faster and reaches a greater peak value than those from the cavity-laden cases. After approximately $t=0.75~{}\mu\mathrm{s}$, the $\left|\bar{\lambda}_{t}\right|$ curves for all of the three cases have reached a plateau or a quasi-plateau with fluctuations. The plateau value for the case of neat NM is noticeably greater than the mean quasi-plateau value for the heterogeneous cases. The amplitudes of the fluctuations on the curve for the case of regularly spaced cavities are greater than those for the case with randomly distributed cavities.

In the experimental studies of SDT phenomena in heterogeneous explosives, the commonly used approach to determine overtake time is via fitting a function to the trajectory of the leading shock front (obtained from the gauge measurement of in-situ particle velocity) and finding the time at which the shock acceleration is maximum, e.g., the approaches used by Hill and Gustavsen Hill and Gustavsen (2002) and Gustavsen et al. Gustavsen, Sheffield, and Alcon (2006). A similar approach of finding the maximum shock acceleration is applied to selected cases of the current study. The measurement of overtake time via this approach closely agrees with that obtained via finding the time of maximum overall reaction rate described in this section. The comparison between these two approaches can be found in Appendix C. The results of overtake time reported in the rest of this paper are obtained via the maximum reaction rate approach only.

IV.5 Overtake time vs. incident shock pressure

The simulation results of the relationship between the overtake time and the input shock pressure for the scenarios with neat NM and regularly and randomly distributed cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ are summarized in Fig. 10 and plotted as blue circles, open green triangles, and solid red triangles, respectively. Note that the data are plotted with linear scales in this figure. For all of these three scenarios, the resulting overtake time decreases as incident shock pressure increases. For an incident shock pressure less than $9.42~{}\mathrm{GPa}$, the overtake time for both cavity-laden cases is less than that of the neat NM case; the overtake time for a random distribution of cavities is slightly, but still noticeably, less than that for the case with regularly spaced cavities. The difference among the resulting overtake times from these three scenarios becomes greater as input shock pressure decreases. For an input shock pressure greater than $10.05~{}\mathrm{GPa}$, the overtake times for all of the three scenarios are very close while those resulting from the cavity-laden NM mixtures are very slightly greater than that for neat NM.

The simulation results of overtake time vs. input shock pressure for the scenarios with neat NM and randomly distributed cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ are replotted in Fig. 11 to compare with experimental data of the SDT overtake time in both neat and GMB-sensitized NM Campbell, Davis, and Travis (1961); Voskoboinikov, Bogomolov, and Apin (1968); Lysne and Hardesty (1973); Sheffield, Engelke, and Alcon (1989); Leal-Crouzet, Baudin, and Presles (2000); Dattelbaum et al. (2009, 2010a). Note that Fig. 11 is plotted with log-log scales–known as a “Pop-plot”–which is how the experimental results of overtake times as a function of input shock pressure are commonly reported in the literature. The simulation results for neat NM are plotted as blue circles. The dashed line is a fitting line to the scatter of experimental data of overtake times reported in several independent studies. Campbell, Davis, and Travis (1961); Voskoboinikov, Bogomolov, and Apin (1968); Lysne and Hardesty (1973); Sheffield, Engelke, and Alcon (1989); Leal-Crouzet, Baudin, and Presles (2000) The simulation results of the overtake times for neat NM qualitatively agree with the experimental data. The best quantitative agreement occurs over a range of input shock pressure from $8.14~{}\mathrm{GPa}$ to $9.42~{}\mathrm{GPa}$. The simulations result in greater overtake times than the experimental data for input shock pressures greater than $9.42~{}\mathrm{GPa}$. For input shock pressure less than $8.14~{}\mathrm{GPa}$, the simulation results of overtake time are slightly smaller than the experimental data.

In Fig. 11, the simulation results of overtake time for the cases with randomly distributed cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ are plotted as red triangles. The black squares and diamonds in Fig. 11 represent the experimentally measured overtake times for SDT processes in guar-gum-gelled NM mixed with GMBs of $d_{\mathrm{c}}=40~{}\mu\mathrm{m}$ and $\phi=0.84\%$ and $2.8\%$, respectively. Note that only weight fractions of GMBs for these experimentally tested mixtures, i.e., $0.36~{}\mathrm{wt}\%$ and $1.2~{}\mathrm{wt}\%$, respectively, are provided in the original publication Dattelbaum et al. (2010a); the corresponding volume fraction can be found in a later technical report Dattelbaum et al. (2010b) via examining the Scanning Electron Microscopy (SEM) images. In these experimental cases, the added GMBs are of lower porosities and a smaller size than those considered in the current simulations. The experimental results for GMB-sensitized NM mixtures demonstrate a greater reduction in overtake time from the results for neat NM than that resulting from the simulations. Possible reasons for the discrepancies between these simulation and experimental results (as shown in Fig. 11) will be discussed in Sect. V.5.

IV.6 Statistical analysis

As revealed in a number of numerical studies, the hot spots formed over the course of a shock-induced collapse of a cavity are of irregular shapes; the distribution of temperature and density within these hot spots is spatially non-uniform. Thus, temperature and size of a hot spot are rather vaguely defined concepts since both of them vary over time and reflect neither how much material is present within such a hot spot nor how fast energy is released from it. To further quantitatively analyze the effect of hot spots based on meso-resolved simulation data, it is of importance to know the amount of material reacting at specific reaction rates throughout the SDT process. To this end, mass-weighted probability density functions (PDFs) over a broad range of volumetric rate of energy release (in the unit of $\mathrm{GW}/\mathrm{m}^{3}$ and denoted as $\dot{q}$) have been obtained via performing statistical analysis on the simulation data. The procedure to obtain these PDFs is described in detail in the following paragraphs.

At each time step, the two-dimensional distribution of all flow and state variables can be obtained directly from the simulations. For each computational cell, the corresponding volumetric rate of energy release can be calculated as follows,

[TABLE]

which is the right-hand side of the energy balance equation in Eq. 1. Via scanning through the simulation data of $\dot{q}$ over the entire domain at specific time $t$, a statistical counting of how many computational cells that are associated with a $\dot{q}$ in a specific interval can then be performed. Since the values of $\dot{q}$ of interest spans over a broad range, i.e., over nearly six orders of magnitude, this range of interest is equally divided on a decadic logarithmic scale. Thus, each specific interval of $\dot{q}$ over which the statistical counting is performed can be expressed as

[TABLE]

where $\Delta n$ is the value of the equal interval on a decadic logarithmic scale of $\dot{q}$. To obtain all of the PDFs reported in this paper, $\Delta n=0.02$ was used. The area under each of these PDF curves can be interpreted as the proportion of the initial total mass of reactive NM releasing energy at volumetric rates within a specific interval of $\dot{q}$ at time $t$.

To take the mass of material reacting at a specific rate into account, the statistical counting is performed on a mass-weighted basis. As one scans through the computational cells one by one, $m(\dot{q},t)$, the probability mass measuring the cumulative occurrence of cells with its local $\dot{q}$ falling within a specific interval at a specific time $t$, can be counted by adding the mass of NM occupying the volume of this cell, i.e., $z_{2}(x,y,t)\rho_{2}(x,y,t)$. After counting the probability mass $m(\dot{q},t)$ over a range of $\dot{q}$ from $10^{2}$ to $10^{10}~{}\mathrm{GW}/\mathrm{m}^{3}$ throughout the entire domain, the mass-weighted probability density function of $\dot{q}$ can be calculated as follows,

[TABLE]

Note that $f(\dot{q},t)$ is based on the normalization with respect to the initial total mass of reactive NM so that the results of simulations with the same domain length but different domain widths can be compared on a common basis; the area under a PDF curve at a specific time is less than one. Since the material in the initially shocked zone and coming into the domain through the left boundary is set to be inert (i.e., product of NM with $\lambda=0$), it does not contribute to the statistical counting of the mass of reacting NM.

The mass-weighted PDFs of volumetric rate of energy release $\dot{q}$ for the cases of neat NM and a heterogeneous NM mixture with randomly distributed of cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ are compared in Fig. 12. Each subfigure of Fig. 12 is at a specific time which is marked at the top-left corner of the plot. The results shown in this figure are over the time from $t=0.4~{}\mu\mathrm{s}$ to $3.2~{}\mu\mathrm{s}$. The PDFs for neat NM and the case with randomly distributed cavities are plotted as blue and red curves, respectively. At $t=0.4~{}\mu\mathrm{s}$, the PDF curve for neat NM exhibits a sharp peak of nearly $f=1$ located at a narrow range slightly above $10^{5}~{}\mathrm{GW}/\mathrm{m}^{3}$ and $f=0$ elsewhere. The PDF curve for the heterogeneous case at $t=0.4~{}\mu\mathrm{s}$ is spread out over a range from approximately $\dot{q}=10^{4}~{}\mathrm{GW}/\mathrm{m}^{3}$ to $10^{6}~{}\mathrm{GW}/\mathrm{m}^{3}$. Without a sharp peak, the PDF curve for the heterogeneous case exhibits a smooth “hump” with its maximum probability density at approximately the same $\dot{q}$ where the high peak of the PDF curve for neat NM is located. As time progresses from $t=0.4~{}\mu\mathrm{s}$ to $2.8~{}\mu\mathrm{s}$, the PDF curve for the neat scenario gradually expands into the region of greater $\dot{q}$ up to above $10^{7}~{}\mathrm{GW}/\mathrm{m}^{3}$; the peak value of $f$ decreases and remains at roughly $\dot{q}=10^{5}~{}\mathrm{GW}/\mathrm{m}^{3}$. For the heterogeneous case, from $t=0.4~{}\mu\mathrm{s}$ to $1.8~{}\mu\mathrm{s}$, the PDF reaches its maximum at roughly the same $\dot{q}$ as where the peak for the neat NM curve is; the area under the PDF curve, as indicated by the red arrows, significantly increases in the range of greater $\dot{q}$ from roughly $10^{6}~{}\mathrm{GW}/\mathrm{m}^{3}$ to $10^{8}~{}\mathrm{GW}/\mathrm{m}^{3}$ over which the PDF for neat NM remains zero. The PDF curve for the heterogeneous case has its hump shifted away from $\dot{q}=10^{5}~{}\mathrm{GW}/\mathrm{m}^{3}$ after $t=2.4~{}\mu\mathrm{s}$ onwards, and eventually at $t=3.2~{}\mu\mathrm{s}$, it evolves to a much smaller hump located at roughly $\dot{q}=10^{6}~{}\mathrm{GW}/\mathrm{m}^{3}$. For the case with neat NM at $t=3.2~{}\mu\mathrm{s}$, a second peak appears in the region where $\dot{q}\geq 10^{8}~{}\mathrm{GW}/\mathrm{m}^{3}$ as indicated by the blue arrow.

In Fig. 13, the mass-weighted PDFs of $\dot{q}$ for the scenarios with an array of regularly spaced cavities (plotted as green curves) and randomly distributed cavities (plotted as red curves) of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ are compared. Note that the results for the case with a random distribution is the same as that shown in Fig. 12 but plotted with a smaller scale of $f(\dot{q},t)$, i.e., from [math] to $0.15$ in this figure instead of from [math] to $1$ in Fig. 12, in order to more clearly demonstrate the difference between the cases with a regular and random distributions. As shown in the two subfigures for $t=0.4~{}\mu\mathrm{s}$ and $0.8~{}\mu\mathrm{s}$, the PDF curves for regular and random cases have roughly the same overall behavior while the curve for the random case appears to be relatively smoother and the one for the regular case exhibits more small-amplitude fluctuations. From $t=1.2~{}\mu\mathrm{s}$ onwards, the area under both curves increases over the range of $\dot{q}$ from roughly $10^{6}~{}\mathrm{GW}/\mathrm{m}^{3}$ to $10^{8}~{}\mathrm{GW}/\mathrm{m}^{3}$. This increase in area over the range of greater $\dot{q}$ is faster for the case with randomly distributed cavities than that for the case with a regular distribution.

V Discussion

The meso-resolved simulation results and post-simulation analysis presented in the previous section reveal various facets of the SDT processes in neat NM and cavity-laden mixtures of NM. In this section, these results are interpreted to elucidate the key mechanisms underlying the shock-initiation behaviors of different scenarios. The findings from the current simulation results are compared with those speculated or experimentally discovered in previous studies. Possible reasons for the discrepancy between the simulation results and experimental data of the SDT overtake times are discussed.

V.1 Neat NM

The simulations for the case of neat NM performed in this study captured an SDT mechanism featuring a sudden thermal explosion after an induction period since the passage of the incident shock and a gradually developed superdetonation eventually overtaking the incident shock wave. All of these key features (as indicated on the $x$-$t$ diagram shown in Fig. 5(a) and the time histories of particle velocity shown in Fig. 7) agree with the SDT mechanism in neat NM elucidated in the literature.Holland, Campbell, and Malin (1957); Cook, Pack, and Gey (1958); Chaiken (1960); Campbell, Davis, and Travis (1961); Sheffield, Engelke, and Alcon (1989). As shown in Fig. 5(a), the thermal explosion occurs at the inflow material interface. The left-propagating blast wave emanating from the explosion propagates through the inflow material and increases the density in this region, which can be observed in the snapshots from $t=2.7~{}\mathrm{\mu s}$ onwards. It can be noticed from Figs. 5(a), 7, and Fig. 9(a) that, following the thermal explosion, the reaction wave develops into a superdetonation over a finite amount of time from approximately $t=2.8~{}\mu\mathrm{s}$ to $3.3~{}\mu\mathrm{s}$. This gradual development of a superdetonation in neat NM better agrees with that experimentally observed by Sheffield et al. via in-situ measurement of particle velocitySheffield, Engelke, and Alcon (1989), rather than Chaiken and Campbell et al.’s speculation Chaiken (1960); Campbell, Davis, and Travis (1961) that a steady superdetonation is immediately initiated by the thermal explosion. Further, it is of interest to note that the history of the global reaction rate for neat NM shown in Fig. 9(a) exhibits a qualitatively similar behavior as that of the temporal profile of the radiance emitted by neat NM undergoing an SDT process, which was recorded by Bouyer et al. for spectral analysis to determine the corresponding temperature profile (as shown in Fig. 4 of this paper) Bouyer et al. (2006).

V.2 Cavity-laden heterogeneous mixtures of NM

As shown by Fig. 5(b) and (c), the $x$-$t$ diagram of $y$-direction averaged pressure contour plot for the case of NM mixed with air-filled cavities, the mechanism governing the SDT is completely different from that for the case of neat NM. The pressure behind the leading shock gradually increases, resulting in a smooth transition from an incident shock to a quasi-steady detonation wave. No evidence of thermal explosion or superdetonation can be found in the $x$-$t$ diagram or the time histories of particle velocities (Fig. 8) for heterogeneous NM subjected to an incident shock pressure of $7.81~{}\mathrm{GPa}$. This difference can also be found by comparing the time histories of the global reaction rate for neat NM and cavity-laden heterogeneous mixtures of NM in Fig. 9(a). The temporal profiles of $\left|\bar{\lambda}_{t}\right|$ for the heterogeneous NM mixtures (i.e., red and green curves) start to gradually increase from a very early stage, at approximately $t=0.4~{}\mu\mathrm{s}$, without undergoing a long induction process. These analyses show that significant energy release is triggered immediately downstream from the leading shock and is gradually accelerated over time, which is consistent with the experimental findings for SDT processes in NM mixed with silica beads and GMBs Dattelbaum et al. (2009, 2010a) and heterogeneous solid explosives Ramsay and Popolato (1965).

This enhancement in SDT sensitivity has long been attributed to the formation and growth of hot spots. As revealed in a number of computational studies, the hot spots formed due to shock–heterogeneity interactions are of irregular shapes and sizes. Ball et al. (2000); Swantek, Shukla, and Austin (2010); Swantek and Austin (2010); Hawker and Ventikos (2012); Lauer et al. (2012); Ozlem et al. (2012); Michael and Nikiforakis (2014); Kapila et al. (2015); Betney et al. (2015); Apazidis (2016); Michael and Nikiforakis (2019a, b); Michall2019Book In addition, as shown in this study (e.g., Fig. 4) and others, the wave interaction due to the collapse of neighboring cavities may further complicate the sensitizing mechanisms of hot spots. Concerning these complexities, using hot-spot number density derived from the volume or weight fraction of the heterogeneous inclusions to quantitatively gauge the effect of hot spot is likely an oversimplification. Taking the advantage of the full flow-field information obtained from the meso-scale simulations, statistical analysis has been performed to reveal the mass spectrum (i.e., mass-weighted PDF) over a broad range in reaction rate at different stages in an SDT process.

For an input shock pressure of $7.81~{}\mathrm{GPa}$, as shown in Fig. 12, the area of the PDF curve for neat NM (blue curve) is initially concentrated at a relatively low reaction rate; the gradual expansion of this curve towards the direction of greater reaction rates indicates the increase in the mass of shocked material. In contrast, the PDF curve for the heterogeneous case (red curve) exhibits a hump-shape and is spread out over a wider range in reaction rate since a very early stage. The amount of mass at greater reaction rates beyond the span of the neat NM curve is likely associated with hot spots. The tail of the heterogeneous PDF curves into the region of lower reaction rates is attributed to the remaining amount of reacting material with smaller mass fractions of reactant $\lambda$. A continuously increasing amount of fast reacting material, as indicated by the red arrows in Fig. 12, further accelerates the initiation process, demonstrating the collective effect of a large number of hot spots regardless of their shapes and sizes.

V.3 Effect of input shock pressure

As shown in Fig. 10, the overtake times for the heterogeneous cases are significantly shorter than those for the cases with neat NM subjected to incident shock pressures less than $9.42~{}\mathrm{GPa}$; for an incident shock pressure greater than this value, the overtake time for neat NM is slightly shorter than those for heterogeneous NM mixtures. For a cavity-laden NM mixture subjected to a relatively weak input shock, as revealed by the simulation results and statistical analysis shown in Figs. 9(a) and 12, the initiation process is dominated by the mechanism of hot spots. The resulting SDT behaviors are thus considered to be in the hot-spot-driven regimes. Dattelbaum et al. (2010b) For increasingly weaker input shock strength, the sensitizing effect of hot spots is more pronounced as shown in Figs. 10 and 11. For further lower input shock pressures, shock focusing and reflection may not be able to induce sufficiently high temperature for fast reaction, and viscous heating is speculated to dominate the formation and growth of hot spots.Field (1992) It is thus of importance to incorporate material viscosity into the model to examine the SDT behavior under low input shock pressures.

If the strength of the input shock is greater than a certain value, given an elevated post-shock temperature, the induction process is very short for both neat NM and heterogeneous mixtures. As shown in Fig. 9(b), the global reaction rate for both neat NM and heterogeneous mixtures start to significantly increase immediately after the incident shock entering the reactive material and remain nearly the same at an early stage of the SDT process. These results suggest that the initiation process following a strong input shock is dominated by the energy release in the bulk material, rather than at localized high-temperature regions. Such SDT behaviors are considered to be in the thermally-driven regimes wherein the energetically diluting effect of adding cavities dominates over the sensitizing effect via producing hot spots. Dattelbaum et al. (2010b) Thus, for greater input shock pressures, the subsequent global reaction rates for the heterogeneous cases are lower than that for neat NM (as shown in Fig. 9(b)); the overtake times resulting from heterogeneous cases are slightly longer than those for neat NM (as shown in Fig. 10). For the selected values of porosity and cavity size (i.e., $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$), an input shock pressure of approximately $9.42~{}\mathrm{GPa}$, as indicated by the gray band in Fig. 10, marks the balance between the sensitizing and energetically diluting effects of the addition of cavities, or the boundary between hot-spot-driven and thermally-driven SDT regimes. The cross-over from a hot-spot-driven regime to a thermally-driven regime for an increasingly greater input shock strength has been experimentally identified from the SDT tests with gelled NM mixed with silica beads of $6\%$ mass fraction and $40~{}\mu\mathrm{m}$ diameter. Dattelbaum et al. (2009, 2010a)

V.4 Regular vs. random distributions of cavities

For the same porosity and size of the cavities, in hot-spot-driven regimes, the overtake times resulting from the cases with randomly distributed cavities are shorter than those from the cases with a regular distribution as shown in Fig. 10(a). Comparing the mass spectrum over reaction rate between the cases with random and regular distributions for an incident shock pressure of $7.81~{}\mathrm{GPa}$ as shown in Fig. 13, it can be noticed that, at $t=1.2~{}\mu\mathrm{s}$, the random distribution starts to result in a greater amount of material reacting at rates above $\dot{q}=10^{6}~{}\mathrm{GW}/\mathrm{m}^{3}$ than that for the regular case. Due to a thermally positive feedback mechanism, the amount of fast reacting material increases faster in the cases with a random distribution, and thus, sooner leads to the detonation overtake.

The further enhancement in shock-initiation sensitivity via randomizing the spatial distribution of cavities might be due to the fact that there are local regions with more closely spaced cavities, i.e., the spacing between neighboring cavities is smaller than the average spacing $\delta_{\mathrm{c}}$ calculated via Eq. 7. As the incident shock passes through such regions, the local number density of hot spots is greater than the average value. Thus, a local cluster of hot spots is formed. The localized explosions triggered by this cluster of hot spots may interact with each other via emanating strong pressure waves. This collective behavior of clustering hot spots may further enhance the energy release rate within the cluster and accelerate the overall SDT process via a thermally positive feedback mechanism.

V.5 Simulation results vs. experimental data

As shown in Fig. 11, the simulation results of overtake times for neat NM agree fairly well with the experimental data for relatively low input shock pressures over a range from $7.55~{}\mathrm{GPa}$ to $9.42~{}\mathrm{GPa}$, but exhibit an increasingly larger deviation as input shock pressure increases beyond $9.42~{}\mathrm{GPa}$. The agreement and discrepancy between simulation results and experimental data for neat NM reflect the fact that, in the reaction rate model adapted for these simulations (Eq. 6), the value of the pre-exponential factor $C$ is calibrated based on the experimental measurement of overtake time of liquid NM for shock pressures ranging from $7.5$ to $9.5~{}\mathrm{GPa}$, and the value for $T_{\mathrm{a}}$ was calibrated based on the in-situ gauge data of particle velocity for an input shock pressure of $9.1~{}\mathrm{GPa}$.

The choice and simplifications made to the EoS of NM and its products also impose major limitations to the quantitative accuracy of the model. As Menikoff and Shaw pointed out, the choice of EoS and reaction rate model is crucial for an accurate modeling of detonation waves in NM.Menikoff and Shaw (2011) Using the same Cochran-Chan EoS (Eqs. 2-4) to govern the reactant and products of liquid NM does not adequately describe the physics in the reaction zone and the subsequent flow expansion. In future efforts, following the approach proposed by Wilkinson et al. Wilkinson et al. (2017), the EoS for NM reactant and reaction rate model can be calibrated using the up-to-date gauge measurement data, and the product EoS can be obtained by fitting to the principal adiabat calculated by an ideal detonation code such as IDex. Simulations treating the reactant and product of an explosive as two miscible materials governed by different EoS can be implemented using the MiNi16 formulation Michael and Nikiforakis (2016, 2018). In addition, considering the constant-volume heat capacity and Grüneisen coefficient for NM to be constant values may also introduce inaccuracies to the simulation results. The thermodynamically consistent EoS developed by Winey et al.Winey et al. (2000) for shocked liquid NM with density- and temperature-dependent heat capacity, thermal pressure coefficient, isothermal bulk modulus, and Grüneisen coefficient can be incorporated in further development of the model.

For the cases under hot-spot-driven regimes, the (perhaps only) available experimentally measured overtake times for SDT in gelled NM mixed with GMBs of $d_{\mathrm{c}}=40~{}\mu\mathrm{m}$ and $\phi=0.84\%$ and $2.8\%$ are significantly less than those resulted from the current simulations with cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$. This difference can be attributed to several factors: (1) The overall volumetric number density of cavities; (2) clustering effect of hot spots due to a spatially non-uniform distribution of cavities; (3) the discrepancy in hot-spot temperature arising from two and three dimensions. Each factor is briefly discussed in the following paragraphs.

In a three-dimensional system, knowing the size and porosity of the cavities in the mixture, the volumetric number density can be calculated as

[TABLE]

For the two experimentally tested mixtures with GMBs of $d_{\mathrm{c}}=40~{}\mu\mathrm{m}$ and $\phi=0.84\%$ and $2.8\%$, the corresponding number densities are $\rho_{\mathrm{N}}\approx 250~{}{\mathrm{mm}}^{-3}$ and $835~{}{\mathrm{mm}}^{-3}$, respectively. For the heterogeneous case considered in the simulation results in Fig. 11, the equivalent three-dimensional number density of cavities is determined to be $\rho_{\mathrm{N}}\approx 153~{}{\mathrm{mm}}^{-3}$. Since the $\rho_{\mathrm{N}}$ of cavities for the simulation cases is smaller than the $\rho_{\mathrm{N}}$ of GMBs in the experimentally tested mixtures, the hot-spot sensitizing effect on the SDT process is likely less significant in the current simulation cases.

As discussed in Sect. V.5, the overall number density of cavities or hot spots may not be the sole parameter that determines the resulting SDT behaviors. The clustering of cavities can likely enhance their sensitizing effect. It is of importance to recall that, in the current simulations, a minimum cavity spacing of $2d_{\mathrm{c}}$ is imposed to the random distribution of cavities. For the cases with cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$, this imposed minimum cavity spacing is $250~{}\mu\mathrm{m}$ while the average spacing is $313~{}\mu\mathrm{m}$, which renders the random distribution in these simulations rather uniform over space. Likely due to this spatial uniformity of the random distribution of cavities, the clustering effect of hot spots is less pronounced in the simulations than that in the experimentally tested mixtures. In future studies, it might be of interest to impose some controllable spatial non-uniformities to the random distribution of meso-scale heterogeneities and probe the effect of such non-uniformities on the SDT behavior of the explosive mixtures.

The shock-induced collapse of a single air-filled cavity in liquid NM has been investigated by Michael and Nikiforakis via both two- and three-dimensional simulations using the same set of governing equations as that in the current work.Michael and Nikiforakis (2019a, b) These authors have demonstrated that the peak temperature reached by a hot spot in a three-dimensional configuration can be nearly $1.5$ times of the peak hot-spot temperature resulting from a two-dimensional scenario (as shown in Fig. $11$ of Ref. [28]). This difference in hot-spot temperature can lead to an order-of-magnitude greater reaction rate at the kernel of a hot spot formed during the shock-induced collapse of a spherical cavity in three dimensions. Hence, in the current two-dimensional simulations, the hot-spot effect is likely less significant than in realistic three-dimensional scenarios.

Further, numerical inaccuracies may be rooted in the diffuse-interface approach to separate the two fluids. A cavity collapses due to the passage of a shock wave, and the initially sharp interface grows into an artificial mixing layer during this process. As pointed out by Ozlem et al.Ozlem et al. (2012) and Michael and NikiforakisMichael and Nikiforakis (2019b), the thermodynamics in this artificial mixing layer are open to questions. Since the thickness of artificial mixing layer decreases for an increasingly finer grid resolution, the numerical convergence of the SDT overtake time results demonstrated in Appendix A indicates that the results are not significantly affected by the numerical uncertainty arising from the artificial mixing layers. To further examine this issue, it is of interest to compare the current results to those obtained from simulations using a sharp-interface approach, such as the level-set methods used to study void collapse in solid explosivesSambasivan, Kapahi, and Udaykumar (2013); Kapahi, Sambasivan, and Udaykumar (2013); Rai and Udaykumar (2015); Rai, Schmidt, and Udaykumar (2017); Rai and Udaykumar (2018).

VI Concluding remarks

The complete SDT process in a heterogeneous NM mixture that contains a statistically significant number of air-filled cavities has been studied via GPU-enabled computational simulations. Without invoking any empirically calibrated, phenomenological reaction models, the characteristic shock-induced initiation process in such a heterogeneous explosive mixture can be captured by these meso-resolved simulations. The SDT mechanisms revealed by the simulations for neat NM and cavity-laden NM mixtures are qualitatively consistent with those hypothesized based upon the experimental evidence. The simulation results also show that the sensitizing effect of the addition of heterogeneities is significant for a relatively weak input shock wave; for a sufficiently strong input shock wave, the presence of cavities slightly impedes the initiation process due to a dilution of the energetic medium. A uniformly random distribution of cavities has a noticeably more pronounced sensitizing effect on the shock-induced initiation of the mixture than a regular array of cavities does. With the detailed flow-field data obtained from the meso-resolved simulations, a framework of statistical analysis has been proposed to quantitatively examine the hot-spot mechanisms underlying the SDT process. Further development towards a quantitatively better agreement between the simulation and experimental results will be carried out in future efforts.

Acknowledgements.

XCM is supported by an NSERC Postdoctoral Fellowship (PDF-502505-2017). Computing resources used in this work were provided by Compute Canada and University of Cambridge. The authors are grateful to J. Loiseau for useful discussions in developing this paper.

Appendix A Numerical convergence study

In this paper, the key result obtained from each SDT simulation is the detonation overtake time. Thus, numerical convergence tests have been performed on this time scale, which characterizes the resulting SDT behavior. Since the cavity diameter is chosen to be $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ for all cavity-laden cases reported in this paper, this diameter is used for the convergence tests.

For the case with an array of regularly spaced cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ subjected to an input shock pressure of $7.81~{}\mathrm{GPa}$, the normalized global reaction rate $\left|\bar{\lambda}_{t}\right|$ (calculated via Eqs. 12 and 13) as a function of time resulting from simulations with various computational grid sizes are compared in Fig. 14. The overtake times obtained from this set of convergence tests for a regular distribution of cavities are summarized in Table 4. The temporal profile of $\left|\bar{\lambda}_{t}\right|$ and the result of overtake time do not change significantly as the computational grid size is further reduced beyond $\mathrm{d}x=1~{}\mu\mathrm{m}$, i.e., $100$ grids across the diameter of a cavity. For this selected case, the length and width of the computational domain are $24~{}\mathrm{mm}$ and $0.313~{}\mathrm{mm}$, respectively. At the highest numerical resolution reported in this study, i.e., $\mathrm{d}x=0.5~{}\mu\mathrm{m}$, there are $48000\times 626$ computational points in the domain. To obtain most of the results of overtake times, the simulations were performed at the second highest resolution, i.e., $\mathrm{d}x=1\mu\mathrm{m}$, with $24000\times 313$ computational points in the domain.

For the case with randomly distributed cavities of $\phi=5\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ subjected to an input shock pressure of $7.81~{}\mathrm{GPa}$, the normalized global reaction rate $\left|\bar{\lambda}_{t}\right|$ as a function of time resulting from simulations with various computational grid sizes are compared in Fig. 15. The overtake times obtained from this set of convergence tests for a regular distribution of cavities are summarized in Table 5. Refining the computational grid from $\mathrm{d}x=2~{}\mu\mathrm{m}$ to $1~{}\mu\mathrm{m}$ results in a much smaller change in overtake time than that resulting from a grid refinement from $\mathrm{d}x=4~{}\mu\mathrm{m}$ to $2~{}\mu\mathrm{m}$. As previously shown for the case with a regular distribution, further reducing $\mathrm{d}x$ below $1~{}\mu\mathrm{m}$ would not lead to a significant change in the global overtake time. Thus, all of the results for random distributions of cavities reported in this paper were obtained from simulations performed at a numerical resolution of $\mathrm{d}x=1~{}\mu\mathrm{m}$. Note that, as the transverse width of the computational domain $W$ has to be sufficiently large so that the overall resulting dynamics for randomly distributed cavities are statistically converged (as discussed in Appendix B). For the current case with an input shock pressure of $7.81~{}\mathrm{GPa}$, the length and width of the computational domain are $24~{}\mathrm{mm}$ and $2~{}\mathrm{mm}$, respectively. There are thus $24000\times 2000$ computational points in the domain. Given the current model formulation, the memory requirement of such a simulation is approximately $5~{}\mathrm{GB}$. Given a $16$-$\mathrm{GB}$ NVIDIA Tesla P100 GPU, it is unlikely feasible to double the numerical resolution for the same domain size. The use of multiple-GPU computing via Message Passing Interface (MPI) will be explored with future efforts to further refine the computational grid so that smaller cavities (on the order of $~{}10~{}\mu\mathrm{m}$) can be resolved.

Appendix B Statistical convergence study

For the cases with randomly distributed cavities, it is of importance to verify whether the results, e.g., the detonation overtake time, are statistically converged so that further increasing the size of the random distribution would not significantly alter the results. To this end, a set of testing simulations has been carried out for the scenario with randomly distributed cavities of $\phi=5\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ subjected to an input shock pressure of $7.81~{}\mathrm{GPa}$ at various widths ($W$) in the $y$-direction of the computational domain. The resulting temporal profiles of $\left|\bar{\lambda}_{t}\right|$ have been compared in Fig. 16. Note that, in order to extend the domain width beyond $W=2~{}\mathrm{mm}$ with a domain length of $24~{}\mathrm{mm}$, the numerical resolution of these simulations is chosen to be $\mathrm{d}x=2~{}\mu\mathrm{m}$. As shown in this figure, the profiles resulting from the simulations with $W=2$, $4$, $6$, and $8~{}\mathrm{mm}$ do not exhibit any significant difference in their overall behaviors and the detonation overtake time (where the value of $\left|\bar{\lambda}_{t}\right|$ peaks). Hence, in order to economize on the use of computing resource, the smallest domain width tested here, i.e., $W=2~{}\mathrm{mm}$, is selected for all of the simulations with randomly distributed cavities reported in this paper.

Appendix C An alternative method to determine detonation overtake time

The trajectory of the leading shock front $x_{\mathrm{s}}(t)$ is tracked via finding the location where pressure first increases to $1.01$ times the initial pressure. A moving-averaged differentiation is performed to the obtained $x_{\mathrm{s}}(t)$ to calculate the time history of leading shock velocity $u_{\mathrm{s}}(t)$ and acceleration $a_{\mathrm{s}}(t)$ as follows, respectively,

[TABLE]

and

[TABLE]

where $\Delta\tau$ is the width of the time window over which the moving-averaged derivative is calculated. In order to filter high-frequency noises, $\Delta\tau=1.5~{}\mathrm{\mu s}$ is chosen to obtain the data of shock acceleration reported in this Appendix.

The shock acceleration histories (thin orange curves against the right $y$-axis) obtained for the cases with neat NM and NM mixed with (b) regularly and (c) randomly distributed cavities of $\phi=8\%$ and $d_{\mathrm{c}}=100~{}\mu\mathrm{m}$ are compared with the corresponding histories of the overall reaction rate (thick blue curves against the left $y$-axis) in Fig. 17. The time at which the shock acceleration reaches its maximum closely agrees with the time of maximum $\left|\bar{\lambda}_{t}\right|$. Since the exact measurement of the time of maximum $a_{\mathrm{s}}$ depends on the time window width $\Delta\tau$ for moving-averaged differentiation (otherwise, influenced by high-frequency fluctuations in the velocity and acceleration histories), it is more convenient and accurate to use the approach of maximum overall reaction rate to determine detonation overtake time for the simulation results in this study.