0
Research Papers

The Escape of High Explosive Products: An Exact-Solution Problem for Verification of Hydrodynamics Codes OPEN ACCESS

[+] Author and Article Information
Scott W. Doebling

Verification & Analysis (XCP-8),
Los Alamos National Laboratory,
Los Alamos, NM 87544
e-mail: doebling@lanl.gov

Manuscript received April 21, 2015; final manuscript received October 18, 2016; published online December 2, 2016. Assoc. Editor: Urmila Ghia.The United States Government retains, and by accepting the article for publication, the publisher acknowledges that the United States Government retains, a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for United States government purposes.

J. Verif. Valid. Uncert 1(4), 041001 (Dec 02, 2016) (13 pages) Paper No: VVUQ-15-1018; doi: 10.1115/1.4035031 History: Received April 21, 2015; Revised October 18, 2016

This paper documents the escape of high explosive (HE) products problem and demonstrates the use of the problem for code verification assessment. The problem, first presented by Fickett and Rivard (1974, “Test Problems for Hydrocodes,” LASL Report, Los Alamos Scientific Laboratory, Los Alamos, NM, Report No. LA-5479), tests the implementation and numerical behavior of a high explosive detonation and energy release model and its interaction with an associated compressible hydrodynamics simulation code. The problem simulates the detonation of a finite-length, one-dimensional piece of HE that is driven by a piston from one end and adjacent to a void at the other end. The HE equation of state is modeled as a polytropic ideal gas. The HE detonation is assumed to be instantaneous with an infinitesimal reaction zone. Via judicious selection of the material specific heat ratio, the problem has an exact solution with linear characteristics, enabling a straightforward calculation of the physical variables as a function of time and space. Implementation of the exact solution in the Python code ExactPack is discussed, as are verification cases for the exact solution code. The problem is used to conduct code verification assessment on a Lagrangian hydrodynamics code.

In computational science and engineering, code verification is “the process of determining that the numerical algorithms are correctly implemented in the computer code and of identifying errors in the software” [1]. The two methods generally used for code verification are: (a) implementation of test problems with exact solutions and (b) the method of manufactured solutions. In this paper, a problem with an exact solution is presented that tests both the correct implementation of a high explosives (HE) detonation and energy release model in a computational physics code, as well as the correct coupling of that model with the compressible fluid dynamics algorithms of the computational physics code. Furthermore, the problem tests the correct calculation of the expansion of the detonation products into a void region.

The paper begins with the definition of the escape of HE products (EHEP) problem, including the history of the problem definition and a review of some past work. Next, the equations for the regions defining the exact solution in x − t space are presented, followed by a derivation of the equations for the physical variables in each region. Following that is a description of a solution algorithm using a point-in-polygon detection algorithm in the x − t plane, including verification cases from Ref. [2]. Next, the evolution of the flow is described, including a physical interpretation of the mathematical solution. Following that, the problem is used to perform code verification assessment on a Lagrangian hydryodynamics code. The paper concludes with a discussion of additional considerations and other cases yet to be explored using the EHEP problem.

The EHEP problem was first published by Fickett and Rivard [2]. The problem concerns a planar one-dimensional rod of HE extending from the origin for a length of x̃, as shown in Fig. 1. The HE is modeled as a polytropic ideal gas with equation of state (EOS) p=ρe(γ1), adiabatic index γ = 3.1 A shock wave, simulating a detonation front, travels through the material at the C–J detonation velocity D. For simplicity, the distinction between HE reactants and products is neglected, and the HE is modeled as a single material. To the left of the HE is a piston moving in the +x direction with velocity up. To the right of the HE is a void, which extends to location x̂.2 At t = 0, the shock wave departs from the origin in the +x direction at the C–J detonation velocity D, and the piston begins to move and isentropically compresses the reaction products. When the detonation wave reaches x̃ at time t̃, the HE has all been consumed, and the material begins to expand isentropically into the void region. At the same time t̃, the arrival of the detonation wave at an interface with lower impedance to the right causes a rarefaction wave to be propagated from x̃ in the −x direction. Fickett and Rivard recommend using the values x̃=1.0cm and up = 0.05 cm μs−1, and a time of tf = 5 μs for comparison of the exact solution to the hydrodynamics code result.

The problem has been analyzed previously to a limited extent. One analysis was documented by Painter [3], who was using this problem to test the implementation of HE burn models in a hydrodynamics code project called Blanca. The test problem helped to expose an error in the programmed burn model and another error in the hydrodynamic algorithm. A more thorough analysis of the problem was documented by Dykema et al. [4]. That report provides a derivation of the formulas for the solution characteristics in the x − t plane (which we rederive in this paper), and offers some physical interpretation of the solution. Those authors applied the problem to analyze HE burn models and hydrodynamic behavior in the 2D hydrodynamics code CALE.

In this paper, we derive the solution formulas for the x − t characteristics, as well as the formulas for the physical variables in each region of the x − t space, following Ref. [4]. Because the problem definition and solution have only been documented in internal reports [2,4], we seek to document the full description of the problem and solution in the archive literature. We also present an original implementation of a solution algorithm for the physical variables, based on a method to determine whether the point is within, or on the boundary of, a polygon in x − t space.

To define the polygons for the solution algorithm, we document the corner points of each region including points that lie on the boundary of the x − t analysis space. The x − t plane and associated polygons are shown in Fig. 2. This solution algorithm has been incorporated into the Python package ExactPack, which is documented and available for public download at the address shown in Ref. [5]. We document and apply verification cases for the exact solution code, as originally proposed in Ref. [2]. These test cases, and others, are included in the code self-testing module of ExactPack. Note that there is a typographical error in the tables in Ref. [2] that is corrected here.

