Verification of a Specialized Hydrodynamic Simulation Code for Modeling Deflagration and Detonation of High Explosives

Unsteady Lagrangian codes are some of the earliest hydrodynamic simulation codes developed for scientific computing [1]. Many modern Lagrangian multiphysics codes exist [2–4], which are able to model multidimensional problems with and many interacting physics processes. Such codes, though very powerful, are complex and difficult to extend to consider new physics models. Specialized codes, though more limited in the types of geometry and physics that they can consider, are generally simpler in their implementation, and it can be easier to prototype new physics models.

This paper considers a specialized hydrodynamic simulation code for use in the study of the deflagration and detonation of high explosives. In future work, this code will be used in the calibration of existing models for the chemical reactions occurring in high explosives and will also be used to examine new physics models for representing this physics process. This paper presents a series of verification exercises which ensure that the code has been implemented correctly to model the relevant physics processes for such experiments.

Modeling the combustion of high explosives requires a numerical scheme well suited to shocked flows. Additionally, problems involving deflagration require models for thermal conduction. Finally, reactive flows require models for the conversion of reactants into one, or more, species of combustion products and methods are required to obtain an equilibrium thermodynamic state within the reacting species.

Before such a code can be used in calibration and model development, a thorough verification exercise is required to ensure that the models have been implemented correctly. Each of the three disciplines, hydrodynamics, thermal conduction, and reactive flows, are considered in three different verification exercises.

Additionally, physics models for the equation of state (EOS) of pure species and mixtures are required as well as physics models for thermal conduction and chemical reaction rates. These models are also verified to ensure that their implementation reflects their intended behavior.

Section 2 briefly describes the underlying numerical scheme, highlighting some of its unique features relevant to the verification exercise. Section 3 describes the various physics models that are available in the code. Section 4 presents the verification exercises. Section 5 provides some conclusions.

The code uses an unsteady Lagrangian scheme formulated in one-dimensional pathline-attached coordinates so that no material flows across the cell interfaces. A typical computational cell is shown in Fig. 1.

where Eq. (6) includes additional terms for the conductive heat flux from adjacent cells, $q″=−k (∂T/∂x)$⁠.

where the n_c reacting species mass fractions, λ, are included in the vector of conserved states. Since the total mass in each cell is conserved in the Lagrangian coordinate system, only the mass fractions need to be evolved in time, not the total mass of each species in the cell.

where $u⋆$ is the velocity at the interface between two adjacent cells. The time-step dt is calculated based on hydrodynamic, thermal conduction, and reaction rate constraints, described in Appendix  A.

The conservation equations in integral form (4)–(6) can be converted into a differential form by performing the integration over the control volume and knowing that, since the mesh follows the pathline, no mass enters or leaves the cell along the edges, though pressure forces and temperature gradients allow momentum and energy to cross the boundary.

where $Ri,j$ is a vector n_c long, giving the rate of chance of each species given the thermodynamic state of the computational cell.

In these source terms, planar geometry corresponds to $ϕ=0$⁠, cylindrical $ϕ=1$⁠, and spherical $ϕ=2$⁠.

where the subscript _r is the cell centered quantity to the right of the interface, and _l is the cell centered quantity to the left of the interface.

The conserved variables can be evolved in time using a simple forward-Euler method (i.e., Eq. (5.1.1) in Ref. [6]) which is $O(dt)$ accurate. A second-order Runge–Kutta method can be used to evolve the conserved variables with $O(dt2)$ temporal accuracy, though this method requires some modifications, presented in Ref. [7], to account for the pathline-attached coordinate system. To achieve second-order spatial accuracy in the solution, a second-order method is also needed to extrapolate the cell centered quantities to the cell edges to solve the approximate Riemann problem (Sec. 2.1.1). Linear min-mod interpolation [8] is used with an additional check that the extrapolation did not change the sign of the internal energy, if it did, the cell centered value was used.

A “ghost” cell is added to the two extremes of the domain to enforce the boundary conditions on the problem. There are five possible hydrodynamic boundary conditions: “free” assumes that the domain continues infinitely in that direction, “vacuum” assumes that the boundary is adjacent to a vacuum, “wall” assumes that the boundary is a fixed wall, “symmetry” assumes that the boundary is an axis of symmetry, and “piston” assumes that a constant velocity piston is forcing the boundary. The state of the ghost cell for each boundary condition is given in Table 1.

