Abstract

A specialized hydrodynamic simulation code has been developed and verified for the simulation of one-dimensional unsteady problems involving the detonation and deflagration of high explosives. To model all the relevant physical processes in these problems, a code is required to simulate compressible hydrodynamics, unsteady thermal conduction, and chemical reactions with complex rate laws. Several verification exercises are presented which test the implementation of these capabilities. The code also requires models for physics processes such as equations of state and conductivity for pure materials and mixtures as well as rate laws for chemical reactions. Additional verification tests are required to ensure that these models are implemented correctly. Though this code is limited in the types of problems it can simulate, its computationally efficient formulation allows it to be used in calibration studies for reactive burn models for high explosives. This study demonstrates how a series of verification tests can be used to ensure that the various physics processes needed to simulate complex phenomenon can be tested to ensure that they are correctly implemented.

1 Introduction

Unsteady Lagrangian codes are some of the earliest hydrodynamic simulation codes developed for scientific computing [1]. Many modern Lagrangian multiphysics codes exist [24], 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.

2 Numerical Formulation

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.

Fig. 1
Schematic of the finite volume region. Also indicated are the “initial” and “final” t levels, tj and tj+1, respectively, as well as the pathline bounds as indicated by the various x locations. The conservative vector state, F, at both t locations is denoted. The deflection of the pathlines in x–t space is the velocity of the interface, u.
Fig. 1
Schematic of the finite volume region. Also indicated are the “initial” and “final” t levels, tj and tj+1, respectively, as well as the pathline bounds as indicated by the various x locations. The conservative vector state, F, at both t locations is denoted. The deflection of the pathlines in x–t space is the velocity of the interface, u.
Close modal
The three conservation laws for unsteady compressible flows in one-dimensional Eulerian coordinates are
ρt+x(ρu)=0
(1)
t(ρu)+x(ρu2+P)=0
(2)
t(ρ(e+12u2))+x(u(ρ(e+12u2)+P))x(kTx)=0
(3)
In streamline attached Lagrangian coordinates, the velocity normal to the cell edges is zero. Using this simplification, the conservation laws can be transformed into discretized equations for conservation of mass
txi12xi+12ρidx=0
(4)
conservation of x-momentum
txi12xi+12ρiuidx=Pi12Pi+12
(5)
and conservation of energy
txi12xi+12ρi(ei+12ui2)dx=Pi12ui12+qi12Pi+12ui+12qi+12
(6)

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

The vector of conserved variables F is
F=Δx(ρ,ρu,ρ(e+12u2)),λ1,,λnc
(7)

where the nc 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.

In this pathline following scheme, the mesh evolves to follow the flow at the interface between the cells. This method of mesh construction ensures the velocity normal to the cell edges is zero, and therefore no mass enters or leaves the cell through these interfaces
xi12,j+1=xi12,j+dtui12,j
(8)

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.

2.1 Conservation Equations for Hydrodynamics.

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.

Integrating Eq. (4) provides an equation for the evolution in time of the first conserved variable
F1i,j+1=F1i,jϕri,jdt(xi+12,jxi12,j)ρi,jui,j
(9)
integrating Eq. (5) yields
F2i,j+1=F2i,jdt(Pi+12Pi12)ϕri,jdt(xi+12,jxi12,j)ui,j2ρi,j
(10)
integrating Eq. (6) yields
F3i,j+1=F3i,jdt(Pi+12,jui+12,jqi+12,j′′Pi12,jui12,j+qi12,j′′)ϕri,jdt(xi+12,jxi12,j)(ui,j(ρi,j(ei,j+12ui,j2)+Pi,j)12(qi+12,j′′+qi12,j′′))
(11)
and finally, the mass fractions are evolved via
F(43+nc)i,j+1=F(43+nc)i,j+dtRi,j
(12)

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

To model one-dimensional problems which have cylindrical and spherical symmetries, Eqs. (9)(11) contain additional source terms to account for the change in the conserved quantity as the radius increases. These assume that the cell represents a finite volume revolved about x = 0
ri,j=12(xi12,j+xi+12,j)
(13)

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

2.1.1 Approximate Riemann Solver.