We also correct a historical misinterpretation of the length scale of the problem as proposed by Fickett and Rivard. The original report cites a length of 10 units for the material, which has been previously interpreted to be 10 cm. We note that a detonation wave moving at D = 0.85 cm μs−1 will have traveled only 4.2 cm during the specified problem duration of 5 μs. In order for the HE products to “escape” into the adjacent void, the detonation wave must reach the end of the material before the final time of the problem. Dykema et al. [4] handled this discrepancy by increasing tf by a factor of 10 or more, which is adequate to generate both the exact solution and hydrodynamics code results as intended. Further review of Ref. [2], however, indicates that the original report uses length units of mm, such that the material extent should be 10 mm or 1.0 cm. By adopting the length x̃=1.0cm into our definition of the problem, we preserve the original scales of the problem in both time and space.

For completeness, the derivation of the exact solution is presented here. The exact solution first appeared in Ref. [2]. The derivation was first documented in Ref. [4], and the derivation here closely follows that one. For brevity, we will apply the constraint γ = 3 early on in the derivation, although it will not be motivated until later.

Thermodynamic Relations for Polytropic Gas Undergoing Chapman–Jouguet (C–J) Detonation.

Assume a polytropic gas, where the adiabatic index is defined as the negative logarithmic slope of the isentrope, using Eq. (2.15) of Ref. [6] Display Formula

(1)γ(lnplnv)S

where the subscript S indicates that the relation holds for constant entropy (i.e., along the isentrope). Integrate Eq. (1) to obtain the isentropic pressure–density relation for a polytropic gas Display Formula

(2)p=kργ=kρ3

where k is an arbitrary constant that will be solved for below by enforcing other thermodynamic relations. We now introduce the sound speed as the square root of the slope of the isentrope in p–ρ space, using Eq. (2.22) of Ref. [6] Display Formula

(3)cs2(pρ)S

Integrate Eq. (3) to obtain the isentropic sounds speed relation for a polytropic gas. Substituting Eq. (2) into Eq. (3), and enforcing γ = 3, yields Display Formula

(4)cs2=γpρ=3pρ

Now, we evaluate Eq. (2) at the C–J conditions to determine the constant k. Start with the C–J state for a polytropic gas as presented in Eq. (2A-22) thru Eq. (2A-25) in Ref. [6], and apply γ = 3, to get Display Formula

(5)pCJ=ρ0D2γ+1=ρ0D24ρCJ=γ+1γρ0=43ρ0uCJ=Dγ+1=D4csCJ=γDγ+1=34D

We evaluate the p–ρ relationship in Eq. (2) at the C–J condition using Eq. (5), and to solve for the constant k as Display Formula

(6)k=27256D2ρ02

Next, we solve Eqs. (2), (4), and (6) to obtain ρ and p. Specific internal energy, e, can be calculated using the ideal gas relation. Thus, the thermodynamic state for the gas can be written as Display Formula

(7)ρ=169ρ0(csD)p=1627ρ0D2(csD)3e=pρ(γ1)

This concludes the derivation of the thermodynamic relations, resulting in expressions for ρ, p, and e as functions of ρ0, D, and cs. The problem definition specifies ρ0 and D. The remainder of the derivation is to generate expressions for cs and u as functions of x and t.

Solution to Euler's Equations Using the Method of Characteristics.

The method of characteristics is a solution technique for one-dimensional unsteady flow, based on computing curves in x − t space, called characteristics. The characteristics have slopes dx/dt which are equal to the velocity at which a small disturbance will propagate in the medium. Along each characteristic there is a quantity that is constant, called the Riemann invariant. The derivation of these quantities is important for the analytic solution of this test problem. For more detail in the derivation, refer to Secs. 96 and 97 of Ref. [7].

Begin with Euler's equations for compressible, isentropic, inviscid flow written as Display Formula

(8)pt+upx+ρcs2ux=0
Display Formula
(9)ut+uux+1ρpx=0

Divide the first equation by ±ρcs and add it to the second to get Display Formula

(10)ut±1ρcspt+(ux±1ρcspx)(u±cs)=0

Define the Riemann Invariants as Display Formula

(11)J±u±dpρcs

Then, Eq. (10) can be rewritten as Display Formula

(12)[t+(u±cs)x]J±=0

Now, write an exact integral of the form Display Formula

(13)dJdt=Jt+Jxdxdt=[t+xdxdt]J

By comparing Eqs. (12) and (13), we conclude that dJ/dt = 0 when dx/dt = u ± cs, and therefore Display Formula

(14)J±=fwhendxdt=u±cs

where f is a constant. The solution curves are called characteristicsDisplay Formula

(15)C±dxdt=u±cs

Thus, we say that along the characteristic curves C±, the Riemann invariants J± are constant. For a polytropic ideal gas, Eq. (11) simplifies to Display Formula

(16)J±=u±2csγ1

To derive an analytic solution for C±, we ask: Under what circumstances are the characteristics straight lines? A sufficient condition for this is when u ± cs are constant. By inspecting Eq. (16), we see that the choice of γ = 3 produces Display Formula

(17)J±=u±cs

implying that, in a given solution region, u ± cs is invariant. Considering Eqs. (15) and (17) together, we conclude that dx/dt is constant within a given solution regime, such that the solutions x(t) are straight lines of the form Display Formula

(18)x=(u±cs)t+f±

where f± are functions of (u ± cs) determined by the initial and boundary conditions.

Solutions for Region Boundaries.

As the problem evolves in space and time, there are five region boundaries that need to be described in the x − t plane. These boundaries are straight lines in x − t space (enabled by choosing γ = 3), with the lines dividing the space into a series of polygons. The boundaries and corresponding regions for the exact solution are shown in Fig. 2. Inside each polygon, formulas for the values of the physical variables (ρ, p, e, u, cs) can be derived.

