A multi-scale simulation of hot spot initiation of detonation utilizing experimental measurements

Empirical and phenomenological hydrodynamic reactive flow models, such as the ignition-and-growth and Johnson–Tang–Forest models, have been effective in predicting the shock initiation and detonation characteristics of various energetic substances. These models utilize the compression and pressure properties of the reacting mixture to quantify its reaction rate. However, it has long been known that the shock initiation of detonation is controlled by local reaction sites called ‘hot spots’. In this study, a hot-spot model based on the temperature-dependent Arrhenius reaction rate is developed. The complex reaction process of the target explosive is addressed by conducting differential scanning calorimetry experiments whereas the reaction rate is determined using the Friedman isoconversional method. The hot spot is approximated by the region of high pressure accumulation due to multiple shock reverberations within the polymer binder, which is surrounded by the bulk explosive. The mesoscale smoothed particle hydrodynamic simulation is adopted to identify the peak temperatures within the hot spots. These peak temperatures obtained from the mesoscale level are then used to initialize the random sites of heat release prior to conducting the full-scale hydrodynamic simulation of the shock-to-detonation transition (SDT). To validate the simulation, the distance to detonation is compared with the reported experimental value to validate the initiation process of the proposed model and an 18-mm-radius rate stick is experimentally tested to confirm the reproducibility of the detonation properties. The comparison shows that the detonation properties and the initiation process of the explosive are well characterized, while no-go conditions are observed if no mesoscale hot-spot model is included in the hydrodynamic simulation. Therefore, the SDT process can be well described by the present model based on multi-scale hot-spot initiation.

Empirical and phenomenological hydrodynamic reactive flow models, such as the ignition-and-growth and Johnson-Tang-Forest models, have been effective in predicting the shock initiation and detonation characteristics of various energetic substances.These models utilize the compression and pressure properties of the reacting mixture to quantify its reaction rate.However, it has long been known that the shock initiation of detonation is controlled by local reaction sites called 'hot spots'.In this study, a hot-spot model based on the temperature-dependent Arrhenius reaction rate is developed.The complex reaction process of the target explosive is addressed by conducting differential scanning calorimetry experiments whereas the reaction rate is determined using the Friedman isoconversional method.The hot spot is approximated by the region of high pressure accumulation due to multiple shock reverberations within the polymer binder, which is surrounded by the bulk explosive.The mesoscale smoothed particle hydrodynamic simulation is adopted to identify the peak temperatures within the hot spots.These peak temperatures obtained from the mesoscale level are then used to initialize the random sites of heat release prior to conducting the full-scale hydrodynamic simulation of the shock-to-detonation transition (SDT).To validate the simulation, the distance to detonation is compared with the reported experimental value to validate the initiation process of the proposed model and an 18-mm-radius rate stick is experimentally tested to confirm the reproducibility of the detonation properties.The comparison shows that the detonation properties and the initiation process of the explosive are well characterized, while no-go conditions are observed if no mesoscale hot-spot model is included in the hydrodynamic simulation.Therefore, the SDT process can be well described by the present model based on multi-scale hot-spot initiation.© 2018 Author(s).All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).https://doi.org/10.1063/1.5041761