There are four additional boundary conditions required for thermal conduction: “adiabatic” assumes that the boundary is perfectly insulated and there is no heat flux, “constant temperature” assumes that the boundary is maintained at a constant temperature, “constant flux” assumes that a constant heat flux passes through the boundary, and “constant convection” assumes that a flow with a constant temperature and convective heat transfer coefficient is impinging on that boundary. The state of the ghost cell for each thermodynamic boundary condition are given in Table 2.

In order to solve the equations given in Sec. 2, additional models are required to represent the physics processes at work in these problems.

By restricting the following discussion of EOS models for fixed-composition mixtures to only EOS models of this form, the analysis can be simplified compared to other presentations.

Seven different equation of state models have been implemented in the code, including a set of complete and thermodynamically consistent EOS models:

• ideal gas [10]

• Tait stiffened gas [11]

• Jones–Wilkins–Lee [12]

• Davis reactants [13]

• Davis products [13]

as well as two incomplete EOS models, generally used to model inert solids.

• linear U_s–U_p [14]

• Murnaghan [15]

These are EOS models typically used when simulating high explosive experiments. Some of these EOS models (Davis reactants, Davis products, ideal gas, and Tait stiffened gas) provide a relationship between temperature, density, and energy. Other EOS models (Murnaghan, Jones–Wilkins–Lee, and linear U_s–U_p) do not provide this relationship so an additional temperature model given by Eqs. (D.16) and (D.17) in Ref. [16] is used. In this model, it is assumed that the specific heat capacity at constant volume (C_v) and the Grüneisen gamma (Γ) are constants.

When material in a cell burns, the reactants are converted into one or more product species. The cell no longer has a homogeneous composition but is composed of a mixture of two or more species, each of which has its own EOS. Hydrodynamic calculations require a single set of thermodynamic properties to represent the state of the cell. Obtaining this state requires some assumptions about how the species in the cell mix. It is assumed that the pressure and temperature of all species are in equilibrium with the specific volume and internal energy of each species varying. Procedures to identify this equilibrium state have been implemented in many different codes, e.g., Appendix D of Ref. [17]. However, a simplified presentation of this procedure is given in Appendix  B which assumes that one reactant species is converted into one or more products and that all EOS models provide the functions described in Sec. 3.1. These assumptions simplify the presentation compared to other sources.

Additionally, given a mixture of two or more species, properties of the mixed cell, such as sound speed and specific heat capacity, are no longer trivially defined. Many simulation codes approximate these properties using a finite difference approach. In this code, these properties are calculated exactly using thermodynamic relations. The derivation of these exact relations for a mixture of two species is given in Sec. 4.6 of Ref. [9]. For the case of more than two species, the calculation is more complex and relies on analysis of the thermodynamic relationships for the mixture. This derivation of these relations for EOS models in the form described in Sec. 3.1 is not commonly described and is thus shown in Appendix  C.

The evolution of the mass fraction of the reacting species is given by an empirical model for the reaction. Four different chemical reaction models have been implemented.

• Arrhenius single- and multistep linear reactions [18]

• Wescott–Stewart–Davis (WSD) [19]

• Arrhenius–Wescott–Stewart–Davis (AWSD) [13]

• scaled unified reactive front (SURF) [20,21]

The Arrhenius-type rate laws provide the most physically correct model of the chemical reaction. The other rate laws are typically used to model detonation in high explosives and are similar to the Arrhenius rate law but have empirical adjustments that are meant to make them better represent the observed behavior of reacting condensed-phase explosives. To model reactive flow, an EOS for the reactants and each products species is required. The AWSD, WSD, and SURF models each have a single reactant and single product phase. The reactants are modeled using the Davis reactants EOS, and the products are modeled by the Davis products model for all calibrations shown in this paper.

Currently, the only models available are constant thermal conduction models, though more complex models that depend on density, temperature, and species concentration have been developed [22] and their implementation will be the focus of future work.

The code was verified using different exact solutions which exercised different sets of physics processes which are needed to model deflagration and detonation of high explosives.

where $hj+1<hj$⁠.

When using the second-order methods described in Sec. 2.1.2 and considering a domain unperturbed by shocks, the expected order of convergence is 2 for each of these problems as they use $second$-order accurate time integration and $second$-order spatial interpolation [6].

The Sedov [24] test problem examines the evolution of a shock generated by the instantaneous introduction of a fixed amount of energy at the origin of the domain. This problem has an exact solution defined in planar, cylindrical, and spherically symmetric coordinates. Solutions for this problem were obtained using the computer code exactpack [25]. This problem tests only the hydrodynamic equations and the approximate Riemann solver, and it does not exercise the thermal condition or chemical reaction models or algorithms.

