Modeling the shock-induced multiple reactions in a random b e d of metallic granules in an energetic material

An investigation of shock–particle interactions in reactive ﬂows is performed using an Eulerian hydrody- namic method with a hybrid particle level-set algorithm to handle the material interface dynamics. The analysis is focused on the meso- to macro-scale numerical modeling of a granular metalized explosive containing randomly distributed metal particles intended to enhance its blast effect. The reactive ﬂow model is used for the cyclotrimethylene-trinitramine (RDX) component, while thermally induced deﬂa- gration kinetics describes the aerobic reaction of the metal particles. The complex interfacial algorithm, which uses aligned level sets to track deforming surface between multi materials and to generate the random shape of granule elements, is described for aluminized and copperized RDX. Then, the shock-induced collapse of metal particles embedded in the condensed phase domain of a high explosive is simulated. Both aluminized and copperized RDX are shown to detonate with a shock wave followed by the burning of the metal particles. The energy release and the afterburning behavior behind the detonat- ing shock wave successfully identiﬁed the precursor that gave rise to the development of deﬂagration of the metal particles. © 2019 The Combustion by Elsevier Inc. reserved.


Introduction
When a shock wave collides with a particle, complex flow structures are generated due to the distortion of the incident pressure wave and the shape deformation of the particles; the diffraction of the rarefaction waves develops in various forms due to the interactions between the shock wave and the downstream particles. The presence of the particles acts as an obstacle, creating distortions in the wave front and causing the overlapping of various types of reflected tensile waves from behind the particles. An additional key feature of this process is that metal particles which are combustible can burn and spherically expand into atmosphere, which is a complex phenomenon not easily understood due to the interactions between a large number of metal particles and the strong shock waves generated from an energetic material [1][2][3] .
Metal particle additives in an energetic material enhance the multiple reaction functionality due to the afterburning characteristic of the particles. Such secondary reactions following the primary detonation of an explosive allow for a longer duration of overpressure, which is an intended blast enhancement effect. To understand the extended burning at high pressure condition of * Corresponding author.
We characterize a multi-functional high explosive that is comprised of 50% RDX (C 3 H 6 N 6 O 6 ) and 35% additive metal powders of aluminum (Al) or copper (Cu) with a 15% HTPB (hydroxylterminated polybutadiene) binder. The overall reaction of RDX is C 3 H 6 N 6 O 6 → 3N 2 + 3H 2 O + 3CO. Because under-oxidized explosives produce free carbon, which can form black smoke, the presence of black smoke is a crude indication of severe underoxidation. Some of the products themselves are fuels, specifically free carbon, C, and carbon monoxide, CO. After the burning or detonation reaction is complete, these products may be free to expand into the air. As this occurs, these products may mix with the oxygen in the air, burst into flame, and burn to CO 2 when the proper mixture with the air is reached. If aluminum particles are involved in this reaction as an additive, then the oxidizing competition with carbon atoms will become more intense. Aluminum is reactive and will react spontaneously with water and/or air to form aluminum oxide. Therefore, the overall reaction of aluminized RDX, in which aluminum is added, is as follows. The second reaction ( Eq. (1.1) ) of the oxidation of aluminum is called the secondary fireball and/or afterburning. Such fireballs can also be fueled by other burnable materials, such as copper, silicon, boron, zirconium, and binders that are mixed with the explosive. Referring to the general formula for a CHNO explosive, C x H y N w O z , we see that, for all the carbon to be burned to CO 2 , we need twice the number of oxygen atoms as we have carbon atoms. Most CHNO-based explosives have a negative oxygen balance (OB%), so they always have a fuel-rich reaction and the remaining unburned fuel is able to burn again when reacting with incoming atmospheric air. Therefore, a metal-fuel-rich energetic material reacts in a combined format of detonation followed by deflagration, occurring in a time-delayed sequence. We assume that aluminum combustion is a stoichiometric process. The afterburning of aluminum consists of two aerobic ( Eqs. (1.1) and ( 1.2 )) reactions and two anaerobic ( Eqs. (1.3) and ( 1.4 )) reactions. However, H 2 O and CO 2 can be considered as final products without additional reactions with Al.
In [7,8] , the cluster particles have been assumed to be a continuum and have been treated at the macro-scale level. However, this approach does not provide the detailed behavior of the individual and collective particles. In addition, because this assumption is only applicable to special situations where the particle size is large ( ∼mm) and the population is high (over 50% wt.), it is not suitable for recent heterogeneous explosives in which micro-to nanometersized particles are primarily used.
Recently, Mehta et al. [2] compared and analyzed the interactions between a single particle of a cylindrical (or spherical) rigid body and a shock wave in a 2D geometry. The flow separation of the cylinder occurred later than that of the sphere, and the superposition of a sound wave was observed downstream. Ling et al. [9] considered the deformation of a particle impacted by a shock wave. As the impact pressure was applied, the rounded particle was gradually deformed into a flat shape and vortex shedding occurred. Boiko et al. [10] experimentally visualized the interaction between a shock wave and a cluster of particles using a shock tube and a high-speed photo camera. It was observed that the coarse particles were dispersed in the atmosphere with faster particle velocities then finer ones. Numerical studies of clouds of particles were performed in Ref. [5] . A number (2 or 11) of aluminum particles were mixed in a composite explosive to explore the complex shock interaction processes for multiple particles. However, the initial arrangement of the aluminum was artificially uniformly distributed, and the size and shape of the aluminum particles were also fixed. Therefore, it is likely that the size and position of the particles are far from the actual geometry of composite powders, which are always randomly distributed.
Modern experimental techniques still lack the resolution necessary to capture these phenomena in extremely precise conditions on a length scale of several micrometers and a time scale of a few microseconds. This leads to a motivation for conducting a series of hydrodynamic simulations for analyzing the interactions between metal particles and RDX in a composite mixture.
Metal is an elastoplastic substance that can deform and flow if thermally activated to burn. The complex process of the shockinduced detonation of a high explosive also requires precisely tuned ignition and growth reaction kinetics to accurately reproduce the detonation process. It is imperative to precisely capture the interface between the metal and the explosive, which is in principle the most difficult task associated with describing the physical response associated with shocking a metalized energetic material, as shown in Fig. 1 . Figure 1 is a schematic of a granular metalized explosive or aluminized RDX combined with a binder. Such heterogeneity in the energetic composition ensures an enhanced blast performance with a longer burning time at an extended blast strength. For simplicity, we have considered only two types of metal granules, i.e., aluminum and copper, together with RDX as a blast enhanced explosive.
The hydrodynamic simulations were performed via two-way coupling of the fluid-structure interaction between the condensed phase flow and the deformation of solid particles at the microscale level. The study aims to accurately simulate the detonation of RDX followed by the later burning of the embedded metal granules. The computational work takes into consideration the randomness of placement, distribution, and shape of the particles along with the chemical reaction and deformation due to strong shock waves.