I. INTRODUCTION
Phenomenological hydrodynamic reactive flow models, such as the ignition-and-growth 1 and Johnson-Tang-Forest 2 models, have been successful in predicting the shock initiation and detonation characteristics of solid explosives.These models utilize the compression and pressure properties of the reacting mixtures to quantify their reaction rates.In addition to these conventional pressure-dependent models, several advanced models describe the shock-to-detonation transition (SDT) phenomenon.The CREST model 3 uses the reaction rate parameters obtained from particle velocity measurements and represents them as a function of the entropy of the unreacted sample.It also describes the porosity of the explosive sample using the snowplough model and the two-temperature model 4 and can calculate the effect of the initial temperature and porosity on the propagation of detonation and ignition.The SURF model is based on the ignition-and-growth model and includes a hot-spot model. 5t is a pressure-dependent rate law in which various physical factors, such as the initial temperature, porosity, and distribution of hot spots, are included in a modular form.
In this work we propose a model that is distinct from the currently available ones, as it separates the mesoscale hot-spot evolution for the SDT from the macroscale Arrhenius-type detonation of a bulk explosive.Thus, the particle-based smoothed particle hydrodynamic (SPH) method identifies the sites of high energy release from the formation of hot spots and then uses the results to initialize the full-scale hydrodynamic simulation to describe the detonation.Certainly, all pressure-based ignition-and-growth models are excellent for attaining a fully developed detonation.However, the ignition part of the SDT process can be handled more precisely because the growth process is dependent on the amount of accumulated energy released from the various sites of hot spots in a bulk sample. 6irst, we describe a method of constructing a temperature-based detonation rate law (Arrhenius) for any unknown explosive composition using a laboratory-scale calorimeter, which necessitates only a small sample mass.According to the existing literature, conventional thermal-based hotspot models, which are generally based on one-step kinetics, fail to fully describe the complex process of hot-spot initiated reactions; therefore, complex decomposition kinetics are preferred for enhancing the simulation. 7In this study, isoconversional decomposition kinetics is experimentally measured using differential scanning calorimetry (DSC).The three typical methods used to determine the reaction kinetics of energetic materials are as follows.The first approach considers a detailed mechanism that covers more than a hundred reactions. 8Most of such high-energy materials have complex reaction paths, which cannot be described by a single reaction step.Nevertheless, kinetics simulations involving multiple specific reactions is time consuming and does not allow efficient chemistry computations.The second choice, a multi-step scheme, assumes only a few reaction steps that determine the global chemical reaction, 9 thus offering a significant reduction in the computation time with reasonable accuracy.Lastly, the method that utilizes Friedman's isoconversional method 10 efficiently determines factors such as the activation energy as well as the pre-exponential factor, which varies with each degree of conversion during the entire reaction.The kinetics obtained through this method does not require multiple chemical steps.Instead, a set of multiple Arrhenius factors-namely, the activation energy and the pre-exponential factor-is constructed based on the local progress of the entire exothermic reaction.Thus, we build the reaction kinetics by following the reaction progress of a designated energetic sample.The kinetic scheme derived in this way features a single step with tabulated Arrhenius terms parametrized with respect to the local reaction progress variable.This provides a computational advantage over conventional multi-step kinetics, as well as precision in the reactive flow simulations. 7fter constructing the Arrhenius-type rate law, we use an SPH-based mesoscale simulation to simulate the temperature evolution during SDT from hot spots.The size of the hot spot is approximately 1 µm and the mesh size is kept smaller than 1 µm to capture the peak temperature evolution of the hot spots.As the diameter of the explosive sample is of the order of 10 mm, the mesh resolution is insufficient for covering the hot-spot domain during the macroscale detonation simulation.Therefore, for a detonation calculation that properly accounts for hot spots, a macroscale hydrocode simulation should be linked to a mesoscale simulation of a hot spot through a Lagrangian-Eulerian multiscale framework.To obtain the peak temperatures in the mesoscale simulation, a Lagrangian SPH is adopted.The hot spot is approximated by the region of high pressure accumulation due to multiple shock reverberations within the lower-density region of a binder (Estane), which is surrounded by the bulk of cyclotetramethylene-tetranitramine (HMX).To describe this region, we consider three different sizes of hot spot bounded by the end walls of the HMX and with one end exposed to the high-velocity impact.The peak temperatures obtained from the mesoscale-level hot spots are then used to initialize the random sites of heat release before conducting the full-scale hydrodynamic simulation of the SDT.In the scanning electron microscope (SEM) image of the explosive sample, we locate the initial sites of designated heat release.To validate the simulation, the distance to detonation is compared with the reported experimental value to confirm the initiation process of the proposed model, and an 18-mm-radius rate stick is designed and tested.
In this study, the target material comprised 95% HMX and 5% Estane binder (by weight) and its initial density after pressing was 1.82 g/cc.Hereafter, we refer to the target material as "95% HMX".