The domain is $1.1 m$ wide and contains an ideal gas material with $γ=7/5$⁠, initially at rest with a density of $1 kg m−3$⁠. A point source of energy is located at the origin, and the amount of energy deposited is $0.0673185 J, 0.311357 J$⁠, and $0.851072 J$ for the planar, cylindrical, and spherical problems, respectively. The left boundary condition is symmetric, and the right boundary condition is free. The problem is initialized with the exact solution after it has evolved for $0.5 s$⁠, the numerical scheme then proceeds for a further $0.5 s$⁠, and the final state is compared to the exact solution.

Figure 2 shows order of convergence for the planar, cylindrical, and spherical problems. Since the solution was initialized using the density and energy of the exact solution, the error metric was the $L1$ norm of pressure. In all cases, as the spatial resolution is increased, the order of convergence approaches the theoretical rate of 2.0 for the second-order numerics. When the entire domain was used to compute the error, only a linear rate of convergence was observed. These results verify the implementation of the underlying numerical methods of Sec. 2.1, the approximate Riemann solver of Sec. 2.1.1, and the second-order methods presented in Sec. 2.1.2.

This verification problem tests unsteady conduction in a solid. The solid is initially at rest and at $300 K$ and is immersed in an environment at $1000 K$⁠, the domain is $5 μm$ wide and has a Biot number of 0.3, and the solution is evolved until the Fourier number is 0.4. The material is copper (⁠$k=401 W m−1 K$⁠, $ρo=8.933 g cm−3$⁠, and $α=117×10−6 m2 s−1$⁠). The left boundary is plane of symmetry and adiabatic, and the right boundary is a fixed wall and subject to a constant convective heat flux.

An exact solution to this problem is given by Incropera and DeWitt [26] (Secs. 5.5 and 5.6) for problems with planar, cylindrical, and spherical symmetries. The exact solution used the first five terms of the series solution, given by Schneider [27]. This tests the thermal conduction model as well as the thermal source term added to Eq. (11). Since the entire domain remains at rest, this problem only tests the thermal conduction models.

Figure 3 shows the order of convergence of the unsteady conduction problem using the second-order methods. The error metric is the dimensionless temperature θ, defined in Ref. [26]. The observed order of convergence approaches the theoretical rate of convergence of 2.0. This verifies the thermal conduction models as well as the additional source terms added to Eq. (11) to account for the cylindrical and spherical symmetries.

Powers [18] provides exact analytic solutions to the Zeldovich, von Neumann, Dohring (ZND) structure for a single- and two-step reaction of a hydrogen–air mixture with an Arrhenius-type rate law. The domain was initialized with the computed exact solution and then allowed to evolve into quiescent material with a constant velocity piston boundary condition on the right boundary. The solution is expected to converge at first order using the first-order numerics. In the region that did not pass through the shock wave, second-order convergence is expected [23] when second-order numerics are used. Though in theory, some energy is transmitted from the reacting region to the quiescent material through conduction, the energy deposited in the quiescent material through shock-heating is much greater than that transported via conduction, so conduction is ignored in this problem. The error metric for these simulations is the $L1$ norm of pressure.

The verification results for the single-step reaction are shown in Fig. 4. The expected rate of convergence for the first-order methods is 1 and for the second-order methods is 2. Both simulations converge at the expected rate verifying the implementation of reactive flow models and multispecies equilibration methods for the case of a single product species. When the second-order numerics were used and the entire domain downstream of the shock was used to assess the error, first-order convergence was observed, as shown in Fig. 5.

The verification results for the two-step reaction are shown in Fig. 6. These simulations had the same expected rates of convergence and, as with the single-step reaction, the solutions converge at close to the expected rates. This verifies the implementation of the multispecies equilibration methods for the more complicated case of a reaction with multiple species of reaction products.

Each EOS was tested at several points and, with the exception of the linear U_s–U_p EOS, each of the five thermodynamically complete EOS models was found to agree so that the relative error between the left and right-hand side of the equations was less than 1 × 10⁻⁷. Not all EOS models are thermodynamically consistent. Though the linear U_s–U_p is appropriate for pure hydrodynamic simulations, it should not be used for problems where the relationship between temperature and pressure is important, such as for reacting flows.

Using known mathematical relations such as Eqs. (21)–(25) to verify the implementation of physics models independent of the numerical methods is not typically part of verification studies. However, in the case of EOS models, the nine functions the EOS model provides are not independent but are all related. Subtle errors in implementation can distort the desired relationships between the functions. This verification exercise was important in both developing the code and ensuring that these models can be used in confidence in combustion problems where the interaction between all the thermodynamic properties of the material is of central importance. In practice, performing these tests revealed several extremely subtle errors in the model implementation.