Governing equations
The rapid and violent reaction from a detonation differs from the classical combustion process in that all the energy transfer is governed by the strong compression waves with a limited contribution from such processes as heat diffusion, typically associated with the slow burning process. The leading part of a detonation front is a strong shock wave propagating into the fresh mixture. This compression wave heats up the material as it triggers the chemical reaction, and a balance is attained such that the reaction effectively supports the shock propagation. We assume that molecular diffusion, thermal conduction, and viscous effects are insignificant since the detonation energy is converted rapidly in the shock to detonation transition process. In practice, the time scales of the chemical reactions are very small compared to the time scale of the fluid dynamics.
The compressible Euler equations in a two-dimensional coordinate system reflect the conservations of mass, momentum, and energy as shown in Eqs. (2) , ( 3 ), and ( 4 ), respectively.
Here, Eq. (3) expresses the compressible equations for an energetic material that undergoes a gas phase transition during the chemical reaction and Eq. (4) expresses the governing equations for the deformation and chemical reaction of the metal particles.
Here r and z are the cylindrical coordinate variables, ρ is the density, and u r and u z are the velocity components in the radial and axial directions, respectively. E = e + (u 2 r + u 2 z ) / 2 is the total energy per unit mass, e is the specific internal energy, and p is the hydrostatic pressure. The temperature is derived from the relationship, e = c v T , where c v is the specific heat capacity at a constant volume. In general, it is assumed that c v and c p are nearly equal for solids, such that γ ≡ c p / c v ∼ 1, but for gaseous detonation products they can be very different. Figure 2 shows c v and c p of RDX according to temperature which shows a similar trend while the ratio ( γ = c p / c v ) follows a decreasing pattern with the increasing temperature. The average value of the ratio is about 1.3, and this value was used in the calculation. The reaction rate, ˙ , is described by the empirical ignition and growth relation obtained from shock to detonation transition data for high explosives, while the Arrhenius law is adapted to calculate the thermally induced reaction of the metal particles [11] . A sharp material interface is guaranteed via the use of a hybrid particle level-set method [12] . Then, the resulting system of hyperbolic equations is solved using third-order Runge-Kutta (RK) and fifth-order essentially non-oscillatory (ENO) methods [13] for the temporal and spatial discretizations, respectively. The code uses a stable high-order explicit RK time integrator with its known stability property. To handle stiffness associated with the limited kinetic schemes used in this study using an explicit integrator, the minimum time step for convective expansion, shock advection, and chemical reaction is chosen efficiently throughout the whole calculation.
Here, the non-spherical stresses in the unreacted high explosive are relatively small in comparison to the dominant hydrostatic pressure of the product gas and therefore are conventionally ignored. The metal particles, however, must be closely monitored for changes associated with deformation, and therefore the Cauchy stress tensor is comprised of the deviatoric and hydrostatic (pressure) stresses as follows: where I 1 and J 1 are the first scalar invariants of the Cauchy stress tensor and the deviator stress, respectively. The deviatoric stress tensor, S ij , and the hydrostatic pressure, p , are taken to be positive in compression. The rate of the deviatoric stress change follows a where each operator is defined as where D ij and D i j are the strain-rate tensor and the deviatoric strain-rate tensor, respectively. The components of the strain tensor are used to derive the yield stress depending on the shear rate in forms of the Johnson-Cook flow stress model. The Johnson-Cook model was applied to obtain the flow stress or the minimum outer force needed for plastic deformation. This model makes use of the equivalent plastic strain, strain rate, and melting temperature [14] as shown in Eq. (13) .
Here, ɛ is the strain tensor, ɛ e and ɛ p are the elastic strain tensor and the plastic strain tensor, respectively, σ Y is the yield stress, and ˙ ε p and ˙ ε p0 are the effective plastic strain rate and the effective plastic strain rate of the quasi-static state, respectively. The normalized temperature is defined according to the reference room temperature ( T 0 ) and the reference melt temperature ( T m ). For conditions where ( T − T 0 ) < 0, we assume that m = 1. A strength model that accounts for the effects of strain hardening, strain-rate hardening, and thermal softening was adopted to describe the dynamic response of the solids. In addition, the strength model constants, i.e., A 0 , B 0 , C 0 , m , and n , are used for the aluminum and copper. As the strain rate approaches zero, the natural log approaches negative infinity and, therefore, the Johnson-Cook model sets C 0 to zero if the strain rate reaches a certain minimum value, usually 1 s − 1 . ˙ ε p0 is commonly set to unity. The material properties and the Johnson-Cook stress model constants are summarized in Table 1 [15] .

