ABSTRACT

The paper aims to theoretically and numerically investigate the confinement effect of inert materials on the detonation of insensitive high explosives. An improved shock polar theory based on the Zeldovich-von Neumann-Döring model of explosive detonation is established and can fully categorize the confinement interactions between insensitive high explosive and inert materials into six types for the inert materials with smaller sonic velocities than the Chapman-Jouguet velocity of explosive detonation. To confirm the theoretical categorization and obtain the flow details, a second-order, cell-centered Lagrangian hydrodynamic method based on the characteristic theory of the two-dimensional first-order hyperbolic partial differential equations with Ignition-Growth chemistry reaction law is proposed and can exactly numerically simulate the confinement interactions. The numerical result confirms the theoretical categorization and can further merge six types of interaction styles into five types for the inert materials with smaller sonic velocity, moreover, the numerical method can give a new type of interaction style existing a precursor wave in the confining inert material with a larger sonic velocity than the Chapman-Jouguet velocity of explosive detonation, in which a shock polar theory is invalid. The numerical method can also give the effect of inert materials on the edge angles of detonation wave front.

Keywords:

Detonation, Confinement Effect, Insensitive High Explosives, Shock Polar Theory, Cell-Centered Lagrangian Method

1. Introduction

Insensitive high explosives (IHEs) are gaining popularity in weapon engineering due to their safety. Because of two main characteristics of IHEs: slower detonation velocity and wider chemical reaction zone, the detonation of IHEs is more easily influenced than the detonation of sensitive explosives by confinement materials. The confinement interaction of different inert materials can change the front shape and propagation speed of detonation shock wave of IHEs, and leads to the IHEs detonation transform to the non-ideal state when its front shape turns into a curve or its propagation speed decreases under the Chapman-Jouguet (CJ) velocity of the explosive detonation.

Earlier studies on the confinement effects focused on the driving capability of IHEs and the accepted work of the driven materials. With the development of numerical simulations and experimental techniques, it is possible to obtain more information on the detonation flowfield to better understand the mechanism of the confinement effect. Aslam and Bdzil [1] [2] [3] divided the confinement effect of inert materials into either a strong or weak confinement based on the shock polar theory, and their two-dimensional Eulerian numerical results verified their division. Hill and Aslam [4] designed a “sandwich” test to measure the detonation front shapes of the PBX-9502 explosive under the confinement of polymethyl methacrylate (PMMA) and SS304 steel. Tarver et al. [5] constructed a three-term ignition-growth chemistry reaction model to examine the confinement effects of bronze and aluminum. Anderson and Jackson [6] found the shock front of IHEs detonation confined by a layered slab became curve in flow measurement. Short and Jackson [7] developed a non-ideal dynamics analysis for high sound-speed metal confiners driven by IHE. Eden and Belcher [8] investigated the effects of brass and beryllium (Be) on an EDC35 explosive in the sandwich test, and observed the presence of a preshocked layer in the undetonated explosive and a precursor wave in Be. Aveille and Carion [9] measured the surface velocity of a copper plate that confined TATB-based explosives and indicated a distinctive difference between the experimental result and the simulation results derived from a first-order Lagrangian method with program-burn detonation model. Balaganskii and Aguereikin [10] experimentally investigated the confinement interaction of porcelain material on TATB-based explosives and also observed the presence of the precursor wave in porcelain. In order to develop much precise theory and numerical method to investigate the confinement effect in IHEs, the author and his colleagues try to propose an improved shock polar theory and a cell-centered Lagrangian hydrodynamic method recently [11] [12] [13] .

The present paper studies the physical mechanism for the confinement effect of different inert materials on PBX-9502, a widely-applied IHE. Part 2 presents the theoretical study, which includes the establishment of an improved shock polar theory and the categorization of the confinement effect for inert materials with a sonic velocity smaller than the CJ velocity of explosive detonation. Part 3 presents the numerical simulation study, which comprises: 1) the full derivation of a second-order two-dimensional cell-centered Lagrangian hydrodynamic method with a phenomenological detonation reaction model; and 2) a detailed analysis of the confinement interaction for typical inert materials with sonic velocities smaller and larger than the CJ velocity of explosive detonation, and a representative comparison and discussion of the detonation edge angles between the improved shock polar theory and the Lagrangian numerical simulation. Part 4 presents the main conclusions.

2. Theoretical Study on the Confinement Effect of Inert Materials

2.1. An Improved Shock Polar Theory

A relatively simple shock polar theory can provide a good leading-order prediction of the confinement effect. The shock polar analysis generally utilizes a reference frame that moves along with the intersection point between the detonation shock wave and the material interface, and considers the matching relation of pressure p behind shock wave to the deflection angle θ of the streamline across the shock wave fronts of the IHEs and the inert materials.

Traditionally, detonation under the shock polar theory is described by the CJ model, which neglects the chemical reaction zone and regards the detonation as a strong discontinuity front with zero thickness [14] . In addition, the confinement interaction is determined by the reflection polar curve of the gas product of explosive detonation and the refraction shock polar curve of the confining inert material. Some evidences [1] [9] indicated that neglecting the chemical reaction zone results in a considerable difference between the confinement effect described by theoretical predictions and the confinement effect observed in experiments.