The procedures to obtain multispecies pressure–temperature equilibrium (Appendix  B) and the properties derived from the equilibrium concentrations of species (Appendix  C) were complex and required some form of verification to ensure that the thermodynamic theory was correct.

The method presented in Appendix  C to calculate specific heats can be tested by comparing the exact specific heat for a mixture calculated via finite differences. The mixture had a fixed composition of two fluids with the properties of HE reactants material and products, modeled using the EOS given in Ref. [13]. The material is at $ρ=2.9343 g cm−3$⁠, $e=7.9684 kJ g−1, and λ=0.54553$⁠. The error should converge at first order, and this is the behavior seen in Table 3. The positive convergence rate shows that the exact analytic models for specific heat capacity are in agreement with a finite difference approximation in the limit of small changes in temperature.

The piston was moving at $0.1 m s−1$⁠, and the fixed-composition mixture was the same as in Sec. 4.4.2. The pressure profile across the domain is shown in Fig. 7 after 0.622 μs of simulated time. The difference between the exact change in pressure and the observed pressure difference was 0.0302 kPa (0.001%). The agreement between the observed pressure change across the shock and Eq. (27) shows that the sound speed for mixtures is calculated correctly.

The verification problems examined in Sec. 4.3 used Arrhenius-type rate laws and ideal gas equations of state. The problems in this section examine the use of complex rate laws typically used to model detonation problems in HE as well as realistic EOS models. When using ideal gas EOS models, the solution to the differential-algebraic equation governing the ZND structure had an analytic solution for each integration time-step. When using realistic EOS models, root finding procedures are required to find the solution to the algebraic problems for each time-step in the differential equation, as described by Menikoff [29]. Direct integration of the underlying equations of the ZND wave provided high accuracy solutions of the ZND profile. Similar to the verification problems in Sec. 4.3, the problem was initialized with the exact solution and allowed to evolve in time. The pressure profile at the end of the solution was compared to the exact solution. Figure 8 shows the convergence of the ZND problem for the AWSD, WSD, and SURF rate laws described in Sec. 3.3 using calibrations for the plastic bonded explosives (PBX) PBX 9501 and PBX 9502. Since the numerical methods governing the evolution of the reaction progressed variables were verified in Sec. 4.3, the metric for success of these verification tests of the reaction rate model's implementation was only to observe a positive convergence rate. All cases converged at a rate close to one, verifying the implementation of these complex reaction rate models in the code.

Numerical methods have been presented and verified for a specialized Lagrangian hydrodynamics code intended to model the deflagration and detonation of condensed-phase high explosives. Several order-of-accuracy verification tests exercised the various physics capabilities that were required in this code. The expected order of convergence was observed in these tests, demonstrating that the numerical methods and models have been correctly implemented in this code. Prior to this verification exercise, the code had been subject to traditional software unit-testing of its various functions. Though this sort of testing is necessary for quality code, these tests are not sufficient to ensure that implementation of the numerical scheme is correct. The process of performing these verification exercises identified errors in the implementation of the numerical methods that unit-testing did not catch and, most likely, could not have been identified any other way. This shows the importance of verification studies in the development of scientific software, alongside traditional software testing methods.

Additionally, choosing a suite of verification problems that tested subsets of the physics models required by the code made it easier to identify where the implementation errors were in the code. An additional feature of the tests shown in this paper is that they have been automated. Every figure and table shown in this paper is regenerated whenever a change is made to the source code. This allows the quality of the software to be tested continuously during development as well as before performing critical calculations. This ensures that the code performs well not only at the snapshot in time when the verification exercise was performed but continues to show the desired rates of convergence as the code is developed.

The authors would like to thank the development teams of the many open source software packages used in the preparation of this paper, including the developers of scipy [31], matplotlib [32], numpy [33], and cython [34]. Additionally, the authors would like to thank their colleagues at Los Alamos National Laboratory for many informative discussions on this topic including Ann E. Wills and C. N. Woods. This work was executed by Los Alamos National Laboratory under Contract No. 89233218CNA000001. Approved for public release LA-UR-20-30207.

• U.S. Department of Energy's Advanced Simulation and Computing program and the Conventional High Explosives Grand Challenge (Funder ID: 10.13039/100000015).

• Los Alamos National Laboratory (Contract No. 89233218CNA000001; Funder ID: 10.13039/100008902).

Variables

c =

isentropic sound speed (m s⁻¹)

c_t =

isothermal sound speed (m s⁻¹)

C_v =

specific heat at constant volume (J kg⁻¹ K⁻¹)