The evolution of the conserved variables from one time-step to the next depends on the current state of the cells and the pressure, velocity, and conductive heat flux at the interface between two cells. The pressure and velocity at the interface are obtained by an approximate Riemann solver, given by Eq. (9.20) in Ref. [5]
u=ρrcrur+ρlclul+(PlPr)ρrcr+ρlcl
(14)
P=ρrcrPl+ρlclPr+ρrcrρlcl(ulur)ρrcr+ρlcl
(15)

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 heat flux at the interface is calculated using Fourier's law
q=ki(Ti1Ti+1)12(xi32xi+32)
(16)

2.1.2 Second-Order Methods.

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.

2.1.3 Hydrodynamic Boundary Conditions.

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.

Table 1

State of the ghost cell for hydrodynamic boundary conditions

BoundaryPghostughost(ρc)ghost
FreeP+u+(ρc)+
Vacuum0u+(ρc)+
Wall, symmetryP+u+(ρc)+
PistonP+2upistonu+(ρc)+
BoundaryPghostughost(ρc)ghost
FreeP+u+(ρc)+
Vacuum0u+(ρc)+
Wall, symmetryP+u+(ρc)+
PistonP+2upistonu+(ρc)+

Note that the plus (+) represents the properties of the penultimate cell.

2.1.4 Thermal Conduction Boundary Conditions.

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.

Table 2

State of the ghost cell for thermal boundary conditions

BoundaryTghostkghost
AdiabaticT+k+
Constant temperatureTBCk+
Constant fluxT++qBCΔxk+k+
Constant convectionT++hBC(T+TinfBC)Δxk+k+
BoundaryTghostkghost
AdiabaticT+k+
Constant temperatureTBCk+
Constant fluxT++qBCΔxk+k+
Constant convectionT++hBC(T+TinfBC)Δxk+k+

Note that the plus (+) represents the properties of the penultimate cell and Δx=(xi(3/2)xi+(3/2)).

3 Physics Models

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

3.1 Equation of State.

An EOS provides a set of functions which relate the thermodynamic properties of a pure material or mixture of materials. The natural thermodynamic variables for constructing an EOS are specific volume and temperature. However, in Lagrangian hydrodynamics, it is convenient to formulate EOS models as functions of density and mass-specific internal energy. Though some equations of state can be easily inverted to provide explicit functions in any set of two thermodynamic variables, the Mie–Grüneisen EOS models commonly used in high explosive (HE) research are typically transcendental functions [9], which cannot be inverted explicitly. These models can only provide the following nine functions:
P=fP(ρ,e),Pρ=fdPdρ(ρ,e),Pe=fdPde(ρ,e)T=fT(ρ,e),Tρ=fdTdρ(ρ,e),Te=fdTde(ρ,e)Cv=fCv(ρ,e)c2=fc2(ρ,e)Γ=fΓ(ρ)

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 UsUp [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 UsUp) 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 (Cv) and the Grüneisen gamma (Γ) are constants.

3.2 Multispecies Pressure–Temperature Equilibrium.

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.

3.3 Chemical Reaction Rate.

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.

3.4 Thermal Conduction.

Converting the temperature distribution in the domain into thermal fluxes at the cell interface requires a model for the thermal conduction in the material
k=fk(ρ,T,λ)
(17)

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.

4 Verification

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.

The error metric was the L1 norm of a property integrated across a domain. This was a subset of the entire computational domain, chosen to exclude regions which had been perturbed by a shock, as such regions could only experience, at best, linear convergence, regardless of the expected order of convergence of the method [23]. Additionally, the presence of discontinuities and boundaries in the solution created disturbances which traveled through the domain at the speed of sound and prevented convergence at rates greater than one, regardless of the numerical scheme. The simulation time and the subset of the domain used to assess the rate of convergence were chosen to exclude regions affected by these disturbances
Lj1=i=1ndxi|uiuexact(xi)|i=1ndxi
(18)
The characteristic grid size of the simulation was
hj=i=1ndxin
(19)
and the observed order of convergence was calculated as
p=ln(Lj+11Lj1)ln(hj+1hj)
(20)

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].

4.1 Hydrodynamic Verification.

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.1m wide and contains an ideal gas material with γ=7/5, initially at rest with a density of 1kgm3. A point source of energy is located at the origin, and the amount of energy deposited is 0.0673185J,0.311357J, and 0.851072J 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.5s, the numerical scheme then proceeds for a further 0.5s, 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.

