## Abstract

For condensed explosives, containing metal particle additives, interaction of the detonation shock and reaction zone with solid inclusions leads to high rates of momentum and heat transfer that consequently introduce non-ideal detonation phenomena. During the time scale of the leading detonation shock crossing a particle, the acceleration and heating of metal particles are shown to depend on the volume fraction of particles, dense packing configuration, material density ratio of explosive to solid particles and ratio of particle diameter to detonation reaction-zone length. Dimensional analysis and physical parameter evaluation are used to formalize the factors affecting particle acceleration and heating. Three-dimensional mesoscale calculations are conducted for matrices of spherical metal particles immersed in a liquid explosive for various particle diameter and solid loading conditions, to determine the velocity and temperature transmission factors resulting from shock compression. Results are incorporated as interphase exchange source terms for macroscopic continuum models that can be applied to practical detonation problems involving multi-phase explosives or shock propagation in dense particle-fluid systems.

## 1. Introduction

Detonation in condensed explosives typically features wave propagation speeds from 6 to 9 mm μs^{−1} and peak pressures from 10 to 50 GPa. The detonation front is often modelled by the von Neumann (VN) shock followed by a reaction zone in the idealized one-dimensional Zel'dovich–von Neumann–Döring (ZND) wave (Fickett & Davis 1979), which travels at a minimum speed for the Chapman–Jouguet (CJ) condition where a sonic point terminates the reaction zone. The addition of heterogeneities such as metal particles into condensed explosives introduces microscopic interaction of the detonation shock and reaction zone with the particles, which produces localized hot spots, transverse waves and high-pressure fluctuations that are manifested in macroscale detonation phenomena.

Particles saturated in liquid explosives form a fundamental system termed wetted powder, or ‘slurry’ (Frost & Zhang 2009). The effect of a large volume fraction of metallic particle additives on the detonation properties of slurry explosives has been investigated by a number of researchers, e.g. Engelke (1979), Baudin *et al.* (1998), Gogulya *et al.* (1998), Haskins *et al.* (2002), Zhang *et al.* (2002) and Kato *et al.* (2006), where the bulk propagation speed has been related to the particle properties, morphology and concentration. In general, adding particles results in a velocity deficit below that of the neat explosive, because some of the chemical energy released goes into heating, and particularly, accelerating the particles.

Detonation failure depends on the competing effects of explosive sensitization owing to the formation of hot spots near the particles, and the increasing detonation reaction-zone length with momentum and energy absorbed into the particles and detonation front curvature caused by lateral expansion. Engelke (1979) and Lee *et al.* (1995) extensively studied the critical diameter of nitromethane (NM) containing silica glass beads; Kurangalina (1969) and Frost *et al.* (2005) studied NM detonation failure with aluminium powders.

Particle ignition delay and reaction times are typically greater than the explosive detonation reaction time scale, and the resulting metal particle combustion heat release occurs predominantly behind the sonic point in the detonation. Baudin *et al.* (1998) suggested that for spherical aluminium particles as small as 100 nm, there is insufficient time for particles to react within the detonation reaction zone. Kato *et al.* (2006) showed that aluminium reaction contributes to the NM detonation only for particles smaller than 2 μm. Furthermore, the results of Yoshinaka *et al.* (2007) for 30 μm aluminium particles shock-compressed to 13–30 GPa in condensed matter showed that the majority of particles did not melt. Thus, micrometric metal particles can usually be considered inert and intact within the condensed detonation reaction zone, and the interaction between the particles and explosive is dominated by shock compression.

Shock interaction and related momentum and heat transfer from the explosive to the particles within the detonation zone are important mechanisms associated with macroscopic detonation initiation, wave propagation, stability and failure phenomena. On the other hand, the detonation transmits a strong shock into the solid particles that rapidly accelerates and heats the metal as the wave passes. Internal wave reflections and external interactions with neighbouring particles dominantly affect the particle velocity and temperature owing to shock compression before viscous flow interaction takes over. While experimental methods can provide information on the bulk detonation response, at present the diagnostics do not have the resolution required to record the individual particle behaviour. Mesoscale numerical simulation is a practical alternative to gain insight into detonation at the grain scale (0.1 μm–1 mm) in condensed heterogeneous matter, where very high resolutions are required to capture the small geometrical features and shock interaction.

Ordered matrices of packed circles, spheres and simple polygons have typically been used to represent layers of heterogeneous condensed matter in mesoscale simulation. Shock interaction with cylindrical ‘particles’ has been simulated in two-dimensions by Mader (1979), Milne (2000) and Zhang *et al.* (2003). Early three-dimensional models of Mader & Kershner (1982) have evolved to modern calculations of Baer (2002) and Cooper *et al.* (2006) that typically employ at least 10^{7} computational cells and contain (100) particles. In general, only a very small material sample (up to millimetre size) can be simulated using mesoscale simulation.

While much of the foregoing literature investigates the influence of heterogeneities on the detonation wave structure and propagation velocity, limited work has addressed the acceleration and heating imparted on metal particles in condensed matter during explosive detonation. Zhang *et al.* (2003) showed that light metal particles achieve 60–100% of the shocked fluid velocity, and demonstrated that the material density ratio of the explosive to solid particle is an important factor affecting particle acceleration. The influence of other parameters affecting particle acceleration and heating have been studied, including particle matrix properties and explosive reactivity (Ripley *et al.* 2005, 2006, 2007).

In the present work, a quantitative description of the resultant momentum and heat transfer is sought in conjunction with determination of the principal shock interaction mechanisms. This is achieved using a theoretical dimensional analysis to identify key parameters, applying three-dimensional mesoscale continuum modelling of packed particle matrices saturated in liquid explosive to understand the mechanisms and behaviour of the key parameters, and by compiling results into transmission factors that quantify the momentum and heat transfer from the explosive to the particles. Finally, the transmission factors are incorporated into momentum and heat exchange source terms developed for a macroscopic fluid dynamics framework, which is suitable for modelling detonation shock compression in metal particle-condensed explosive systems.