The derivation of the formulas for the straight-line boundaries is presented here, following Dykema et al. [4]. We add to that the derivation of the corner points where the straight-line boundaries intersect. We document the equations for and physical significance of the straight-line characteristics and the regions defined by those boundaries. We also define a “boundary point” for each boundary where it crosses the outer edge of the problem domain, as defined by xmax and tmax.

Boundary A.

At t = 0, the detonation wave begins moving from the origin in the +x direction at a velocity of D. The detonation wave can be described as Display Formula

(19)x=Dt

Past the point where the detonation ceases, (x̃,t̃), the slope of the C+ characteristic emanating from (x̃,t̃) must still satisfy the C–J conditions as stated in Eq. (5), so u = D/4 and cs = 3D/4. Then, the C+ characteristic has the form, from Eq. (18), of Display Formula

(20)x=(u+cs)t+f+=Dt+f+

Since this line must pass through (x̃,t̃), the constant f+ can be evaluated, and yields the equation for boundary A beyond (x̃,t̃) as Display Formula

(21)xx̃=D(tt̃)

which, since x̃=Dt̃, reduces to Display Formula

(22)x=Dt

The boundary point for boundary A is: Display Formula

(23)(xbdyA,tbdyA)=(xmax,xmaxD)

Boundary B.

Boundary B is the surface of the piston. The x − t relationship that describes this boundary is Display Formula

(24)x=upt

The boundary point for boundary B is: Display Formula

(25)(xbdyB,tbdyB)=(uptmax,tmax)

Boundary C.

Boundary C is the leading piston characteristic. At early times, it describes the point where the rarefaction region behind the shock levels out into the region of constant physical properties. This C+ characteristic has the form, from Eq. (18), of Display Formula

(26)x=(u+cs)t+f+

Since boundary C passes through the origin, f+ = 0. Boundary C must be consistent with boundary B, where u = up, and the C characteristic of boundary A, where u − cs = −D/2. Solving these two constraints together yields the expression for boundary C as Display Formula

(27)x=(2up+D2)t

The boundary point for boundary C is Display Formula

(28)(xbdyC,tbdyC)=((2up+D/2)tmax,tmax)

Boundary D.

Boundary D is the C characteristic of the wave that emanates from (x̃,t̃). Thus, u = D/4 and cs = 3D/4 (as with boundary A), and it has the form, from Eq. (18) of Display Formula

(29)x=(ucs)t+f=D2tf

Since this line must pass through (x̃,t̃), the constant f can be evaluated to be 3x̃/2, and yields the equation for boundary D as Display Formula

(30)x=12(3x̃Dt)

Boundary D does not have a boundary point as it does not extend to the edges of the analysis domain.

Intersection Points.

For the solution method we are using, it is necessary to define the locations of the points where the characteristic boundary lines intersect. The first point of interest is the intersection of boundaries A and D, which is simply the length of the HE medium x̃ and the arrival time of the shock at that location, t̃Display Formula

(31)(xAD,tAD)=(x̃,t̃)

where t̃=x̃/D.

The second point of interest is the intersection of boundaries C and D, which can be computed by simultaneously solving Eqs. (27) and (30) to get Display Formula

(32)(xCD,tCD)=3x̃2(1D4up+2D,12up+D)

The third point of interest is the intersection of boundaries B and D, which can be computed by simultaneously solving Eqs. (24) and (30) to get Display Formula

(33)(xBD,tBD)=3x̃2up+D(up,1)

Boundary E.

Boundary E describes the result of the interaction of the rarefaction wave of boundary D with the piston. This C+ characteristic has the form, from Eq. (18) of Display Formula

(34)x=(u+cs)t+f+

Observing that boundary E emanates from (xBD, tBD), immediately we can rewrite the expression as Display Formula

(35)xxBD=(u+cs)(ttBD)

Because boundary E intersects the piston, u = up, and because boundary D is a C characteristic emanating from (x̃,t̃), we know that J = u − cs = −D/2. Thus, cs = up + D/2, and boundary E has the form Display Formula

(36)xxBD=(2up+D/2)(ttBD)

and we can observe immediately that boundary E is parallel to boundary C. Invoking the expression for intersection point (xBD, tBD) from Eq. (33) yields the solution for boundary E Display Formula

(37)x=32x̃+(2up+D2)t

The boundary point for boundary E is Display Formula

(38)(xbdyE,tbdyE)=(32x̃+(2up+D/2)tmax,tmax)

Axis Boundary Points.

Additional boundary points are necessary in order to provide polygon vertices along the x and t axes: Display Formula

(39)(xbdyx,tbdyx)=(xmax,0)
Display Formula
(40)(xbdyt,tbdyt)=(0,tmax)

The boundaries defined above divide the x − t plane into several regions, as shown in Fig. 2. Each region has a distinct formula for the particle velocity u and the sound speed cs (these are derived for each region below). Once cs is known, the density ρ, pressure p, and specific internal energy e can be calculated as defined in Eq. (7). Region 00 is behind the piston and therefore outside of the problem domain. The remainder of the regions are described below. The solution strategy for each region is to identify C± and/or J± in that region, and then solve that system of equations for u and cs.

Regions 0V and 0H: Ahead of Detonation Wave.

Regions 0V and 0H are ahead of the detonation wave, so the pressure p and particle velocity u are both zero. In Region 0V, where x>x̃, which is ahead of the detonation wave in the vacuum, the density ρ is also zero. The vertices of the polygon that describes Region 0V are Display Formula