Equations of state
The correct and suitable constitutive relationship that relates the pressure as a function of the other thermodynamic properties is used to address the mathematical closure of the governing laws of conservation. The equations of state (EOSs) for different elements involved in the problem are identified and combined to reflect unburned reaction states as opposed to the reacted hot product states of the high explosive or reactive metals. For metals, the Mie-Grüneisen EOS [16] is adopted where the corresponding pressure is related to the internal energy, such that where p 0 and e 0 are the pressure and internal energy of a reference state, respectively. The shock relations and Hugoniot equations for the conservation laws are as follows.
u shock = c 0 + s u particle (17) ρ 0 u shock = ρ( u shock − u particle ) (18) p = ρ 0 ( c 0 u particle + su 2 particle ) = ρ 0 u particle u shock (19) Here, the shock velocity is u shock and the material particle velocity is u particle . c 0 and s are the bulk sound speed and the linear Hugoniot slope coefficient, respectively. The shock velocity and the particle velocity follow a linear relationship, and ρ is considered to be a constant. (20) Here, p H is the pressure on the Hugoniot and e H is the internal energy per unit mass on the Hugoniot, such that As a result, the formally defined Mie-Grüneisen EOS becomes Since the sound speed is defined as and setting ρ = ρ 0 0 , the sound speed of the Mie-Grüneisen EOS is given as follows.
To describe the unreacted state of the high explosive, the Jones-Wilkins-Lee (JWL) form of Eq. [17] is used: Here, A, B, C, R 1 , and R 2 are the material-dependent JWL parameters with ω being the Grüneisen coefficient of the explosive and e 0 = ρ 0 C v T .
The reacted product state of either the metals or the high explosive is also given by the JWL product form, such that The JWL EOSs are empirical, and the parameters are obtained from a fitting of the cylinder expansion test results, as well as the use of a thermo-chemical equilibrium code such as CHEETAH [18] . In the present study, the optimal parameterizations of the JWL EOSs were incorporated using multiple CHEETAH runs to satisfy the empirical fitting constraint.
The subsequent sound speeds for both the unreacted and reacted EOSs are given as follows.
The unreacted and reacted EOSs were combined into the single expression shown in Eq. (30) by means of the product mass fraction ( λ) and reactant depletion (1 − λ).
The combined sound speed is then calculated using Eq. (31) .
These partial equations calculate the pressures and sound speed of the reactant and product according to the reaction progress variable, λ. If λ = 0, the pressure EOS for unreacted is used, and λ = 1 corresponds to a completed reaction. The equations are useful for calculating the pressure and subsequent sound velocity for the shock to detonation transition problems. Table 2 summarizes all of the EOSs used in the present simulation [19] .