Because the chemical reaction zone of IHEs is relatively wide, the Zeldovich-von Neumann-Döring (ZND) model is more precise for the analysis of the confinement effect as it considers the chemical reaction structure. Under the ZND detonation model, a leading shock wave and a closely following chemical reaction zone move at CJ velocity, which is the intrinsic velocity of a explosive, such that the leading shock wave and the chemical reaction zone form an indivisible whole structure called the detonation wave, and the detonation wave can self-sustainingly propagate by means of the releasing energy of the chemical reaction. When a detonation wave interacts with a confining material, its leading shock wave first refracts into the confinement material. As a result, it is more reasonable to establish a shock polar at the leading shock wave of the unreacted explosives. Figure 1 presents the typical flow state of detonation in IHEs that are confined by a metal with a smaller sonic velocity than the CJ velocity of explosives detonation. The flowfield will become steady state when the reference frame moves at constant speed along with intersection point “O” between the leading shock front of detonation and the interface of the confinement material. Conversely, the detonation will transform to a non-ideal state due to the lateral expansion, which results in a curved detonation shock front. Under the moving reference frame, the flow before the leading shock wave is supersonic and the flow behind the leading shock wave is subsonic. It is the subsonic velocity that

makes the leading shock wave not generate any reflection wave. Therefore, the confinement interaction is determined by the polar curve of the leading shock wave within the unreacted explosive and the polar curve of the refraction shock wave within the confining material. Evidently, the main difference between the traditional shock polar theory and the improved shock polar theory lies in the difference between neglecting and considering the chemical reaction zone while adopting different detonation models. When the propagation speed of detonation shock and the equations of state (EOS) of the explosives and the confining inert materials are known, the expression of the improved shock polar curve can be analytically obtained.

For the unreacted explosives, a Jones-Wilkins-Lee (JWL)-type of EOS is generally adopted, of which the form of the shock polar curve of the leading shock wave can be expressed as follows:

$\{\begin{array}{l}p=\frac{b\left(1\right)a\left(\overline{v}\right)-b\left(\overline{v}\right)a\left(1\right)}{b\left(1\right)-b\left(1\right)\left(1-\overline{v}\right)b\left(\overline{v}\right)/\left(2{\rho }_0\right)}\hfill \\ \text{tg}\theta =\frac{\sqrt{{\rho }_0{D}_{CJ}^{2}p\left(1-\overline{v}\right)-p^2}}{{\rho }_0{D}_{CJ}^2-p}\hfill \end{array}$ (1)

where $\overline{v}=\frac{v}{v_0}$ ; $a\left(\overline{v}\right)=A\left(1-\frac{\omega }{R_1\overline{v}}\right)\exp \left(-R_1\overline{v}\right)+B\left(1-\frac{\omega }{R_2\overline{v}}\right)\exp \left(-R_2\overline{v}\right)$ ;

$b\left(\overline{v}\right)=\frac{\omega {\rho }_0}{\overline{v}}$ ; ${\rho }_0$ is the standard density and v is the specific volume, and $D_{CJ}$

is CJ velocity of the explosives. In addition, $A,B,R_1,R_2$ and ω are the constants related to the unreacted explosives [5] .

When the detonation product expands, under the JWL-type of EOS, the form of the polar curve of the rarefaction wave can be expressed as follows:

$\{\begin{array}{l}\frac{\text{d}e}{\text{d}\overline{v}}=-v_{0}a\left(\overline{v}\right)-v_{0}b\left(\overline{v}\right)e\hfill \\ \frac{\text{d}{q}^2}{\text{d}\overline{v}}=\pm \frac{2c^2}{\overline{v}}\hfill \\ \frac{\text{d}\theta }{\text{d}\overline{v}}=\mp \frac{\sqrt{c^2\left(q^2-c^2\right)}}{\overline{v}{q}^2}\hfill \\ p=A\exp \left(-R_1\overline{v}\right)+B\exp \left(-R_2\overline{v}\right)+C{\overline{v}}^{-\left(1+\omega \right)}\hfill \end{array}$ (2)

where e is the specific internal energy, q is the local velocity of detonation, and c is the sonic velocity of detonation, $c^2=v_0{\overline{v}}^{2}A{R}_1\exp \left(-R_1\overline{v}\right)+v_0{\overline{v}}^{2}B{R}_2\exp \left(-R_2\overline{v}\right)+v_{0}C\left(1+\omega \right){\overline{v}}^{-\omega }$ . The initial conditions of the above ordinary differential equations are the values of the explosives at sonic state.

For the inert materials, a $p\left(\rho ,T\right)-e\left(\rho ,T\right)$ -type of EOS is generally adopted, wherein the form of the shock polar curve of the refraction shock wave can be expressed as follows:

$\{\begin{array}{l}e-e_0=\frac{1}{2}\left(p+p_0\right)\left(v_0-v\right)\hfill \\ \text{tg}\theta =\pm \frac{\sqrt{{\left({\rho }_{0}D\right)}^2\left(p-p_0\right)\left(v_0-v\right)-{\left(p-p_0\right)}^2}}{{\rho }_0{D}^2-\left(p-p_0\right)}\hfill \\ p=p\left(\rho ,T\right)=p_c\left(\rho \right)+p_n\left(\rho ,T\right)+p_e\left(\rho ,T\right)\hfill \\ e=e\left(\rho ,T\right)=e_c\left(\rho \right)+e_n\left(\rho ,T\right)+e_e\left(\rho ,T\right)\hfill \end{array}$ (3)

where $v_0,{\rho }_0,p_0$ and D are the initial specific volume, initial density, initial pressure, and shock wave speed, respectively, and the three terms on the right hand side of the equality sign for $p\left(\rho ,T\right)$ and $e\left(\rho ,T\right)$ denote the cold pressure and energy, the thermal pressure and energy of ions, and the thermal pressure and energy of electrons, respectively [15] .

Under the shock polar theory, the edge angle of detonation (the angle between the leading shock wave and the interface of the undisturbed confinement

material at their intersection point) is $\alpha =\arcsin \left(\frac{1}{D}\sqrt{\frac{p-p_0}{{\rho }_0\left(1-{\rho }_0/\rho \right)}}\right)$ , which is

an important variable to characterize non-ideal state. Table 1 presents the edge angle of PBX-9502 detonation under copper confinement from the improved shock polar theory and the traditional shock polar theory, and the improved shock polar theory has much nearer result to the experiment, which validates the appropriateness of the improved shock polar theory.

2.2. Categorization of the Confinement Effect Based on the Improved Shock Polar Theory

If the sonic velocity of a confinement material is smaller than the CJ velocity of the explosive detonation, the intersection point “O” in Figure 1 will exist, thereby allowing the adoption of the shock polar theory. Figure 2 presents six types of matching ways for the pressure behind shock and deflection angle of the

streamline, wherein the solid line represents the polar curve of the leading shock wave of detonation and the dashed line represents the polar curve of the refraction wave for the confinement material, and the points “Ce” and “Cm” represent the sonic state of polar curve. The six types of inert materials are steel (Fe), phenolic resin, PMMA, silicone rubber foam, lithium hydride (LiH), and lithium deuteride (LiD). The six presented types of matches correspond to the six types of confinement effects.

Figure 2(a) presents an intersection point at the strong branch of the polar curve of the unreacted explosive, which defines this type of confinement as a “strong confinement”. In this case, the flow behind the leading shock wave is subsonic, and no reflected wave travels back into the IHEs. At the same time, the intersection point is located at the weak branch of the polar curve of the confinement material which indicates that the flow behind the refraction shock wave is supersonic. Figure 2(b) presents an intersection point located at the rarefaction branch originating from the sonic point of the polar curve of the unreacted explosive, which defines this type of confinement as a “weak confinement”. In this case, the flow behind the leading shock wave expands to the sonic state through a chemical reaction, from which the leading shock wave is then transmitted into the confinement material by means of a Prandtl-Meyer rarefaction fan. At the same time, the intersection point is located at the weak branch of the polar curve of the confinement material which indicates that the flow behind the refraction shock wave is supersonic. Figure 2(c) presents a polar curve for the confinement material that intersects with both the strong and rarefaction branches of the polar curve of the unreacted explosive, which defines this type of confinement as a “one strong and one weak solution confinement”. Additionally, two intersection points were observed at the strong and weak branches of the polar curve of the confinement material. Figure 2(d) presents three intersection points between the polar curves, for explosive one interaction point is located at the supersonic weak branch and the remaining two intersection points are located at the rarefaction weak branch originated from the sonic state of the unreacted explosive. Conversely, the confinement material exhibits one intersection point located at the supersonic weak branch and the remaining two intersection points located at the strong branch, which defines this type of confinement as a “three weak solutions confinement”. Figure 2(e) presents an intersection point located at the supersonic weak branch of the polar curve of the unreacted explosive and at the strong branch of the polar curve of the confinement material, thereby defining this type of confinement as a “one supersonic weak solution confinement”. Figure 2(f) indicates the absence of an intersection point between the polar curves of the unreacted explosive and the confinement material, thereby defining this type of confinement as a “zero solution confinement”, however, a Prandtl-Meyer rarefaction fan originating from the sonic point in the confinement material could intersect the supersonic weak branch of the polar curve of the unreacted explosive, which means that a local expansion acceleration could appear in the confinement material.

3. Numerical Study on the Confinement Effect of the Inert Materials

The shock polar theory only presents the flow state near the intersection point between the leading shock wave of detonation and the interface of the confinement material, however, numerical simulations can provide more details on the entirety of the flowfield for both the detonated explosives and the confinement materials. Additionally, the shock polar theory can only analyze the confinement effect of inert materials that have smaller sonic velocities than the CJ velocity of the explosives. For the inert materials having larger sonic velocity than the CJ velocity of the explosives, a precursor wave in the confinement material may appear ahead of the detonation wave, thereby the flow near the interaction point become unsteady and the shock polar theory become invalid. Therefore, numerical simulations are essential. The confinement interaction of inert materials on IHEs can be regarded as a compressible multi-material flow problem and be suitable to adopt a second-order, cell-centered Lagrangian hydrodynamic method. It is well known that the key point of the cell-centered Lagrangian method is its determination of the velocity at the cell vertex, and the presented study obtains the velocity at the cell vertex by means of constructing a two-dimensional Riemann solver at the cell vertex based on the characteristic theory of multidimensional reactive flow equations.

3.1. A Second-Order, Cell-Centered Lagrangian Hydrodynamic Method for Chemistry Reaction Flows

3.1.1. Governing Equations for the Chemistry Reaction Flows

The governing equations of the chemical reaction flows in the Lagrangian formulation are presented as follows:

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega \left(t\right)}\text{d}\Omega }={\displaystyle {\int }_{\partial \Omega \left(t\right)}u\cdot \text{d}l}$ (4.1)

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega \left(t\right)}\rho \text{d}\Omega }=0$ (4.2)

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega \left(t\right)}\rho u\text{d}\Omega }=-{\displaystyle {\int }_{\partial \Omega \left(t\right)}p\text{d}l}$ (4.3)

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega \left(t\right)}\rho E\text{d}\Omega }=-{\displaystyle {\int }_{\partial \Omega \left(t\right)}pu\cdot \text{d}l}$ (4.4)

$\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{\Omega \left(t\right)}\rho \lambda \text{d}\Omega }={\displaystyle {\int }_{\Omega \left(t\right)}\rho r\text{d}\Omega }$ (4.5)