A. DSC experimental setup
The DSC experiments were performed with a Mettler Toledo DSC 3 instrument utilizing a 40-µL aluminum crucible.Sealed aluminum pans were used to endure the pressure generated by the exothermic reaction of the energetic materials.Additionally, because evaporation is faster than decomposition when the sample is in an open-pan DSC, 11 a closed-pan DSC was employed to observe the exothermic chemical reaction.Furthermore, a purge flow of nitrogen at 80 mL/min was maintained throughout the experiment.The energetic materials were tested under non-isothermal experimental conditions, where the product of the sample mass and the heating rate was maintained below 1 mg • C/min to prevent self-heating of the sample. 12The sample masses were in the range 0.2-1.7 mg.Four constant heating rates of 0.5, 1.0, 2.0, and 4.0 • C/min were utilized between 35 and 600 • C.

B. Isoconversional method
The DSC is a type of calorimetry performed using a very small sample.A sample pan and a reference pan are placed in the furnace where the heat is applied, and both are heated to the same temperature.In this situation, the difference between the exothermic energy values of the sample pan and the reference pan is caused by the chemical reaction of the sample, and the difference of the generated energy is transferred between the two pans by thermal equilibrium.The transferred heat flow is measured as a DSC signal.The DSC traces have the form shown in Fig. 1, where the mass fraction of the product λ and the reaction rate dλ/dt are obtained from equations ( 1)-(3). 13 dt = S(t) -B(t) Here, S(t) is the DSC signal as a function of time t.B(t) is called the baseline, which is a reference line for determining the magnitude of the heat flow.The construction of the baseline involves the superposition of tangents at each side of the exothermic signal peak.Each tangent is linked each other through the product mass fraction in the baseline function.The reaction rate at time t is the instantaneous heat flow divided by the sum of the instantaneous energy releases during the full chemical reaction process.The product mass fraction at the time is given by the sum of the instantaneous energy releases divided by the total energy release.
The DSC signals of 95% HMX are plotted in Fig. 2. In this figure, the peak value of the signal increases with the heating rate.Moreover, the reaction begins and terminates at a higher temperature as heating rate increase.
The heat of the reaction (Q) was calculated from the DSC signal at each heating rate using Eq. ( 4).The average value of all heating rates was assumed to be the effective heat of the reaction for the given sample.The maximum deviation in the experimental values was found to be 8%, which is considered acceptable.
Using equation ( 4), the obtained heat of the reaction for the 95% HMX sample is Q = 1293.6J/g.This value indicates the energy released per unit mass of the mixture.
The reaction rate was calculated using a relationship expressed in the form of the Arrhenius equation, as shown in Eq. (5).
Here, R, t, T, A λ , and E λ indicate the universal gas constant, time, temperature, frequency factor, and activation energy for the instantaneous product mass fraction λ, respectively.The function f (λ) expresses the dependence of the reaction rate on the reaction progress.Applying a logarithm on both sides of Eq. ( 5) yields Eq. ( 6).
In the logarithmic plot obtained from Eq. ( 6), -E λ /R is the slope of the linear curve and ln[A λ f (λ)] corresponds to its intercept with the vertical axis, ln[dλ/dt].The reaction rate and reaction progress were obtained from the DSC signals shown in Fig. 2 to conduct Friedman analysis from the Arrhenius plot.
Friedman analysis was used to extract the kinetics parameters of the target material and the analysis was conducted on the Arrhenius plot (reaction rate-inverse of temperature) using the DSC experimental data.The Friedman analysis results for 95% HMX are presented in Fig. 3, where each measured DSC signal corresponds to a reaction progress curve on the Arrhenius plot.The straight line connects equal reaction progress states for each signal.The slope of these lines and their intercept with the vertical axis represent the activation energy and pre-exponential factor at that particular reaction progress, λ, respectively.For instance, the straight line in Fig. 3  to a reaction progress from 0 (unreacted) to 1 (completely reacted).With this approach, we obtain a high correlation factor of 0.99, which represents the linear relationship between the DSC signals and the heating rates in the Arrhenius plot.
Friedman analysis assumes that the straight line connecting equal reaction progress states for each DSC signal passes through the same reaction progress state of other DSC signals of higher or lower heating rates.This means that a straight line also passes through the other temperature regions, where no experiment was conducted.This is because the straight line is obtained statistically from four different signals and, for a high correlation factor, the DSC signals at all other heating rates will also be distributed along this straight line.This is why we can apply the extracted kinetics in the higher-temperature regions.
Regarding the feasibility of the extracted kinetics for describing the SDT process, another important point is the pressure conditions during the DSC experiment.During the DSC experiment that was performed for kinetics extraction, the reaction heat of the samples was measured using the closed pan.As mentioned in Section II A, the closed pan is a sealed pan with a lid attached using a sealing machine, which can withstand the pressure of up to approximately 100 atm produced by the exothermic reaction of the energetic material.This means that the DSC experiment was performed in an environment that represents the explosion pressure of an amount of target energetic material in the milligram scale, not at conditions of 1 atm.However, the sample used was very small, so it did not reach the detonation pressure that would be produced from a larger sample.Nevertheless, the DSC experiment was conducted at elevated pressure conditions, which coincided with the applied sample mass.
The activation energy and pre-exponential factor as a function of λ are shown in Fig. 4. The figure indicates that the extracted kinetics describes the complete process of the chemical reaction through a set of Arrhenius parameters obtained at an instantaneous state of the reaction progress.