Pressure-induced "fast" chemical reaction: detonation of high explosives
The rate of production of the burned mass is governed by the chemical species equation:  Table 3 The C-J conditions and detonation model parameters for RDX.
where w i is the reaction rate and λ i is the reaction progress variable or product mass fraction. The reactive flow model consists of ignition and growth steps [7] , as shown in Eq. (33) .
Here, the constants I, a, G, and b are the unknown parameters while λ = 0 and λ = 1 specify the unreacted and reacted states, respectively. The degree of compression due to a shock is defined as μ = ρ/ ρ 0 − 1. The procedure for defining these four unknowns is discussed in Ref. [20] , where a series of standard rate stick tests is used. For RDX, the constants on the ignition I and growth G were set to 5.8 × 10 7 s -1 and 2.4 × 10 6 s -1 GPa -b , respectively. The pressure sensitivity b was set to 1.1, and the compression sensitivity a was set to 4.0. The C-J conditions and detonation parameters for RDX are shown in Table 3 . Figure 3 compares the size effect curves of the aluminized RDX with the experimental and numerical values from Kim et al. [7] . The hydrodynamic simulation is shown to reproduce the detonation velocities of unconfined rate sticks with five different radii (i.e., 0.025 mm −1 , 0.050 mm −1 , 0.075 mm −1 , 0.100 mm −1 , and 0.125 mm −1 ). The speed of the detonation wave was approximately 650 0-750 0 m s −1 , and the error bars of the present calculation are calculated from five trials for each radius. The size effect data using the parameterization and the present method of coupling RDX and aluminum into their respective components are in good agreement with both the referenced experimental data and the numerical data.

Temperature-induced "slow" chemical reaction: deflagration of metal particles
The ignition of particles exposed to the high-temperature environment resulting from the hot product gasses of a detonated explosive depends on the amount of heat added to the particles and the activation energy threshold. Inside the hot gas environment  where the RDX detonation has already passed, the convective heating of the metal particles rapidly increases the enthalpy of the drifting particles toward the activation of metal ignition. However, if the particles expand too rapidly, the energy will be dissipated, and heat loss will occur before reaching the critical temperature for the onset of metal particle deflagration. If only a small portion of the cloud particles is successfully ignited during an unsteady expansion, the combustion wave will propagate and nearby particles will also become reactive. Therefore, a local heterogeneous surface reaction can trigger a sympathetic reaction among the particles and develop into metal flames in a sequence of multiple explosions. This is why it is necessary to understand the energy transfer between the particles and ambience via the use of well-described detonation (for RDX) and deflagration (for Al/Cu) kinetics. The chemical reaction of the metals is governed by a temperature-based Arrhenius law to describe the deflagration, such that The parameters of the reactions, i.e., the activation energy and the pre-exponential factor, are determined by applying differential scanning calorimetry (DSC). The samples are heated at varying rates, and the peak reaction temperatures are recorded for each rate.
Rearranging Eq. (35) and taking the logarithm yields ln The two basic parameters ( d λ/dt and λ) are determined from the DSC exotherm, and Eq. (36) can be solved using multiple linear regression. The activation energy ( E a ) and the pre-exponential factor ( Z ) are obtained from the slope and intercept of the plot, respectively. These parameters are summarized in Table 4 , and the Arrhenius rates are plotted as a function of the temperature in Fig. 4 [21,22] .

Hybrid particle level-set algorithm with alignment fix
To obtain a sharp interface between two different materials, a hybrid particle level-set method [12] was developed. The motion of a level set follows an equation that describes the time evolution of the material interface, ∂φ ∂t Here, the interface of each material is a zero level set, φ = 0 . φ < 0 indicates the inside and φ > 0 indicates the outside of a material. This equation is integrated using a fifth-order scheme in space and a third-order Runge-Kutta method in time [23] . While calculating the interface level-set function, a drastic change in the material properties may give rise to an undesired distortion of the interface. To remedy this well-known weakness of any Eulerian fixed-mesh method, a periodic re-initialization is adapted by solving the following equation until steady state is reached: where d is the grid size.
Two nonlinear characteristics intersecting at the interface for system are given as  Then u I and p I can be calculated directly as p I = ρ l C l p r + ρ r C r p l + ρ l C l ρ r C r ( u l − u r ) ρ l C l + ρ r C r = w l p r + w r p l + w l w r ( u l − u r ) The relevant cases, namely gas-solid conditions are as follows.
At the interface, if ρ r = ρ l , u r = −u l , p r = p l , then u I = 0 , p I = 2 p l . It is necessary to note the dissipation characteristics of any ENO scheme; the repeated re-initialization of distance function level sets leads to round off errors in the actual interface and often violates the required mass conservation. Accordingly, the present in-terfacial algorithm is also subjected to meeting and sufficiently addressing these concerns.
The handling of drastic interactions between two distinct phases or materials can result in level-set warping, as illustrated in Fig. 5 . In particular, the definition of the interface normal vector suffers difficulties under such conditions, causing nearby local variables to be incorrectly defined or to converge to completely nonphysical values. The proposed strategy is to foresee such warping incidents associated with harsh shock conditions or high strainrate deformations. In particular, the following case studies are summarized to further describe the strategy adopted herein. A useful alignment of a level set, where there is the kink or warping occurs, is considered. In principle, the correction procedure is intended to avoid such situations as illustrated in Fig. 6 (a)-(c) and to restart the calculation using the aligned level sets. In these situations, the velocity normal is not fully and correctly defined; therefore, the ghost nodes cannot be defined. A Laplacian average process is used to smooth the kink, and the smoothing operation is described per-vertex as where the odd number (2 n + 1) is the filter width. The level-set values are artificially designated so that the configuration of the level set has the form shown in Fig. 7 . This type of smoothing does not affect the overall accuracy if the grid size becomes sufficiently small. The conservative variables at the center grid point are re-calculated via a distance-based interpolation using the same material grid points around the center grid point.