where ρ is the density, u is the velocity, p is the pressure, E is the specific total energy, $E=e+u^2/2$ , e is the specific internal energy, λ is the chemical mass fraction, r is the chemical reaction rate about the three-term ignition-growth reaction law [5] , and $\Omega \left(t\right)$ is a moving control volume with the boundary surface $\partial \Omega \left(t\right)$ and the boundary differential length $\text{d}l$ . For two-dimensional flows, $u=ui+vj$ , where u and v are the component velocities.

Equations (4.1)-(4.5) are the conservation laws for geometry, mass, momentum, total energy, and chemical mass fraction, respectively. In particular, the geometry conservation law is essential for the Lagrangian finite volume method as it ensures the compatibility between the grid motion and the area variation of a control volume. The EOS of the JWL-type for explosives and the $p\left(\rho ,T\right)-e\left(\rho ,T\right)$ -type for inert materials were adopted.

3.1.2. The Finite Volume Scheme

For a discretized control volume ${\Omega }_c$ with a mass $m_c={\displaystyle {\int }_{{\Omega }_c}\rho \text{d}\Omega }$ and an area $A_c={\displaystyle {\int }_{{\Omega }_c}\text{d}\Omega }$ , the average value of any physical variable f can be defined as

${\overline{f}}_c=\frac{1}{m_c}{\displaystyle {\int }_{{\Omega }_c}\rho f\text{d}\Omega }$ . Therefore, Equation (4.2) becomes an algebraic equation

${\overline{\rho }}_c{A}_c=m_c$ , and Equations (4.1), (4.3), (4.4), and (4.5) can be expressed as follows:

$\frac{\text{d}{\overline{U}}_c}{\text{d}t}=-\frac{1}{m_c}{\displaystyle {\int }_{\partial {\Omega }_c}H\cdot n\text{d}l}+{\overline{r}}_c$ (5)

where ${\overline{U}}_c={\left(-{\overline{\tau }}_c,{\overline{u}}_c,{\overline{v}}_c,{\overline{E}}_c,{\overline{\lambda }}_c\right)}^{\text{T}}$ , ${\overline{r}}_c={\left(0,0,0,0,{\overline{r}}_c\right)}^{\text{T}}$ , ${\overline{\tau }}_c=A_c/m_c$ , n is the outward unit vector normal to the boundary surface of the control volume, and $H={\left(u,p,0,pu,0\right)}^{\text{T}}i+{\left(v,0,p,pv,0\right)}^{\text{T}}j$ is the tensor of the flux.

For any non-overlapping structured quadrilateral mesh with sides denoted by $I_k$ (k =1, 2, 3, 4), the semi-discrete finite volume discretization of Equation (5) can be written as

$\frac{\text{d}{\overline{U}}_c}{\text{d}t}=-\frac{1}{m_c}{\displaystyle \sum _{k=1}^4{\displaystyle {\int }_{I_k}H\cdot n\text{d}l}}+{\overline{r}}_c$ (6)

Based on the trapezium rule to approximate the integration ${\displaystyle {\int }_{I_k}H\cdot n\text{d}l}$ about the numerical interface flux and the mathematical meaning of the semi-discrete model to describe the instantaneous behavior of a dynamics system at some initial time, a temporal-spatial second-order and implicit-explicit Runge-Kutta scheme can be used to solve Equation (6) as follows (the stiffness of chemical reaction leads to the implicit scheme to solve the chemical mass fraction equation):

$\begin{array}{l}{\overline{U}}_c^{(\ast )}={\overline{U}}_c^{\left(n\right)}-\frac{\Delta t}{2m_c}{\displaystyle \sum _{r=1}^4\left[H_r\left(E_0{R}_c{\overline{U}}_c^{\left(n\right)}\right)+H_{r+1}\left(E_0{R}_c{\overline{U}}_c^{\left(n\right)}\right)\right]\cdot n_{r,r+1}{L}_{r,r+1}}+\frac{\Delta t}{2}{\overline{r}}_c^{(\ast )}\\ {\overline{U}}_c^{\left(n+1\right)}=\frac{{\overline{U}}_c^{\left(n\right)}+{\overline{U}}_c^{(\ast )}}{2}-\frac{\Delta t}{4m_c}{\displaystyle \sum _{r=1}^4\left[H_r\left(E_0{R}_c{\overline{U}}_c^{(\ast )}\right)+H_{r+1}\left(E_0{R}_c{\overline{U}}_c^{(\ast )}\right)\right]\cdot n_{r,r+1}{L}_{r,r+1}}\\ +\frac{\Delta t}{2}{\overline{r}}_c^{\left(n+1\right)}\end{array}$ (7)

where E₀ is the multidimensional Riemann solver to compute the instantaneous evolution solution of a mesh vertex at time $t_n^{+}=t_n+0$ , namely $\overline{U}\left(t_n^{+}\right)=E_0\overline{U}\left(t_n\right)$ , R_c is a reconstruction operator to transform the average value of a mesh cell to a spatially linear distribution within the mesh cell [19] , r and r + 1 is the counterclockwise numbering of the neighboring vertices of the quadrilateral mesh with periodic way such that the fifth vertex is the same as the first vertex, and $L_{r,r+1}$ is the length of a side with vertices numbered as r and r + 1 for the quadrilateral mesh.

In the present notation, $t_n^{+}$ denotes the infinitely small time interval in terms of the initial time t_n, expressed as

$t_n^{+}=t_n+0=\underset{\tau \to 0}{\lim }\left(t_n+\tau \right)$

where τ is a small time interval. Here, E₀ may be also named as the vertex solver to generate the solution of the mesh vertex, of which the following subsection will present its expression.

3.1.3. Vertex Solver E₀