## 2. Regimes for detonation interaction with particles

A natural starting point for condensed explosives is to consider a system of a liquid explosive containing a high volume fraction of metal particles. Liquid explosives are considered homogeneous in the absence of additives, and detonation of liquid explosives has a reaction length (*L*_{R}) ranging from a few microns to a few millimetres. Sheffield *et al.* (2002) estimated the reaction zone in NM to be about 300 μm, whereas Engelke & Bdzil (1983) predicted a NM reaction length of 36 μm. The detonation reaction length can be reduced to a few micrometres when NM is sensitized by liquid diethylenetriamine or triethylamine. Nitroglycerin is a liquid above 13^{°}C, with a reaction-zone length of 210 μm (Dobratz & Crawford 1985). Isopropyl nitrate is another liquid explosive with an estimated reaction-zone length on the order of 1 mm (Zhang *et al.* 2002). Trinitrotoluene (TNT) is a liquid above 81^{°}C, where the reaction-zone length is 0.9–1.1 mm (Dobratz & Crawford 1985). Thus for liquid explosives, a range for detonation reaction-zone length of 10^{−6}<*L*_{R}<10^{−3} m is assumed.

Various metal particles have been added to liquid explosives to form a dense heterogeneous system. This includes magnesium, aluminium, titanium, steel, copper and tungsten in a size range of 10^{−8}<*d*_{p}<10^{−3} m (Zhang *et al.* 2001; Frost *et al.* 2002, 2005; Haskins *et al.* 2002; Kato *et al.* 2006), where *d*_{p} denotes a spherical particle diameter. A parameter relating the particle diameter and reaction length scale can therefore be defined as *δ*=*d*_{p}/*L*_{R}. A potential range of 10^{−5}<*δ*<10^{3} can be derived from the ranges of *d*_{p} and *L*_{R} as justified earlier. Three regimes for detonation interaction with particles can be identified according to the value of *δ*, as illustrated schematically in figure 1. In reality, the leading shock structure is three-dimensional with transverse waves and transmission into particles, followed by a multitude of shocks reflected from particles.

### (a) Small particle limit (*d*_{p}/*L*_{R}≪1)

At the limit of , the detonation shock front is considered inert (i.e. the von Neumann shock) during the early interaction, which can then be represented by a Heaviside step function. Let us define a shock interaction time scale as the time required for the detonation front to cross a particle: *τ*_{S}=*d*_{p}/*D*_{0}, where *D*_{0} is the shock velocity. Within the shock interaction time, the detonation reaction-zone length is no longer a parameter and the response is represented by a single length scale of the particle diameter. For an inert planar shock crossing a particle, the dynamic response of the particle and the surrounding flow field at any given time can be scaled by the particle diameter using geometrical similarity when employing inviscid governing equations and rate-independent material models (Zhang *et al.* 2003). Thus, the computational results for a system of liquid explosive containing particles of a given size can be scaled to systems of the same liquid explosive with any diameter particles within the small particle limit.

### (b) Large particle limit (*d*_{p}/*L*_{R}≫1)

For , the reaction length becomes negligibly small and the detonation wave can therefore be considered as a discontinuity of the CJ front, separating the fresh explosive from its detonation products. The interaction consists of diffraction of a thin CJ detonation front dictated by the curved boundary of the particle, followed by unsteady expanding products flow controlled by the rear boundary. The particle acceleration and heating are then characterized by the single length scale of the particle diameter, similar to the small particle limit (*d*_{p}/*L*_{R}≪1), whereas the Heaviside VN shock is replaced with the CJ shock and expanding products rear flow.

### (c) Intermediate regime (*d*_{p}/*L*_{R}∼1)

The case of *δ*=*d*_{p}/*L*_{R}∼1 lies between the above two limits and is therefore most complex because both the particle diameter and reaction-zone length scales play a role in the particle acceleration and heating. In this regime, the particle interacts with both the VN shock front and the expanding reacting flow in the detonation reaction zone that terminates at a sonic point. Locally, the reaction zone is affected by the particle presence, resulting in a decreased reaction-zone length at hot spots and an increased reaction-zone length in the expansion flow around the particle.

### (d) Definition of transmission factors

In order to describe the effect of the acceleration of solid particles in shock and detonation of a condensed explosive, a velocity transmission factor, *α*, is defined as the ratio of the particle mass-averaged velocity, *u*_{p}, after an interaction time, *τ*, over the shocked fluid velocity, *u*_{f1}:
2.1

In general, the velocity transmission factor varies between 0 for perfect reflection off a rigid body and 1 for perfect transmission into a particle with the same material properties as those of the host fluid. Similarly, a temperature transmission factor, *β*, is defined as the ratio of the particle mass-averaged temperature, *T*_{p}, after an interaction time, over the shocked fluid temperature, *T*_{f1}:
2.2where 0≤*β*≤1 for most metals. In (2.1) and (2.2), with *τ*_{S}=*d*_{p}/*D*_{0}. The characteristic shock interaction time, *τ*_{S}, is used for a single particle in condensed matter. For a dense solid particles-fluid system, transmission factors are measured after *τ*=2*τ*_{S} such that the immediate effect of wave reflections both within particles and in the voids between neighbouring particles is included. This comprises the majority of the acceleration and heating owing to primary shock transmission, while subsequent internal waves further influence the final velocity and temperature achieved during shock compression. This time frame also accounts for the influence of transverse wave reflections and upstream/downstream particle reflections in a densely packed matrix.

## 3. Factors affecting particle acceleration and heating

### (a) Dimensional analysis

The force acting on a particle during detonation of a liquid explosive containing dense solid particles is assumed to be a function of relevant dimensional parameters, e.g. *F*_{d}=*f*(*d*_{p},*D*_{0},*ρ*_{f0},*ρ*_{s0},*ϕ*_{s0},*p*_{0},*μ*_{f0},*L*_{R}), where *d*_{p} is the particle diameter, *D*_{0} is the detonation velocity, *ρ*_{f0} is the initial density of the explosive fluid, *ρ*_{s0} is the initial particle material density, *ϕ*_{s0} is the initial solid particle volume fraction, *p*_{0} is the ambient pressure, *μ*_{f0} is the molecular viscosity of the explosive and *L*_{R} is the detonation reaction-zone length. The Buckingham *π* theorem (Buckingham 1914) considers primary dimensions (length, mass, time and temperature) to reduce the number of parameters describing a physical behaviour into a minimum set of dimensionless groups. Assuming *d*_{p}, *D*_{0} and *ρ*_{f0} are independent variables among the nine parameters in the force expression, *π* theorem yields that there are six dimensionless groups, which are as follows:
3.1

In (3.1), *C*_{d} is called an ‘effective’ drag coefficient, *Re* is the particle Reynolds number and *M*_{0} is the Mach number of the detonation shock (*M*_{0}=*D*_{0}/*a*_{f0}, where *a*_{f0} is the sound speed). Note that the flow compressibility is represented by the shock Mach number, instead of the flow Mach number commonly used. Thus, the force expression can be rewritten in a non-dimensional form:
3.2

Similarly, heat transfer to a particle during detonation in a solid particle-explosive system is assumed to be *Q*_{c}=*f*(*d*_{p},*D*_{0},*ρ*_{f0},*ρ*_{s0},*ϕ*_{s0},*T*_{0},*μ*_{f0},*k*_{f0},*c*_{p},*L*_{R}), where *T*_{0} is the ambient temperature, *k*_{f0} is the thermal conductivity of the explosive and *c*_{p} is the explosive fluid heat capacity. Using *d*_{p}, *D*_{0}, *ρ*_{f0} and *T*_{0} as the independent variables, the *π* theorem gives seven dimensionless groups:
3.3

In (3.3), *Nu* is an ‘effective’ Nusselt number (*Nu*=*hd*_{p}/*k*_{f}, where *h* is the convective heat transfer coefficient) and *Pr*=*μ*_{f}*c*_{p}/*k*_{f} is the Prandtl number. Finally, the non-dimensional form of the heat transfer equation can be expressed by
3.4

The same non-dimensional parameters can also be obtained from the mass, momentum and energy conservation equations that govern the flow. Defining dimensionless variables , *t**=*t*/*τ*_{S}=*tD*_{0}/*d*_{p}, , *ρ**_{f0}=*ρ*_{f}/*ρ*_{f0}, *ρ**_{s}=*ρ*_{s}/*ρ*_{s0} and *ϕ**_{s}=*ϕ*_{s}/*ϕ*_{s0}, and substituting them into the continuity equation for the solid particles-explosive system, one obtains
3.5which shows the non-dimensional parameters of *ϕ*_{s0} and *ρ*_{f0}/*ρ*_{s0}.

Let additional dimensionless variables be defined as , *σ**=*σd*_{p}/*μ*_{f0}*D*_{0}, *μ**=*μ*/*μ*_{f0}, , , , *q**=*qd*_{p}/*k*_{f0}*T*_{0}, where *F*_{0} is a reference force, *E* is the total specific energy and *Q*_{R} represents a chemical reaction energy source. Substituting all the dimensionless variables into the three-dimensional Navier–Stokes equations for the fluid phase of the particles-explosive system, one derives the linear momentum equations and energy equation of the fluid phase in dimensionless form. For example, the energy equation becomes
3.6

Thus, (3.5) and (3.6) provide the non-dimensional parameters of *Re*, *Pr*, *M*_{0}, *ϕ*_{s0}, *δ*, *C*_{d} and *Nu*. The same *Re*, *M*_{0}, *ϕ*_{s0} and *C*_{d} are also derived for the momentum equation, which has been omitted for brevity.

For dilute particles-gas flow (*ϕ*_{s0}≪1; *ρ*_{f0}/*ρ*_{s0}≪1) with a small particle-limiting scale (*d*_{p}/*L*_{R}≪1), the drag coefficient function (3.2) and Nusselt number function (3.4) approach the well-established classic forms of *C*_{d}=*f*(*Re*,*M*_{0}) and *Nu*=*f*(*Re*,*Pr*,*M*_{0}). The dimensional analysis for dense solid particles-reactive fluid systems shows the additional parameters, including the density ratio of fluid to solid particle, solid volume fraction and ratio of particle diameter to detonation reaction-zone length.

### (b) Discussion of the parameters

This dimensional analysis is aimed at studying the particle acceleration and heating during detonation shock compression in a dense solid particle-condensed explosive system. There are two major forces acting on a solid particle during the detonation process: shock compression and viscous drag, with corresponding time scales characterized by the shock interaction time, *τ*_{S}, and velocity relaxation time, *τ*_{V} (given in Rudinger (1980)), respectively. The ratio of shock interaction time to viscous relaxation time is
3.7where *C*_{d} denotes a classic drag coefficient. For particles larger than 0.1 μm, the viscous relaxation time is an order of magnitude greater than the shock interaction time for a typical example of aluminium particles-liquid NM system (*ρ*_{s}=2785 kg m^{−3} for inert aluminium, *ρ*_{f,CJ}=1539 kg m^{−3}, *u*_{f,CJ}=1.764 mm μs^{−1}, *D*_{CJ}=6.612 mm μs^{−1} for NM detonation, and assuming *C*_{d}∼1 for an array of spheres in compressible flow). Thus, shock compression is the dominant force for the particle acceleration during the detonation process and the viscous drag can be assumed to be negligible. Consequently, the Reynolds number can be considered less important in the effective drag equation (3.2) and the flow can be assumed to be inviscid.

Similarly for heat transfer, the ratio of shock interaction time to thermal relaxation (convection) time, *τ*_{T} (given in Rudinger (1980)), is
3.8

For detonation of the inert aluminium particles-liquid NM system, the thermal relaxation time is two orders of magnitude greater than the shock interaction for *d*_{p}=0.1 μm and at least three orders of magnitude greater than the shock interaction time for *d*_{p}>1 μm. Therefore, shock compression heating controls the particle temperature during the detonation process and convective heat transfer can be assumed negligible. Consequently, the Reynolds number and Prandtl number can therefore be considered less important in the Nusselt number equation (3.4), and the flow can be assumed to be completely inviscid and non-heat-conducting.

The detonation shock Mach number of solid and liquid explosives at their theoretical maximum density has a narrow range of 2.5<*M*_{0}<4 in general. Shock velocity and pressure effects were investigated in a previous work (Zhang *et al.* 2003) by varying the inert shock pressure from 5 to 20 GPa. While the resulting momentum transfer to particles remained proportional to the shocked fluid velocity, the variation in velocity transmission after the shock interaction was less than 10 per cent for a shock velocity range of 4–9 mm μs^{−1} for metal particles in cyclotrimethylenetrinitramine (RDX) explosives.

For a step shock wave passing a spherical metal particle in condensed matter, the velocity transmission factor, *α*, was studied for magnesium, beryllium, aluminium, nickel, uranium and tungsten in liquid NM and various solid RDX densities subjected to a shock pressure range of 5–20 GPa (Zhang *et al.* 2003). The particle velocity after the shock interaction time, *τ*_{S}, was found to strongly depend on the initial density ratio of explosive to metal and can be expressed by
3.9where *a* and *b* are constants independent of the particle and explosive matter.

The volume fraction of solid particles is defined by *ϕ*_{s}=*n*_{p}*m*_{p}/*ρ*_{s}, where *n*_{p} is the number density of particles and *m*_{p} is the mass of an individual particle, and has a range of 0≤*ϕ*_{s}≤1. For example, in Tritonal (80 wt% TNT, 20 wt% Al), the solid volume fraction is about *ϕ*_{s0}≈0.13. Loose polydisperse powders typically have *ϕ*_{s0}≈0.5; random packing in particle beds produces higher volume fractions, for instance, *ϕ*_{s0}=0.58–0.62 in slurries containing Al particles saturated in liquid NM (Haskins *et al.* 2002; Frost *et al.* 2005). For ordered packing of mono-sized spheres, *ϕ*_{s0}=0.74 is a theoretical geometrical maximum for a close-packed matrix. The particle acceleration and heating depends on the solid volume fraction and packing configuration.

For inert shock interaction, there is only a single characteristic length scale, i.e. the particle diameter, *d*_{p}. The detonation process takes place in a reaction-zone length, *L*_{R}, which is an additional length scale that influences the particle acceleration and heating. Therefore, the ratio *δ*=*d*_{p}/*L*_{R} appears in equations (3.2) and (3.4) and has been used to characterize three regimes of the detonation interaction (figure 1). The effect of the volume fraction of solid particles and the ratio of particle diameter to detonation reaction-zone length is investigated in this paper using a mesoscale modelling approach.

## 4. Mesoscale computational approach

A system of spherical Al particles in NM forms the prototype heterogeneous mixture for the study of detonation interaction with particles. Nitromethane (CH_{3}NO_{2}) is a uniform, low-viscosity liquid explosive and is therefore assumed to follow ZND detonation theory. The incident shock pressures encountered in the NM detonation are 13–23 GPa, far exceeding the yield strength of aluminium, which is 0.035 GPa for 99 per cent commercially pure aluminium (Callister 1997). Therefore, the aluminium material strength has been neglected by assuming that only volumetric strain occurs in the particles. Hence, neglecting viscosity and thermal conductivity as justified earlier, the hydrodynamic response of both liquid NM and solid aluminium particles can be computed using the three-dimensional inviscid Euler equations.

The governing equations are solved using the Harten, Lax, van Leer with contact correction (HLLC) approximate Riemann solver scheme (Batten *et al.* 1997). The three materials are tracked using mass fraction continuity equations for the liquid NM explosive, NM gaseous detonation products and solid Al particles treated as inert, with equations of state for each described below. The mixture of unreacted explosive and detonation products within the detonation zone and the material boundaries at the metal particle surface are treated using continuum mixture methods following Benson (1992).

### (a) Equations of state

The NM and the aluminium particles are modelled using the Mie-Grüneisen equation of state (EOS):
4.1which is an expansion approximation around the shock Hugoniot. In (4.1), *ν* is the specific volume, *e* is the specific internal energy and *Γ*_{S} is the so-called Grüneisen gamma, which is the negative of the log slope along an isentrope. The Hugoniot states (subscript H) are determined using the pressure relation:
4.2in which the coefficients *C* and *S* are determined from experiments (Marsh 1980). The shock Hugoniot parameters are fit to the linear approximation *D*_{0}=*C*+*Su*_{f1}, where *D*_{0} is the shock velocity and *u*_{f1} is the fluid velocity.

The temperature-density behaviour is modelled using fitting to the data of Walsh & Christian (1955): 4.3

The temperature fitting coefficients (*F*_{S}, *G*_{S}, *H*_{S}, *I*_{S} and *J*_{S}) and Hugoniot parameters used for NM and aluminium are summarized in table 1, which were selected from those available for several common materials (Mader 1998).

The expansion of the gaseous detonation products is represented by the Jones–Wilkins–Lee EOS (Lee *et al.* 1968):
4.4

The fitting coefficients *A*=277.2 GPa, *B*=4.934 GPa, *C*=1.223 *GPa*, *R*_{1}=4.617, *R*_{2}=1.073 and *ω*=0.379 are calculated using the Cheetah thermochemical code and Becker–Kistiakowsky–Wilson–Sandia library (Fried *et al.* 1998).

### (b) Nitromethane reaction model

The NM detonation model follows Mader (1998), who simulated reaction-zone lengths ranging from 0.24 to 70.5 μm. For this work, the length of the reaction zone has arbitrarily been fixed and the particle diameter is varied to study a range of *δ*=*d*_{p}/*L*_{R}. A single-step Arrhenius reaction model was employed for the NM detonation,
4.5with *Z*=8.0×10^{12} s^{−1} and *E*_{a}=1.672×10^{5} J mol^{−1} K^{−1}. A heat of detonation of Δ*H*_{det}=5.725 kJ cm^{−3} was determined using the Cheetah chemical equilibrium code (Fried *et al.* 1998). The resulting reaction-zone length is 2 μm when measured from the VN spike to the sonic point (location where *D*=*u*_{f}−*a*_{f}). In this model, the reaction is 99.9 per cent complete at the CJ point. A comparison of the numerical CJ point with the Cheetah equilibrium results shows reasonable agreement. Figure 2 shows the resulting detonation wave pressure profiles for a closed boundary condition at *x*/*d*_{p}=0. Other more sophisticated kinetics schemes and modern equations of state are established in the literature, e.g. Lysne & Hardesty (1973), Hardesty (1976) and Winey *et al.* (2000); however, the present approach is still adequate to study the mechanical and thermal interaction of particles with shock and detonation flow.