Initialization of embedded granules in two-phase domain
Circular particles can be generated by defining their center point and a radius; conversely, for a polygon with a random shape, the number of vertices ( N v ), the distance from the center to each vertex, and the angle between the center and the vertices must be randomly determined. The coordinates of the center point ( x, y ) of each random granule are set at random, and the numbering of each vertex is sorted counterclockwise while the distance r i from the center and the angle θ i are determined.
At this time, the angle between any two vertices must not exceed 180 °to ensure an acute angle. The position coordinate of each vertex is given by the following formulas: And the distance r i ' from the center of gravity and the angle θ i ' are calculated by Figure 8 shows the random generation of a sorted pentagonal granule and re-assigning to the center of gravity. With the known coordinates of each vertex of a polygon, the first-order linear equation for each side of the polygon can be determined. Then, the subsequent signed shortest distance from each node in the Eulerian grid to the polygon is constructed, which is the desired level set.
To retain the intended accuracy associated with extrapolating the material properties across the interface, the minimum number of grids inside the granule must be met. A ghost band of grids that is centered at the zero level set and populated in opposite directions from an interface is required. The thickness of the band is determined according to the accuracy of the spatial discretization. Here, the distance from the center of mass to any side is set to 6 x . Therefore, the ghost-band thickness is set to be greater than twice this requirement, i.e., 12 x .
In addition, when generating the initial random granules, overlapping granules are removed so as to maintain the minimum band width during the simulation. Therefore, initialization with a randomly structured level-set domain strictly adheres to tracking the interface between a cloud of particles with a high-order level-set tracking technique. Figure 9 shows the initial particle distribution implemented by random allocation.
The percent weight of each metal granule within a composite sample (or a given computational domain) is estimated as follows. The surface area of the particles is obtained using partial area segmentation and the discretized curve from a computer-aided image.  For aluminum, which is spherically shaped, the following equations are used. Here, r i is the radius of each particle and N is the number of particles. R and h are the radius and height of the explosive, respectively. For copper particles, which are random polygons, the area is obtained using the formula for a polygon: Here, n i indicates the number of each i th particle ( n i angular). The initial geometries satisfying the composite fractions of aluminized and copperized RDX, shown in Table 5 , are used in all calculations in this study.

Strong shock collapse of a single particle
Before considering a randomly distributed cloud of metal particles, a single granule subject to a strong shock wave is considered. Such a simple consideration has previously been reported by Ripley et al. [24] and Lieberthal et al. [25] , both of which are useful for comparison with the present calculation. Figure 10 shows the initial geometry for the computation of the interaction between a shock wave and a single particle. Figure 11 shows timed images of the shock and single aluminum particle interaction via contours of the shadowgraph in the 1st row, pressure in the 2nd row, Al species in the 3rd row, and Fig. 11. Timed images of the shock and single aluminum particle interaction according to the shadowgraph (1st row), pressure (2nd row), Al species (3rd row), and longitudinal velocity (4th row) contours. longitudinal velocity in the 4th row. Here, shadowgraph field was calculated by the Laplacian of density, ∇ 2 ρ defined as the divergence ( ∇• ) of the gradient ( ∇ρ).
The initial shock wave propagates in the upward direction from the bottom of the domain. It starts to collide with the bottom of the aluminum particle at approximately 0.25 μs. The shadowgraph and pressure contours from 0.5 μs to 0.75 μs show that the transmitted wave passing through the aluminum is faster than the shock wave outside the aluminum. The density difference between RDX and aluminum is approximately 1.68 times. The sound velocity of aluminum is 5500 m s −1 , while that of the RDX is much lower. The propagating speed of the main shock is approximately 7800 m s −1 . Therefore, as it progresses inside the aluminum, it is already propagating faster than the speed of sound. Consequently, the shock does not get faster in the aluminum. However, because the expansion wave is reflected according to the sound velocities of RDX and aluminum, it can be seen that the wave propagates at different speeds in each medium. Therefore, one can see that the expansion wave in aluminum is faster in this case. Another interesting feature is the shape change of the aluminum particle. The shock impact causes the particle to flatten and the ends of its sides to protrude slightly, resulting in a high value (red) in the longitudinal velocity contour. The evolution of the shape change is consistent with those reported in Refs. [24,25] .
The evolution of the level set is only concerned with the fluid velocity and has nothing to do with the conservation of mass. Therefore, it is imperative to check whether the mass of the metal particles remains constant from the beginning to the end of the computation. If ideal, the total numbers of particles prescribed by the zero level sets would maintain the total mass throughout the simulation. However, some loss of mass is expected because the time-evolved level-set boundary is continuously reconstructed in the discrete domain. Figure 12 shows the time histories of the mass changes of the aluminum and copper particles during the simulation. The shape of the particle defined by the level set is traced at every time step in the calculation, and the mass is calculated using the following equation.
Here, H is (52) Fig. 13. Computational schematic of a rate stick of metalized RDX and probe locations.
After using the level-set alignment, the simulation is shown to retain approximately 95% of its original mass. Mass conservation for the suggested aligned level-set method is confirmed as illustrated.