The interface flux of a mesh cell is only correlated to three physical variables, namely u, v, and p from Equation (5). Therefore, the solutions of the three presented variables, specifically u, v, and p can be obtained for the vertex solver. The central idea of the vertex solver is to compute the theoretical solution along every bi-characteristic direction for a small time interval from the given initial conditions by means of the characteristic theory about the hyperbolic partial differential equations, of which certain approximation operations and limit operations can be derived to obtain the analytical solution at the infinitely small time interval [20] [21] . This can be regarded as a generalization of the original idea of the Godunov scheme for multidimensional hyperbolic conservation laws.

To obtain the theoretical solution about the nonlinear hyperbolic system, a suitable linearized treatment in terms of the primitive variables is necessary to reduce the bi-characteristics to straight lines. For this purpose, it is convenient to start with the following flow equations about the primitive variables:

$\frac{\text{d}q}{\text{d}t}+A_1\left(q\right)\frac{\partial q}{\partial x}+A_2\left(q\right)\frac{\partial q}{\partial y}=r_q$ (8)

where $q=\left[\begin{array}{c}\rho \\ u\\ v\\ p\\ \lambda \end{array}\right]$ , $A_1\left(q\right)=\left[\begin{array}{ccccc}0& \rho & 0& 0& 0\\ 0& 0& 0& \frac{1}{\rho }& 0\\ 0& 0& 0& 0& 0\\ 0& \rho c^2& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$ , $A_2\left(q\right)=\left[\begin{array}{ccccc}0& 0& \rho & 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{\rho }& 0\\ 0& 0& \rho c^2& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$ , and $r_q=\left[\begin{array}{c}0\\ 0\\ 0\\ -r{\frac{\partial p}{\partial \lambda }|}_{\rho ,e}\\ r\end{array}\right]$ , wherein c is the sonic velocity under the chemical reaction.

The flow equations are linearized by freezing the Jacobian matrices about the reference state $\tilde{q}={\left(\rho ,u,v,p,\lambda \right)}_{t=t_n}^{\text{T}}$ at the initial point $\tilde{P}={\left(x,y,t\right)}_{t=t_n}$ . The linearized system with frozen Jacobian matrices can be written as follows:

$\frac{\text{d}q}{\text{d}t}+A_1\left(\tilde{q}\right)\frac{\partial q}{\partial x}+A_2\left(\tilde{q}\right)\frac{\partial q}{\partial y}=r_q$ (9)

The presented special characteristic surface, namely the so-called characteristic cone or Mach cone, generated by all the bi-characteristic lines passing an evolution point $P={\left(x,y,t\right)}_{t=t_n+\tau }$ (τ is the evolution time) is considered. This solution at the evolution point $P={\left(x,y,t\right)}_{t=t_n+\tau }$ for Equation (9) can be given in terms of the conditions at the initial time $t_n$ .

Consideration of any unit vector denoted by $n\left(\theta \right)={\left(\cos \theta ,\sin \theta \right)}^T$ , $\theta \in \left[0,2\text{π}\right]$ , presents a matrix pencil $A\left(\tilde{q},\theta \right)=\cos \theta A_1\left(\tilde{q},\theta \right)+\sin \theta A_2\left(\tilde{q},\theta \right)$ , which has five real eigenvalues, namely ${\lambda }_1=\tilde{c}\cos \theta$ , ${\lambda }_{2,3,4}=0$ , and ${\lambda }_5=-\tilde{c}\cos \theta$ , and five corresponding linearly independent right eigenvectors. The eigenmatrix R is constructed by the five right eigenvectors, and the characteristic variables can be defined as $w=R^{-1}q$ .

Multiplying the system in Equation (9) by $R^{-1}$ from the left generates an eigen-system as follows:

$\frac{\text{d}w}{\text{d}t}+B_1\left(\tilde{q}\right)\frac{\partial w}{\partial x}+B_2\left(\tilde{q}\right)\frac{\partial w}{\partial y}=r_w$ (10)

where $B_k\left(\tilde{q}\right)=R^{-1}{A}_k$ and $r_w=R^{-1}{r}_q$ , where k = 1, 2.

Equation (10) can then be transformed into the following quasi-diagonalized system:

$\frac{\text{d}w}{\text{d}t}+{\Lambda }_1\frac{\partial w}{\partial x}+{\Lambda }_2\frac{\partial w}{\partial y}=s+r_w$ (11)

where ${\Lambda }_k=\left({\lambda }_{k,1},{\lambda }_{k,2},\cdots ,{\lambda }_{k,5}\right)$ (k = 1, 2) is a diagonal matrix.

Given the initial condition at time $t_n$ , the solution at the evolution point $P={\left(x,y,t\right)}_{t=t_n+\tau }$ of the component $w_l$ (l = 1, 2, 3, 4, 5) of the characteristic variables w for Equation (11) can be obtained from the characteristic theory of the two-dimensional linear hyperbolic partial differential equations as follows:

$w_l\left(x,y,t_n+\tau ,\theta \right)=w_l\left(x-{\lambda }_{1,l}\tau ,y-{\lambda }_{2,l}\tau ,t_n\right)+{\widehat{s}}_l+{\widehat{r}}_{w,l}$ (12)

where ${\widehat{s}}_l={\displaystyle {\int }_{t_n}^{t_n+\tau }{s}_l\left[x-{\lambda }_{1,l}\left(t_n+\tau -\xi \right),y-{\lambda }_{2,l}\left(t_n+\tau -\xi \right),\xi \right]\text{d}\xi }$ and ${\widehat{r}}_{w,l}={\displaystyle {\int }_{t_n}^{t_n+\tau }{r}_{w,l}\left[x-{\lambda }_{1,l}\left(t_n+\tau -\xi \right),y-{\lambda }_{2,l}\left(t_n+\tau -\xi \right),\xi \right]\text{d}\xi }$ .

Therefore, the solution of Equation (9) may be obtained by multiplying Equation (12) with the right eigenmatrix R from the left, $q=Rw$ , and then integrating about θ from 0 to 2π [20] . The detailed expressions of the primitive variables are as follows:

$\begin{array}{l}u\left(P\right)=\frac{1}{2}u\left(\tilde{P}\right)+\frac{1}{\text{2π}}{\displaystyle \underset{0}{\overset{\text{2π}}{\int }}\left[-\frac{p\left(Q\right)}{\tilde{\rho }\tilde{c}}\cos \theta +u\left(Q\right){\cos }^2\theta +v\left(Q\right)\sin \theta \cos \theta \right]\text{d}\theta }\\ +\frac{1}{\text{2π}}{\displaystyle \underset{0}{\overset{\text{2π}}{\int }}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}S\left[z+\tilde{c}\left(t_n+\tau -\xi \right)n\left(\theta \right),\xi ,\theta \right]\cos \theta \text{d}\xi \text{d}\theta }}\\ -\frac{1}{2\tilde{\rho }}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}\frac{\partial p\left(z,\xi \right)}{\partial x}\text{d}\xi }+\frac{1}{\tilde{\rho }\tilde{c}}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}r\frac{\partial p\left(z,\xi \right)}{\partial \lambda }\text{d}\xi }\end{array}$ (13)