Fig. 2
Verification of Sedov problem with planar, cylindrical, and spherical symmetries, second-order numerics. The solution at each resolution is shown as solid lines, and the error as dashed colored lines. The dashed-dotted black line shows the exact solution used to initialize the problem, and the dashed black line shows the exact solution used to compute the error. (a) Planar, (b) planar, (c) cylindrical, (d) cylindrical, (e) spherical, and (f) spherical.
Fig. 2
Verification of Sedov problem with planar, cylindrical, and spherical symmetries, second-order numerics. The solution at each resolution is shown as solid lines, and the error as dashed colored lines. The dashed-dotted black line shows the exact solution used to initialize the problem, and the dashed black line shows the exact solution used to compute the error. (a) Planar, (b) planar, (c) cylindrical, (d) cylindrical, (e) spherical, and (f) spherical.
Close modal

4.2 Thermal Conduction Verification.

This verification problem tests unsteady conduction in a solid. The solid is initially at rest and at 300K and is immersed in an environment at 1000K, 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=401Wm1K, ρo=8.933gcm3, and α=117×106m2s1). 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.

Fig. 3
Verification of thermal conduction with planar, cylindrical, and spherical symmetries, second-order numerics. The solution at each resolution is shown as solid lines, and the error as dashed colored lines. The dashed black line is the exact solution at the final time. (a) Planar, (b) planar, (c) cylindrical, (d) cylindrical, (e) spherical, and (f) spherical.
Fig. 3
Verification of thermal conduction with planar, cylindrical, and spherical symmetries, second-order numerics. The solution at each resolution is shown as solid lines, and the error as dashed colored lines. The dashed black line is the exact solution at the final time. (a) Planar, (b) planar, (c) cylindrical, (d) cylindrical, (e) spherical, and (f) spherical.
Close modal

4.3 Reactive Flow Verification.

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.

Fig. 4
Convergence for the single-step reaction after 1.5 μs of simulated time. Solid colored lines are the solution at each resolution, and dashed colored lines are the error. The dotted line is the exact solution at the final time. The vertical dotted line is the position of the shock at t = 0, and the vertical dashed lines show the region where the error is computed. (a) First order, (b) first order, (c) second order, and (d) second order.
Fig. 4
Convergence for the single-step reaction after 1.5 μs of simulated time. Solid colored lines are the solution at each resolution, and dashed colored lines are the error. The dotted line is the exact solution at the final time. The vertical dotted line is the position of the shock at t = 0, and the vertical dashed lines show the region where the error is computed. (a) First order, (b) first order, (c) second order, and (d) second order.
Close modal
Fig. 5
Convergence for the two-step reaction after 1.25 ns of simulated time. Solid colored lines are the solution at each resolution, and dashed colored lines are the error. The dotted line is the exact solution at the final time. The vertical dotted line is the position of the shock at t = 0. The error is computed across the entire domain. (a) Second order and (b) second order.
Fig. 5
Convergence for the two-step reaction after 1.25 ns of simulated time. Solid colored lines are the solution at each resolution, and dashed colored lines are the error. The dotted line is the exact solution at the final time. The vertical dotted line is the position of the shock at t = 0. The error is computed across the entire domain. (a) Second order and (b) second order.
Close modal

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.

Fig. 6
Convergence for the two-step reaction after 1.25 ns of simulated time. Solid colored lines are the solution at each resolution, and dashed colored lines are the error. The dotted line is the exact solution at the final time. The vertical dotted line is the position of the shock at t = 0, and the vertical dashed lines show the region where the error is computed. (a) First order, (b) first order, (c) second order, and (d) second order.
Fig. 6
Convergence for the two-step reaction after 1.25 ns of simulated time. Solid colored lines are the solution at each resolution, and dashed colored lines are the error. The dotted line is the exact solution at the final time. The vertical dotted line is the position of the shock at t = 0, and the vertical dashed lines show the region where the error is computed. (a) First order, (b) first order, (c) second order, and (d) second order.
Close modal

4.4 Model Verification

4.4.1 Equation of State Models for Pure Species.

All EOS models define a function relating pressure to density and energy. For almost all models, the derivatives of this function are not given in the published sources and must be derived; often requiring the assistance of a computer-aided algebra system, such as maxima [28]. Verification tests ensure that the integral of these derivative functions agrees with the function itself. The following integrals were performed, and the left-hand side was observed to agree with the right-hand side within a relative tolerance of 1 × 10−7:
ρ1ρ2dρPρ(ρ,e)=P(ρ2,e)P(ρ1,e)
(21)
e1e2dePe(ρ,e)=P(ρ,e2)P(ρ,e1)
(22)
ρ1ρ2dρTρ(ρ,e)=T(ρ2,e)T(ρ1,e)
(23)
e1e2deTe(ρ,e)=T(ρ,e2)T(ρ,e1)
(24)
Additionally, the relationship between the functions for temperature and for pressure is tested using thermodynamic relations. For a thermodynamically consistent EOS, the following relationship should be true:
1T(PE|VPTTE|V)=1T2TV|E
(25)