Prior to the detonation wave entering a packed particle matrix, the detonation is run out to a distance of 1 mm (five times the running distance depicted in figure 2). This was performed to reduce the gas expansion rate in the Taylor wave during the interaction with the metal particles. Numerical simulation of the NM detonation was performed on a high-resolution mesh with a 10^{−8} m cell size, giving three to five cells across the shock and 200 cells in the reaction zone for *L*_{R}=2 μm. A steady detonation was obtained after a running distance of 0.1 mm. A selected spatial wave distribution was used to initialize three-dimensional meshes of the same resolution containing the particle bed. For *d*_{p}=1 μm, the particles are resolved with 100 cells per *d*_{p}, which is sufficient to provide a grid-independent solution.

### (c) Initial particle packing configurations

Realistic heterogeneous explosive mixtures feature random packing and non-spherical particles. To fundamentally address the acceleration and heating of metal particles during detonation, ordered matrices of uniform spherical particles are considered for various inter-particle spacings, covering a range of volume fractions from 0.093≤*ϕ*_{s0}≤0.740. The regular spacing and monodisperse size assumptions avoid dense particle–particle collisions because all particles are subjected to the same acceleration processes. Three idealized packing configurations were considered as illustrated in figure 3; however, the face-centred cubic lattice arrangement was selected for most of the calculations because the closest packing represents the saturation limit, whereas the solid volume fraction can be varied by adjusting the spacing between particles. Candidate particles are studied six to eight layers into the matrix, where the interaction has become quasi-steady in the absence of starting and end effects (Ripley *et al.* 2006, 2007). The dilute limit () is simulated using a two-dimensional axi-symmetric model of a single spherical particle, whereas the dense limit () is simulated using a one-dimensional model of a semi-infinite solid slab. In three-dimensions, plane-symmetric domains containing 20–40 particles were resolved using up to 32 million mesh points and 100 central processing units in parallel. Reflective boundary conditions were used to represent a periodic solution in a semi-infinite array of ordered particles.