The strong shock ignition of a metalized explosive stick (Al-RDX, Cu-RDX) of infinite diameter
The computational domain of a finite rate stick of metalized RDX is shown in Fig. 13 . Three distinct materials, atmospheric air, RDX, and metal granules (Al and Cu), were initially brought into contact via a zero level set.
Granular particles must undergo elastic-plastic deformation due to spherical shock compression. Once stresses exceeding the yield strength of the aluminum particles are applied, the particles experience bending and deformation. A detonation wave with a peak pressure of up to ∼35 GPa is generated by the reaction of RDX.
Numerical tracking of the transient interfacial interactions between granular metals and a high explosive is very challenging. Figure 14 shows an interpretation of the interactions between RDX and randomly distributed aluminum particles. Each particle acts as an obstacle to a propagating detonation wave of RDX, and as such reflected waves are generated in the reverse direction. Reflections spread out along a circle and overlap each other in the wake. The metal granules do not collapse immediately because of their strength and stiffness. Instead, they undergo severe deformation, generating a form of tensile wave. This wave is formed by the reflection of the compressive stress while the detonation wave of RDX impacts the backside interface of the highly dense aluminum particles. These waves collide with neighboring aluminum particles and form reflection waves. With time, the tensile waves develop into a very diverse and complex pattern overlapping and coalescing into the resulting flow field.
Each aluminum particle advances in the forward direction with the particle velocity as its shape is deformed by the pressure of the shock wave. The particle velocity is approximately 30 0-70 0 m s − 1 . As a result, there are slight differences in the velocities and shapes of each particle. Because the pressure and velocity of the shock wave acting on each particle may differ slightly, the geometrical arrangement, shape deformation, and particle velocity may be locally different. Considering such a realistic stochastic distribution, it is possible to directly analyze the behaviors of all the particles and their interactions with the flow field, which develops in a very complicated form, therefore providing a meaningful result.
In the case of copperized RDX, as shown in Fig. 15 , unlike the aluminum particles, the copper particles have more complicated irregular shapes. This geometric factor causes irregularities of refraction and diffraction of detonation waves propagating in the longitudinal direction. A shock wave that passes through a circular particle surrounds the particle and is superimposed on the back surface. Conversely, the polygonal copper particles have sharp corners that cause rapid refraction in the direction of the shock. This leads to a more complicated reactive flow field than in the case of the reacting aluminized RDX. In addition, in the interaction between the shock and the particle, the reflected waves formed by the tensile wave from the initial collision and its repetitive collisions between the neighboring particles appear downstream in a very irregular pattern without an isodirectional tendency.
These results are attributed to only the geometric factors, and it can be seen that the polygonal obstacles make the progression of the shock waves more erratic than do the spherical particles. In Fig. 16 , it is interesting to compare the reactive flow fields of aluminized RDX and copperized RDX behind the detonation wave front. The propagating speed of the detonation wave, as visualized through the profiles of the product mass fraction, is the same in both cases, which means that the detonation speed is ultimately determined by the RDX component. In other words, metal particles that have a slow burning rate relative to RDX do not substantially affect the propagation of the initial detonation wave. Note, however, that the magnitude of the pressure perturbation is different in the two cases. The pressure fluctuation is greater in the copperized RDX than in the aluminized RDX. Therefore, the deviation in the pressure perturbation is calculated to be larger because the reactive flow of the copperized RDX is reflected along the shape of the embedded particle to form various fluctuations. Because the    values of the initial yield stress and the plastic hardening modulus of aluminum are higher than those of copper, their mechanical responses to the shock impact of the detonation wave appear to be different. Aluminum does not show a noticeable change in shape immediately after colliding with a shock pressure of several tens of GPa; however, copper gradually changes in shape over time.
The chemical reaction of the metal particles does not reach the critical energy point at which deflagration starts; therefore, the characteristics of the metal particles do not become prominent at short time scales. However, because the pre-exponential factor of aluminum is much higher than that of copper, the chemical reaction of aluminum proceeds more rapidly to thermal runaway.