(41)(x̃,0),(x̃,t̃),(xmax,xmaxD),(xmax,0)

and the values of the thermodynamic variables in Region 0V are Display Formula

(42)cs=0u=0p=0ρ=0e=0

In Region 0H, where x<x̃, the density is equal to ρ0. The vertices of the polygon that describes Region 0H are Display Formula

(43)(0,0),(x̃,t̃),(x̃,0)

and the values of the thermodynamic variables in Region 0H are Display Formula

(44)cs=0u=0p=0ρ=ρ0e=0

Region I: Behind Detonation Wave.

Region I represents the rarefaction region just behind the shock wave, before it is impacted by either the rarefaction wave from the interface (boundary D) or the characteristic from the piston (boundary C). The polygon vertices for region I are Display Formula

(45)(0,0),(x̃,t̃),3x̃2(1D4up+2D,12up+D)

To derive the values for the thermodynamic variables in Region I, we first recognize that the C+ characteristics for this region emanate from the origin, thus they take the form Display Formula

(46)x=(u+cs)t

The C− characteristics in this region emanate from the detonation wave (boundary A); thus from the C–J conditions of Eq. (5)Display Formula

(47)J=ucs=D2

Solving these two equations simultaneously yields Display Formula

(48)cs=12(xt+D2)u=12(xtD2)

Region II: Isentropic Expansion Into Void.

Region II represents the isentropic expansion of the reaction products into the void. Thus, the density and pressure will drop, and the velocity will increase as the expansion proceeds in space and time. The polygon vertices for region II are Display Formula

(49)(x̃,t̃),3x̃2(1D4up+2D,12up+D)(xmax,xmax2up+D/2),(xmax,xmaxD)

The thermodynamic state in Region II can be determined by recognizing that the C+ characteristics emanate from the origin and thus have the form Display Formula

(50)x=(u+cs)t

The C characteristics emanate from (x̃,t̃), and thus have the form Display Formula

(51)xx̃=(ucs)(tt̃)

Solving these two equations simultaneously yields Display Formula

(52)cs=12(xtxx̃tt̃)u=12(xt+xx̃tt̃)

Region III: Fluid Pushed by Piston.

Region III has constant density and pressure, as the fluid is moved to the right by the action of the piston. It extends from the piston to the trailing wave behind the shock. The polygon vertices for Region III are Display Formula

(53)(0,0),3x̃2(1D4up+2D,12up+D),3x̃2up+D(up,1)

The thermodynamic properties in Region III are computed by recognizing that the C+ characteristics in this region emanate from the origin and thus have the form Display Formula

(54)x=(u+cs)t

The C characteristics are from boundary A, and thus satisfy Display Formula

(55)J=ucs=D2

Where the C characteristics intersection boundary B (the piston), it is known that u = up. Thus, from J, the solution in Region III is Display Formula

(56)cs=up+D2u=up

Region IV: Fluid Expansion and Rarefaction.

Region IV represents the expansion of the fluid from Region III after it interacts with the rarefaction wave coming off of the original fluid/void interface. The polygon vertices for Region IV are Display Formula

(57)3x̃2(1D4up+2D,12up+D),3x̃2up+D(up,1)(32x̃+(2up+D2)tmax,tmax),(xmax,xmax2up+D/2)

The values of the thermodynamic properties in Region IV are computed by first recognizing that the C+ characteristics are the same as in Region III Display Formula

(58)u+cs=2up+D2

and the C characteristics emanate from (x̃,t̃) and have the form Display Formula

(59)xx̃=(ucs)(tt̃)

Solving these equations simultaneously leads to the solution for Region IV Display Formula

(60)cs=up+D2(12xx̃Dtx̃)u=up+D2(12+xx̃Dtx̃)

Region V: Rarefaction Interaction With Piston.

Region V is the result of the original rarefaction wave from the fluid/void interface reflecting off of the piston. The polygon vertices for Region V are Display Formula

(61)3x̃2up+D(up,1),(uptmax,tmax),(32x̃+(2up+D2)tmax,tmax)

The values of the thermodynamic properties in Region V are computed by first recognizing that the C characteristics are the same as in Region IV Display Formula

(62)xx̃=(ucs)(tt̃)

Where these characteristics intersect the boundary B (the piston), we define a point (xp, tp). Along the surface of the piston we know that u = up and x = upt. Combining these constraints, along with x̃=Dt̃ into the expression for C, we get Display Formula

(63)uptpDt̃=(upcs(tp))(tpt̃)

We solve this equation for cs(tp) to get Display Formula

(64)cs(tp)=(Dup)t̃tpt̃

This expression is the sound speed along boundary B in the Region V at time tp. The C+ characteristics emanate from the piston surface at a point that changes with time, which we denote (xp, tp). The C+ characteristics have the form Display Formula

(65)xxp=(u+cs)(ttp)

Now, given a point (x, t) within Region V, we can use Eqs. (64) and (65) to compute the time tp at which C+ intersect boundary B. The result is Display Formula

(66)tp=t̃xupt+(Dup)txupt+(Dup)t̃

Now, we substitute Eq. (66) into Eq. (64) to get an expression for cs(t) along boundary B in Region V. The result is Display Formula

(67)cs(t)=xupt+(Dup)t̃tt̃

Now, recalling that u = up along boundary B, we can write the expression for the J+ Riemann invariant, u + cs, in region V as Display Formula

(68)J+=u+cs=x+(D2up)t̃tt̃

Solving Eqs. (68) and (62) simultaneously yields Display Formula

(69)cs=(Dup)t̃tt̃u=xupt̃tt̃