C. Kinetics verification
To verify the kinetics extracted by following the method described, the DSC measurements were reproduced using the proposed kinetic parameters.The equations that describe the evolution of the temperature and the chemical progress are as follows.
Here, ẇ is the heating rate.During the DSC experiment, samples of a few milligrams were utilized, and the reaction heat was released to the reference material. 13Thus, the course of the reaction was not influenced by the generated heat, which is neglected in Eq. (7).The Arrhenius parameters obtained from Fig. 4 were used in the calculation.
A comparison of the reaction progress obtained from the DSC experiments and from the simulations is shown in Fig. 5.The simulation results agree well with the experimental measurement, suggesting that the kinetics extraction is suitable for describing the reaction progress of 95% HMX.

III. MESOSCALE HOT-SPOT SIMULATION VIA SPH
The mesoscale simulation was conducted using a fully Lagrangian SPH method.We constructed an SPH code and performed validation simulations.The temperature evolution data calculated from the SPH simulation were used in the macroscale hydrocode simulation to initialize random regions of peak heat release during the SDT process.

A. Numerical method
Any function f (x) can be represented through its convolution with a kernel function, W, as described in Ref. 14. Kernel interpolants are the basis of the SPH method.
where h is the smoothing length defining the influence domain of the kernel function.The kernel function has the following two properties: where δ(x − x ) is the delta function.
where m is a constant that defines the effective (non-zero) area of the smoothing function, called 'support domain' (mh).At the interpolation point x j , the support domain defines the volume V j .V j is defined using the mass m j and the density ρ j as follows for a particle located within mh.Thus, any interpolation point x j can be associated with the particle of mass m j and density ρ j .With the help of the integral representation ( 9), the discrete or particle approximation of any function at the location of particle x i can be written as Here, W ij denotes W (x i − x j , h), N is the number of interpolation points (hereafter referred to as the number of particles), and f j is the value of the function associated with particle j.If a kernel with compact support is chosen, then the size of the support domain defines the number of neighboring particles.
With the help of (13), the particle approximation for a gradient of the field function at the location of particle i can be written as where ∇W ij = ∂W ij /∂x i and the gradient is calculated with respect to x i .In this equation, the term f i has been added to ensure that the derivatives of a constant function vanish.If W ij is a function of the distance between particles, r ij = x i − x j , then and The momentum equation is written as follows By applying the particle approximation of Eqs. ( 9)-( 16), Eq. ( 17) is converted into an SPH form.
Equation ( 18) represents the SPH form of the momentum equation as solved in the simulation.We consider the spring force between particles to describe the solid behavior of HMX and Estane.The spring force is f spring,α = kx α with an effective spring constant of k = E l 0 (x i , x j , t) .α represents the spatial dimension, E is Young's modulus, and l 0 (x i , x j , t) is the direction vector in the undeformed state between particles i and j.The value of Young's modulus is listed in Table I.
The energy equation is written as follows By applying the particle approximation of Eqs. ( 9)-( 16), Eq. ( 19) is converted into an SPH form.v ij represents the velocity difference between particles i and j and e i represents the internal energy of particle i.The temperature of particle i, T i , is calculated from the following equation where c v, i represents the specific heat of particle i at constant volume.
For a kernel function, we adopt the quadratic spike function to avoid particle clustering issues.The kernel function is written as follows where q = r/h, r is the distance between two particles, d is the number of spatial dimensions, and α d is the normalization constant.
The density is calculated from the updated volume.In the SPH simulation, an initial mass is assigned to each particle.The initial mass m ini is calculated by m ini = ρ 0 A/N, where ρ 0 is the initial density of the material and A is the considered area.In our simulation, because each particle is bounded by the spring force, it moves with the spring velocity at each time step.Subsequently, as the position of the particle is changed, the relevant volume also changes and the density is updated accordingly.
In addition to the governing equations given above, the Mie-Gruneisen equation of state (EOS) was used to calculate the pressure of the material that satisfies the mathematical closure condition. where where η = 1 − (ν/ν 0 ) and ν = 1/ρ.
The Mie-Gruneisen factors and material properties are summarized in Table I. 15,16 B. SPH code validation