## 5. Mesoscale results

### (a) Results for the small particle limit (*d*_{p}/*L*_{R}≪1)

Figure 4 illustrates the particle deformation in a three-dimensional particle matrix (*d*_{p}=10 μm, ) resulting from a 10.1 GPa inert Heaviside step shock travelling from left to right. The deformed particles resemble a saddle shape and, owing to the inviscid hydrodynamics, are strongly influenced by the complex shock reflections from neighbouring particles. The deformation of the first layer is different because the reflected shock wave is not subsequently re-reflected from upstream particles. The severe deformation in the rear flow indicates that the aluminium particles will likely be damaged or fragmented if material shear stress and failure are considered.

Figure 5*a* shows the pressure histories at the leading edges of the first eight layers of packed particles. The pressure reaches a peak at the perfect reflection value of 45.4 GPa and rapidly expands to a sustained quasi-steady pressure plateau centred below the one-dimensional wave transmission value of 20.9 GPa. The reverberating pressure oscillates with a period of 0.5*τ*_{S} corresponding to the wave transit time within the interstitial pores contained between packed particles. These resonant fluctuations would be dampened with random packing and non-spherical particles as shown by Baer (2002).

Figure 5*b* shows the mass-averaged particle velocity for three packing configurations. Within 1*τ*_{S}, the velocity in all three matrices exceeds the one-dimensional wave transmission value of 1.11 mm μs^{−1} owing to the internal rarefaction as the wave exits the trailing edge of the particle. Subsequent shock reflection and interaction with neighbouring particles causes velocity fluctuations proportional to the particle spacing that vary with packing configuration. A time scale of 2*τ*_{S} is sufficient to capture one full interaction cycle that includes the internal wave reverberation and successive expansions and compressions from upstream, downstream and neighbouring particles. Because severe particle deformation and resonant fluctuations occur after the shock interaction time scale in the rear flow, they do not influence the evaluation of the velocity and temperature transmission factors.