$\begin{array}{l}v\left(P\right)=\frac{1}{2}v\left(\tilde{P}\right)+\frac{1}{2\text{π}}{\displaystyle \underset{0}{\overset{\text{2π}}{\int }}\left[-\frac{p\left(Q\right)}{\tilde{\rho }\tilde{c}}\sin \theta +u\left(Q\right)\cos \theta \sin \theta +v\left(Q\right){\sin }^2\theta \right]\text{d}\theta }\\ +\frac{1}{\text{2π}}{\displaystyle \underset{0}{\overset{\text{2π}}{\int }}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}S\left[z+\tilde{c}\left(t_n+\tau -\xi \right)n\left(\theta \right),\xi ,\theta \right]\sin \theta \text{d}\xi \text{d}\theta }}\\ -\frac{1}{2\tilde{\rho }}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}\frac{\partial p\left(z,\xi \right)}{\partial y}\text{d}\xi }+\frac{1}{\tilde{\rho }\tilde{c}}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}r\frac{\partial p\left(z,\xi \right)}{\partial \lambda }\text{d}\xi }\end{array}$ (14)

$\begin{array}{l}p\left(P\right)=\frac{1}{\text{2π}}{\displaystyle \underset{0}{\overset{\text{2π}}{\int }}\left[p\left(Q\right)-\tilde{\rho }\tilde{c}u\left(Q\right)\cos \theta -\tilde{\rho }\tilde{c}v\left(Q\right)\sin \theta \right]\text{d}\theta }\\ -\frac{1}{\text{2π}}\tilde{\rho }\tilde{c}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}S\left[z+\tilde{c}\left(t_n+\tau -\xi \right)n\left(\theta \right),\xi ,\theta \right]\text{d}\xi \text{d}\theta }}+{\displaystyle \underset{t_n}{\overset{t_n+\tau }{\int }}r\frac{\partial p\left(z,\xi \right)}{\partial \lambda }\text{d}\xi }\end{array}$ (15)

where$\begin{array}{c}S\left(z,t,\theta \right)=\tilde{c}\left[\frac{\partial u\left(z,t,\theta \right)}{\partial x}{\sin }^2\theta +\frac{\partial v\left(z,t,\theta \right)}{\partial y}{\cos }^2\theta \right]\\ -\frac{\tilde{c}}{2}\left[\frac{\partial u\left(z,t,\theta \right)}{\partial y}+\frac{\partial v\left(z,t,\theta \right)}{\partial x}\right]\sin 2\theta \end{array}$ , $z=\left(x,y\right)$ ,

$z+\tilde{c}\left(t_n+\tau -\xi \right)n\left(\theta \right)=\left[x+\tilde{c}\left(t_n+\tau -\xi \right)\cos \theta ,y+\tilde{c}\left(t_n+\tau -\xi \right)\sin \theta \right]$ , andQ denotes the initial position ${\left(x+\tilde{c}\tau \cos \theta ,y+\tilde{c}\tau \sin \theta ,t\right)}_{t=t_n}$ .

For discretized grids, ${\theta }_{ib}$ and ${\theta }_{ie}$ are denoted as the starting and ending angles of a mesh with common vertices to update coordinates. Some approximation operations with similar procedures [19] are required to simplify the computation of the integrals in Equations (13)-(15) involving the pressure gradient terms and source terms, of which the limit operations are presented in terms of $t_n^{+}=\underset{\tau \to 0}{\lim }\left(t_n+\tau \right)$ . The ensuing analytical expressions of the vertex solver E₀ are as follows:

$\begin{array}{c}u\left(t_n^{+}\right)=\frac{1}{\text{π}}{\displaystyle \sum _{i=1}^4[-\frac{p_i}{\tilde{\rho }\tilde{c}}\left(\sin {\theta }_{ie}-\sin {\theta }_{ib}\right)}+u_i\left(\frac{{\theta }_{ie}-{\theta }_{ib}}{2}+\frac{\sin 2{\theta }_{ie}-\sin 2{\theta }_{ib}}{4}\right)\\ -v_i\frac{\cos 2{\theta }_{ie}-\cos 2{\theta }_{ib}}{4}]\end{array}$ (16)