Riemann problem at shock condition
To validate our SPH code, we calculated the Riemann problem at shock conditions.Through this validation process, we confirmed the accuracy of the description of the material interface and the shock interface.We initialized the problem in the following manner.The length of the shock tube was 1 mm and the material interface was initially located at x = 0.5 mm.For the left side of material interface, the properties of HMX at the shock state were assigned as ρ left = 2400 kg/m 3 , P left = 9.5×10 9  The results are shown in Fig. 6.The solid line represents the exact solution of the Riemann problem and the dotted line marks the calculation shown at 0.04 µs.Upon initiation, the shock wave is propagated to the right, while rarefaction is propagated to the left with the contact discontinuity shown.The result is in good agreement with the exact solution, suggesting that the SPH formulation was properly implemented to describe the interface motions.

Square-patch problem
The second validation problem involved the rotation of a square patch.This problem was chosen to test the ability of the code to handle the free-surface motion precisely and establish the robustness of the constructed SPH method. 17Although this validation problem was conducted for the fluid state in particular, it is directly applicable to the particle clustering of the hot-spot problem.
A square patch of a viscous fluid was initially subjected to the following velocity field: v x = ωy, v y = −ωx.Here, ω is the angular velocity, set to 100 s −1 .At the start of the simulation, the square patch underwent large deformation of its free boundary due to a centrifugal force.We investigated the capability of our method to handle such complex instabilities without particle clustering.The initial length of each side was 2 m, and the particle number was 2601.The Tait EOS was applied, and the initial pressure field was set to zero everywhere.
The simulation results are plotted in Fig. 7.The dotted line shows the exact solution for the location of the instability. 17If the particle clustering problem was not resolved, then the location of the instability would differ from the exact solution because the particle approximation would be interrupted by particle clustering.As shown in the results, the SPH code was able to match the exact solution with high precision.