Figure 6 provides a summary of the velocity and temperature transmission results for the small particle limit, using a 10.1 GPa inert shock. Comparison of the packing configurations shows that the transmission factors for the simple cubic packing were considerably different. This is primarily due to this particular packing matrix that provides two propagation channels: one through the linear array of stacked particles, and the other uninhibited though the column of liquid in the void space. The simple cubic packing is an unrealistic configuration in practice and is not given further consideration. The present three-dimensional results are compared with two-dimensional cylindrical results of Zhang *et al.* (2003). The agreement is improved for a higher number of two-dimensional cylinders, as this better approximates an infinite array of particles as in the three-dimensional configuration.

Both the velocity and temperature transmission factors for the close-packed and body-centred matrices follow an inverted U-shaped function. The maximum transmission factors occurred for 0.25≤*ϕ*_{s0}≤0.40, whereas the minimum transmission factors resulted at the dilute limit (*ϕ*_{s0}=0) and high volume fraction limit (*ϕ*_{s0}=1). In between the volume fraction limits, the spacing of the particles affects the arrival time and magnitude of reflected shocks from upstream and downstream particles, in addition to affecting the local flow diffraction. The maximum transmission factors are a result of superposition of transmitted shock waves that act in a coherent manner.

### (b) Results for the large particle limit (*d*_{p}/*L*_{R}≫1)