$\begin{array}{c}v\left(t_n^{+}\right)=\frac{1}{\text{π}}{\displaystyle \sum _{i=1}^4[\frac{p_i}{\tilde{\rho }\tilde{c}}\left(\cos {\theta }_{ie}-\cos {\theta }_{ib}\right)}-u_i\frac{\cos 2{\theta }_{ie}-\cos 2{\theta }_{ib}}{4}\\ +v_i\left(\frac{{\theta }_{ie}-{\theta }_{ib}}{2}-\frac{\sin 2{\theta }_{ie}-\sin 2{\theta }_{ib}}{4}\right)]\end{array}$ (17)

$p\left(t_n^{+}\right)=\frac{1}{\text{2π}}{\displaystyle \sum _{i=1}^4[p_i\left({\theta }_{ie}-{\theta }_{ib}\right)}-\tilde{\rho }\tilde{c}{u}_i\left(\sin {\theta }_{ie}-\sin {\theta }_{ib}\right)+\tilde{\rho }\tilde{c}{v}_i\left(\cos {\theta }_{ie}-\cos {\theta }_{ib}\right)]$ (18)

where i is the counterclockwise numbering of the mesh cells with common vertices.

3.2. Validation of the Cell-Centered Lagrangian Hydrodynamic Method

1) The steady structure of a one-dimensional planar detonation wave

The Von Neumann spike values of the PBX-9502 explosive are (in cm-g-μs units) p_N = 0.375 and u_N = 0.253; and the values for the CJ state are p_CJ = 0.285, u_CJ = 0.192, and D_CJ = 0.7655. The analyzed explosive is 5.0 cm long and is initiated by the CJ condition on its left-hand side. Figure 3 presents the pressure and velocity distributions in the chemical reaction zone where the mesh sizes are Δx = 1/100, 1/200, 1/500, and 1/1000 cm respectively. According to Figure 3, the shock front of detonation wave is well-resolved and the spurious oscillation does not appear in the vicinity of the shock discontinuity. Meanwhile, mesh sizes that were less than 1/500 cm (about 50 meshes in the reaction zone) generated solutions that agreed well with the exact solutions, which indicates the good resolution and convergence of the presented method. Figure 4 presents changes in the pressure and velocity at several typical times during the course of the unsteady propagation of the detonation, in which the discretized mesh is Δx = 1/500 cm and the corresponding times are t = 0.06, 0.12, 0.24, 0.48, 0.96, 1.44, 1.92, 2.40, 2.88, 3.36, 3.84, 4.32, 4.80, and 5.28 μs. From the results, the pressure rapidly increased to reach steady state about 3.84 μs after the initiation of the CJ conditions, and the propagation velocity of the detonation shock was about 0.7670 cm/μs under steady state. The course change was almost identical with the experimental results [16] .

2) The front shape of the detonation shock wave in an explosive confined by steel

In an explosive-confiner experiment [18] , the high-speed streak camera was

utilized to record the front shape of the detonation shock wave in the insensitive explosive and the front shape of the refraction shock wave in steel. Figure 5 presents the confinement configuration, wherein PBX-9502 and steel were 7.0 cm in length, and the PBX-9502 and steel were 2.0 cm and 0.5 cm in thickness, respectively. In the simulation, the upper boundary was set as a solid wall and the explosive was initiated at the CJ conditions by a plane detonator on the right end, of which the discretized mesh had a scale of 500 cells/cm. Figure 6 presents a comparison of the shock fronts between the experimental fit by the author himself through experiments [18] (black solid line) and the numerical extraction from the pressure contours. It can be found that the simulating result were close to the experimental result.

3.3. Analysis of the Numerical Result on the Confinement Effect

3.3.1. Properties of the Flow Fields under Different Material Compressibility

The confinement materials were characterized into two groups based on their respective compressibility, of which one group exhibited a smaller sonic velocity than the CJ velocity of the detonated explosive, which exhibited a steady flow near the intersection point between the leading shock wave of detonation and the interface of the IHE/inert material. In comparison, the other group exhibited a larger sonic velocity than the CJ velocity of the detonated explosive, which exhibited an unsteady flow near the intersection point.

To compare the above confinement materials with the aforementioned theoretical analysis, confinement materials with smaller sonic velocities were also selected, specifically steel, phenolic resin, PMMA, silicone rubber foam, LiH, and LiD. The simulation results presented in Figure 7 show the leading shock wave, the sonic line, and the end line of the chemical reaction zone in the detonated explosive, as well as the pressure contours in the confinement material. Figure 7(a) presents the “strong confinement”, wherein a leading shock wave in the IHE was transmitted into the inert material with no reflected wave traveling back into the IHE. In addition, a subsonic region was observed behind the detonation shock and a supersonic region was observed behind the confinement shock. The leading shock wave curved slightly backward and exhibited a large edge angle, a thin reaction zone width, and a long subsonic distance, thereby producing alternating compression and expansion regions in the confinement material. Figure 7(b) presents the “weak confinement”, in which the leading shock wave exhibited a significant backward curve such that the edge angle was smaller and the width of the reaction zone was thicker than that of the “strong confinement” case. In addition, the sonic line and the leading shock wave intersected at the IHE/inert material interface, thereby generating a supersonic flow behind the refraction shock wave and a subsequent Prandtl-Meyer rarefaction in the confinement material. From the shock polar theory, the Prandtl-Meyer rarefaction was generated from the match between the expansion pressure of the leading shock in the IHE and the lower pressure of the refracted shock in the confinement material. Figure 7(c) presents a “one strong and one weak confinement” case, which exhibited the same flow structure as the “weak confinement” case. Figure 7(d) presents the “three-weak confinement” case, which exhibited a similar detonation flow structure as the “weak confinement” case with the exception of the existence of a subsonic region in the confinement material. The subsonic flow was derived from the match between the expansion pressure of the leading shock in the IHE and the higher refracted shock pressure in the confinement material. Figure 7(e) presents the “one-supersonic weak confinement”, which exhibited a leading shock wave that was generally curved backwards but curved forward locally near the interaction point between the leading shock wave and the interface of confinement material, thereby generating a