Each EOS was tested at several points and, with the exception of the linear UsUp 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−7. Not all EOS models are thermodynamically consistent. Though the linear UsUp 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.

4.4.2 Fixed Composition Mixture Specific Heat.

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.

Specific heat at constant volume is defined as
CveT|Ve2e1T(ρ,e2,λ)T(ρ,e1,λ)
(26)

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.9343gcm3, e=7.9684kJg1,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.

Table 3

Convergence of specific heat to the analytically computed value of 1.6493 × 10−3 kJ g−1 K−1

ΔeRelative errorConvergence rate
J kg−1
3.20 × 10−28.89 × 10−3
1.60 × 10−24.49 × 10−30.99
8.00 × 10−32.26 × 10−30.99
4.00 × 10−31.13 × 10−31.00
ΔeRelative errorConvergence rate
J kg−1
3.20 × 10−28.89 × 10−3
1.60 × 10−24.49 × 10−30.99
8.00 × 10−32.26 × 10−30.99
4.00 × 10−31.13 × 10−31.00

4.4.3 Fixed-Composition Mixture Sound Speed.

The sound speed calculation defined in Appendix  C was verified by comparing the pressure change in a domain of uniform composition perturbed by a piston that sent a sonic wave through the domain. The pressure change across a shock is obtained from the Hugoniot jump conditions (see Chap. 16 in Ref. [14]). In the limit of a weak disturbance, the shock speed approaches the sound speed, and the pressure change across this wave is given by
ΔP=ρcupiston
(27)

The piston was moving at 0.1ms1, 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.

Fig. 7
Pressure profile in a fixed-composition mixture disturbed by a piston
Fig. 7
Pressure profile in a fixed-composition mixture disturbed by a piston
Close modal

4.4.4 Complex Reaction Rates.

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.

Fig. 8
Convergence of the ZND solution for realistic EOS models and rate laws. (a) PBX 9502 using AWSD rate law [13]. t = 0.1 μs. (b) PBX 9502 using WSD rate law [19]. t = 0.1 μs. (c) PBX 9501 using SURF rate law [30]. t = 25 ns.
Fig. 8
Convergence of the ZND solution for realistic EOS models and rate laws. (a) PBX 9502 using AWSD rate law [13]. t = 0.1 μs. (b) PBX 9502 using WSD rate law [19]. t = 0.1 μs. (c) PBX 9501 using SURF rate law [30]. t = 25 ns.
Close modal

5 Conclusions

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.

Acknowledgment

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.

Funding Data

  • 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).

Nomenclature

Variables
c =

isentropic sound speed (m s−1)

ct =

isothermal sound speed (m s−1)

Cv =

specific heat at constant volume (J kg−1 K−1)

e =

mass-specific internal energy (J kg−1)

F =

vector of conserved quantities

CFL =

Courant–Friedrichs–Lewy number (dimensionless)

h =

characteristic grid size (m)

hBC =

convective heat transfer coefficient on boundary (W m−2 K−1)

k =

thermal conductivity (W m−1 K−1)

L1 =

the absolute value norm of a quantity (dimensionless)

nc =

the number of chemically reacting species (dimensionless)

p =

rate of convergence (dimensionless)

P =

pressure (Pa)

q =

heat flux (W m−2)

r =

radius (m)

R =

chemical reaction rate (s−1)

t =

temporal coordinate (s)

T =

temperature (K)

u =

velocity (m s−1)

upiston =

velocity of piston boundary condition (m s−1)

V =

specific volume (Vρ1) (m3 kg−1)

x =

spatial coordinate (m)

Symbols
α =

thermal diffusivity (m s−2)

Γ =

Grüneisen gamma (dimensionless)

θ =

dimensionless temperature (TTinf)/(ToTinf) (-)

λ =

species mass fraction (dimensionless)

μ =

Mach angle (rad)

ρ =

density (kg m−3)

ϕ =

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

Appendix A: Time-Step Size

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.