e =

mass-specific internal energy (J kg⁻¹)

F =

vector of conserved quantities

CFL =

Courant–Friedrichs–Lewy number (dimensionless)

h =

characteristic grid size (m)

h_BC =

convective heat transfer coefficient on boundary (W m⁻² K⁻¹)

k =

thermal conductivity (W m⁻¹ K⁻¹)

$L1$ =

the absolute value norm of a quantity (dimensionless)

n_c =

the number of chemically reacting species (dimensionless)

p =

rate of convergence (dimensionless)

P =

pressure (Pa)

$q″$ =

heat flux (W m⁻²)

r =

radius (m)

R =

chemical reaction rate (s⁻¹)

t =

temporal coordinate (s)

T =

temperature (K)

u =

velocity (m s⁻¹)

u_piston =

velocity of piston boundary condition (m s⁻¹)

V =

specific volume (⁠$V≡ρ−1$⁠) (m³ kg⁻¹)

x =

spatial coordinate (m)

Symbols

α =

thermal diffusivity (m s⁻²)

Γ =

Grüneisen gamma (dimensionless)

θ =

dimensionless temperature $(T−Tinf)/(To−Tinf)$ (-)

λ =

species mass fraction (dimensionless)

μ =

Mach angle (rad)

ρ =

density (kg m⁻³)

$ϕ$ =

symmetry term (dimensionless)

Modifiers

chem =

property of chemical equilibrium models

cond =

property of thermal conduction models

crit =

time-step limiting value

ghost =

property of a ghost cell

hyrdo =

property of hydrodynamic models

_i =

spatial cell index

_j =

temporal cell index

^⋆ =

property at an interface

The time-step size controlled the stability of the temporal evolution of the solution. Hydrodynamics, thermal conduction, and reaction rates all affected the stability of the solution. The largest stable time-step for each of these limiting factors was determined, and then the smallest of these critical time-steps was used for the temporal integration described in Sec. 2.1.2. This method ensures that the critical time-step for each physical process is not exceeded. In some cases, where multiple physical processes are controlling the time-step, a CFL number closer to $1/2$ may be required for numerical stability.

As reactions proceed in the material and reactants are converted into $nprod≥1$ product species, the thermodynamic state of the cell is no longer defined by the properties of a single species but from a fixed-composition mixture of multiple species. For materials defined by complex equations of state, there is no analytic solution to the thermodynamic state of this mixture; it must be solved for in an iterative manner. This is a common problem in many hydrodynamic codes, the presentation in this section simplifies the derivation by assuming that the reaction has one reactant and one or more products and that all materials are defined by equations of state which provide the functions and derivatives listed in Sec. 3.1.

The set of Eqs. (B1)–(B4) give a set of $2nprod+2$ equations for the $2nprod+2$ unknowns $x={ρreact,ereact,{ρprodn,eprodn∀n∈nprod}$⁠. This set of equations can be solved using Newton–Raphson iteration [6].

Let $εi=0−f(xi)$⁠, where f returns a vector $2nprod+2$ long, evaluating the right-hand side of Eqs. (B1)–(B4).

this process is iterated until the Euclidean norm of the error vector is less than an absolute tolerance, in this case $1×10−6$⁠.

and the equations which make up the error vector ε are in the same order as Eqs. (B1)–(B4).

An initial vector, x₀, is required to start the Newton–Raphson iteration. This vector is obtained from a guessed pressure for the mixture and a guessed density for each reacting species so that $er0=e(ρ̂r,P̂)$⁠, $epn0=e(ρ̂pn,P̂)∀n∈nprod$ and $ρr0=ρ̂r$⁠, and $ρpn0=ρ̂pn∀n∈nprod$⁠. For the first iteration of the hydrocode, the guessed pressure is the initial pressure, and the guessed densities for reactants and products are each the mixture density. The pressure and species densities of the equilibrium mixtures are stored and used as the initial guess for the next iteration of the hydrocode.

Once an equilibrium mixture of reactants and products is obtained from the pressure–temperature equilibration routine, additional properties such as the specific heat capacity and sound speed of the mixture are needed. For two component mixtures, Menikoff [9] provides simple formulas for these properties. For mixtures of more than two components, more laborious investigations of the thermodynamic derivatives of the mixture are required.

where n_spec is the total number of reacting species including reactants and products. (Note that specific volume is the reciprocal of density so $∂V/∂ρ=−1/ρ2$⁠.) Evaluation this expression requires knowledge of the other thermodynamic derivatives for the mixture. Each of these derivatives can be constructed from the properties of the constituent species as follows.