1 Introduction
The theory of geometrical shock dynamics (GSD) was developed in the mid 1950s by Chisnell and Whitham, and it was originally applied to the development of supersonic aircraft and hypersonic re-entry from orbit and near orbit (Chisnell Reference Chisnell1955; Whitham Reference Whitham1956, Reference Whitham1957). At the heart of GSD is the premise that a shock in an ideal gas propagates in the direction normal to itself, and the shock is sensitive to conditions only immediately behind the shock, which can be estimated from the Rankine–Hugoniot conditions. GSD theory is based on the concept of weak shock curvature in the sense that adjacent normal unit vectors on the shock are nearly parallel. These are rational approximations that assume that particle flow streamlines through the shock vary slowly in the direction perpendicular to the shock normal and that the flow is quasi-steady. Whitham’s extensive developments on GSD are well documented in his classic text, ‘Linear and Non-linear waves’ (Whitham Reference Whitham2011).
The original theory was developed in the limit of weak shock waves, which means that behind the shock the flow is nearly sonic. Hence the closure condition of GSD, also known as Whitham’s characteristic rule, uses flow states estimated from the Rankine–Hugoniot conditions and evaluates the forward characteristic to obtain the shock surface motion rule. This step generates an intrinsic, shock-surface-based evolution law, a relationship between flow area and Mach number, that is used to compute the motion of the shock surface. The original formulas were developed for an ideal gas.
Whitham demonstrated that the motion rule could be used to calculate the location of the shock in cases of shock diffraction over regular and irregular surfaces. Even though Whitham’s relationship between area and Mach number is valid only in the acoustic limit of weak shock waves, it is reasonably predictive for significantly large Mach numbers, well outside the domain of its asymptotic validity. Thus as an engineering tool GSD has been very useful for analysing nonlinear shock diffraction in many applications. In a sense, GSD is an early example of shock front tracking, whereby the leading shock in the fluid is computed a priori, without recourse to a numerical solution of the Euler equations.
Taylor blast wave theory (TBW) is another theory of shock dynamics that is based on the opposite asymptotic limit – an infinitely strong shock disturbance. TBW theory is based on a similarity solution of the Euler equations discovered independently by Taylor, Sedov and von Neumann during World War II. Taylor used it to analyse the yield of the first atomic blast, and his first unclassified paper on this problem was published in 1950 (Taylor Reference Taylor1950). Similar to GSD, the analysis was carried out for an ideal gas with applications to air shocks. Since the blast wave solution is in fact an asymptotic solution of the full Euler equations, one not only finds an asymptotic approximation to the motion of the lead shock, but also a self-similar solution for the flow field behind the lead shock. Therefore, this theory makes a prediction of the observed wind velocity at a fixed distance from the blast that is related to the total energy released in the blast. By measuring the wind velocity, one estimates the blast yield. Most importantly, since the radial shock displacement
$R(t)$
is determined explicitly as being proportional to a power of time,
$t^{2/5}$
, the radial normal shock velocity and shock acceleration are also determined with relation to the shock curvature
$\unicode[STIX]{x1D705}=2/R$
. TBW theory is valid in the limit of a very strong shock, and like GSD, it generates an intrinsic shock motion rule that relates the shock curvature to its normal speed and acceleration. Also like GSD, TBW theory is predictive and an effective engineering tool for blast calculation in regimes that are outside the asymptotic limits of its derivation.
As an interesting exercise in the theory of compressible flow, we investigate if it is possible to link the theory of GSD that is valid in the acoustic limit of weak shock waves, and the TBW theory that is valid for infinitely strong shock waves. Both theories, as shown below, have distinct motion laws for the lead shock that relate the normal shock velocity, the normal shock acceleration and the local curvature. We ask if one can link them together in the regime of intermediate shock strength, in which the shock Mach number is significantly greater than unity but not infinite. If this is possible and shown to be predictive, then one has an extended form of Whitham’s GSD theory that can be used to make improved predictions in engineering applications. There is precedent to using empirical models to extend GSD theory, such as shock propagation in water (Cates & Sturtevant Reference Cates and Sturtevant1997) and in gases with complex geometries (Aslam & Stewart Reference Aslam and Stewart1999).
A key application of interest is the design of explosive systems, which are any multi-part or multi-component systems that have both pure molecular explosives, such as the nitramines HMX or RDX, that are embedded with confining materials such as metals or plastics. Such explosive systems are used in defence, aerospace and various mining and demolition applications. They function by initiating explosive charges in such a manner that detonation waves sweep through the explosive domain, becoming bound by the inert particles. As a practical matter, modelling and numerically simulating explosive systems is a very difficult, multi-scale problem. The detonation front in the explosive materials has a thin chemical reaction zone that lies behind the leading shock. In almost all technological applications the thickness of the reaction zone is extremely thin, often
$O(10^{-3})$
or less compared to the characteristic flow dimensions (Bdzil & Stewart Reference Bdzil and Stewart2012).
Standard engineering practice to model explosive systems relies on a reduced model called ‘program burn’ (PB), originally invented by Wilkins et al. (Reference Wilkins1963), that uses a simple motion rule, equivalent to the Huygens construction of geometric optics, to propagate a detonation in the explosive material regions of the device. Self-propagating detonation shocks mainly travel at speeds close to the Chapman–Jouguet detonation velocity
$D_{CJ}$
, such that the flow at the end of the supporting reaction zone is sonic in the frame of the lead normal shock. The model simulation specifies the initial location of the detonation shock, then computes the motion of the shock according to the detonation shock propagation law. At each point in the explosive domain of the device, a time of arrival (TOA) of the detonation shock can be computed from the transient solution of the shock motion, and in practice once the shock crosses a computational cell its TOA value is stored in that cell. Strictly speaking, the TOA is a spatial field of the form TOA
$(x,y,z)$
. The TOA field depends on the geometry of the explosive, the placement of the inerts in the device and the location of the initial shock front or detonation points.
Early PB algorithms used the simple rule that the shock normal velocity is equal to the constant
$D_{CJ}$
for each explosive region, and this computed the TOA field a priori. The PB model simulation used the stored TOA values to ‘light’ the explosive at the appropriate times, when the detonation shock is supposed to cross each unreacted material node. A change in the equation of state is made at that material point, that corresponds to the amount of energy that would be released due to the reaction of the material. The action of material volume expansion then generates fluid motion in the ‘burnt’ regions where energy has been deposited to increase the fluid pressure and change the density. Computational models that use PB are effectively treating the pre-computed shock front as a propagating surface with a delta function forcing of the states across the propagating lead shock (Bdzil, Stewart & Jackson Reference Bdzil, Stewart and Jackson2001; Kapila, Bdzil & Stewart Reference Kapila, Bdzil and Stewart2006). The rational asymptotic theory of detonation shock dynamics (DSD) replaces the simple constant velocity Huygens construction with a motion rule that relates the normal shock velocity to its total curvature and leads to a modified PB algorithm that uses the DSD motion rules for the explosive shock (Bdzil & Stewart Reference Bdzil and Stewart2012).
Even though DSD modifications lead to dramatic improvements, the PB algorithm fails to be predictive when the sweeping detonation shock encounters an embedded inert. The transmitted shock passage time through that inert, after collision with the detonation shock, is shorter than the pre-calculated TOA time that only considers propagation in the explosive regions. The resolution is to compute the shock in the inert and the explosive both, to account for this discrepancy in the TOA times. This issue presents the practical need to compute both the dynamics of the detonation shock in the explosive and the dynamics of the lead shock in the inerts. This motivates the development of a theory of GSD for inerts that can be used for condensed phase materials. From this we can develop new hybrid algorithms to account for the motion of the lead shock wave in the explosive system, whether it be a detonation shock or an inert shock. This can lead to a vastly improved algorithm that can accurately approximate the states at the lead front of the system.
In § 2 we offer a mathematical summary of TBW and GSD theories, showing how we can describe both theories with a
${\dot{D}}_{n}$
–
$D_{n}$
–
$\unicode[STIX]{x1D705}$
equation of the same form. We follow up in § 3 with a numerical simulation of a radially expanding blast shock in air, occupying the intermediate region between TBW and GSD. Using a model fit we describe the shock transition from the blast wave to the acoustic state. In § 4 we describe the multi-dimensional Bryson and Gross experiment, using the experimental design to develop a composite solution for the intermediate blast range.
In the second part of this paper, beginning in § 5, we adapt the principles of GSD to metals that follow the Mie–Gruneisen equation of state. Finally, we describe a hybrid shock dynamics theory that combines GSD and DSD to simulate a hybrid explosive with embedded inert material. We test this theory with a slab explosive in § 6 and a unit cell explosive in § 7.
2 The shock dynamics relations for blast waves and geometrical shock dynamics
2.1 The shock dynamics of high intensity blast waves
Following the original blast wave analysis by (Taylor Reference Taylor1950), we describe the motion of the shock in an ideal gas as intrinsic to the shape and velocity of the shock. We assume that there is a blast energy
$E$
, initially concentrated at a single point, that is released into an ideal gas. The blast expands radially, and we designate the radius of the blast shock as
$R(t)$
. We also define the shock speed as
$\text{d}R/\text{d}t\equiv D_{n}$
, and the shock acceleration as
$\text{d}^{2}R/\text{d}t^{2}\equiv \text{d}D_{n}/\text{d}t\equiv \dot{D_{n}}$
. We assume that the shock velocity is dependent on the blast energy
$E$
and the inert gas density
$\unicode[STIX]{x1D70C}_{0}$
. Making a simple dimensional argument, the relation between
$R$
and
$t$
takes the form
$$\begin{eqnarray}\displaystyle R(t)=k\left(\frac{E}{\unicode[STIX]{x1D70C}_{0}}\right)^{1/5}t^{2/5}, & & \displaystyle\end{eqnarray}$$
where
$k$
is a constant related to the energy in the flow. Importantly, one finds that
$R$
is proportional to
$t^{2/5}$
. We represent the shock curvature as
$\unicode[STIX]{x1D705}=j/R$
, where
$j=1$
in the case of a two-dimensional (2-D) circular blast and
$j=2$
in the case of a 3-D spherical blast. With this in mind, we derive a universal relation for the shock motion as
$$\begin{eqnarray}\displaystyle \frac{{\dot{D}}_{n}R}{D_{n}^{2}}=-\frac{3}{2}\quad \text{or }{\dot{D}}_{n}=-\left(\frac{3}{2j}\right)D_{n}^{2}\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$
Although this relation is specifically derived under the assumption of radial symmetry, we note that the second form of (2.2) can be interpreted as an intrinsic motion rule for the propagation of the shock surface. This rule holds for all ideal gases, regardless of their physical properties, and in § 5.3 we prove that this rule also holds for metals that follow the Mie–Gruneisen equation of state (EOS). For a shock surface that is not radially symmetric, the TBW relation is a partial differential equation with a hyperbolic character, and it is amazingly independent of the properties of the material. We note that Taylor’s analysis also developed the spatial structure of the flow behind the blast wave, which is determined by the solutions of the self-similar Euler equations.
2.2 Geometrical shock dynamics derived in the acoustic limit
Whitham (Reference Whitham2011) gives a full account of the acoustic limit of a shock propagating down a non-uniform tube with a slowly varying cross-sectional area. By deriving and validating a characteristic rule, he derives the GSD motion rule as a relationship between the Mach number of the shock and the area of the tube. Aslam (Reference Aslam1996) describes the GSD motion rule in terms of the intrinsic properties of the shock, relating its velocity, acceleration and curvature in the form shown in (2.13). We present a quick derivation of the GSD shock motion rule for an ideal gas.
Whitham shows that in the acoustic limit of a weak shock, a characteristic forms behind the shock. With this in mind, the derivation proceeds as follows. First, we consider the intrinsic properties of the fluid: the density
$\unicode[STIX]{x1D70C}$
, the particle velocity
$u$
and the pressure
$p$
. We represent the sound speed squared as
$c^{2}$
, and for an ideal gas with constant specific heat,
$c^{2}=\unicode[STIX]{x1D6FE}p/\unicode[STIX]{x1D70C}$
, where
$\unicode[STIX]{x1D6FE}$
is the ideal gas adiabatic index. Some texts use the variable
$a$
instead of
$c$
.
Whitham expresses the GSD motion rule in the form,
$$\begin{eqnarray}\displaystyle \frac{\text{d}p}{\text{d}x}+\unicode[STIX]{x1D70C}c\frac{\text{d}u}{\text{d}x}+\frac{\unicode[STIX]{x1D70C}c^{2}u}{u+c}\frac{1}{A}\frac{\text{d}A}{\text{d}x}=0, & & \displaystyle\end{eqnarray}$$
where
$A(x)$
is the tube’s cross-sectional area. We obtain the intrinsic form of the GSD motion rule simply by replacing
$A^{\prime }(x)/A(x)$
, by the shock curvature
$\unicode[STIX]{x1D705}$
. In addition, we use the formula for the forward characteristic speed
$$\begin{eqnarray}\displaystyle \frac{\text{d}x}{\text{d}t}=u+c & & \displaystyle\end{eqnarray}$$
to substitute the total time derivative for the total
$x$
derivative along the characteristic. With this in mind, we obtain the alternate, equivalent formula
$$\begin{eqnarray}\displaystyle \frac{\text{d}p}{\text{d}t}+\unicode[STIX]{x1D70C}c\frac{\text{d}u}{\text{d}t}+\unicode[STIX]{x1D70C}c^{2}u\unicode[STIX]{x1D705}=0, & & \displaystyle\end{eqnarray}$$
where the total time derivative is interpreted as the time rate of change of the fluid states evaluated along the
$C+$
characteristics. Equation (2.5) is known as Whitham’s characteristic rule.
The Rankine–Hugoniot conditions govern the change in fluid states across the shock. They require that the mass flux, momentum flux and energy are continuous. The equations are
$$\begin{eqnarray}\displaystyle & \displaystyle u_{s}=D_{n}\left(1-\frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}_{s}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p_{s}-p_{0}=\unicode[STIX]{x1D70C}_{0}u_{s}D_{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle e(p_{s},\unicode[STIX]{x1D70C}_{s})-e(p_{0},\unicode[STIX]{x1D70C}_{0})=\frac{1}{2}(p_{s}+p_{0})\left(\frac{1}{\unicode[STIX]{x1D70C}_{0}}-\frac{1}{\unicode[STIX]{x1D70C}_{s}}\right), & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}_{s},u_{s}$
and
$p_{s}$
are the density, particle velocity and pressure behind the shock,
$\unicode[STIX]{x1D70C}_{0},u_{0}$
and
$p_{0}$
are the density, particle velocity and pressure ahead of the shock and
$D_{n}$
is the normal shock velocity. We solve these equations for
$\unicode[STIX]{x1D70C}_{s}$
,
$u_{s}$
and
$p_{s}$
in terms of
$D_{n}$
. We also express the sound speed on the shock
$c_{s}$
in terms of
$D_{n}$
.
Since
$D_{n}$
is a function of time, we compute the time derivatives of the fluid states behind the shock as
$$\begin{eqnarray}\displaystyle \frac{\text{d}p}{\text{d}t}=\frac{\text{d}p_{s}}{\text{d}D_{n}}{\dot{D}}_{n},\quad \frac{\text{d}u}{\text{d}t}=\frac{\text{d}u}{\text{d}D_{n}}{\dot{D}}_{n}. & & \displaystyle\end{eqnarray}$$
Making these substitutions into (2.5) and solving for
${\dot{D}}_{n}$
, we derive the general form of the
${\dot{D}}_{n}-D_{n}-\unicode[STIX]{x1D705}$
relation
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-\frac{p_{s}c_{s}^{2}u_{s}}{\displaystyle \frac{\text{d}p_{s}}{\text{d}D_{n}}+p_{s}a_{s}\displaystyle \frac{\text{d}u_{s}}{\text{d}D_{n}}}\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$
When we specifically evaluate for an ideal gas, we obtain the intrinsic form of the result obtained by Whitham:
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-c_{0}^{2}\frac{M^{2}-1}{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D705}, & & \displaystyle\end{eqnarray}$$
where
$M=D_{n}/c_{0}$
is the shock Mach number,
$c_{0}=\sqrt{\unicode[STIX]{x1D6FE}p_{0}/\unicode[STIX]{x1D70C}_{0}}$
is the inert sound speed and
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}=\left(1+\frac{2}{\unicode[STIX]{x1D6FE}+1}\frac{1-\unicode[STIX]{x1D707}^{2}}{\unicode[STIX]{x1D707}}\right)\left(1+2\unicode[STIX]{x1D707}+\frac{1}{M^{2}}\right),\quad \unicode[STIX]{x1D707}^{2}=\frac{(\unicode[STIX]{x1D6FE}-1)M^{2}+2}{2\unicode[STIX]{x1D6FE}M^{2}-(\unicode[STIX]{x1D6FE}-1)}.\quad & & \displaystyle\end{eqnarray}$$
This GSD motion rule applies when
$M$
is greater than but close to
$1$
. We note that equations (2.2) and (2.12) are both of the form
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-f(D_{n})\unicode[STIX]{x1D705}, & & \displaystyle\end{eqnarray}$$
and since
$f(D_{n})$
is strictly positive in both the TBW and GSD equations, the dynamics of the wave front evolution is hyperbolic in character.
2.3 Summary
In the two cases discussed, the blast wave limit and the acoustic limit, we have suggested that the dynamics of the lead shock behind a disturbance propagates according to a hyperbolic front evolution of the general form (2.13). While we can derive or postulate the two limiting cases of very strong and very weak disturbances, we ask if it is possible to find a similar evolution equation that is valid for any intermediate Mach shock number. And since one can envision a strong disturbance that starts out as a blast and decays to a weak wave, we also ask if the decaying transient of a blast wave to an acoustic wave determines a universal shock evolution equation on the shock speed.
We test this conjecture directly via simulation of the Euler equations and by comparing the numerically generated transients with the TBW and GSD theoretical solutions. This test is a simulation of a blast singularity in ambient air that expands outwards according to the Euler equations. We track the lead shock motion
$R(t)$
, then we compute the dimensionless ratio
${\dot{D}}_{n}/D_{n}^{2}\unicode[STIX]{x1D705}$
.
In the blast wave limit, according to equation (2.2), we expect this value for a spherical wave to be equal to
$-3/4$
. If we simulate the flow until the shock becomes extremely weak, this ratio transitions to
$0$
. In the next section we examine how such simulations of the relaxation from a blast wave to an acoustic disturbance compare with the two limits, for the well-studied case of an ideal gas.
(a) A radially expanding blast wave front of radius
$R$
propagating at velocity
$D_{n}$
. The variables
$\unicode[STIX]{x1D70C},u$
and
$p$
represent the density, material velocity and pressure inside the blast, and they are radially symmetric around the centre. The material outside of the blast is in its cold, inert state and it has the constant properties
$\unicode[STIX]{x1D70C}_{0},u_{0}$
and
$p_{0}$
. (b) The radius of the blast shock
$R(t)$
versus time. The TBW, GSD and intermediate regions are shown, although the boundaries between them are not explicitly defined.
3 Comparisons of simulations with blast wave and geometric shock dynamics radial transients for an ideal gas
3.1 Ideal gas test
We conduct this test using the NEWCODE software, our in-house code that numerically solves Euler equation problems. We thoroughly discuss NEWCODE and its numerical method in the appendix A. We first test the accuracy of NEWCODE by performing a simulation of a radially expanding blast wave in an ideal gas, specifically atmospheric air. Figure 1 shows a portrayal of the simulation domain, but the simulation is processed in one dimension with the assumption of three-dimensional axisymmetry. The physical properties of air are
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{0}=0.0013~\text{g}~\text{cm}^{-3},\quad p_{0}=1.02\times 10^{-4}~\text{Mbar},\quad \unicode[STIX]{x1D6FE}=1.4,\quad c_{0}=0.331~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The initial size and velocity of the hotspot is
$R=0.1$
cm and
$D_{n}=1000c_{0}=331~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
. The simulation runs up to
$t=100~\unicode[STIX]{x03BC}\text{s}$
, which is enough time for the hotspot to decay into an acoustic steady state.
(a) A simulation of a spherical blast wave expansion in atmospheric air. The TBW limit, GSD limit and simulation data are shown. (b) A dimensionless graph that shows the transition from the TBW state to the GSD state for a blast wave expansion in aluminium. The shape of the curve was found by fitting the simulation data to a model.
We track the shock location as
$R(t)$
over the course of the simulation, and from these data we derive
$\unicode[STIX]{x1D705}$
,
$D_{n}$
and
${\dot{D}}_{n}$
. In figure 2(a) we plot
${\dot{D}}_{n}/\unicode[STIX]{x1D705}$
in terms of
$D_{n}$
, up to
$D_{n}=1.5~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
. We also plot the TBW and GSD predictions of
${\dot{D}}_{n}/\unicode[STIX]{x1D705}$
, which according to equations (2.2) and (2.11), are determined uniquely by the value of
$D_{n}$
. Note that in theory, the data presented in figure 2(a) are independent of the initial values of
$R$
and
$D_{n}$
.
The blast shock initially has a high velocity, high curvature and extremely negative acceleration, as shown in the lower right corner of this figure. Note that the simulation data at this point closely match the TBW prediction. The shock front expands rapidly and loses velocity, and the shock transitions to the data shown in the upper left corner of this figure. When the shock velocity becomes very small, less than
$0.5~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, the GSD model becomes a better prediction. It is important to be clear that in terms of time, the shock lies in the TBW state very briefly. It takes only approximately
$0.13~\unicode[STIX]{x03BC}\text{s}$
for the shock velocity to decay from
$331$
to
$10~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, then another
$7.0~\unicode[STIX]{x03BC}\text{s}$
for the velocity to decay to
$1~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
. At this point the shock will gradually approach but never completely reach its ambient state of
$D_{n}=c_{0}$
.
3.2 Intermediate shock transient in dimensionless variables
Our goal is to quantitatively describe the transition from the strong shock limit to the weak shock limit by portraying the data in such a way that the TBW and GSD predictions appear as single points, with the transient appearing as a curve between them. Define the dimensionless variables
$$\begin{eqnarray}\displaystyle X=t\unicode[STIX]{x1D705}D_{n},\quad Y=\frac{{\dot{D}}_{n}}{\unicode[STIX]{x1D705}D_{n}^{2}}. & & \displaystyle\end{eqnarray}$$
In the TBW limit, the radius of the blast wave follows the equation
$R(t)=(C_{1}t)^{2/5}$
, where
$C_{1}$
is a constant, according to equation (2.1). Therefore, we compute
$X=2j/5$
and
$Y=-3/(2j)$
, where
$j=1$
for the 2-D model and
$2$
for the 3-D model. Far into the acoustic limit, when the shock velocity has decayed all the way down to the bulk sound speed
$c_{0}$
, we can model the radius of the blast wave by
$R(t)\sim C_{2}+c_{0}t$
, where
$C_{2}$
is another constant. Given enough time that constant becomes negligible, and the radius of the shock is just
$R(t)\sim c_{0}t$
. Therefore, as
$t\rightarrow \infty$
we compute
$X=j$
and
$Y=0$
.
As the shock transitions from the TBW state to the GSD state, we expect
$X$
to transition from
$2j/5$
to
$j$
, and
$Y$
to transition from
$-3/(2j)$
to
$0$
. Since both the strong shock limit and the weak shock limit have functions for the shock radius in the form of
$R(t)=Ct^{k}$
, in the intermediate transient phase we assume that
$R(t)\sim C(t)t^{k(t)}$
, where
$C(t)$
and
$k(t)$
are functions slowly varying in time. The value of
$C(t)$
is initially dependent on the total energy of the blast wave, and over time it decreases to
$c_{0}$
. The value of
$k(t)$
begins at
$2/5$
in the TBW limit and increases to
$1$
in the GSD limit.
Figure 2(b) shows the simulation data from figure 2(a) in terms of
$X$
and
$Y$
. The TBW and GSD limits are depicted as individual points whereas the curve between them is the intermediate transient. The simulation begins in the lower left corner and decays to the upper right corner, but note that the two marked points are asymptotic limits – the simulation can never completely reach either one. Also, the two asymptotic limits are independent of the properties of the gas, but the transition curve is unique to atmospheric air.
The NEWCODE simulation results show that we can fit the relationship between
$Y$
and
$X$
to a function of the form
$$\begin{eqnarray}\displaystyle Y=a-b\text{e}^{-cX}, & & \displaystyle\end{eqnarray}$$
where
$a,b,c$
are positive constants that we find through a nonlinear model fit algorithm. This particular model form was found empirically, as it well describes the convergent nature of this system. For atmospheric air, the best model fit is
$a=0.00643,b=18.2,c=3.97$
. This curve is shown superimposed on the simulation data in figure 2(b).
In summary, using NEWCODE we simulated an air blast shock that originated as a small, high energy hotspot, which rapidly expanded and decelerated, eventually decaying to a large but weak air disturbance. We plotted the radius of this shock over time, and we demonstrated that the relationship between the curvature, shock velocity and shock acceleration initially follows TBW theory, then transitions to GSD theory. In order to describe this transition quantitatively, we described the same data using dimensionless variables, then fit the simulation data to an exponential function. This model fit function provides a description of the blast shock evolution that encompasses both the strong and weak shock limits.
4 Example 1: simulation of two-dimensional shock diffraction in air
Our next goal is to determine a
${\dot{D}}_{n}$
–
$D_{n}$
–
$\unicode[STIX]{x1D705}$
relation that we can apply to shock surfaces of any shape and with intermediate shock strength. First, we replicate the Bryson and Gross experiment in NEWCODE over a range of Mach numbers. We measure the shock–shock diffraction velocity in each of these simulations, then we use these data to fit a composite model that bridges the TBW and GSD limits. We use this composite model to simulate the Bryson and Gross experiment again using a level set shock code, then we compare our two simulation methods.
The Bryson and Gross experiment, first conducted in the early 1960s, is one of the first experiments to validate GSD theory. The experiment consists of a solid object, either a cylinder, sphere or cone, placed in ambient air (Bryson & Gross Reference Bryson and Gross1961). A piston loaded shock wave, with an initial Mach number of 2.85, passes over the solid object, and the resulting shocks are observed using a schlieren light screen. The interaction of the shock front and the solid object causes a Mach shock diffraction that expands outwards. The discontinuity between the diffracted and undisturbed regions of the shock is called the shock–shock, and GSD theory provides a prediction for the diffraction rate of the shock–shock. Bryson and Gross measured the location of the shock–shock in their experiments and showed that they matched the GSD predictions within experimental error. Subsequent simulation research has verified these results and discussed a further analysis of the shock wave (Drikakis et al. Reference Drikakis, Ofengeim, Timofeev and Voionovich1997; Ofengeim & Drikakis Reference Ofengeim and Drikakis1997).
We start by simulating the Bryson and Gross experiment in NEWCODE. The simulation domain is two-dimensional, with a solid cylinder of radius 1 cm embedded in the centre. A high pressure region of air creates a blast that propagates upwards at constant speed
$D_{0}=M_{0}c_{0}$
, where
$M_{0}$
is the Mach number and
$c_{0}$
is the sound speed of air. When the shock approaches the cylinder, it becomes orthogonal to the cylindrical boundary. As shown in figure 4(a), this orthogonal boundary condition causes the Mach shock to start travelling around the cylinder at a faster velocity than the incident shock, which remains flat with velocity
$D_{0}$
. The intersection of the Mach shock and the incident shock, indicated by a dotted line, is the shock–shock. The shock–shock travels outwards at speed
$v_{ss}$
, which experiments show is approximately constant.
Suppose we can model the shock propagation with an equation of the form
${\dot{D}}_{n}=-f(D_{n})\unicode[STIX]{x1D705}$
. Then according to (A 17) in the appendix A, the shock–shock should diffract at a rate approximately equal to
$v_{ss}=\sqrt{f(D_{0})}$
. Both the TBW and GSD theories have shock propagation laws of this form, where
$f(D_{n})$
comes from (2.2) for TBW theory and equation (2.11) for GSD theory. We refer to the TBW and GSD forms of
$f(D_{n})$
as
$f_{TBW}(D_{n})$
and
$f_{GSD}(D_{n})$
respectively.
Using an asymptotic bridging method, we seek a composite
${\dot{D}}_{n}$
–
$D_{n}$
–
$\unicode[STIX]{x1D705}$
relation that is accurate in the intermediate shock region of the form
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-f_{comp}(D_{n})\unicode[STIX]{x1D705},\quad f_{comp}(D_{n})=a\text{ Erf}(x-c_{0})f_{TBW}(D_{n})+b\text{ Erfc}(x-c_{0})f_{GSD}(D_{n}), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where Erf is the error function, Erfc is the complimentary error function and
$a$
and
$b$
are model fit constants. By measuring
$v_{ss}$
in NEWCODE and fitting it to the equation
$\sqrt{f_{comp}(D_{0})}$
or
$\sqrt{f_{comp}(M_{0}c_{0})}$
, we can plot a relation between
$M_{0}$
and
$v_{ss}$
and find a model fit for the values of
$a$
and
$b$
.
Figure 3 shows a photograph of the original Bryson and Gross experiment, along with a schlieren diagram of the equivalent NEWCODE simulation, with
$M_{0}=2.85$
(Hernández, Lieberthal, & Stewart Reference Hernández, Lieberthal and Stewart2017). NEWCODE successfully tracks the location of the Mach shock, incident shock, and shock–shock, closely agreeing with the experimental data.
Schlieren diagram of shock diffraction over a cylinder: (a) the experimental density photograph from Bryson & Gross (Reference Bryson and Gross1961), and (b) the equivalent NEWCODE simulation.
We run the NEWCODE simulation at a variety of Mach numbers, ranging from
$1$
to
$3$
in increments of
$0.25$
, and from
$3$
to
$10$
in increments of
$1$
. In each simulation we measured the value of
$v_{ss}$
from the time of arrival field. Figure 4(b) shows the NEWCODE simulation values of
$M_{0}$
and
$v_{ss}$
, along with the GSD and TBW predictions of
$v_{ss}$
. As we expect, the simulation results for
$v_{ss}$
lie between the GSD and the TBW predictions. GSD theory is a better predictor for low Mach numbers, but as the Mach number increases, the TBW prediction becomes more accurate. We fit these data to (4.1) and find the model fit values
$a=0.7174,b=3.561$
. The composite solution is also shown in figure 4(b). Note that the Mach number needs to be as high as
$30$
until the NEWCODE solution agrees with the TBW prediction, however few Bryson and Gross experiments have been conducted at that velocity.
4.1 GSD/TBW composite simulation
We now use the composite model to simulate the Bryson and Gross experiment by substituting (4.1) into equation (A 14). The simulation code, written in C++, solves (4.1) using a method of lines approach, and the software is similar to what we used in Lieberthal, Bdzil & Stewart (Reference Lieberthal, Bdzil and Stewart2014) at comparable resolutions. While the shock is either below or above the cylinder, the simulation computes the shock using a Cartesian arrangement of nodes. While the shock is attached to the cylinder, the simulation uses a polar arrangement of nodes centred on the cylinder. This allows the software to process the orthogonal boundary condition without any dependence on the simulation resolution.
Figure 5 shows the results of the Bryson and Gross simulation, with initial Mach numbers of
$2.85$
and
$10$
. The two top figures show the results of the GSD/TBW composite simulation, and the two bottom figures show the equivalent simulations computed in NEWCODE. The
$M_{0}=2.85$
simulation is printed in time increments of
$\unicode[STIX]{x0394}t=0.2~\unicode[STIX]{x03BC}\text{s}$
and the
$M_{0}=10$
simulation is printed in time increments of
$\unicode[STIX]{x0394}t=0.05~\unicode[STIX]{x03BC}\text{s}$
. Comparing the time of arrival profile for both simulations results in precision within 3 % relative error, which serves as a validation of the GSD/TBW composite model. The shock–shock diffraction rate in each composite model simulation is
$0.707~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
and
$2.577~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, respectively. In comparison, the respective NEWCODE diffraction rates are
$0.719~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
and
$2.431~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
.
The two simulations disagree slightly in the region close to the surface of the sphere. This discrepancy is due to NEWCODE’s use of a Cartesian coordinate frame. NEWCODE simulates the reflective boundary condition fairly well, but it is impossible to perfectly compute a radial shock front in rectangular coordinates. This error, while slight, builds up over time to result in approximately a 10 % relative difference between the two simulation codes along the particle boundary. Running the NEWCODE simulation at higher resolutions alleviates this problem somewhat, at the cost of computational efficiency.
(a) A representation of the Bryson and Gross simulation, with the Mach shock, the incident shock and the shock–shock pictured. (b) The relation between the Mach number
$M_{0}$
and the shock–shock diffraction rate
$v_{ss}$
, as measured in NEWCODE. The GSD and TBW predictions for
$v_{ss}$
are also shown, along with the composite solution.
Shock data for the Bryson and Gross experiment with initial Mach numbers of 2.85 and 10, conducted with the composite model and with NEWCODE. Shock–shocks are indicated by dotted lines.
5 GSD composite models for the Mie–Gruneisen equations of state
5.1 Mie–Gruneisen equation of state
In this section, we apply the principles of GSD theory to non-ideal materials, specifically metals that follow the Mie–Gruneisen (MG) EOS. The MG EOS is traditionally written in the form
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}(p,v)=v\left(\frac{\text{d}p}{\text{d}e}\right)_{v}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6E4}$
is the Gruneisen parameter,
$p$
is the pressure,
$e$
is the internal energy per unit mass and
$v$
is the specific volume. We assume that
$\unicode[STIX]{x1D6E4}$
is dependent only on
$v$
, and we rewrite the equation in terms of the density
$\unicode[STIX]{x1D70C}$
, to find the equation
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D70C})=\unicode[STIX]{x1D6E4}_{0}\left(\frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}}\right)^{q}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6E4}_{0}$
is the Gruneisen constant. We further assume that
$q=1$
, which is a reasonable assumption for most metals (Miller & Puckett Reference Miller and Puckett1996; Mitchell Reference Mitchell1981).
With this in mind, we define the temperature corrected version of the MG EOS as
$$\begin{eqnarray}\displaystyle p(e,v)=\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D70C})\unicode[STIX]{x1D70C}e+f(\unicode[STIX]{x1D70C})=\unicode[STIX]{x1D6E4}_{0}\unicode[STIX]{x1D70C}_{0}e+f(\unicode[STIX]{x1D70C}). & & \displaystyle\end{eqnarray}$$
The function
$f(\unicode[STIX]{x1D70C})$
is defined conditionally depending on whether the density is greater than or less than the inert density (Arienti, Morano & Shepherd Reference Arienti, Morano and Shepherd2004).
$$\begin{eqnarray}\displaystyle f(\unicode[STIX]{x1D70C})=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{\unicode[STIX]{x1D70C}_{0}c_{0}^{2}\unicode[STIX]{x1D712}}{(1-s\unicode[STIX]{x1D712})^{2}}\left[1-\frac{\unicode[STIX]{x1D6E4}_{0}}{2}\unicode[STIX]{x1D712}\right]\quad & \text{if }\unicode[STIX]{x1D70C}\geqslant \unicode[STIX]{x1D70C}_{0}\\ c_{0}^{2}(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})\quad & \text{if }\unicode[STIX]{x1D70C}<\unicode[STIX]{x1D70C}_{0},\end{array}\right. & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D712}=1-(\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x1D70C}),c_{0}$
is the ambient bulk speed of sound and
$s$
is the Rankine–Hugoniot slope coefficient.
The constant
$s$
is defined as the slope of a
$U_{p}-U_{s}$
curve, found through piston driven shock experiments. Noting that the piston velocity
$U_{p}$
and the shock velocity
$U_{s}$
are analogous to the particle velocity
$u$
and the shock velocity
$D_{n}$
in this paper, this particular form of the MG EOS is defined on the assumption that
$u$
and
$D_{n}$
are linearly related:
$$\begin{eqnarray}\displaystyle D_{n}=c_{0}+su. & & \displaystyle\end{eqnarray}$$
For our problem of interest, we normally only expect
$\unicode[STIX]{x1D70C}$
to be greater than
$\unicode[STIX]{x1D70C}_{0}$
. We can also write the MG equation in the form
$$\begin{eqnarray}\displaystyle e(p,v)=\frac{1}{\unicode[STIX]{x1D6E4}_{0}\unicode[STIX]{x1D70C}_{0}}(p-f(\unicode[STIX]{x1D70C})). & & \displaystyle\end{eqnarray}$$
The upper limit of
$\unicode[STIX]{x1D70C}$
for which this EOS holds is
$\unicode[STIX]{x1D70C}_{max}=s\unicode[STIX]{x1D70C}_{0}/(s-1)$
, and
$f(\unicode[STIX]{x1D70C})$
is not defined above this limit.
We define the sound speed as
$$\begin{eqnarray}\displaystyle c^{2}=\frac{p+\unicode[STIX]{x2202}e/\unicode[STIX]{x2202}v}{\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x2202}e/\unicode[STIX]{x2202}p}=\frac{p-\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x2202}e/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x2202}e/\unicode[STIX]{x2202}p}. & & \displaystyle\end{eqnarray}$$
For this particular EOS, the sound speed simplifies to
$$\begin{eqnarray}\displaystyle c^{2}=\frac{p+\unicode[STIX]{x1D70C}^{2}\displaystyle \frac{1}{\unicode[STIX]{x1D6E4}_{0}\unicode[STIX]{x1D70C}_{0}}f^{\prime }(\unicode[STIX]{x1D70C})}{\unicode[STIX]{x1D70C}^{2}\displaystyle \frac{1}{\unicode[STIX]{x1D6E4}_{0}\unicode[STIX]{x1D70C}_{0}}}=\unicode[STIX]{x1D6E4}_{0}\unicode[STIX]{x1D70C}_{0}\frac{p}{\unicode[STIX]{x1D70C}^{2}}+f^{\prime }(\unicode[STIX]{x1D70C}). & & \displaystyle\end{eqnarray}$$
We also compute
$$\begin{eqnarray}\displaystyle f^{\prime }(\unicode[STIX]{x1D70C})=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{c_{0}^{2}\unicode[STIX]{x1D70C}_{0}^{2}[(\unicode[STIX]{x1D6E4}_{0}+s)(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{0})+\unicode[STIX]{x1D70C}]}{[s(\unicode[STIX]{x1D70C}_{0}-\unicode[STIX]{x1D70C})+\unicode[STIX]{x1D70C}]^{3}}\quad & \text{if }\unicode[STIX]{x1D70C}\geqslant \unicode[STIX]{x1D70C}_{0}\\ c_{0}^{2}\quad & \text{if }\unicode[STIX]{x1D70C}<\unicode[STIX]{x1D70C}_{0}.\end{array}\right. & & \displaystyle\end{eqnarray}$$
5.2 Rankine–Hugoniot relations
Define
$\unicode[STIX]{x1D70C}_{s},u_{s},p_{s}$
as the density, particle velocity and pressure on the surface of the blast wave shock front. Then the Rankine–Hugoniot laws govern the jump conditions across the wave front, ensuring that the mass flux, momentum flux and energy are continuous along the shock:
$$\begin{eqnarray}\displaystyle & \displaystyle u_{s}=D_{n}\left(1-\frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}_{s}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p_{s}-p_{0}=\unicode[STIX]{x1D70C}_{0}u_{s}D_{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle e(p_{s},\unicode[STIX]{x1D70C}_{s})-e(p_{0},\unicode[STIX]{x1D70C}_{0})=\frac{1}{2}(p_{s}+p_{0})\left(\frac{1}{\unicode[STIX]{x1D70C}_{0}}-\frac{1}{\unicode[STIX]{x1D70C}_{s}}\right). & \displaystyle\end{eqnarray}$$
In general we specify that
$u_{0}=0$
, and in the case of an MG metal we assume that
$p_{0}=0$
.
By combining the Rankine–Hugoniot equations and the
$e(p,\unicode[STIX]{x1D70C})$
function for the MG EOS, we can solve for
$\unicode[STIX]{x1D70C}_{s},u_{s},p_{s}$
in terms of
$D_{n}$
. The solution is
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}(D_{n})=\frac{D_{n}s\unicode[STIX]{x1D70C}_{0}}{c_{0}+D_{n}(s-1)},\quad u_{s}(D_{n})=\frac{D_{n}-c_{0}}{s},\quad p_{s}(D_{n})=\frac{D_{n}(D_{n}-c_{0})\unicode[STIX]{x1D70C}_{0}}{s}\qquad \quad & & \displaystyle\end{eqnarray}$$
and the equation for the energy per unit mass is
$$\begin{eqnarray}\displaystyle e_{s}=e(p_{s},\unicode[STIX]{x1D70C}_{s})=\frac{(D_{n}-c_{0})^{2}}{2s^{2}}. & & \displaystyle\end{eqnarray}$$
Note that the equation for
$u_{s}(D_{n})$
is consistent with the linear shock Hugoniot equation (5.5). Also note that
$\unicode[STIX]{x1D70C}_{s}$
approaches
$\unicode[STIX]{x1D70C}_{0}$
and
$u_{s},p_{s}$
and
$e_{s}$
approach
$0$
as
$D_{n}\rightarrow c_{0}$
, as we would expect for a nearly inert material. We also substitute these solutions into our equation for the sound speed to find
$$\begin{eqnarray}\displaystyle c_{s}^{2}=\frac{[c_{0}+D_{n}(s-1)]^{2}[D_{n}^{2}(\unicode[STIX]{x1D6E4}_{0}+2s)-D_{n}sc_{0}-\unicode[STIX]{x1D6E4}_{0}c_{0}^{2}]}{D_{n}c_{0}s^{3}}. & & \displaystyle\end{eqnarray}$$
5.3 Taylor blast wave limit
It has been established that the non-corrected version of the MG EOS (with
$f(\unicode[STIX]{x1D70C})=0$
) is self-similar, but the corrected version is not (Holm & Logan Reference Holm and Logan1983). However, in the asymptotic limit of high
$D_{n}$
, we may take our solutions to the Rankine–Hugoniot relations and approximate to the highest order of
$D_{n}$
:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{s}(D_{n})\rightarrow \frac{s\unicode[STIX]{x1D70C}_{0}}{s-1},\quad u_{s}(D_{n})\rightarrow \frac{D_{n}}{s},\quad p_{s}(D_{n})\rightarrow \frac{D_{n}^{2}\unicode[STIX]{x1D70C}_{0}}{s}. & & \displaystyle\end{eqnarray}$$
We also find the asymptotic limit of the energy:
$$\begin{eqnarray}\displaystyle e_{s}(D_{n})\rightarrow \frac{D_{n}^{2}}{2s^{2}}. & & \displaystyle\end{eqnarray}$$
With this in mind, the following equation holds in the high velocity limit:
$$\begin{eqnarray}\displaystyle e_{s}\sim \frac{p_{s}}{2(s-1)\unicode[STIX]{x1D70C}_{s}}. & & \displaystyle\end{eqnarray}$$
Note that this is the same equation as the ideal EOS with
$2(s-1)=\unicode[STIX]{x1D6FE}-1$
. In other words, under a high velocity blast wave, a Mie–Gruneisen metal behaves like an ideal gas. The traditional TBW theory applies, and the radius of the blast wave front is given by the equation
$R(t)=(Ct)^{2/5}$
. Therefore, the
${\dot{D}}_{n}-D_{n}-\unicode[STIX]{x1D705}$
relation is
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-\left(\frac{3}{2j}\right)D_{n}^{2}\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$
For this equation to be reasonably precise, we require the Mach number,
$M=D_{n}/c_{0}$
to be at least approximately 20. In a heterogeneous explosive, the Chapman–Jouguet velocity which drives the experiment is of the order of magnitude of
$1~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, so the Mach number rarely exceeds
$2$
or
$3$
. Therefore, the TBW limit is not generally relevant in modelling hybrid explosives.
5.4 Geometrical shock dynamics limit
In § 2.2 we derived the general form for the GSD shock propagation law:
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-\frac{\unicode[STIX]{x1D70C}_{s}c_{s}^{2}u_{s}}{\displaystyle \frac{\text{d}p_{s}}{\text{d}D_{n}}+\unicode[STIX]{x1D70C}_{s}c_{s}\frac{\text{d}u_{s}}{\text{d}D_{n}}}\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$
Equation (5.20) is valid for any equation of state in the acoustic limit, and this result agrees with the equation found in Saenz, Taylor & Stewart (Reference Saenz, Taylor and Stewart2012) for an inert material. We substitute our Rankine–Hugoniot solutions for the MG EOS to find
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{D}}_{n}=-\frac{D_{n}c_{s}^{2}(D_{n}-c_{0})s}{2D_{n}^{2}(s-1)+D_{n}c_{s}s+D_{n}c_{0}(3-s)-c_{0}^{2}}\unicode[STIX]{x1D705}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{s}^{2}=\frac{[c_{0}+D_{n}(s-1)]^{2}[D_{n}^{2}(\unicode[STIX]{x1D6E4}_{0}+2s)-D_{n}sc_{0}-\unicode[STIX]{x1D6E4}_{0}c_{0}^{2}]}{D_{n}c_{0}s^{3}}. & \displaystyle\end{eqnarray}$$
We now proceed to test our TBW and GSD theories with a simulated blast in aluminium, analogous to the ideal air simulations we presented in § 3.1.
5.5 Aluminium test
Since we have validated NEWCODE’s simulation method in § 3.1, we proceed to run the same test with aluminium. The Mie–Gruneisen parameters for aluminium (Steinberg Reference Steinberg1996) are
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{0}=2.785~\text{g}~\text{cm}^{-3},\quad c_{0}=0.5328~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1},\quad s=1.338,\quad \unicode[STIX]{x1D6E4}_{0}=2.00.\qquad \quad & & \displaystyle\end{eqnarray}$$
We set up the model identically to the ideal gas test but with an initial shock velocity of only
$D_{n}=10c_{0}=5.328~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
. The initial blast radius is still 0.1 cm. As in the ideal gas test, the blast wave state lasts for an extremely short period of time. After about
$0.32~\unicode[STIX]{x03BC}\text{s}$
, the shock front has slowed down to 20 % of its original velocity. Around this time the wave starts to gradually decay to its acoustic limit. It takes approximately
$4.02~\unicode[STIX]{x03BC}\text{s}$
until the blast wave velocity is within 10 % of
$c_{0}$
. As
$t\rightarrow \infty$
, the system will continue to approach but never quite reach its acoustic limit.
A plot of
${\dot{D}}_{n}/\unicode[STIX]{x1D705}$
versus
$D_{n}$
of the simulation data, the GSD prediction, and the TBW prediction up to
$D_{n}=2~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
are given in figure 6(a), which is analogous to figure 2(a). In the low velocity limit,
$D_{n}<1~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, the GSD prediction almost perfectly describes the wave decay behaviour. At higher velocities, the simulation graph begins to run parallel to the TBW prediction, but the TBW prediction does not serve as a good approximation unless the velocity is much higher. Since the typical Chapman–Jouguet velocity
$D_{CJ}$
of a high explosive is usually less than
$1~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, the GSD relation alone is an adequate estimation for most heterogeneous explosives and there is no practical need to find a composite solution.
The transition from the TBW state to the GSD state can be modelled using the same dimensionless variables as in the ideal gas test, and equation (3.3) still holds. For aluminium, the best model fit is
$a=0.0166,b=9.86,c=3.19$
.
Since we have derived a GSD shock motion law for metals and shown that the law agrees with the results of an Euler equation simulation, we now seek to simulate heterogeneous materials consisting of high explosive fluids and inert metals.
(a) A simulation of a spherical blast wave expansion in aluminium. The TBW limit, GSD limit and simulation data are shown. (b) A dimensionless graph that shows the transition from the TBW state to the GSD state for a blast wave expansion in aluminium. The shape of the curve was found by fitting the simulation data to a model.
6 Example 2: simulation of an infinite two-dimensional slab with a DSD/GSD hybrid model
In this section, we seek to combine the GSD shock motion law for metals with the DSD shock motion law for high explosives (HE) in order to simulate the shock propagation through a hybrid explosive/inert material. In Lieberthal et al. (Reference Lieberthal, Bdzil and Stewart2014) we used DSD to simulate a detonation shock propagating around a solid spherical particle embedded in an HE. Our goal is to use GSD to expand upon this simulation, incorporating the shock propagation through the particle itself.
A two-dimensional infinite slab. The HE is modelled using detonation shock dynamics and the metal is modelled using geometrical shock dynamics.
To briefly summarize DSD, we typically model high explosives with a shock motion rule that directly relates the shock normal velocity and the curvature. The simplest such law is the linear
$D_{n}-\unicode[STIX]{x1D705}$
equation,
$$\begin{eqnarray}\displaystyle D_{n}=D_{CJ}(1-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D705}), & & \displaystyle\end{eqnarray}$$
where
$D_{CJ}$
is the Chapman–Jouguet velocity of the HE and
$\unicode[STIX]{x1D6FC}$
is a curvature correction constant related to the reaction zone length (Bdzil, Lieberthal & Stewart Reference Bdzil, Lieberthal and Stewart2010). Substituting into (A 14), this DSD law is a parabolic partial differential equation for
$g$
that is first order in time and second order in space. Since this law is only dependent on
$D_{n}$
and not
${\dot{D}}_{n}$
, unlike the GSD laws discussed in the previous section, the DSD law is self-sustaining. The high energy of the detonation shock causes the unreacted fluid to transition into a product phase, releasing more energy into the shock. Any shock front that is left undisturbed will eventually revert back into a flat wave with zero curvature and constant speed
$D_{CJ}$
.
For the simplest example of an HE/inert hybrid, we consider the two-dimensional infinite slab shown in figure 7. The material in the bottom and left halves of the domain (
$z<0$
or
$x<0$
) is an HE fluid, and the material in the upper right quadrant (
$z>0$
and
$x>0$
) is a metal. Initially, before
$t=0$
, the shock front is flat, propagating through the HE at a constant speed. At time
$t=0$
, the domain splits into HE on the left and metal on the right.
We model the incident shock propagation in metal using the GSD law in equation (5.21), as for most explosive fluids, we do not expect
$D_{CJ}$
to be high enough for the TBW limit to apply. Then similarly to § 4, when the shock front is only slightly perturbed, we expect to see a travelling wave with a characteristic velocity equal to
$\sqrt{f_{GSD}(D_{0})}$
.
6.1 Shock polar analysis
The values of the confinement angles
$\unicode[STIX]{x1D714}_{c}$
and
$\unicode[STIX]{x1D714}_{i}$
depend on the material properties of the HE and the metal and the speed of the shock front, and we find these values through shock polar analysis. Figure 8 shows a shock propagating through an HE/metal interface. The point where the shock and interface meet travels vertically at speed
$D_{0}$
, and the shock normal velocities in the HE and the metal are equal to
$D_{0}\sin (\unicode[STIX]{x1D714}_{c})$
and
$D_{0}\sin (\unicode[STIX]{x1D714}_{i})$
respectively.
The shock causes the interface behind it to deflect on the inert side. Using the Rankine–Hugoniot relations, we can find values for the pressure behind the shock and the deflection of the interface
$\unicode[STIX]{x1D703}$
, independently on the HE and the inert sides. Define
$p_{1}$
and
$\unicode[STIX]{x1D70C}_{1}$
as the pressure and density behind the shock on the HE side and
$p_{2}$
and
$\unicode[STIX]{x1D70C}_{2}$
as the pressure and density behind the shock on the inert side.
If the HE is modelled using a Jones–Wilkins–Lee EOS, described in § 6.2, then
$p_{1}$
and
$\unicode[STIX]{x1D70C}_{1}$
must be determined by solving the Rankine–Hugoniot relations numerically for different values of
$\unicode[STIX]{x1D714}_{c}$
. The pressure on the inert side is given by
$$\begin{eqnarray}\displaystyle p_{2}=\frac{D_{0}\sin (\unicode[STIX]{x1D714}_{i})(D_{0}\sin (\unicode[STIX]{x1D714}_{i})-c_{0})\unicode[STIX]{x1D70C}_{0}}{s}, & & \displaystyle\end{eqnarray}$$
and the density is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}_{2}=\frac{D_{0}\sin (\unicode[STIX]{x1D714}_{i})s\unicode[STIX]{x1D70C}_{0}}{c_{0}+D_{0}\sin (\unicode[STIX]{x1D714}_{i})(s-1)}. & & \displaystyle\end{eqnarray}$$
The deflection of the interface can be computed from either the HE or the inert side:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703} & = & \displaystyle \unicode[STIX]{x1D714}_{c}-\tan ^{-1}\left[\displaystyle \frac{\unicode[STIX]{x1D70C}_{0,HE}}{\unicode[STIX]{x1D70C}_{1}}\tan (\unicode[STIX]{x1D714}_{c})\right]\quad \text{or}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x03C0}-\unicode[STIX]{x1D714}_{i}+\tan ^{-1}\left[\displaystyle \frac{\unicode[STIX]{x1D70C}_{0,inert}}{\unicode[STIX]{x1D70C}_{2}}\tan (\unicode[STIX]{x1D714}_{i})\right].\end{eqnarray}$$
By setting
$p_{1}=p_{2}$
and balancing the two equations for
$\unicode[STIX]{x1D703}$
, we can create a shock polar parametric plot and find values for
$\unicode[STIX]{x1D714}_{c}$
and
$\unicode[STIX]{x1D714}_{i}$
simultaneously. A more thorough mathematical explanation can be found in Brown & Ravichandran (Reference Brown and Ravichandran2013).
A representation of a shock passing through an HE/inert interface. We can determine the values of
$\unicode[STIX]{x1D714}_{c}$
and
$\unicode[STIX]{x1D714}_{i}$
by balancing the shock pressure and the deflection angle.
6.2 Slab simulation
For our simulation we use PBXN-9 as the high explosive and aluminium as the metal. We describe the detonation of PBXN-9 with an ignition and growth model. The energy equation is the Jones–Wilkins–Lee (JWL) EOS:
$$\begin{eqnarray}\displaystyle e=\frac{1}{\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}}\left[p-A\left(1-\frac{\unicode[STIX]{x1D714}}{R_{1}}\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}}\right)\exp \left(-R_{1}\frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}}\right)-B\left(1-\frac{\unicode[STIX]{x1D714}}{R_{2}}\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}}\right)\exp \left(-R_{2}\frac{\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}}\right)\right],\quad & & \displaystyle\end{eqnarray}$$
where
$A,B,R_{1},R_{2}$
and
$\unicode[STIX]{x1D714}$
are physical constants. The reactant and product PBXN-9 each have a different set of constants. The reaction rate law is
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}\displaystyle \frac{\text{d}\unicode[STIX]{x1D706}}{\text{d}t}\; & =\; & \displaystyle \frac{\text{d}\unicode[STIX]{x1D706}_{ignition}}{\text{d}t}+\frac{\text{d}\unicode[STIX]{x1D706}_{growth}}{\text{d}t}+\frac{\text{d}\unicode[STIX]{x1D706}_{completion}}{\text{d}t},\\ \displaystyle \frac{\text{d}\unicode[STIX]{x1D706}_{ignition}}{\text{d}t}\; & =\; & \left\{\begin{array}{@{}ll@{}}I(1-\unicode[STIX]{x1D706})^{b}\left(\displaystyle \frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}}-1-a\right)^{x}\quad & \text{if }0<\unicode[STIX]{x1D706}<F_{igmax}\\ 0\quad & \text{otherwise},\end{array}\right.\\ \displaystyle \frac{\text{d}\unicode[STIX]{x1D706}_{growth}}{\text{d}t}\; & =\; & \left\{\begin{array}{@{}ll@{}}G_{1}(1-\unicode[STIX]{x1D706})^{c}\unicode[STIX]{x1D706}^{d}p^{y}\quad & \text{if }0<\unicode[STIX]{x1D706}<F_{G1max}\\ 0\quad & \text{otherwise},\end{array}\right.\\ \displaystyle \frac{\text{d}\unicode[STIX]{x1D706}_{completion}}{\text{d}t}\; & =\; & \left\{\begin{array}{@{}ll@{}}G_{2}(1-\unicode[STIX]{x1D706})^{e}\unicode[STIX]{x1D706}^{g}p\quad & \text{if }F_{G2min}<\unicode[STIX]{x1D706}<1\\ 0\quad & \text{otherwise},\end{array}\right.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D706}$
is the fraction reacted and
$I,G_{1},G_{2},F_{igmax},F_{G1max},F_{G2min},a,b,c,d,e,g,x,y$
and
$z$
are physical constants (Tarver & Urtiew Reference Tarver and Urtiew2010).
From these equations or from experiment, we can compute the DSD constants for PBXN-9 as
$D_{CJ}=0.8559~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
and
$\unicode[STIX]{x1D6FC}=0.07948$
cm (Hull Reference Hull1997). PBXN-9 also has a sonic angle of
$\unicode[STIX]{x1D714}_{s}=0.887$
(Holman Reference Holman2010). For aluminium, we use the MG EOS and the GSD law described in § 5.5. Through shock polar analysis, we find the confinement angles at speed
$D_{CJ}$
to be
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{c}=1.373,\quad \unicode[STIX]{x1D714}_{i}=1.768. & & \displaystyle\end{eqnarray}$$
Figure 9 shows our simulation results up to
$t=3.5~\unicode[STIX]{x03BC}\text{s}$
, in increments of
$0.1~\unicode[STIX]{x03BC}\text{s}$
. The shock front on the boundary converges almost immediately to a speed of
$0.84~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, or approximately
$98\,\%$
of
$D_{CJ}$
. The angle boundary conditions on the interface create a disturbance that propagates laterally through both mediums. The characteristic velocity of this propagation is
$0.43~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
in the HE and
$0.47~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
in the metal. Equation (A 17) predicted the characteristic velocity in the metal to be
$0.47~\text{cm}~\unicode[STIX]{x03BC}\text{s}^{-1}$
, so this validates our prediction that the shock front perturbation can be modelled with the wave equation.
The shock front in the HE is self-similar. The shape of the front is a parabola that stretches to the left at the characteristic speed. The shock front in the metal, on the other hand, transitions from a concave to a convex shape, then becomes flat, except in the region immediately adjacent to the boundary where the angle condition holds. A comparison of the shock time of arrival with the equivalent ALE3D simulation gives very similar results, with an average relative error of approximately
$2\,\%$
.
Simulation results of a shock passing through a PBXN-9/aluminium interface. The dotted lines represent the characteristics in the HE and the metal.
7 Example 3: high explosive unit cell array embedded with inert spherical particles
We can use the same principles to consider a series of unit cells of PBXN-9, each embedded with an aluminium sphere (Lieberthal et al.
Reference Lieberthal, Bdzil and Stewart2014). The simulation domain is now axisymmetric instead of two-dimensional, so the level set equations change slightly. We define
$z$
as the height of the shock front within the unit cell and
$r$
as the distance from the centre, so that the level set function and (A 14) described in the appendix A are now:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}(r,z,t)=z-g(r,t)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D705}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|}\right)=-\frac{g_{rr}}{(1+g_{r}^{2})^{3/2}}-\frac{g_{r}}{r\sqrt{1+g_{r}^{2}}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle D_{n}=-\frac{\unicode[STIX]{x1D713}_{t}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|}=\frac{g_{t}}{(1+g_{r}^{2})^{1/2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{D}}_{n}=\frac{\text{d}D_{n}}{\text{d}t}=\frac{g_{tt}}{(1+g_{r}^{2})^{1/2}}-\frac{g_{r}g_{t}g_{rt}}{(1+g_{r}^{2})^{3/2}}. & \displaystyle\end{eqnarray}$$
Our model consists of a series of cylindrical unit cells with a diameter and height of
$1$
cm. An inert spherical particle is embedded in the axial centre of each unit cell. Because this model is axisymmetric, we may consider it as a two-dimensional simulation, and the shock front as a space curve embedded in two dimensions.
As in the two-dimensional slab problem, the shock propagation in the HE is governed by a linear
$D_{n}-\unicode[STIX]{x1D705}$
law:
$$\begin{eqnarray}\displaystyle D_{n}=D_{CJ}(1-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D705}) & & \displaystyle\end{eqnarray}$$
and in the metal sphere is governed by the GSD law:
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-f_{GSD}(D_{n})\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$
Since shock polar analysis tends to assume a flat interface, we explore two schools of thought when dealing with the curved boundary interface. The first method assumes that once the shock front reaches its sonic state, we can treat the sphere boundary as if it is locally flat in the region of the point of attachment. This method is shown in figure 10 and is modelled as follows.
A flat detonation wave in an HE medium will normally travel upwards at constant speed
$D_{CJ}$
. The interaction between the detonation wave and an HE/inert interface is governed by the angle between the shock front and the particle surface
$\unicode[STIX]{x1D714}$
, as shown in figure 10. The standard DSD angle boundaries are such that when
$\unicode[STIX]{x1D714}<\unicode[STIX]{x1D714}_{s}$
, we assume a continuous outflow condition in which the wave effectively behaves as though the particle is not present. When
$\unicode[STIX]{x1D714}$
approaches the sonic attachment angle
$\unicode[STIX]{x1D714}_{s}$
,
$\unicode[STIX]{x1D714}$
is fixed to the value of the confinement angle
$\unicode[STIX]{x1D714}_{c}$
. In the case of a weak confinement (plastic) material, we set
$\unicode[STIX]{x1D714}_{c}=\unicode[STIX]{x1D714}_{s}$
, so the attachment angle is constrained to the value of
$\unicode[STIX]{x1D714}_{s}$
. In the case of a strong confinement (metal) material, we have a value of
$\unicode[STIX]{x1D714}_{c}$
that is greater than
$\unicode[STIX]{x1D714}_{s}$
, so the angle will experience a sharp increase at the point of attachment (Bdzil et al.
Reference Bdzil, Lieberthal and Stewart2010).
The shock front in the inert material has a similar outflow interface condition while the wave in the HE is in its detached state, as shown in figure 10. When the wave becomes attached, that is, when
$\unicode[STIX]{x1D714}$
approaches
$\unicode[STIX]{x1D714}_{s}$
, we confine the boundary angle on the inert side to be
$\unicode[STIX]{x1D714}_{i}$
. The value of
$\unicode[STIX]{x1D714}_{i}$
depends on the shock propagation velocity tangent to the surface of the sphere, but this value does not tend to vary much. This angle condition gets relaxed slightly when the shock front reaches the top half of the particle to prevent the shock front from becoming completely vertical.
However, when we compare this simulation to a traditional direct numerical simulation (DNS) method, as described in § 7.1, we find that the attachment angles are not constant. For our second method, we ran the DNS model first using the ALE3D code and manually compute the attachment angles on the HE and metal sides. We found that the values of
$\unicode[STIX]{x1D714}_{c}$
and
$\unicode[STIX]{x1D714}_{i}$
can be modelled as piecewise linear functions of
$\unicode[STIX]{x1D719}$
, the angle on the sphere of the point of attachment. For the PBXN-9/aluminium model, these functions are shown in figure 11. There is not currently a theory that can predict these functions, but by artificially enforcing these conditions in the DSD/GSD hybrid model, we can find a much better match with the DNS model.
Regardless of the method we use, the wave travels around the sphere, experiencing a net decrease in velocity, until it reaches the top of the sphere and reattaches to itself. The wave then travels upwards until it approaches the next unit cell and interacts with another particle. Further explanation of this simulation model can be found in Lieberthal et al. (Reference Lieberthal, Bdzil and Stewart2014).
At time
$t_{1}$
, the interface angle
$\unicode[STIX]{x1D714}$
is less than the sonic angle, so no boundary condition is enforced on the interface. When
$\unicode[STIX]{x1D714}$
approaches the sonic angle
$\unicode[STIX]{x1D714}_{s}$
at time
$t_{2}$
, we enforce a confinement condition that sets the external angle to
$\unicode[STIX]{x1D714}_{c}$
and the internal angle to
$\unicode[STIX]{x1D714}_{i}$
. The shock front propagates up through both media, as shown at time
$t_{3}$
, until it reaches the top of the particle. The wave attaches with itself at time
$t_{4}$
and becomes detached from the sphere. It then travels upwards into the next unit cell.
7.1 Simulation test
To form a basis of comparison, we first ran the unit cell simulation with ALE3D using a DNS method. we use the reactive flow model described in equations (6.5) and (6.6) to model the PBXN-9, and the Mie–Gruneisen model for aluminium. The simulation was performed in a two-dimensional domain under the assumption of axisymmetry. The shock front results of this simulation are shown in figure 12(c).
At each time step, we manually computed the angle between the shock front and the sphere boundary, on both the interior and exterior of the sphere. We plotted these data in terms of
$\unicode[STIX]{x1D719}$
, the angle on the sphere at which the shock is attached. Note that
$\unicode[STIX]{x1D719}=-\unicode[STIX]{x03C0}/2,0,\unicode[STIX]{x03C0}/2$
represent the south pole, equator and north pole of the sphere, respectively. These data were modelled to a series of piecewise linear functions, as shown in figure 11. The exterior attachment angle is originally at 0, meaning the shock is parallel to the sphere boundary. The angle becomes obtuse at exactly
$\unicode[STIX]{x1D719}=0$
, then climbs up to an attachment angle of
$\unicode[STIX]{x03C0}$
, meaning the shock is again parallel but in the opposite direction. The interior attachment angle begins at
$\unicode[STIX]{x03C0}$
, becomes acute at
$\unicode[STIX]{x1D719}=0.53$
then drops to
$0$
at the top of the sphere. The angles nearly but not quite complement each other, that is
$\unicode[STIX]{x1D714}_{c}+\unicode[STIX]{x1D714}_{i}$
is close to but not exactly
$\unicode[STIX]{x03C0}$
at all times. The subtle differences between the two angles are what result in the shock propagation patterns.
The measured attachment angles on the HE (
$\unicode[STIX]{x1D714}_{c}$
) and aluminium side (
$\unicode[STIX]{x1D714}_{i}$
) as a function of
$\unicode[STIX]{x1D719}$
for the unit cell simulation performed in ALE3D.
For the DSD/GSD hybrid model, we used the same physical parameters as in the two-dimensional slab simulation in § 6.2. We ran the simulation using the two angle boundary methods discussed in § 7, first using the traditional sonic state confinement condition and then by artificially matching the boundary angles produced by the ALE3D simulation.
Figure 12 shows a comparison of our simulations, with the wave front displayed in increments of
$0.1~\unicode[STIX]{x03BC}\text{s}$
. It is clear by inspection that the second angle boundary method produces a much closer match to the ALE3D simulation. However, all three models capture the transition from a convex to concave wave front shape. What is most significant is the difference in computation time. Since the DSD/GSD model only needs to compute the propagation of the shock front, a 1-D curve, the simulation only takes approximately five minutes to compute. Since ALE3D solves the Euler conservation equations at every node at each time step, the simulation requires several hours of computational time. Therefore, the DSD/GSD model presents a massive increase in computational efficiency.
Figure 13 shows the relative error in time of arrival between the DSD/GSD hybrid simulation and the ALE3D simulation, using both angle boundary methods. With both methods, the largest error occurs inside the particle when the shock front becomes concave. The simulation using shock polar analysis has a maximum relative error of
$25\,\%$
and an average error of
$6.8\,\%$
. The method using artificially enforced boundary angles has a maximum relative error of
$8.6\,\%$
and an average error of
$1.9\,\%$
. As we would expect, the second method produces much more accurate data.
A comparison of a unit cell simulation performed with (a) the DSD/GSD hybrid model using the sonic state confinement conditions, (b) the DSD/GSD hybrid model with artificially enforced confinement angles and (c) the reactive flow model in ALE3D. Each wave represents a time difference of
$0.1~\unicode[STIX]{x03BC}\text{s}$
.
Contour plots of the relative error in time of arrival between the ALE3D unit cell simulation and the DSD/GSD hybrid simulation using (a) shock polar boundary angles and (b) artificially enforced boundary angles.
Finally, we run the simulation over 100 unit cells and capture the shock pressure at each time step using the method with artificially enforced boundary angles. We average these data over each unit cell to create a representative pressure profile, as shown in figure 14. The shock pressure is always higher in the particle than in the HE, as expected. The pressure in the HE is nearly constant, with a slight increase in pressure when the wave front is travelling around the particle. The particle undergoes the strongest shock loading in the bottom of the sphere, near the centre. By the time the detonation shock front reaches the top of the sphere, the pressure of the shock front has decayed by half.
The pressure profile of an aluminium sphere under shock loading, averaged over 100 unit cells.
8 Conclusion
In this paper we demonstrate the following. First, we show the trajectory of the transition of an expanding shock in atmospheric air from a blast wave state to an acoustic state. For our first case study, we use these data to simulate shock diffraction around a sphere at a wide range of Mach numbers, providing an extension to the classical Bryson and Gross theory.
Next, we derive and validate a theory that governs the propagation of shock fronts in inert materials that follow the Mie–Gruneisen equation of state, which includes most metals and plastics. The principles in this paper could be applied to any inert material that follows other equations of state, although more complicated materials might require a numerical solution. We present a guide to simulating hybrid explosives that consist of HE fluids and inert materials, although we restrict ourselves to relatively simple models. Our second and third case studies involve, respectively, a two-dimensional slab interface between an HE and a metal and an array of metal particles embedded in HE unit cells.
When making use of this theory, one must take care to estimate the range of shock velocities used in their simulation model. For most hybrid DSD/GSD models, the GSD theory alone is sufficient, but if one expects to see extremely high shock velocities in the simulation, it is necessary to develop a composite model that includes the GSD and TBW asymptotic limits. This must be done numerically for each individual material.
The method described in this paper to produce the most accurate unit cell simulation requires that one runs the same simulation using a DNS method first. This may seem to negate the point of using the hybrid method in the first place, however, one can run the DNS method over a single unit cell, record the boundary angles over time and then use those angles to run the hybrid method over hundreds of unit cells. Still, there is a great need for a theory of shock polar analysis that can accurately describe curved boundary interfaces.
The principles in this paper could inform the design of a more complicated simulation software that uses a nodal level set model to simulate mixed material explosives with arbitrary geometries, using the methods discussed in Aslam (Reference Aslam1996) and Hernández, Bdzil & Stewart (Reference Hernández, Bdzil and Stewart2013). This would require a separate level set propagation algorithm for the HE and the inert, but it would give considerable insight into studying the propagation of detonation shocks in hybrid explosives.
Acknowledgements
We would like to thank Dr J. Saenz, Dr S. Yoo, our coworkers, and the University of Illinois in Urbana-Champaign. This research was funded by Eglin Air Force Base (FA8651-10-1-0004; Advanced Modeling and Simulation Technologies for Micro-Munitions) and the Air Force Office of Scientific Research (FA9550-06-1-0044; Analysis of Multi-Scale Phenomena and Transients in Explosive and Complex Energetic Systems, FA9550-12-1-0422; Computational and Analytical Modeling of Advanced Energetic Materials).
Appendix A
A.1 NEWCODE hydrodynamic solver
The NEWCODE is our in-house hydrodynamic software that solves the three-dimensional reactive, compressible Euler equations on a uniform structured grid (Xu, Aslam & Stewart Reference Xu, Aslam and Stewart1997; Hernández & Stewart Reference Hernández and Stewart2013). The code is fully parallel and is based on a domain-decomposition model which uses the message-passing interface paradigm. The underlining numerical solver uses the method-of-lines to reduce the Euler equations to a system of ordinary differential equations, allowing us to solve the temporal and spatial problem independently. We use a higher-order weighted essential non-oscillatory (WENO) interpolation with fifth-order convergence and a positive preserving scheme to solve the spatial problem, and we solve the temporal problem via a third-order total variation diminishing (TVD) Runge–Kutta scheme. We use level sets to represent internal boundaries embedded in the flow fluid. We then apply reflective boundary conditions, with ghost nodes sorted by their connectivity to other internal boundary nodes, allowing the enforcement of the boundary conditions via a local and explicit approach. For a multi-component EOS mixture we assume that all components are uniformly well mixed, sharing one common pressure and temperature. We also assume mechanical equilibrium, that all components have the same particle velocity. Given a chemical reaction with
$N$
components and under the above assumptions we define the total mixture internal energy and saturation condition (mixture volume) as,
$$\begin{eqnarray}\displaystyle & \displaystyle e=\displaystyle \mathop{\sum }_{i=1}^{N}\unicode[STIX]{x1D706}_{i}e_{i}(p,v), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle v=\displaystyle \mathop{\sum }_{i=1}^{N}\unicode[STIX]{x1D706}_{i}v_{i}(p,T). & \displaystyle\end{eqnarray}$$
The solution of the above nonlinear equations involves an equilibration algorithm to solve for the mixture pressure and temperature. As part of the nonlinear solution process we define the following equation which we obtain by setting (A 1)–(A 2) equal to zero:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}F_{1}=e-\displaystyle \mathop{\sum }_{i=1}^{N}\unicode[STIX]{x1D706}_{i}e_{i}(p,v)=0,\\ F_{2}=v-\displaystyle \mathop{\sum }_{i=1}^{N}\unicode[STIX]{x1D706}_{i}v_{i}(p,T)=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
As part of the solution process we also need to invert
$v_{i}(p,T)$
which is unknown. Given a fixed pressure and temperature, we can solve
$v_{i}(p,T)$
using the EOS thermal and/or mechanical form.
$$\begin{eqnarray}\displaystyle f_{i}=\unicode[STIX]{x1D703}(v_{i})=0. & & \displaystyle\end{eqnarray}$$
We use a globally convergent Newton–Raphson method to equilibrate the pressure and temperature, and for the
$v_{i}(p,T)$
inversion, we use a hybrid Newton–Bisection solver.
A.2 NEWCODE numerical methods
To numerically solve the reactive Euler equations, we discretize the computational domain into a structured finite difference grid in which
$\text{d}x=\text{d}y=\text{d}z$
. Node points are labelled in the computational grid by
$i$
,
$j$
,
$k$
indices, which represent
$x$
,
$y$
and
$z$
coordinates along the physical domain. We use the method-of-lines approach to reduce the Euler equations to a system of ordinary differential equations, allowing us to solve the temporal and spatial problem independently, according to chapter 9.2 of LeVeque (Reference LeVeque2007). The Euler equations can then be rewritten as
$$\begin{eqnarray}\displaystyle \frac{\text{d}u}{\text{d}t}=L(u), & & \displaystyle\end{eqnarray}$$
where
$L(u)$
is the spatial operator on
$\boldsymbol{u}$
and is given as
$$\begin{eqnarray}\displaystyle L(u) & = & \displaystyle S_{i,j,k}-\left[\frac{1}{\text{d}x}(\hat{f}_{i+1/2,j,k}-\hat{f}_{i-1/2,j,k})+\frac{1}{\text{d}y}({\hat{g}}_{i,j+1/2,k}-{\hat{g}}_{i,j-1/2,k})\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\frac{1}{\text{d}z}({\hat{h}}_{i,j,k+1/2}-{\hat{h}}_{i,j,k-1/2})\right],\end{eqnarray}$$
where
$S$
is the source term,
$\hat{f}$
,
${\hat{g}}$
and
${\hat{h}}$
are the interface fluxes in the
$x$
,
$y$
and
$z$
coordinate directions respectively. Fluxes are ‘split’ into right running and left running waves by the following Lax–Friedrich flux splitting scheme:
$$\begin{eqnarray}\displaystyle f^{\pm }(u)={\textstyle \frac{1}{2}}(f(u)\pm \unicode[STIX]{x1D6FC}u), & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FC}$
is the dissipation coefficient defined as
$\unicode[STIX]{x1D6FC}=\max (f^{\prime }(u))$
. We use the local Lax–Friedrich formulation for the dissipation coefficient, which in the context of the Euler equations is defined as:
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}=\max (|u_{j}|+c,|u_{j+1}|+c), & & \displaystyle\end{eqnarray}$$
where
$c$
is the sound speed. We use higher-order WENO interpolation of Jiang & Shu (Reference Jiang and Shu1995) and the positive preserving scheme of Hu, Adams & Shu (Reference Hu, Adams and Shu2013) to approximate the interface fluxes. For our numerical simulations we use a third-order WENO scheme with fifth-order convergence. Generally, higher-order schemes like WENO are not guaranteed positive preserving which could lead to non-physical values for density and pressure. In these cases we use the positive preserving method presented in Hu et al. (Reference Hu, Adams and Shu2013). This scheme uses a flux limiter to detect critical numerical fluxes that may lead to negative values for density and pressure. It then limits these fluxes by combining the higher-order WENO flux with a first-order positive preserving Lax–Friedrich flux formulation.
Numerical integration of the system of ordinary differential equations (ODEs) is performed using a TVD third-order Runge–Kutta solver,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}u^{1}=u^{n}+\unicode[STIX]{x0394}tL(u^{n}),\\ u^{2}=\frac{3}{4}u^{n}+\frac{1}{4}u^{1}+\frac{1}{4}\unicode[STIX]{x0394}tL(u^{1}),\\ u^{n+1}=\frac{1}{3}u^{n}+\frac{2}{3}u^{2}+\frac{2}{3}\unicode[STIX]{x0394}tL(u^{2}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$
A.3 Intrinsic descriptions of shock dynamics
Description of shock dynamics requires a discussion of surface propagation and its intrinsic description in surface coordinates. One can find extended discussion of intrinsic shock attached coordinates in both detonation (Yao & Stewart Reference Yao and Stewart1996) and flame theory (Matalon, Cui & Bechtold Reference Matalon, Cui and Bechtold2003).
We can represent the shock front in general terms by a level set function
$\unicode[STIX]{x1D713}(x,y,z,t)=0$
. Figure 15 shows a sketch of a shock propagating into an unshocked region. We choose the normal vector to the shock surface to point towards the direction of propagation into unshocked material. The shock normal and the total shock curvature at that point on the surface is given by
$$\begin{eqnarray}\displaystyle \hat{n}=\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|}\quad \text{and}\quad \unicode[STIX]{x1D705}\equiv \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{n}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|}\right). & & \displaystyle\end{eqnarray}$$
The normal shock velocity and normal shock acceleration are given by
$$\begin{eqnarray}\displaystyle D_{n}=-\frac{\unicode[STIX]{x1D713}_{t}}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|}\quad \text{and}\quad {\dot{D}}_{n}=\frac{\unicode[STIX]{x2202}D_{n}}{\unicode[STIX]{x2202}t}+D_{n}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}D_{n}. & & \displaystyle\end{eqnarray}$$
We define the time of arrival field as
$T(x,y,z)$
, the time the shock first intersects the point
$(x,y,z)$
in space. If we assume the TOA field is uniquely defined, then we can define the shock location by the surface function
$\unicode[STIX]{x1D713}(x,y,z,t)=T(x,y,z)-t$
. For any specified time
$t_{0}$
, such that
$t_{0}=T(x,y,z)$
, the level set function
$\unicode[STIX]{x1D713}(x,y,z,t_{0})=0$
corresponds to the shock surface at that time. From this it follows that
$$\begin{eqnarray}\displaystyle D_{n}=\frac{1}{|\unicode[STIX]{x1D735}T|}\quad \text{and}\quad {\dot{D}}_{n}=\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}+D_{n}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\right)D_{n}. & & \displaystyle\end{eqnarray}$$
For a two-dimensional shock, it is useful to express the shock surface in terms of a preferred direction such as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(x,y,t)=y-g(x,t)=0, & & \displaystyle\end{eqnarray}$$
where
$g(x,t)$
represents the
$y$
coordinate of the shock front for a given
$x$
and
$t$
.
A 3-D shock surface. The outward shock normal points towards the unshocked region of the material. The shock location corresponds to the level set equation
$\unicode[STIX]{x1D713}(x,y,z,t)=0$
. The shock curvature is equal to
$1/R$
and the shock normal velocity is equal to
$\text{d}n/\text{d}t$
.
We can use (A 13) to express
$\unicode[STIX]{x1D705}$
,
$D_{n}$
and
${\dot{D}}_{n}$
in terms of
$g(x,t)$
(Lieberthal et al.
Reference Lieberthal, Bdzil and Stewart2014). Note that the notation
$g_{x}$
represents the partial derivative
$\unicode[STIX]{x2202}g/\unicode[STIX]{x2202}x$
, etc. It follows that
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D705}=-\frac{g_{xx}}{(1+g_{x}^{2})^{3/2}},\quad D_{n}=\frac{g_{t}}{(1+g_{x}^{2})^{1/2}},\quad \text{and}\quad {\dot{D}}_{n}=\frac{g_{tt}}{(1+g_{x}^{2})^{1/2}}-\frac{g_{x}g_{t}g_{xt}}{(1+g_{x}^{2})^{3/2}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
In this paper we discuss three basic forms of the intrinsic surface propagation law. In the case of inert shock dynamics, both the blast wave and acoustic limit surface propagation laws have the form
$$\begin{eqnarray}\displaystyle {\dot{D}}_{n}=-f(D_{n})\unicode[STIX]{x1D705}. & & \displaystyle\end{eqnarray}$$
We can illustrate the equation type by a local frozen coefficient analysis, where we assume that locally the shock is travelling at speed
$D_{0}$
in the
$y$
-direction. Then we write
$$\begin{eqnarray}\displaystyle g(x,t)=D_{0}t+\unicode[STIX]{x1D716}g_{1}(x,t), & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D716}$
is asymptotically small. Substitution of equation (A 16) into (A 14) then into equation (A 15) leads to a local evolution equation for the shock displacement
$g_{1}$
:
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}^{2}g_{1}}{\unicode[STIX]{x2202}t^{2}}-f(D_{0})\frac{\unicode[STIX]{x2202}^{2}g_{1}}{\unicode[STIX]{x2202}x^{2}}=0, & & \displaystyle\end{eqnarray}$$
which is a second-order linear wave equation with the wave speed for lateral disturbances
$\sqrt{f(D_{0})}$
.
Detonation shock dynamics, which we use for high explosive fluids, involves shock motion rules that we can use with the PB algorithm. The most common is the
$D_{n}$
–
$\unicode[STIX]{x1D705}$
relation that is written as
$\unicode[STIX]{x1D705}=F(D_{n})$
, with the property that when the shock is flat, or
$\unicode[STIX]{x1D705}=0$
,
$D_{n}=D_{CJ}$
for that explosive. Again, if we assume locally the shock is travelling at speed
$D_{0}$
in the
$y$
-direction, then using (A 14) and (A 16) we obtain the local evolution equation
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}g_{1}}{\unicode[STIX]{x2202}t}=-F^{\prime }(D_{0})\frac{\unicode[STIX]{x2202}^{2}g_{1}}{\unicode[STIX]{x2202}x^{2}}, & & \displaystyle\end{eqnarray}$$
which is a parabolic evolution equation, provided that
$F^{\prime }(D_{0})>0$
, which is a requirement for stable propagation.