The strong shock ignition of a metalized explosive stick (Al-RDX, Cu-RDX) of finite diameter intended for multiple reactions
To understand the afterburning process of sufficiently metalized explosives, the detonation and evolution of the post-detonation flow need to be numerically simulated. We simulated the explosion of oxygen-deficient components containing spherical aluminum and polygonal copper metal particles. The afterburning effect was obtained due to the prolonged overpressure and heating of the explosive products in the air. The afterburning process can be completely irregular or random because it is governed by explosive mixing and the atmospheric explosion. However, the main physics of the multiple reactions of such a blast-enhanced energetic mixture usually follow two stages: the initial reaction (or first detonation of RDX) followed by sporadic afterburning (or metal particle combustion). In the detonation process, which is often oxygen deficient, the RDX transforms into hot gaseous products consisting of carbon dust and carbon monoxide because it has a negative oxygen balance ( −21%). These carbonic gasses combine with oxygen in the atmosphere during subsequent reactions and are converted into carbon dioxide. The heat of combustion of the secondary deflagration of the metal is much higher than the RDX detonation energy. The afterburning occurs over a much larger area compared to the first explosion. Accordingly, it is a very complex hydrodynamic process involving detonation propagation, shock reflection, and particle interactions in a very high-pressure and high-temperature environment. Figure 17 shows the simulation results for the pressure (upper) and total energy (lower) evolution for copperized RDX. A total of 50 polygonal copper particles were randomly considered in the RDX. The calculation results are depicted in half using axisymmetry. The detonation wave propagating in the longitudinal direction of the RDX component is observed up to approximately 40 μs. In this case, as observed in Fig. 15 , cellular-like structures develop in the wake while the detonation wave collides with the randomly distributed particles in the RDX. Once the detonation reaches the outline surface of the explosive charge, a blast wave is formed and propagates into the surrounding air. A flow velocity occurs in the exhaust gas flow so that the particles move in the longitudinal and radial directions. This high-temperature and high-pressure environment is maintained for a certain period of time, and heat energy is supplied to the particles. After 70 μs, ignition starts for some particles. This begins to occur when the applied energy exceeds the activation energy. The particle burning progresses gradually, and the second stage of afterburning develops. The energy contour at the bottom of Fig. 17 is useful to understand the energy transfer process in the reactive exhaust gas flow. During the propagation of the initial detonation in the RDX, an energy gradient appears between the particles (at ∼40 μs). In the transient period, chemical energy is transferred from the exhaust gas to the particles (at ∼70 μs). Eventually, after 90 μs, afterburning occurs from the sufficiently heated particles and these particles and their neighboring particles are burned together.
The ignition process for copper particles is elaborated in Fig. 18 by showing additional later time images. The resulting deformation and interface evolution are shown for the reaction progress variable of copper from 90 μs to 100 μs with 1-μs intervals. Before the particles are ignited, only the mechanical deformation due to the external shock impact is observed. Once the particles are ignited, their size and shape are deformed dramatically in the form of a metal flame. Looking at individual particles, the ignition appears to be evenly distributed across the particles, with the result that the highest part of the species is the surface. In other words, the surfaces of the particles react with the atmospheric oxygen first and, therefore, surface reaction characteristics are observed.
The dispersion and combustion of aluminum particles in the post-detonation flow are discussed in Fig. 19 . This figure shows timed images of the pressure (upper), total energy (middle), and reaction progress variable (lower) contours for aluminized RDX. The shock propagation of the RDX component and the subsequent afterburning of the Al particles are clearly captured. A blast wave followed by the hot detonation product gasses is the primary constituent of the post-detonation flow. The detonation front in the explosive charge can be described as a high-pressure and hightemperature reaction zone separating the unburned aluminum particles and the detonation product gasses. The detonation in the RDX develops a strong shock wave with pressure on the order of 10 9 ∼10 10 Pa. This high pressure provides a trigger for generating the required enthalpy on the aluminum particles to result in afterburning. The variation in the pressure of the exhaust gas and the acceleration of the condensed phase flow generates energy transfer and induces the afterburning of aluminum particles. The combustion of the aluminum particles and therefore the energy release due to afterburning generates the second peak pressure. When ignited, these particles react rapidly and generate high-pressure and high-temperature flows leading to a blast wave. Afterburning occurs in the form of the gaseous thermal expansion due to the reaction progress, and the interface of the flame develops very irregularly due to the dispersal and mixing of two or more different phases with different densities in the atmosphere.
It is vital to understand the detonation wave propagation and the development of the afterburning to characterize the postdetonation flow and blast enhancement. Followed by the analysis of the condensed phase flow ensuing homogeneous explosion, Fig. 20 show the timed histories of pressure evolution for aluminized RDX measured at 6 probes along the centerline and 6 more probes along the void region. In the left figure, the traveling speed of the first peak pressure was about 60 0 0 m/s. This is consistent with the size effect behavior of the aluminized RDX rate stick considered above. The important feature is the appearance of the second peak after the first peak and transient period. The second pressure is observed in the probe P 3 located at the middle of the explosive at first. The location where afterburning begins first is determined by the interaction with the RDX reaction, which is influenced by the initial distribution, size and density of particles. Therefore, given the randomly configured real situation, it will be very difficult to specify the first location of the afterburning. Nevertheless, since the occurrence time of the second peak is about 70 ∼ 100 μs, the approximate transient time can be secured. As shown in the right figure, the pressure of the blast wave was measured to be lower than the pressure in the centerline of the explosive because the interface developed in the multiphase flow with different densities. However, since the afterburning of the particles appears in the form of a metal flame, the pressure drop does not occur rapidly in the downstream, but rather the pressure is continuously increased or maintained to 300 μs. After that, it can be seen that the pressure profiles are maintained for the calculated time ( ∼500 μs).
To understand the triggering mechanism of the multiple reactions in the metalized RDX, the strain evolution of the aluminum is analyzed. The effective plastic strain is a monotonically  increasing scalar value that is calculated incrementally as a function of the plastic component of the rate of deformation tensor. The value of the effective plastic strain is the integral of the stepwise increments of the plastic deformation for a calculated time period. Figure 21 shows timed contours of the effective plastic strain. This figure indicates the level of plastic strain during the detonation phase ( t 1 -t 4 ) and subsequent aluminum afterburning ( t 5 -t 10 ). The strain value starts to increase rapidly after the afterburn spreads as a flame starting at t 4 .
In reality, aluminum particles exist in the state of alumina (Al 2 O 3 ), and the melting point of aluminum or the ignition temperature is 926 K while the melting point of alumina is 2300 K. By considering alumina, the aluminum will not begin to react at 926 K, and the burning will begin after oxide layer has been stripped off at 2300 K or higher. However, this study assumes that pure aluminum particles are used without the coated alumina. Thus, combustion starts with a melting temperature at 926 K. Figure 22 shows timed profiles along the centerline for effective plastic strain at selective times from t 1 to t 10 . The strain does not change much because the aluminum particles maintain their stiffness. However, afterburning takes place after the transition period of t 4 ∼ t 5 ( ∼ 926 K), and relatively large deformation is observed from t 6 . Therefore, the phase change due to the aluminum reaction gives rise to a large deformation in each granule from the randomized bed of aluminum in RDX.