A.1 Hydrodynamic Time-Step.
The hydrodynamic time-step was calculated so that an acoustic wave emanating from a pathline would not cross either of the adjacent pathlines before the end of the time-step. A schematic of two adjacent computational cells with an acoustic wave emanating from their shared interface is shown in Fig. 9 
Δtcrit1=CFLdx1c1+u1uc
(A1)
Δtcrit2=CFLdx2c2u2+uc
(A2)
Fig. 9
Diagram of two cells showing how the critical step in Δt was calculated. μ is the Mach angle, μ=arctan(c).
Fig. 9
Diagram of two cells showing how the critical step in Δt was calculated. μ is the Mach angle, μ=arctan(c).
Close modal
A.2 Thermal Time-Step.
The critical time-step for thermal conduction was given by Sec. 5.9 in Ref. [26]
Δttherm=dx2CFLtherm2α
(A3)
A.3 Chemical Reaction Time-Step.
In order to prevent very fast reactions from completing in less than one time-step, the reaction rate was used to limit the time-step size
Δtchem=CFLchemmax(R)
(A4)

Appendix B: Multispecies Pressure-Temperature (P–T) Equilibrium

As reactions proceed in the material and reactants are converted into nprod1 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.

All the reacting species were in pressure equilibrium
0=Preact(ρreact,ereact)Pprod1(ρprod1,eprod1)0=0=Preact(ρreact,ereact)Pprodnprod(ρprodnprod,eprodnprod)
(B1)
as well as temperature equilibrium
0=Treact(ρreact,ereact)Tprod1(ρprod1,eprod1)0=0=Treact(ρreact,ereact)Tprodnprod(ρprodnprod,eprodnprod)
(B2)
Additionally, the total specific volume of the mixture must be completely divided between the constituents
0=1ρλrρreactn=1nprodλpnρprodn
(B3)
Similarity, the specific internal energy of the mixture must be completely divided between the constituents
0=eλrereactn=1nprodλpneprodn
(B4)

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

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

A new value of x is obtained which minimizes the error
xi+1=xi+εixiJ·εi
(B5)

this process is iterated until the Euclidean norm of the error vector is less than an absolute tolerance, in this case 1×106.

The Jacobian matrix J can be evaluated analytically given an analytic EOS model which provides the functions specified in Sec. 3.1
J=[fρrferfρ1fe100fρrfer00fρnprodfenprodgρrgergρ1ge100gρrger00gρnprodgenprodλrρr20λp1ρp120λpnρpn200λr0λp10λpn]
(B6)
where
fρr=Preact(ρreact,ereact)ρreactgρr=Treact(ρreact,ereact)ρreactfer=Preact(ρreact,ereact)ereactger=Treact(ρreact,ereact)intereactfρn=Pprodn(ρprodn,eprodn)ρprodngρn=Tprodn(ρprodn,eprodn)ρprodnfen=Pprodn(ρprodn,eprodn)eprodngen=Tprodn(ρprodn,eprodn)eprodn
Recall that J=ε/x, where the vector x is
x=[vrerv1e1vnprodenprod]
(B7)

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

B.1 Initial Guess.

An initial vector, x0, 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̂)nnprod and ρr0=ρ̂r, and ρpn0=ρ̂pnnnprod. 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.

Appendix C: Multispecies Mixture Properties

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.

C.1 Specific Heat Capacity.
The specific heat at constant volume (Cv) for the mixture is calculated using Eq. (3.2) of Ref. [35]
Cvmix=PT|VmixeV|TmixPV|Tmix+n=1nspecλn(CvnPT|VnenVn|TPVn|T)
(C1)

where nspec 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.

Using Eq. (3.26) of Ref. [35], the mixture isothermal sound speed is
PV|Tmix=(n=1nspecλn1PVn|T)1
(C2)
Using Eq. (3.30) of Ref. [35]
PT|Vmix=PV|Tmixn=1nspecλnPT|VnPVn|T
(C3)
and using Eq. (3.31) of Ref. [35]
eV|Tmix=PV|Tmixn=1nspecλnenVn|TPVn|T
(C4)
The thermodynamic derivatives of the reacting species are calculated as
PVn|T=ctnρn2
(C5)
where the isothermal sound speed can be related to known properties of the EOS using Eq. (4.9) of Ref. [9]
ctn2=cn2Γn2CvnT
(C6)
Additionally, using Eq. (3.2) of Ref. [35]
enVn|T=TρnΓnCvnP
(C7)
and using Eq. (3.8) of Ref. [35]
PT|Vn=ρnΓnCvn
(C8)
C.2 Sound Speed.
Knowing the specific heat of the mixture and other important thermodynamic derivatives, the isothermal sound speed of the mixture can be calculated using Eq. (3.38) of Ref. [35]
(ρmixcmix)2=PV|TmixCvmixn=1nspecλncn2CvnVnPVn|T
(C9)