C. Temperature evolution within a polymer binder (Estane) upon impact
The SPH was used to calculate the peak temperature of a hot spot that was approximated by the region of high pressure accumulation due to multiple shock reverberations within the Estane, which was surrounded by the bulk of HMX.The hot spot was generated in the Estane, instead of the HMX, because it has lower wave impedance.As the pressure wave was propagated from the HMX to the  Estane, the end walls prevented the compression wave from escaping and allowed pressure building, leading to a peak temperature.
To recreate this process via the SPH technique, the initial simulation domain was defined as in Fig. 8.The left end of the HMX is impacted at a speed of 700 m/s.Here, the red symbols represent HMX and the gray ones represent Estane.In the simulation, the gap size was defined as the size of the binder, namely, Estane.Gaps of 1, 3, and 5 µm were considered for the distance between two bounding HMX walls along the impact direction, and the maximum temperature within such gaps was examined.The entire length of the simulation domain was 10 µm, with the Estane region located in the middle.Thus, the variation in gap size corresponded to a specific hot-spot size, and the corresponding HMX wall thickness also varied.For all gap sizes considered, the maximum temperature occurred near the center of the Estane upon the application of pressure reverberations in the Estane.The SPH simulation required 2500 particles with a smoothing length of 1.3∆x for the corresponding quadratic spike function.
Figure 9 shows the temperature distribution at the instant when a peak temperature is achieved.The corresponding interface between HMX and Estane begins at 4.5, 3.5, and 2.5 µm for a gap distance of 1, 3, and 5 µm, respectively.Accordingly, the peak values occur within the low-density region of the sample owing to shock-wave coalescence inside the gap.The peak temperatures and the corresponding gap sizes are summarized in Table II.

IV. SDT EXPERIMENT AND HYDROCODE SIMULATION
The data obtained from the SPH-based mesoscale simulation can be utilized to initialize a localized hot spot capable of achieving instantaneous heat release.Subsequently, the mesoscale SPH simulation identifies the peak temperatures at the site of the hot spot, and these peak temperatures are then used to initialize random sites of heat release prior to the full-scale hydrodynamic SDT simulation.Here, the SEM image of the explosive sample is used to obtain the distribution of the localized sites of heat release.
To validate the simulation, an 18-mm-radius rate stick is experimentally tested, and its measured and calculated detonation properties are compared.

A. Rate-stick test
The r-ate-stick test is an experiment aimed at characterizing the detonation velocity according to the radius of the target explosive.In our research, a rate stick with a diameter of 18 mm was tested for the validation of the proposed model.Figure 10 depicts a schematic of the test.A booster charge is initiated to allow the development of a full detonation in the stick to measure the velocity of the detonation.The steady detonation wave is observed by utilizing a cylindrical charge of sufficient length to inhibit any booster overdrives or initiation transients, which can persist for up to five times the initial charge diameter.A charge length-to-diameter ratio of 7 was chosen to adequately capture the fully developed detonation velocities.Electronic pins were utilized to measure the propagation velocity of the detonation wave.We observed that if the installed pin was exposed to a high-temperature pressure shock wave, breakage occurred, and a time-to-digital converter (TDC) was used to measure the arrival time of the shock wave by sensing an ascending (or descending) waveform of input pulses.We used the pin CA-1041 (DYNASEN, Inc.) and employed the high-speed counter Lec-4208 (LeCroy, Inc.) as a wide-range real-time TDC.The TDC was designed to perform time measurements in real time and required a wide dynamic range with high resolution.The propagation of the detonation front along the length of each charge was measured using five equally spaced sensor probes.In each test, the TDC was used with a bandwidth of 1 GHz or a frequency resolution of 1 ns (109 s).

B. Numerical simulation setup
The conservation laws of mass, momentum, and energy in axisymmetric coordinates are described by the equation Here, the vectors represent the conservation variables U, the spatial fluxes in the axial and radial directions are denoted as E and F, respectively, and the source term S represents a variety of multimaterial loading conditions.In the Eulerian framework, each vector is described as follows.
FIG. 10.Schematic of the rate-stick test.