To evaluate the exact solution for a given value of x and t, the solution algorithm used is as follows:

  1. (1)Define a set of polygons corresponding to the regions defined in Eqs. (41), (43), (45), (49), (53), (57), and (61).
  2. (2)Given a point (x, t), determine in which polygon the point resides
  3. (3)Compute the values of cs and u using Eqs. (42), (44), (48), (52), (56), (60), or (69).
  4. (4)Compute the values of ρ, p, and e using Eq. (7).

The x − t plane and associated polygons are shown in Fig. 2. For step 2 above, an efficient method is needed to determine within which polygon the given point (x, t) lies. To accomplish this, we apply two functions to the list of polygon vertices assembled in step 1. The first function determines whether the point lies within the polygon. This function is built into the standard Python toolbox matplotlib and is called as contains_point,3 where poly is the list of vertices of the polygon, and point is the point being queried. This function returns True only if the point lies within the boundary of the polygon. We also need to query, however, whether the point lies on the boundary of the polygon. To accomplish that robustly, we need a second function. The second function is a custom function called point-on-boundary that determines whether the point lies on the boundary of the polygon. For this function, we loop over all pairs of vertices in the polygon, and query whether the point lies on the line formed by the pair of vertices to within a given tolerance, using a simple triangle inequality.

Fickett and Rivard [2] provide solution tables for the EHEP problem that can be used to verify the exact solution code, using the following generic values for Comp B explosive: Display Formula

(70)ρ0=1.6gcm3D=0.85cmμs1

These solution tables are shown in Tables 16. The implementation of the exact solution in ExactPack, found in the module4 reproduces the results in these tables, to the level of precision provided in the tables. These test cases are in the module5 within ExactPack. Note that there is a typographical error in one of the original tables in Ref. [2] that is corrected here.

We evaluate the exact solution using the values for Comp B explosive from Ref. [2], as listed in Eq. (70). The solutions of the physical variables in x − t space with (xmax, tmax) = (4.5 cm, 6.0 μs) are shown in Fig. 3. Snapshots of each of the physical variables at t = 0.5, 1.0, 2.0, 3.1, 4.0, 5.0 μs are shown in Fig. 4 for density, Fig. 5 for pressure, Fig. 6 for specific internal energy, and Fig. 7 for particle velocity.

We now describe the evolution of the physical solution at some representative times, and the reader is referred to Figs. 47 for illustration. The first phase of the problem involves a simple detonation wave traveling through the HE, while the piston begins its motion (at a velocity much lower than the detonation velocity). The physical variable snapshots at t = 0.5 μs and t = 1.0 μs are typical of this phase. The shock profile is readily apparent in the plots of each physical variable. The peak values at the shock are determined by the Rankine–Hugoniot shock-jump conditions, so they remain constant as the shock traverses the HE material. However, it can be seen that the extent of the rarefaction wave behind the shock grows during this phase. The region of constant conditions between the piston and the tail of the rarefaction wave is also apparent in the plot of each variable.

The next phase begins when the detonation wave reaches the end of the HE material. According to classical shock physics, when a shock wave impacts an interface with a material of lower impedance to the right, a shock wave is transmitted to the right, and a rarefaction wave is transmitted to the left. In this case, the void has, by definition, zero impedance and therefore no shock can be transmitted to the right. There is still a rarefaction wave generated that traverses the (now reacted) HE material back toward the piston. At the same time, the reacted HE material begins to expand isentropically into the void. This phase is fully developed at t = 2.0 μs. There is still a region of constant conditions between the piston and the head of the left-moving rarefaction wave. Perhaps the most startling development of this phase is the extreme increase in particle velocity that occurs as the gas expands into the void. The gas, with high pressure and energy, is no longer confined once the detonation wave attenuates, and so expands rapidly. The leading edge of the gas will accelerate until it reaches the detonation velocity D, for sufficiently large values of x and t. It can be seen in Fig. 7 that the leading-edge velocity has already reached D at t = 2.0 μs. The observation that u reaches a maximum value of D as the gas expands can be proven mathematically by taking the limit of u from Eq. (52) as xx̃ and tt̃, and by considering that x = Dt along boundary A Display Formula

(71)limxx̃,tt̃12(xt+xx̃tt̃)=xt=D

The next phase begins at about t = 3.1 μs, when the left-moving rarefaction wave impacts the piston. This occurs when boundary D intersects boundary B in Fig. 2. Region III disappears as the region of constant properties evident at t = 2.0 μs shrinks to a point against the piston wall, and then re-emerges as a rarefaction traveling to the right. The physical variable snapshots at t = 4.0 μs are typical of this phase. The region of constant properties starts to grow again. The gas continues to expand into the void with leading edge velocity D.

The final state of the problem is considered by Fickett and Rivard to be at t = 5.0 μs. This is their recommended time to compare the hydrodynamics code results to the exact solution. All of the physical regimes of the problem have been exercised, including shock compression, rarefaction, reflection, and isentropic expansion. The piston continues its slow traverse to the right, the gas continues to expand into the void with leading-edge velocity D, and the region of constant properties continues to increase in extent even as the magnitudes of the density, pressure, and specific internal energy decrease. Eventually, the density, pressure, and specific internal energy all tend toward smaller and smaller constant values distributed over a larger and larger extent, and the leading edge continues to move into the void at velocity D.

In this section, the EHEP problem will be employed to perform a code verification assessment on the Los Alamos National Laboratory code FLAG, a 1D/2D/3D unstructured-grid Arbitrary Lagrangian–Eulerian (ALE) multiphysics code [8]. For the purposes of this example, the code will be used in staggered-grid Lagrangian hydrodynamics mode with a 1D Cartesian mesh.