relatively thick reaction zone width and a subsonic region in the confinement material. The subsonic flow resulted in a higher refracted shock pressure in the confinement material. Figure 7(f) presents a “no solution confinement” case, wherein the leading shock wave and the sonic line both curved forward to generate an intersection point between the sonic line and the end line of the chemical reaction inside the detonation flow. In addition, a local subsonic region was observed in the confinement material, and the refraction shock in the confinement material “pulled” the leading shock in the IHE and moved ahead. From the shock polar theory, the ambient sonic velocity of the confinement material was lower than the detonation speed of the IHE, and the pressure induced by a normal shock traveling at a speed equal to the detonation speed was lower than the von Neumann spike pressure of the IHE. The “pulling” behavior from the Prandtl-Meyer rarefaction fan in the confinement material intersected with the supersonic weak branch of the IHE and generated obtuse edge angle.

The above numerical result of the “one strong and one weak confinement” case exhibited almost identical flow structures with the “weak confinement” case. Therefore, the two types of confinement interactions may be regarded as the same type, such that the six types of confinement interactions categorized by the improved shock polar theory can be merged into five types. The flow state of the strong branch in the “one strong and one weak confinement” case is a type of unphysical solution that cannot be conducted in reality due to its higher entropy production.

The confinement material with a larger sonic velocity is selected as Beryllium (Be), whose sonic velocity is about 0.799 cm/μs at standard state. Figure 8 presents the simulation results, in which material property is described by the Johnson-Cook constitutive model to characterize the elastic-plastic dynamics. Figure 8(a) indicates a wholly curved forward leading shock wave, an obtuse edge angle, and a sonic line that intersects with the end line of the chemical reaction zone inside the detonated explosive. Figure 8(b) is a zoomed-in version of the results presented in Figure 8(a), wherein the refraction wave in Be moved faster than the detonation wave and turned into a precursor wave such that the distance between the precursor wave and the detonation wave increased with time, thereby defining this case as unsteady state flow. Meanwhile, the precursor wave in Be generated refraction shock into the unreacted explosive to precompress the unreacted explosive. Figure 9 indicates the presence of a double-wave structure in the precursor stress wave in Be that consisted of an elastic-plastic wave with a velocity of about 1.353 cm/μs and a shock wave with a velocity of about 0.856 cm/μs. The double-wave structure is frequently observed in mineral applications [22] .

3.3.2. Edge Angle of Detonation

The edge angle of detonation is a very important parameter based on the detonation shock dynamics and is generally applied to IHE engineering calculations [23] . Table 2 presents the edge angles of detonation in PBX-9502 for several types of confinement materials. The results derived from the improved shock polar theory agreed well with the numerically simulated results for both the strong and weak confinements, in which the difference was restricted to 5%.

A notable discrepancy between the improved shock polar theory and the numerical simulation was observed in the “one supersonic weak solution” case for LiH confinement. In this case, the match of the shock polar was observed through a Prandtl-Meyer rarefaction fan in LiH into the explosive, wherein the refraction shock wave in LiH “pulled” the leading shock wave of detonation. The leading shock front exhibited a backward curved shape for the most part and a locally forward curved shape only within a very small region near the explosive interface, in which the edge angle is obtuse. The flow structure in this case is presented in Figure 10, and the edge angle in Table 2 was obtained from the fitting dashed line in the circle in Figure 10.

4. Conclusions

The presented paper demonstrated the significant confinement effect of the inert materials on the flow behavior of the system for IHEs and inert materials. The following conclusions were drawn:

1) An improved shock polar theory was established on the basis of the ZND model for IHEs detonation. The theory yielded seven categories of confinement effects. Six categories were for the inert materials that exhibit a smaller sonic velocity than the CJ velocity of the explosive detonation, which have steady flow structures and can be analyzed by the improved shock polar theory, and the remaining category was for the inert material that has a larger sonic velocity than

the CJ velocity of the explosive detonation, which has unsteady flow structure and cannot be analyzed by the improved shock polar theory.

2) A second-order, cell-centered Lagrangian hydrodynamic method is proposed on the basis of the characteristic theory of the two-dimensional first-order hyperbolic partial differential equations with chemical reaction law. The numerical method confirmed the theoretical categorization and presented the detailed flowfield structures.

a) For the inert materials that had a smaller sonic velocity than the CJ velocity of the explosive detonation, the six types of confinement effects divided by the improved shock polar theory might be merged into five types, and the edge angles of detonation could be cheaply and conveniently obtained from the improved shock polar theory with enough accuracy for the strong and weak confinement cases that are most frequently applied in engineering.

b) For the inert material that had a larger sonic velocity than the CJ velocity of the explosives detonation, the precursor wave was produced in the confinement material and generally exhibited a double-wave structure, including an elastic-plastic wave and a following shock wave. In addition, the refraction of the precursor wave could precompress the unreacted explosive.

Acknowledgements

This research is supported by Natural Science Foundation of China with Grant No. 11772066 and 11272064, Science Challenge Project with Grant No. TZZT2017-A1-HT001-F, and Foundation of Chinese Key Laboratory of Computational Physics with Grant No. 9140C 690101150C 69300.

Cite this paper

Yu, M. (2018) The Confinement Effect of Inert Materials on the Detonation of Insensitive High Explosives. World Journal of Mechanics, 8, 161-181. https://doi.org/10.4236/wjm.2018.85012

References