105217-12
Kim, Lee, and Yoh AIP Advances 8, 105217 (2018) Here, ϕ takes the values 0 and 1 for rectangular and cylindrical coordinates, respectively, ρ represents the density, and u z and u r denote the velocity components in axial and radial coordinates, correspondingly.The equation E = e + (u 2 z + u 2 r )/2 expresses the total energy per unit mass, where e is the specific internal energy p is the hydrostatic pressure, and Q = 1293.6J/g, as extracted in Section II B.
The parameterization of the JWL EOS is performed by multiple CHEETAH runs for use in the empirical fitting procedure.The obtained JWL EOS parameters are summarized in Table III.For the temperature calculation, the following simplified form of the JWL temperature equation is adopted, 18 where we assume that the specific heat c v is constant.
When describing the detonation process induced by the reaction of the energetic material, we consider the energetic material to be comprised of the solid-state reactant and the gaseous product.In this situation, the properties of the material at any point are determined by the mixing of the material properties of these two states.According to Lee and Tarver, 1 equilibrium pressure conditions are typically assumed when using the JWL EOS.This means that the temperature and internal energy of the reactant and the product are different.Based on this concept, we adopted the method proposed by Clutter et al. to calculate the equilibrium pressure. 19Using Eqs. ( 23), (27), and (28), the internal energy in the unreacted and the reacted state can be expressed as follows.
Then, the total internal energy is calculated from the product mass fraction λ as follows.
As a result, the equilibrium pressure becomes In Eq. ( 26), ρ λQ denotes the source energy associated with the reaction of the HMX and ρ λHS Q HS represents the reaction energy released owing to hot spots similar to those described in Ref. 20: Here, P ac denotes the criterion for temperature increase in the Estane.If the pressure increases above P ac , the Estane starts to exhibit pressure reverberations.In our simulation, P ac = 0.75 GPa and λ HS represents the fraction of hot spots that react, starting from 1 and reaching 0 when a hot spot finally reacts.In the above equation, λHS represents the reaction rate and ξ denotes the induction progress variable, which has an initial value of 0. This variable is associated with the induction time and can be tracked using the following model.
where the induction rate is The induction time τ is the time required to create the hot spot following the shock passage and represents the time between the delivery of the shock and the attainment of peak temperature.Its values for the different hot-spot sizes are τ 1µm = 1.5ns, τ 3µm = 2.4ns, and τ 5µm = 3.4ns.The heat of the reaction of a hot spot can be defined as follows.
Here, c v represents the constant-volume specific heat of the HMX, whereas T peak is the maximum temperature achieved when the shock has passed through the gap.The reference temperature is denoted by T s .We utilized T peak data obtained from the mesoscale simulation.The computational domain of the rate-stick test is shown in or between the binder walls is examined using SEM images. 21The average distance between binders and its deviation are and are then entered into the random distribution function of the hydrocode.For instance, if five meshes represent the average distance between the binders surrounded by the HMX, the distance with the deviation is expressed as 4 or 6 meshes.Through this process, the hot-spot elements shown in the SEM image can be represented realistically in the hydrocode.Simulations are performed with the following conditions: input pressure of 1.62, 3.04, 3.86, 4.74, 7.80, and 10.10 GPa and corresponding impact velocity of 300, 500, 600, 700, 1000, and 1200 m/s, respectively.The left, right, and upper boundaries have outflow boundary conditions.