The mesh consists of equal-length axial segments, extending from the origin to the recommended distance of x̃=1.0cm. (Note that it is not necessary to choose an extent for the void material because in this Lagrangian analysis, everything outside of the defined mesh is inherently assumed to be a void.) The material is assigned an ideal gas equation of state with γ = 3. The HE is initialized using the recommended properties of Comp B as described in Eq. (70) (ρ0 = 1.6 g cm−3, D = 0.85 cm μs−1). The heat of reaction q is required as an input to the code, and is computed using Eq. (2.34) in Ref. [6] Display Formula

(72)q=D2/2(γ21)

which yields a value of q = 0.04515625 Mbar cm3 g−1. At time t = 0, the HE is detonated from the left end. The detonation model used for this analysis is a Huygens construction (known internal to the code as the Lund model [9]), forming a wavefront that moves at a constant velocity through the mesh. It is a programmed burn model, meaning that the detonation is calculated based on the expected arrival times of waves emanating from the specified detonation points, and is not reactive to the values of the physical variables. The energy release of the HE is sourced into the energy equation at each zone as the detonation wave reaches it. The boundary condition on the left end of the mesh is a constant velocity moving to the right at up = 0.05 cm μs−1. The right end is free to expand, modeling expansion into the void. The calculation is run to a final time of t = 5.0 μs. The NZ = 50 mesh and velocity profile for the calculation at several simulation times are shown in Fig. 8.

The basic code verification process used here, following Refs. [1,10] is:

  1. (1)Run calculation using a series of successively finer grids, NZ = 50, 100, 200, 400 corresponding to zone sizes of h = 0.02, 0.01, 0.005, 0.0025.
  2. (2)For each grid, generate the exact solution at the center point of each zone and at each nodal point.
  3. (3)Compute weighted L1 error norm between exact solution and code output at a series of time snapshots.
  4. (4)Compute convergence rate of the weighted L1 error norm at each time snapshot.
  5. (5)Make determination about correctness of code behavior based on results.

The quantities of interest for code verification are the zonal physical variables (density, pressure, and specific internal energy), and the nodal velocity. The verification error metric will be the weighted L1 error norm between the exact and calculated solution values for each quantity of interest: Display Formula

(73)εn=i=1NZwi|ficalcfiexact|i=1NZwi

For the zonal quantities, wi is the length of the zone, and for nodal quantities, wi = 1.

A key element of the verification analysis is the calculation of the convergence rate of the L1 error norm as a function of the initial grid length hi. The standard error ansatz for verification analysis is assumed Display Formula

(74)ε=Ahpobs

where ε is the vector of errors εn as computed using Eq. (73), and h is the vector of initial zone length values hn. The value of pobs is determined from Eq. (74) using nonlinear optimization [11].

The formal order of convergence for the staggered-grid hydrodynamic algorithm being studied is palg = 2. As documented in Ref. [12], the expected error convergence rate of the L1 norm for this type of problem is palg/(palg + 1), which for this algorithm would yield an expected convergence rate of 2/3. This value will be used as the fiducial for this verification analysis.

The verification analysis is conducted at three values of simulation time t. The value of t = 5.0 μs is recommended by Fickett and Rivard [2] in the original definition of the problem. This time is after the shock wave has attenuated, the rarefaction wave has reflected off of the piston, and the gas has expanded a significant distance into the void. Conducting a verification analysis at this time is valuable as all regimes of the problem have been exercised. In particular, the response at this time will be sensitive to the code correctly representing the expansion of the high-velocity gas into the void, which the author believes was the main focus of the problem at the time of the original problem definition. However, it is also valuable to assess the correctness of the solution at earlier times, in order to isolate the physics of interest in particular regimes of the problem. Therefore, the verification analysis is also conducted at t = 0.5 μs, when the detonation wave is moving through the gas, and t = 3.1 μs, when the initial rarefaction wave interacts with the piston. To simplify the illustration of the results, one zonal quantity (density) and one nodal quantity (velocity) are shown for each time snapshot. The remaining zonal quantities (pressure and specific internal energy) exhibit convergence behavior similar to that of density.

The code verification results at t = 0.5 μs are shown for density in Fig. 9 and for velocity in Fig. 10.6 The density and velocity both converge at a rate closed to the fiducial value. At the finest resolution, the velocity exhibits some overshoot of the exact solution at the shock.

The code verification results at t = 3.1 μs are shown for density in Fig. 11 and for velocity in Fig. 12. At this time, the density converges at a slightly higher rate than at t = 0.5 μs, with about the same magnitude of accuracy. However, the velocity fails to converge at this time. It is apparent from Fig. 12 that the velocity is inaccurate at the expansion front. This inaccuracy is driving the lack of global convergence. This situation may indicate that the code, as currently used, is not going to accurately calculate the velocity of a gas expansion into a void region (although it does accurately calculate the zonal variables). To further explore the situation, consider Fig. 13 that shows the convergence analysis at t = 3.1 μs, but restricted to the spatial domain between x = 0 and x = 1.0 cm. The velocity in this case is convergent at a rate comparable to the fiducial value. This case supports the hypothesis that the source of the nonconvergence of velocity at t = 3.1 μs is the region of the expansion front.

The code verification results at t = 5.0 μs are shown for density in Fig. 14 and for velocity in Fig. 15. The calculations at this time demonstrate a trend similar to that observed at t = 3.1 μs; specifically, the density converges, the velocity over the entire spatial domain does not converge, and the velocity over a domain that excludes the expansion front (see Fig. 16) converges. Once again, inaccuracy in the velocity at the expansion front is driving the nonconvergence, although the zonal variables are accurate and convergent.