References

1.
Morgan
,
N.
, and
Archer
,
B. J.
,
2021
, “
On the Origins of Lagrangian Hydrodynamic Methods
,”
ANS J. Nucl. Technol.
,
207
(
Suppl. 1
), pp.
S147
S175
.10.1080/00295450.2021.1913034
2.
Burton
,
D. E.
,
1994
, “
Multidimensional Discretization of Conservation Laws for Unstructured Polyhedral Grids
,”
Second International Workshop on Analytical Methods and Process Optimization in Fluid and Gas Mechanics (SAMGOP-94)
, Holiday Base, Arzamas, Russian Federation.
3.
Buren
,
K. L. V.
,
Canfield
,
J. M.
,
Hemez
,
F. M.
, and
Sauer
,
J. A.
,
2012
, “
Code Verification of the HIGRAD Computational Fluid Dynamics Solver
,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-UR
12
21165
.
4.
Noble
,
C. R.
,
Anderson
,
A. T.
,
Barton
,
N. R.
,
Bramwell
,
J. A.
,
Capps
,
A.
,
Chang
,
M. H.
,
Chou
,
J. J.
,
2017
, “
Ale3d: An Arbitrary Lagrangian-Eulerian Multi-Physics Code
,” Lawrence Livermore National Laboratory, Livermore, CA, Report No. LLNL-TR 732040.
5.
Toro
,
E. F.
,
1993
,
Riemann Solvers and Numerical Methods for Fluid Dynamics
,
Springer
, Berlin.
6.
Press
,
W. H.
,
Flannery
,
B. P.
,
Teukolsky
,
S. A.
, and
Vetterling
,
W. T.
,
1989
,
Numerical Recipes: The Art of Scientific Computing (FORTRAN)
,
Cambridge University Press
,
Cambridge, UK
.
7.
Andrews
,
S.
, and
Aslam
,
T.
,
2019
, “
On the Direct Construction of the Steady Traveling Solution to High Explosive Sandwich, Cylinder and Aquarium Test Via a Streamline Finite Volume Approximation
,”
J. Comput. Phys.
,
395
, pp.
653
670
.10.1016/j.jcp.2019.06.029
8.
Leveque
,
R. J.
,
2002
,
Finite Volume Methods for Hyperbolic Problems
(Cambridge Texts in Applied Mathematics),
Cambridge University Press
,
Cambridge, UK
.
9.
Horie
,
Y.
, ed.,
2007
,
Shock Wave Science and Technology Reference Library, Vol. 2: Solids I
,
Springer
,
Berlin
.
10.
Anderson
,
J. D.
,
1990
,
Modern Compressible Flow With Historical Perspective
,
McGraw-Hill
, New York.
11.
Tait
,
P.
,
1888
,
Physics and Chemistry of the Voyage of HMS Challenger
, Vol.
II
,
Neil and Company
,
Edinburgh, UK
.
12.
Lee
,
E. L.
,
Hornig
,
H. C.
, and
Kury
,
J. W.
,
1968
, “
Adiabatic Expansion of High Explosive Detonation Products
,” Lawrence Radiation Laboratory, Livermore, CA, Report No. UCRL 50422.
13.
Aslam
,
T. D.
,
2018
, “
Shock Temperature Dependent Rate Law for Plastic Bonded Explosives
,”
J. Appl. Phys.
,
123
(
14
), p.
145901
.10.1063/1.5020172
14.
Cooper
,
P. W.
,
1996
,
Explosives Engineering
,
Wiley
, Hoboken, NJ.
15.
Murnaghan
,
F. D.
,
1944
, “
The Compressibility of Media Under Extreme Pressures
,”
Proc. Natl. Acad. Sci.
,
30
(
9
), pp.
244
247
.10.1073/pnas.30.9.244
16.
Menikoff
,
R.
,
2012
, “
Deflagration Wave Profiles
,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-UR
12
20353
.
17.
Cook
,
A. W.
,
2009
, “
Enthalpy Diffusion in Multicomponent Flows
,”
Phys. Fluids
,
21
(
5
), p.
055109
.10.1063/1.3139305
18.
Powers
,
J. M.
,
2016
,
Combustion Thermodynamics and Dynamics
,
Cambridge University Press
,
Cambridge, UK
.
19.
Wescott
,
B. L.
,
Stewart
,
D. S.
, and
Davis
,
W. C.
,
2005
, “
Equation of State and Reaction Rate for Condensed-Phase Explosives
,”
J. Appl. Phys.
,
98
(
5
), p.
053514
.10.1063/1.2035310
20.
Menikoff
,
R.
, and
Shaw
,
M. S.
,
2012
, “
The SURF Model and the Curvature Effect for PBX 9502
,”
Combust. Theory Modell.
,
16
(
6
), pp.
1140
1169
.10.1080/13647830.2012.713994
21.
Menikoff
,
R.
,
2017
, “
SURF Model Calibration Strategy
,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-UR
17
22073
.
22.
Leiding
,
J. A.
,
Velizhanin
,
K. A.
,
Aslam
,
T. D.
,
Andrews
,
S. A.
,
Perriot
,
R. T.
, and
Cawkwell
,
M. J.
,
2021
, “
A New Parametrization of a 1-Step Thermal Decomposition Model of PBX-9501
,”
J. Appl. Phys.
(in preparation).
23.
Banks
,
J. W.
,
Aslam
,
T.
, and
Rider
,
W. J.
,
2008
, “
On Sub-Linear Convergence for Linearly Degenerate Waves in Capturing Schemes
,”
J. Comput. Phys.
,
227
(
14
), pp.
6985
7002
.10.1016/j.jcp.2008.04.002
24.
Sedov
,
L. I.
,
1953
, “
On Integration of One Dimensional Motion of a Gas
,”
Proc. USSR Acad. Sci.
,
XC
(
5
), p.
735
.
25.
Singleton
,
R. L.
,
Israel
,
D. M.
,
Doebling
,
S. W.
,
Woods
,
C. N.
,
Kaul
,
A.
,
Walter
,
J. W.
, and
Rogers
,
M. L.
,
2017
, “
Exactpack Documentation
,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-UR
16
23260
.
26.
Incropera
,
F. P.
, and
DeWitt
,
D. P.
,
2002
,
Fundamentals of Heat and Mass Transfer
, 5th ed.,
Wiley
, Hoboken, NJ.
27.
Schneider
,
P. J.
,
1955
,
Conduction Heat Transfer
,
Addison-Wesley
,
Cambridge, MA
.
28.
Maxima
,
2014
, “
Maxima, a Computer Algebra System. Version 5.34.1
,” accessed Jan. 4, 2022, http://maxima.sourceforge.net location is online
29.
Menikoff
,
R.
,
2015
, “
Detonation Wave Profiles
,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-UR
15
29498
.
30.
Lambert
,
D. E.
,
Stewart
,
D. S.
,
Yoo
,
S.
, and
Wescott
,
B. L.
,
2005
, “
Experimental Validation of Detonation Shock Dynamics in Condensed Explosives
,”
J. Fluid Mech.
,
546
(
1
), pp.
227
253
.10.1017/S0022112005007160
31.
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
T. E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
,
2020
, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,”
Nature Methods
, 17, pp. 261–272.10.1038/s41592-019-0686-2
32.
Hunter
,
J. D.
,
2007
, “
Matplotlib: A 2D Graphics Environment
,”
Comput. Sci. Eng.
,
9
(
3
), pp.
90
95
.10.1109/MCSE.2007.55
33.
Oliphant
,
T.
,
2006
,
A Guide to NumPy
,
Trelgol Publishing
.
34.
Behnel
,
S.
,
Bradshaw
,
R.
,
Citro
,
C.
,
Dalcin
,
L.
,
Seljebotn
,
D. S.
, and
Smith
,
K.
,
2011
, “
Cython: The Best of Both Worlds
,”
Comput. Sci. Eng.
,
13
(
2
), pp.
31
39
.10.1109/MCSE.2010.118
35.
Mattsson
,
A. E.
,
2016
, “
Short Introduction to Relations Between Thermodynamic Quantities
,” Sandia National Laboratories, Albuquerque, NM, Report No. SAND
2016
2112
.