Conclusions
This study considers a full-scale hydrodynamic process that includes a step by step description of how such detonation of high explosive with embedded metal particles of various shapes must be modeled and calculated. The analysis is focused on the meso-to micro-scale simulations of a metalized RDX with randomly populated particles for the intended afterburning effect. We have developed a pseudo randomly crystallizing algorithm and the aligned level-set method for tracking the instantaneously deforming material boundaries and collapse of individual metal particles subject to a detonating shock impact. The hydrodynamic simulations were performed via the two-way coupling of the fluidstructure interaction between the condensed phase flow and the deformation of the solid particles at the microscale level. The study is aimed to accurately simulate the detonation of a basis explosive (RDX) followed by the later burning of embedded metal granules, namely aluminum or copper. An initial detonation shock pressure on the order of 10 9 -10 10 Pa caused the deformation of the shape of the particles. Because the metal particles burned at later times, they effectively gave rise to a prolonged afterburning following the primary detonation. The precise simulation of deforming material interfaces through which the energy transfer as well as the thermo-chemical reaction of the metallic granules occurs has been a challenging task; consequently, no earlier attempts have adequately reported simulating the shock ignition of a metalized explosive. The detailed numerical simulation reveals that the intrinsic mechanism of spontaneous afterburning of metal particles is strongly related to the energy transfer from the detonation of RDX and the ignition sensitivity to the required activation energy. The energy release and the expansion rate behind the detonation wave give rise to the developing metal flame associated with the afterburning of the enhanced blast energetic material. The additional complexities to a three-dimensional domain that includes the granular structures of metals, other oxidizers, explosives, and binders could further advance the current state of the art for the calculation of interactions between a shock and granular metalized energetic materials.