This article documents the exact solution to the EHEP problem and presents a discussion of the physical interpretation of the solution. Application of the problem to a simple verification analysis of a hydrodynamics code indicates that the calculated values of zonal variables (density, pressure, and specific internal energy) demonstrate convergent behavior at multiple times in the simulation over the entire problem domain. The nodal velocities, however, are not convergent near the expansion front. This result can help the code user to interpret the credibility of more complex calculations from this code and should be a topic of further study.

The EHEP problem, as originally published by Fickett and Rivard, is quite useful for code verification but has some limitations. One limitation is the choice of up allowed under this exact solution. This article and the exact solution code assume properties of up also assumed by Fickett and Rivard, specifically that the piston is moving slowly to the right, with velocity bounded as 0 < up < uCJ. This is the case of an unsupported detonation wave contained by a slow-moving piston.7 If the piston has value up = uCJ, then, the theory states that the flow transitions from unsupported to supported, the rarefaction zone behind the shock (Region I in Fig. 2) vanishes and the region of constant properties (Region III in Fig. 2) extends from the piston to the shock front (boundary B to boundary A in Fig. 2). The value of uCJ under the assumption of instantaneous detonation with infinitesimal reaction zones, and polytropic EOS, is given in Eq. (2.19) of Ref. [6] as uCJ = D/(γ + 1). In this problem γ = 3, thus, the value up = D/4 is where the flow transitions from unsupported to supported. If the piston exceeds this velocity, the flow is considered to be overdriven. Referring to Fig. 2, graphically the situation of the rarefaction region (Region I) disappearing occurs when boundary A and boundary C are coincident. Solving Eqs. (22) and (27) simultaneously yields up = D/4, identical to what is predicted from the theory above. Thus, the flow is unsupported when up < D/4, supported when up = D/4, and overdriven when up > D/4. In this article, only the forward-moving, unsupported case is considered (0 < up < D/4). Another limitation of this problem is that the HE EOS is limited to a polytropic gas formulation. The assumption of a polytropic gas with γ = 3 is fundamental to the derivation of the exact solution, and it is unclear whether the exact solution could be extended to other cases.

This work was sponsored by the Advanced Simulation and Computing program within the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396. This document has been approved for unlimited release by Los Alamos National Laboratory under release number LA-UR-15-22547. The author wishes to thank Los Alamos National Laboratory colleagues James Kamm, Nathaniel Morgan, James Quirk, Diane Vaughan, and Nathan Woods for review of this manuscript and helpful recommendations.

  • cs =

    material sound speed

  • C± =

    characteristic curves

  • C–J =

    pertaining to Chapman–Jouguet detonation theory

  • D =

    detonation velocity

  • e =

    specific internal energy

  • ficalc,fiexact =

    quantity of interest, calculated and exact, for zone i

  • hn =

    initial zone length for grid n

  • J± =

    Riemann invariants

  • NZ =

    number of zones in hydrodynamics calculation

  • p =

    pressure

  • palg =

    formal order of convergence of the computational algorithm

  • pobs =

    observed convergence rate of the error in quantity of interest

  • q =

    heat of reaction (per unit mass)

  • std =

    measure of goodness of fit of error ansatz to convergence data

  • t =

    time

  • tf =

    final time

  • tmax =

    maximum value of time domain

  • u =

    particle velocity

  • up =

    piston velocity

  • wi =

    weighting for norm calculation for zone i

  • x =

    position

  • xmax =

    maximum value of position domain

  • x̃ =

    linear extent of HE material

  • x̂ =

    linear extent of void material

  • γ =

    adiabatic index (specific heat ratio)

  • ρ =

    density

  • εn =

    error in quantity of interest for grid n

  • The default units for this paper are cm, g, μs.

ASME 2006 , “ Performance Test Code Committee 60: Verification and Validation in Computational Solid Mechanics,” Guide for Verification and Validation in Computational Solid Mechanics, American Society for Mechanical Engineers, New York, Standard No. ASME V&V 10-2006.
Fickett, W. , and Rivard, C. , 1974, “ Test Problems for Hydrocodes,” Los Alamos Scientific Laboratory, Los Alamos, NM, Report No. LA-5479.
Painter, J. W. , 2001, “ BLANCA Project: More on the Escape of HE Products Problem,” Los Alamos National Laboratory, Los Alamos, NM, Repot No. LA-UR-01-6549.
Dykema, P. , Brandon, S. , Bolstad, J. , Woods, T. , and Klein, R. , 2002, “ Level 1 V. & V. Test Problem 10: Escape of High Explosive Products,” Lawrence Livermore National Laboratory, Livermore, CA, Report No. UCRL-ID-150418.
Doebling, S. W. , Israel, D. M. , Singleton, R. L. , Woods, C. N. , Kaul, A. , and Walter, J. , 2016, “ ExactPack v1.4,” Technical Report, Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-CC-14-047, accessed on Sept. 26, 2016, http://github.com/losalamos/exactpack
Fickett, W. , and Davis, W. C. , 1979, Detonation: Theory and Experiment, University of California Press, Berkeley, CA.
Landau, L. D. , and Lifshitz, E. M. , 1959, Fluid Mechanics, Pergamon Press, Oxford.
Burton, D. , 1990, “ Conservation of Energy, Momentum, and Angular Momentum in Lagrangian Staggered-Grid Hydrodynamics,” Lawrence Livermore National Laboratory, Livermore, CA, Report No. UCRL-JC-105926.
Mandell, D. , Burton, D. , and Lund, C. , 1998, “ High Explosive Programmed Burn in the Flag Code,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-13406.
V&V20 Committee, 2009, “ Verification and Validation in Computational Fluid Dynamics and Heat Transfer,” American Society for Mechanical Engineers, New York, Standard No. ASME V&V 20-2009.
Srinivasan, G. , Vaughan, D. , and Doebling, S. , 2016, “ Verification of Code Convergence Order Using Regression and Optimization Methods,” ASME J. Verif. Valid. Uncert. (submitted).
Banks, J. , Aslam, T. , and Rider, W. , 2008, “ On Sub-Linear Convergence for Linearly Degenerate Waves in Capturing Schemes,” J. Comput. Phys., 227(14), pp. 6985–7002. [CrossRef]
Copyright © 2016 by ASME
View article in PDF format.