In the *d*_{p}/*L*_{R}≫1 case, early simulation results (Ripley *et al.* 2005) for the detonation wave over a single aluminium particle indicated that the particle velocity is first increased as the detonation front crosses the particle and is then reduced, subject to the Taylor expansion flow conditions. The Taylor expansion after the thin detonation front reduces the detonation flow velocity, thus changing the direction of the drag and reversing the momentum transfer direction from the particles to the detonation expansion flow. The large particle limit was studied preliminarily using *d*_{p}=30 μm (Ripley *et al.* 2006), which corresponds to *d*_{p}/*L*_{R}=15. The unsteady Taylor expansion effect was minimized by running the detonation sufficiently far from the initiation location.

For *d*_{p}/*L*_{R}≫1, the VN shock can be neglected and the detonation interaction with particle matrices can be better represented by a CJ shock. The jump in the host liquid does not follow the Hugoniot because the shock is reactive, where the chemical reaction rate is essentially infinite and the post-shock flow contains hot expansion products. Figure 7 illustrates a CJ shock (*D*_{0}=6.69 mm μs^{−1}, *P*_{CJ}=13.8 *GPa*, *u*_{f1,CJ}=1.827 mm μs^{−1} and *ρ*_{f1,CJ}=1.551 *g* *cm*^{−3}) travelling through a NM–Al matrix. For a particle spacing of 1*d*_{p}, the volume fraction is *ϕ*_{s0}=0.093, and the CJ shock front profile approaches that of a planar wave prior to arrival at successive particle leading edges; this leads to a flattening of the particles during deformation.

The particle velocity and temperature histories are shown in figure 8. Within 2*τ*_{S}, both the incident shock and internal rarefaction accelerate the particle, while an increase in particle temperature owing to shock compression is followed by a decrease from the rarefaction expansion before lateral compression continues the particle heating. Figure 8*a* illustrates a decrease in particle velocity over 3–4 *τ*_{S} as a result of a reflected wave returning from the downstream particles. Further particle acceleration must be influenced by viscous drag, which is beyond the scope of the present paper. Similarly, the particle heating is only affected by shock compression in the non-heat-conducting assumption valid within the time frame considered, although the downstream NM products are hot (*T*_{f1,CJ}=3657 *K*).

### (c) Results for the intermediate regime (*d*_{p}/*L*_{R}∼1)

Computation of diffraction of a reactive shock over a single particle in a three-dimensional mesh showed that the shock is both transmitted into the metal at the particle forward stagnation point, and also reflected back into the reaction zone, thereby increasing the reaction rate. A Mach stem forms as the shock in the explosive diffracts over the particle. Depending on the impedance ratio, the shock inside the metal particle may travel ahead of the incident shock, transmitting an oblique shock into fresh explosive ahead of the detonation front. Behind the diffracting detonation shock, expansion of the flow on the backside of the particle causes the local reaction-zone length to increase. At the particle trailing edge, either the diffracted shock arrives first, or the converging shock inside the particle is retransmitted back into the fresh explosive. This can initiate the explosive locally ahead of the detonation shock. Subsequently, a strong rarefaction forms within the metal particle, contributing to both a large acceleration and a decrease in temperature.

For packed particle matrices, the situation is more complex with reverberating waves and lateral interactions that further influence the particle and surrounding flow. Figure 9 illustrates the irregular propagation pattern of the detonation front in an Al–NM matrix (*d*_{p}=1 μm, *d*_{p}/*L*_{R}=0.5), where the detonation travels within the explosive contained in the voids and narrow channels between particles in addition to shock propagation within the metal. In close-packed matrices, transmission of the shock through the particle can subsequently pre-compress and initiate detonation in the fresh explosive on the far side of the particle prior to the diffracted shock arrival (hot-spot mechanism). For the conditions depicted in figure 9, the NM behind the particle trailing edge reaches 800 K owing to shock transmission and prior to reaction.

Figure 10 shows numerical pressure gauge results for detonation in a close-packed matrix, with the peak reflected VN shock pressure followed by oscillations at a frequency proportional to particle diameter and later by the Taylor wave expansion. Figure 10*a* shows a precursor shock with magnitude of about 10 GPa (temperature approx. 1300 K) that resulted from the reflection of the shock transmitted through the particle trailing edge into the void prior to arrival of the diffracted VN shock. Figure 10*b* shows a smaller precursor shock exiting the trailing edge and a larger peak pressure owing to collision of the diffracted shocks behind the particle trailing edge.

Figure 11 shows the particle velocity and temperature histories in a packed particle matrix for *d*_{p}/*L*_{R}=0.2, where the detonation reaction zone is much larger than the particle size. This is evident in the particle velocity history result, where the velocity first increases during interaction with the VN shock, and then decreases mainly over 5*τ*_{S} during the expansion inside the reaction zone. The 2*τ*_{S} assumption is especially justified in the temperature history where the first peak is reached between 1 and 2*τ*_{S}.

Figure 12*a* shows the detonation shock velocity through the particle matrices immersed in liquid explosive. The bulk propagation velocity was measured using wave time of arrival between consecutive numerical gauge stations located in an array perpendicular to the averaged detonation wavefront. Under the detonation conditions studied here, i.e. aluminium particles in NM, shock waves reflected from the leading edge of the particles cause an increase in the local NM density and pressure, thereby increasing the detonation velocity. Furthermore, shock waves travelling within the metal particle are transmitted into the fresh liquid explosive ahead of the detonation wave diffracting around the curved particle surface, thereby pre-compressing the explosive and increasing the local detonation velocity. These two local hot-spot factors both contribute to bulk propagation speeds in excess of the CJ value in fresh NM. In the present numerical calculations, the greatest bulk propagation velocities (up to 7.4 mm μs^{−1}) were observed for the highest metal mass fraction condition in combination with the smallest particle diameter (or longest reaction zone) within the region of the first particle layer. Afterwards, the momentum and energy transferred into the particles within the detonation zone competes with the local hot-spot factors, resulting in a quasi-steady propagation velocity with a deficit. The velocity deficit increases with an increase in solid volume fraction and a decrease in particle diameter. The maximum velocity deficit corresponded to a wave speed of 5.3 mm μs^{−1} and occurred for the small particle limit.

Figure 12*a* shows increasing instability in the detonation front velocity for higher solid volume fractions, which is an expected feature due to higher momentum and heat loss. The initially transient shock velocity becomes quasi-steady after travelling a distance of 6*d*_{p} into the matrix. Figure 12*b* illustrates the quasi-steady detonation shock velocity as a function of increasing metal mass fraction. In comparison with Cheetah chemical equilibrium predictions using inert aluminium, where velocity and temperature equilibrium is assumed for all phases (Fried *et al.* 1998) and essentially , the mesoscale results for detonation shock velocity are consistently higher because the relative velocity of the solid phase remains below the flow velocity following the shock interaction. Similarly, the Cheetah equilibrium temperature is also lower than the mesoscale detonation results owing to the same phase-non-equilibrium nature.

### (d) Particle acceleration

The mass-averaged particle velocity history achieved during detonation interaction with a single particle is illustrated in figure 13*a*. Because a fixed numerical mesh resolution was maintained (100 cells μm^{−1}), such that the number of cells in the ideal detonation reaction zone (*L*_{R}=2 μm) was unchanged, the number of cells across the particle diameter increased with *δ* and the corresponding computational effort increased exponentially. Thus for a large *δ*, limited duration particle velocity histories are available. For *δ*<0.5, the peak particle velocity exceeds the CJ value, illustrating a significant influence of the VN spike on the particle acceleration. For *t*/*τ*_{S}≫1, the particle velocity equilibrates below the CJ value (*u*_{CJ}=1.827 mm μs^{−1}) in the absence of Taylor expansion.