C. Results and discussion
Figure 12 shows the grid resolution test results for the detonation peak pressure at t = 10.5 µs.
For the low-resolution case, we observe an irregular tail structure due to the insufficiently resolved computational mesh.The figure also shows how numerical artifacts or tail oscillations disappear with the refined grid; as the grid size increases, the position and pressure of the detonation peak are converged.
Figure 13 displays simulation results that show the difference in the shock-initiated process when hot spots are considered and not considered.It is observed that the shock wave is propagated, following a 700 m/s impact at the lower end of the explosive.At the initial stage of propagation, the SDT does not occur immediately.Noticeably, an empirical model, such as the ignition-and-growth model, would consider this transition to occur automatically at the beginning of the process.Moreover, our approach uses the Arrhenius rate with the mesoscale hot-spot model to calculate the entire SDT process with the method described earlier.As shown in Fig. 13(b), an initiation failure is observed if no hot-spot model is considered into the hydrocode simulation.As the shock passes through the hot spots, their reaction heat is released, and a vast energy release is surmounted at the shock front instantaneously.This process of hot-spot-initiated SDT is shown in Fig. 13(a), which presents the development of a full detonation wave with randomly driven heat-release sites created by hot spots.We also investigate whether the transition to detonation occurs for various input pressures.Six input pressures were considered and the results are summarized in Table IV.The results confirmed that the transition to detonation does not occur for an input pressure lower than 3.86 GPa.The input pressure considered is higher than P ac , which can cause local reaction at the potential hotspot site, resulting in an energy increase.However, this increase in energy is insufficient to cause the transition to detonation through the superposition of the energies emitted from the hot spots, resulting in shock-pressure dissipation.For input pressures greater than 4.74 GPa, SDT did occur.
Regarding the initiation characteristics of these SDT cases, the measured distance to detonation is shown in Fig. 14 and Table V. Figure 14 shows that the higher the input pressure is, the faster the initiation occurs.For an input pressure of 4.74 GPa, the initiation occurs only after the shock has passed for a considerable time.After initiation, as shown in Fig. 13, the hot-spot reactions are superimposed, and the resulting energy is transferred to a shock front.Subsequently, the transition to detonation occurs.Moreover, the number of hot spots in the initiation stage also varies greatly depending on the magnitude of the input pressure.
Our distance-to-detonation results are compared with literature values (experimental data for 95% HMX 22 ) for the same loading pressure in Table V.In the simulation, the distance to detonation was considered as the front point of the shock, where the detonation pressure was measured (Fig. 14).The distance to detonation is reproduced according to the loading pressure.
Figure 15 shows the fully developed detonation wave using pressure and temperature contours.The present simulation shows structural complexity with irregular flow structures in the temperature profile due to the existence of multiple energy release sites associated with hot spots.The experimental, theoretical, and calculated detonation properties are also compared.The detonation velocity was measured from the rate-stick pins, while both the temperature and the pressure were obtained from the thermo-chemical equilibrium code. 23The detonation velocity was found to be 8720 m/s.The development of the pressure profile along the vertical direction is plotted in Fig. 16.The irregular pressure structure generated in the ignition and transition stages changed into a smooth pressure structure in the fully developed detonation.The detonation pressure was 35.1 GPa, while the corresponding temperature was 3284.2 • K and the detonation density was 2563 kg/m 3 .A comparison between the experimental and calculated values is shown in Table VI.Besides a slight underestimation of the velocity, the comparison suggests that the method proposed in this paper is adequate for describing the complex SDT phenomenon using the hot-spot model.

V. CONCLUSION
We have proposed a model that is distinct from currently available SDT models, in which we separate the mesoscale hot-spot evolution for the SDT from the macroscale Arrhenius-type detonation of a bulk explosive.Thus, the particle-based SPH method identifies the sites of high energy release from the formation of hot spots, and the results are then used to initialize the full-scale hydrodynamic simulation to describe the detonation.Although all pressure-based ignition-and-growth models are excellent for describing a fully developed detonation, the ignition part of the SDT process can be handled more precisely by our model, as the growth process depends on the amount of accumulated energy released from the various hot-spot sites in a bulk sample.
Additionally, we describe a method of constructing a temperature-based detonation rate law (Arrhenius) for any unknown explosive composition using a laboratory-scale calorimeter, which necessitates only a small amount of sample mass.The obtained rate law accompanies the mesoscale particle hydrodynamic simulation of the shock-induced ignition, providing go/no-go criteria for given impact stimuli.As a result, we were able to accurately reproduce the velocity of the detonation for the tested explosive sample.
FIG. 5. Comparison between reaction progress obtained from the experiment and from the simulation.

FIG. 8 .
FIG. 8. Schematic of a mesoscale impact simulation via SPH.The gap distance is measured in Estane binder bounded by two HMX walls.

Fig. 11 .
FIG. 11.Simulation domain of the rate stick (half domain) with random hot spots.

TABLE I .
Material properties and Mie-Gruneisen factor for 95% HMX.

TABLE II .
Maximum temperature for each Estane gap size, calculated with SPH.

TABLE IV .
SDT results according to input pressure conditions.

TABLE V .
Comparison of distance to detonation with previous experimental results for various input pressures.

TABLE VI .
Comparison between experimental (rate-stick test) and simulated detonation properties.