References

ASME 2006 , “ Performance Test Code Committee 60: Verification and Validation in Computational Solid Mechanics,” Guide for Verification and Validation in Computational Solid Mechanics, American Society for Mechanical Engineers, New York, Standard No. ASME V&V 10-2006.
Fickett, W. , and Rivard, C. , 1974, “ Test Problems for Hydrocodes,” Los Alamos Scientific Laboratory, Los Alamos, NM, Report No. LA-5479.
Painter, J. W. , 2001, “ BLANCA Project: More on the Escape of HE Products Problem,” Los Alamos National Laboratory, Los Alamos, NM, Repot No. LA-UR-01-6549.
Dykema, P. , Brandon, S. , Bolstad, J. , Woods, T. , and Klein, R. , 2002, “ Level 1 V. & V. Test Problem 10: Escape of High Explosive Products,” Lawrence Livermore National Laboratory, Livermore, CA, Report No. UCRL-ID-150418.
Doebling, S. W. , Israel, D. M. , Singleton, R. L. , Woods, C. N. , Kaul, A. , and Walter, J. , 2016, “ ExactPack v1.4,” Technical Report, Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-CC-14-047, accessed on Sept. 26, 2016, http://github.com/losalamos/exactpack
Fickett, W. , and Davis, W. C. , 1979, Detonation: Theory and Experiment, University of California Press, Berkeley, CA.
Landau, L. D. , and Lifshitz, E. M. , 1959, Fluid Mechanics, Pergamon Press, Oxford.
Burton, D. , 1990, “ Conservation of Energy, Momentum, and Angular Momentum in Lagrangian Staggered-Grid Hydrodynamics,” Lawrence Livermore National Laboratory, Livermore, CA, Report No. UCRL-JC-105926.
Mandell, D. , Burton, D. , and Lund, C. , 1998, “ High Explosive Programmed Burn in the Flag Code,” Los Alamos National Laboratory, Los Alamos, NM, Report No. LA-13406.
V&V20 Committee, 2009, “ Verification and Validation in Computational Fluid Dynamics and Heat Transfer,” American Society for Mechanical Engineers, New York, Standard No. ASME V&V 20-2009.
Srinivasan, G. , Vaughan, D. , and Doebling, S. , 2016, “ Verification of Code Convergence Order Using Regression and Optimization Methods,” ASME J. Verif. Valid. Uncert. (submitted).
Banks, J. , Aslam, T. , and Rider, W. , 2008, “ On Sub-Linear Convergence for Linearly Degenerate Waves in Capturing Schemes,” J. Comput. Phys., 227(14), pp. 6985–7002. [CrossRef]

Figures

Grahic Jump Location
Fig. 1

The escape of HE products problem

Grahic Jump Location
Fig. 2

The x − t diagram of the exact solution regions

Grahic Jump Location
Fig. 3

The exact solution of the physical variables in the x − t plane

Grahic Jump Location
Fig. 4

Snapshots of the exact solution for density at various times

Grahic Jump Location
Fig. 5

Snapshots of the exact solution for pressure at various times

Grahic Jump Location
Fig. 6

Snapshots of the exact solution for specific internal energy at various times

Grahic Jump Location
Fig. 7

Snapshots of the exact solution for particle velocity at various times

Grahic Jump Location
Fig. 8

Mesh and velocity values for FLAG 1D calculation at four times

Grahic Jump Location
Fig. 9

Convergence of calculated density at t = 0.5 μs

Grahic Jump Location
Fig. 10

Convergence of calculated velocity at t = 0.5 μs

Grahic Jump Location
Fig. 11

Convergence of calculated density at t = 3.1 μs

Grahic Jump Location
Fig. 12

Convergence of calculated velocity at t = 3.1 μs

Grahic Jump Location
Fig. 13

Convergence of calculated velocity at t = 3.1 μs over domain x = [0, 1.0] cm

Grahic Jump Location
Fig. 14

Convergence of calculated density at t = 5.0 μs

Grahic Jump Location
Fig. 15

Convergence of calculated velocity at t = 5.0 μs

Grahic Jump Location
Fig. 16

Convergence of calculated velocity at t = 5.0 μs over domain x = [0, 3.0] cm

Tables

Table Grahic Jump Location
Table 1 Verification data for EHEP exact solution code at t = 1.25 μs, reproduced from Table 6.1 of Ref. [2]
Table Grahic Jump Location
Table 2 Verification data for EHEP exact solution code at t = 2.0 μs, reproduced from Table 6.1 of Ref. [2]
Table Grahic Jump Location
Table 3 Verification data for EHEP exact solution code at t = 5.0 μs, reproduced from Table 6.1 of Ref. [2] (boldface indicates that typographical error in original document has been corrected).
Table Grahic Jump Location
Table 4 Verification data for EHEP exact solution code at example times and positions, reproduced from Table 6.2 of Ref. [2]
Table Grahic Jump Location
Table 5 Verification data for EHEP exact solution code at example times and positions, reproduced from Table 6.2 of Ref. [2]
Table Grahic Jump Location
Table 6 Verification data for EHEP exact solution code at example times and positions, reproduced from Table 6.2 of Ref. [2]

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In