Figure 13*b* shows the mesoscale results of the single particle velocity and corresponding transmission factor as a function of *δ*=*d*_{p}/*L*_{R}. Both the velocity and the velocity transmission factor decrease from VN to CJ mainly over the interval from 0.1≤*δ*≤1. The single particle acceleration results are bound by the small particle limit () and the large particle limit (). The resulting velocity transmission factor was fit to the sigmoidal function:
5.1where *δ*_{0}=0.38 and *w*=0.25. In equation (5.1), and .

Figure 14*a* illustrates the particle velocity transmission in an Al–NM matrix as a function of volume fraction. Close-packed spheres are used, and the volume fraction is adjusted by changing the inter-particle spacing. For the various *δ* considered, there is a twofold range in particle velocity between the interaction limits. There is weak similarity in the results for various *δ*, which are bound by limiting cases of for VN shock interaction and for CJ shocked flow conditions. The maximum *α* generally occurs at *ϕ*_{s0}=0.2, above which the transmission decreases linearly for increasing *ϕ*_{s0}. The regime of primary interest for dense granular flow in condensed explosives is from 0.2<*ϕ*_{s0}<0.6.

The velocity transmission factor results were fit to a function of volume fraction and detonation reaction-zone length: *α*=*α*(*ϕ*_{s0},*δ*). The volume fraction effect was represented by a second-order polynomial function of *ϕ*_{s0}. The reaction-zone-length influence is exhibited as a shift in the velocity transmission factor, which was represented using an exponential function of *δ* and a fixed offset. The resulting function contains five fitting constants, as follows:
5.2where *c*_{1}=0.2, *c*_{2}=0.1, *c*_{3}=0.4, *c*_{4}=3.3 and *c*_{5}=0.36. A comparison of the fitting function with the mesoscale results is given in figure 14*b*. The agreement is reasonable for *δ*<1; the largest differences occurred for high volume fractions in the large particle limit.

### (e) Particle heating

Figure 15*a* shows the mass-averaged shock compression temperature for a single particle for various *δ*. Smaller particles achieve higher temperatures, although the maximum temperature during shock interaction is much less than the CJ flow temperature (*T*_{CJ}=3657 *K*). Oscillations in the temperature histories are caused by successive compressions and expansions owing to wave reverberation inside the particle.

Figure 15*b* shows the mesoscale results of the single particle temperature and corresponding transmission factors, *β*, which are monotonic decreasing functions for increasing *δ*=*d*_{p}/*L*_{R}. Both the temperature and the temperature transmission factor decrease from the VN to CJ limiting values mainly over the interval from 0.1≤*δ*≤10. Similar to the velocity transmission factor, the temperature transmission factor for a single particle is fit to the sigmoidal function:
5.3where *δ*_{0}=0.70, *w*=0.45, *β*_{CJ}=0.170 and *β*_{VN}=0.337.

Figure 16 shows the particle temperature and corresponding *β* transmission factors in the Al particle-liquid explosive matrix as a function of volume fraction. The temperature transmission factors are scaled using the VN shocked fluid temperature. Similar to particle acceleration, the particle heating is bound between the small particle and large particle-limiting cases, and the effect of solid volume fraction on *β* is reduced for larger particles (). For a given volume fraction, *β* increases with the reaction-zone length (decreasing *δ*). For volume fractions relevant to dense granular flow in condensed explosives, i.e. 0.2<*ϕ*_{s0}<0.6, the range of *β* is limited to 0.32–0.40 across two orders of magnitude for the ratio of *d*_{p}/*L*_{R}.

## 6. Application to macroscopic models

Having gained an understanding of the particle interactions with the explosive fluid and neighbouring particles during shock and detonation in an idealized dense or packed particle-condensed matter system, let us try to apply the mesoscale results to formulate macroscopic momentum and energy transfer functions. These functions can then be incorporated into interphase exchange source terms that are used during the shock compression time scale. The resulting approach can be applied to macroscopic continuum modelling of practical problems such as detonation in a multi-phase explosive or shock propagation in a dense particle-fluid system.

During the detonation shock interaction time, the momentum transfer rate, i.e. the force, *F*_{p}, and the heat transfer rate, *Q*_{p}, acting on a macroscale solid particle-condensed explosive control volume containing *n*_{p} particles each with mass, *m*_{p}, can be approximated by
6.1and
6.2

Because the mesoscale results showed that the increase in mass-averaged particle velocity was mostly linear within the interaction time scale, a simple acceleration is assumed: *u*_{p}(*t*)=*u*_{p0}+(*u*_{p1}−*u*_{p0})(*t*−*t*_{0})/*τ* for *t*_{0}<*t*<*t*_{0}+*τ*, where *t*_{0} is the shock arrival time at the particle leading edge. Combining the standard definition for drag force on spherical particles,
6.3with (6.1) and assuming *u*_{p0}=0 and *u*_{f}≈*u*_{f1}, an effective drag coefficient is obtained for the shock compression interaction:
6.4

In (6.4), *α*=*α*(*ρ*_{f0}/*ρ*_{s0},*ϕ*_{s0},*δ*) can be obtained from the mesoscale results in §5, keeping in mind the limitations of uniform spherical particles in ordered matrices.

Similarly, the particle heating rate in the mesoscale results can be assumed constant within the interaction time, and the particle temperature can therefore be approximated by a linear function: *T*_{p}(*t*)=*T*_{p0}+(*T*_{p1}−*T*_{p0})(*t*−*t*_{0})/*τ*. Combining the standard definition for convective heat transfer on spherical particles,
6.5with (6.2) and assuming *T*_{f}≈*T*_{f1}, an effective Nusselt number is obtained for the shock compression interaction:
6.6

Note that the physical parameters *C*_{d} (equation (6.4)) and *Nu* (equation (6.6)) feature a time-dependence in order to facilitate their numerical implementation. The effective drag force and heat transfer models are suitable for two-phase macroscopic models of dense particles such as Baer & Nunziato (1986).

## 7. Conclusions

Dimensional analysis showed that the particle acceleration force and heat transfer during detonation shock compression in a dense solid particle-condensed explosive system are a function of the material density ratio of explosive to particle, *ρ*_{f0}/*ρ*_{s0}, the volume fraction, *ϕ*_{s0} and the ratio of particle diameter to detonation reaction-zone length, *δ*=*d*_{p}/*L*_{R}, apart from Reynolds number, *Re*, Prandtl number, *Pr* and Mach number, *M*_{0}. While viscosity and heat conduction are important in later time, they can be neglected when compared with the other parameters during the shock compression time scale. Thus, the acceleration force and heat transfer can be described by an effective drag coefficient, *C*_{d}=*f*(*ρ*_{f0}/*ρ*_{s0},*ϕ*_{s0},*δ*,*M*_{0}), and the heat transfer is represented by an effective Nusselt number, *Nu*=*f*(*ρ*_{f0}/*ρ*_{s0},*ϕ*_{s0},*δ*,*M*_{0}). Mesoscale simulations of spherical aluminium particles immersed in liquid NM were conducted by varying *ϕ*_{s0} and *δ* at a given *M*_{0} and *ρ*_{f0}/*ρ*_{s0}, which were known to be important parameters. The full range of *ϕ*_{s0} and *δ* was studied, where *δ* ranged between the small particle limit, with essentially inert shock interaction, to the large particle limit, with infinitely thin detonation front interaction followed by detonation products expansion flow.

Features of heterogeneous detonation were explored including: detonation instability and velocity deficit; pressure front fluctuations with peaks up to four times the CJ detonation pressure and periods proportional to the particle size; and, transverse waves and hot spot characteristics of locally enhanced pressure and temperature fronts. These physics are consistent with macroscopic phenomena observed in published experiments. Detonation failure was not considered in this study because the mesoscale calculations were infinite diameter (no charge edge effects).

From the mesoscale simulation, a shock compression velocity transmission factor, *α*=*f*(*ρ*_{f0}/*ρ*_{s0},*ϕ*_{s0},*δ*,*M*_{0}), and temperature transmission factor, *β*=*f*(*ρ*_{f0}/*ρ*_{s0},*ϕ*_{s0},*δ*,*M*_{0}) were obtained to summarize the acceleration and heating behaviour within the detonation shock interaction time. The maximum particle acceleration occurred at *ϕ*_{s0}=0.25, whereas the maximum shock compression heating occurred over a wider range of solid volume fraction *ϕ*_{s0}=0.43–0.74. Scaling of the transmitted velocity and temperature using the von Neumann state is convenient because it is easily obtained from analytical shock relationships; however, this scaling showed a strong dependence on the ratio of particle diameter to reaction-zone length.

Overall, velocity and temperature transmission factors obtained using mesoscale calculations can be tabulated and fit to simple functions of density ratio, volume fraction and the ratio of particle diameter to detonation reaction-zone length, for a given Mach number. Fitting in this fashion allows practical models to be developed for engineering calculations. As an example, the results were applied to formulate functions for macroscopic momentum and energy transfer between the two phases during detonation shock compression. These functions can then be used as the interphase exchange source terms applied to macroscopic continuum modelling of practical problems such as detonation of a multi-phase explosive or shock propagation in a dense particle-fluid system.

- Received October 2, 2011.
- Accepted January 13, 2012.

- This journal is © 2012 The Royal Society