Verification Studies for the Noh Problem using Non-ideal Equations of State and Finite Strength Shocks

TL;DR

This paper extends the Noh verification problem to realistic equations of state and uses self-similarity solutions for code verification of shock simulations in various geometries.

Contribution

It develops exact solutions for the Noh problem with non-ideal EOSs, enabling more realistic and comprehensive code verification for shock physics.

Findings

Exact solutions for non-ideal EOSs in Noh problem

Applicable to planar, cylindrical, and spherical shocks

Used for verification of the FLAG hydrocode

Abstract

The Noh verification test problem is extended beyond the commonly studied ideal gamma-law gas to more realistic equations of state (EOSs) including the stiff gas, the Noble-Abel gas, and the Carnahan-Starling EOS for hard-sphere fluids. Self-similarity methods are used to solve the Euler compressible flow equations, which in combination with the Rankine-Hugoniot jump conditions provide a tractable general solution. This solution can be applied to fluids with EOSs that meet criterion such as it being a convex function and having a corresponding bulk modulus. For the planar case the solution can be applied to shocks of arbitrary strength, but for cylindrical and spherical geometries it is required that the analysis be restricted to strong shocks. The exact solutions are used to perform a variety of quantitative code verification studies of the Los Alamos National Laboratory Lagrangian hydrocode FLAG.

Tables

Table 1. Table 1: Table of convergence orders for planar stiff gas. The units for c s subscript 𝑐 𝑠 c_{s} are c ​ m / μ ​ s 𝑐 𝑚 𝜇 𝑠 cm/\mu s and for p 0 subscript 𝑝 0 p_{0} the units are M ​ b ​ a ​ r 𝑀 𝑏 𝑎 𝑟 Mbar .

	0.00	0.01	0.10	1.00
0.0	0.998	0.987	1.017	0.997
0.1	0.993	1.062	0.989	0.974
0.5	0.998	0.978	1.001	0.988
0.6	1.018	0.988	1.017	0.987
1.0	1.011	0.975	1.008	1.039
1.5	0.973	1.087	0.975	0.989
2.0	0.987	1.019	1.019	0.985

Equations

$$\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial r}+\rho \left(\frac{\partial u}{\partial r}+\frac{mu}{r}\right)$$

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}+\frac{1}{\rho }\frac{\partial p}{\partial r}$$

$$\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial r}+K\left(\frac{\partial u}{\partial r}+\frac{mu}{r}\right)$$

$$K\left(p,\rho \right)\equiv \rho {\left(\frac{\partial p}{\partial \rho }\right)}_S=\rho {\left(\frac{\partial p}{\partial \rho }\right)}_e+\frac{p}{\rho }{\left(\frac{\partial p}{\partial e}\right)}_{\rho }$$

$$\xi =\frac{ar}{t}$$

$$-\xi \frac{d\rho }{d\xi }+au\frac{d\rho }{d\xi }+\rho \left(a\frac{du}{d\xi }+\frac{mau}{\xi }\right)$$

$$-\xi \frac{du}{d\xi }+au\frac{du}{d\xi }+\frac{a}{\rho }\frac{dp}{d\xi }$$

$$-\xi \frac{dp}{d\xi }+au\frac{dp}{d\xi }+K\left(a\frac{du}{d\xi }+\frac{mau}{\xi }\right)$$

$$(u_2-D){\rho }_2$$

$$p_2+(u_2-D){\rho }_2{u}_2$$

$$e_2+\frac{p_2}{{\rho }_2}+\frac{1}{2}(u_2-D{)}^2$$

$$r_s(t)=\frac{t}{a}=D\cdot t.$$

$${\rho }_1$$

$$p_1$$

$$u_1$$

$${\rho }_2$$

$$p_2$$

$$e_2$$

$$u_2$$

$$D=\frac{p_1-p_2}{{\rho }_2{u}_1}.$$

$$r_s(t)=D\cdot t=\left(\frac{p_1-p_2}{{\rho }_2{u}_1}\right)t.$$

$$K\left(\frac{mau}{\xi }\right)=0.$$

$${\left(\frac{\partial e}{\partial \rho }\right)}_p=\frac{p}{{\rho }^2}$$

$$e(\rho ,p)+\frac{p}{\rho }=f(p)$$

$${\rho }_1=\text{constant}$$

$$e_1+\frac{p_1}{{\rho }_1}=\text{constant}.$$

$$e_1+\frac{p_0}{{\rho }_0}{\left(1-\frac{u_{0}t}{r}\right)}^{-m}=\text{constant}$$

$$p(\rho ,e)=\rho e(\gamma -1)$$

$$p(\rho ,e)=p_0+\left[c_s^2-(\gamma -1)\frac{p_0}{{\rho }_0}\right](\rho -{\rho }_0)+\rho (\gamma -1)(e-e_0)$$

$$p=\rho e(\gamma -1)+c_s^2(\rho -{\rho }_0).$$

$$p_{\infty }$$

$$e_{\infty }$$

$$p+p_{\infty }=\rho (\gamma -1)(e+e_{\infty })$$

$$p(\rho ,e)=\frac{\rho e(\gamma -1)}{1-b\rho }$$

$$p(\rho ,e)=\rho e(\gamma -1)Z(\eta )$$

$$Z(\eta )=\frac{1+\eta +{\eta }^2-{\eta }^3}{(1-\eta {)}^3}$$

Full text

1 Introduction

The Noh problem [1, 2, 3] is a well-studied and widely used verification test problem in the field of computational hydrodynamics of compressible fluids. Typically, it is applied to an ideal polytropic gas [2, 4] and serves as a tool for testing the accuracy of sophisticated hydrocodes. Axford first obtained analytic solutions to the Noh problem for the non-ideal stiff gas EOS and for a specialized form of the Mie-Grüneisen EOS [5]. Recently the topic was revisited and generalized by Ramsey et al.[3], whose formal solutions we employ here. As the Noh problem is also an example of the broader class of Riemann problems, the one-dimensional planar studies of Menikoff and Plohr, Banks, Whitham, and Toro are also relevant to the development of exact, non-ideal solutions in this context [6, 7, 8, 9].

In this work we review the general solution of the Noh problem, which is most rigorously derived using Lie group methods [3], and then apply the general solution to three more realistic, non-ideal EOSs: the stiffened gas [10, 11, 12, 13, 14, 15, 16] , the Noble-Abel gas [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], and the Carnahan-Starling EOS for a model fluid composed of hard spheres [29, 30, 31]. The exact solutions for these non-ideal EOSs are then compared to the numerical results obtained using the Los Alamos National Laboratory (LANL) compressible-flow code FLAG [32, 33, 34, 35].

The Noh problem consists of a strong shock forming as a compressible gas moves at a constant velocity towards a rigid wall (in planar geometry) or as the gas is compressed with constant radial velocity towards a central axis (cylindrical) or a point (spherical) in the curvilinear cases. (See Figs. 1 and 2.) Alternatively, and equivalently, in planar geometry the problem can be visualized as shock wave driven by a piston moving at a constant velocity traveling into a quiescent ideal gas [2]. The problem is defined by specifying the geometry; the EOS of the fluid; and the initial velocity, density, and specific internal energy of the fluid. It is typically applied to an ideal gas initialized at zero pressure and energy (a convenient, but rather unrealistic, initial condition corresponding to a temperature of absolute zero), giving rise to an infinitely strong shock since the speed of sound at the initial conditions is zero. In this paper, we explore extensions of the Noh problem beyond these idealized conditions.

The purpose of this investigation is threefold: to demonstrate the feasibility of extending hydrodynamic verification studies to more physically realistic materials and initial conditions; to use the results to assess a particular hydrocode’s ability to accurately capture the shock formation in the Noh problem; and to elucidate the spatial convergence order of the hydrocode under non-ideal gas conditions. In previous studies, acceptable accuracy has been found to be first order in spatial-convergence tests for hydrocodes including any discontinuities [36, 37, 38, 39].

We begin by reviewing the general one-dimensional solution of the Noh problem derived using the conservation equations, self-similarity methods, and Rankine-Hugoniot jump conditions. Following this, we describe the EOSs and the motivation for choosing these particular forms for the verification studies. Specific analytic solutions for each EOS are then compared with numerical results obtained using the Lagrangian hydrocode FLAG, so that the spatial-convergence order can be calculated. Outcomes for the relative error between numerical and exact results and the observed spatial-convergence order are presented for all three non-ideal EOSs, with a particular focus on planar geometry, where finite strength shocks can also be considered. Results for infinitely strong shocks in cylindrical and spherical geometries (still one-dimensional) using the Noble-Abel EOS are also presented. We conclude with a summary of our findings and ideas for future work.

2 Theory

An analytic solution to the Noh problem can be constructed from piecewise isentropic solutions of the inviscid one-dimensional Euler compressible flow equations in combination with the Noh initial and boundary conditions and the Rankine-Hugoniot jump conditions [15] applied across the shock interface. In the absence of an external field (i.e. no gravity or outside forces), the one-dimensional Euler equations for mass, momentum and internal energy are[11, 40]

[TABLE]

where $t$ denotes time and $r$ represents distance from the rigid wall, center axis, or center point depending on whether the geometry is planar, cylindrical, or spherical (distinguished by $m=0$, $1$, or $2$, respectively). The quantities $\rho$, $u$, $p$, and $K$ represent density, velocity, pressure, and adiabatic bulk modulus, respectively. The density and pressure state variables are related to the specific internal energy, $e$, through an EOS of the form $e(\rho,p)$ or $p(\rho,e)$. The EOS enters into the analysis via the adiabatic bulk modulus, $K$[5], defined as

[TABLE]

where $S$ denotes entropy. The second identity conveniently expresses the bulk modulus as a function of $e$ and $\rho$, common independent thermodynamics variables in hydrocodes, making it straightforward to derive $K$ from an EOS written in the form $p(e,\rho)$.

Lie group methods [40, 41] leverage the symmetry properties of the one-dimensional Euler equations to construct a similarity variable of the form

[TABLE]

where $a$ is a constant parameter with units of inverse velocity. Substitution of Eq. (5) into Eqs. (1)-(3) enables the reduction of the PDEs in question to ODEs in terms of $\xi$,

[TABLE]

The Noh solution is piecewise solution, with the discontinuity in the flow interpreted as a shock wave. As a result, we apply the Rankine-Hugoniot jump conditions [15] at the position of the shock, $r_{s}(t)$, to ensure conservation of mass, momentum, and energy across the discontinuity

[TABLE]

where the subscripts $1$ and $2$ represent the unshocked and shocked regions, respectively, and where $D=dr_{s}/dt$ is the shock velocity. For Eqs. (9)-(11) to be invariant under the same group transformations as the PDEs (i.e., for Eqs. (9)-(11) to be expressible solely in terms of $\xi$), $D=\text{const.}$ and thus,

[TABLE]

The final contribution to the general solution is the conditions of the shocked and unshocked regions characteristic of the Noh problem. Note that $u_{1}$ and $u_{2}$ have pre-defined values as consequences of the assumption on the velocity field in the unshocked and shocked regions. The assumed conditions are constant and inwardly directed velocity in the unshocked region, $u_{1}=u_{0}=\text{const.}<0$, and zero velocity in the shocked region, $u_{2}=0$. The effect is the shock wave propagating outward as seen in Fig. 2.

2.1 General Noh Solution

The Noh solution can be found by solving the one-dimensional system of ODEs, Eqs. (6)-(8), and applying the Rankine-Hugoniot jump conditions and the conditions unique to the Noh problem. The following procedure may be used to solve the Noh problem:

By applying the constant velocity condition in the unshocked region, Eqs. (6)-(8) can be reduced and Eq. (6) integrated with the initial condition of density to get the general solution for first the unshocked density, then the unshocked pressure from Eq. (7), and lastly, the unshocked energy from Eq. (8) by selecting either $K=0$ or $m=0$. 2. 2.

The condition of zero velocity can be applied to Eqs. (6) and (8) within the shocked region. From this we notice that it must be the case that both pressure and density in the shocked region are constant. 3. 3.

The density in the shocked region can be calculated by taking the conservation of mass equation from the jump conditions, Eq. (9), and substituting $u_{2}=0$, $u_{1}=u_{0}$, $D=1/a$ and $\rho_{1}$. 4. 4.

The pressure in the shocked region is found by taking the conservation of momentum equation from the jump conditions, Eq. (10), substituting for $u_{2}$, $u_{1}$, $p_{1}$ and $\rho_{2}$ then solving for pressure. 5. 5.

The energy in the shocked region can also be calculated by taking the conservation of energy equation from the jump conditions, Eq. (11), substituting for $u_{2}$, $u_{1}$, $p_{1}$, $\rho_{2}$ and $p_{2}$ and then solving for energy.

In this fashion, the general solution for the unshocked region is found to be

[TABLE]

The general solution for the shocked region is

[TABLE]

Lastly, we can apply the condition $u_{2}=0$ to the conservation of momentum Rankine-Hugoniot jump condition, Eq. (10), and solve for the shock speed to get

[TABLE]

Given Eq. (12), the shock location is determined to be

[TABLE]

The last element required is the EOS in the form $e_{2}(\rho_{2},p_{2})$ and $e_{1}(\rho_{1},p_{1})$ to create a system of equations with a closed-form solution [3]. If a closed-form solution cannot be found, as with the Carnahan-Starling EOS in Section 2.6, a highly accurate implicit solver can be used. If we relax the definition of verification and allow quasi-analytic solutions to be the exact solution for verification studies, then we can perform verification studies on the Noh problem with any convex EOS.

2.1.1 Constraints on the Initial Conditions for the Curvilinear cases

The constant velocity Noh condition imposes a severe constraint on the initial pressure in cylindrical and spherical geometries. Since $u_{1}=u_{0}$ and $p_{1}=p_{0}$ are both constant, the energy Euler equation, Eq. (8), reduces to

[TABLE]

Since $u$ is a negative constant in the unshocked region and $a$ must be a non-zero constant to have a nontrivial solution then either $m=0$ or $K=0$. In the planar case ($m=0$), this analysis shows that the general solutions can be applied to an EOS that is a convex function and has a corresponding bulk modulus form or with any initial pressure since $K$ need not equal zero. However, if implementing the curvilinear cases ($m\neq 0$), $K$ must be equal to zero by Eq. (22). Since $K=\rho c_{s}^{2}$, where $c_{s}$ is the is the sound speed in the fluid, this implies that $c_{s}=0$ and that the shock will be infinitely strong.

The requirement of a zero initial sound speed in cylindrical and spherical geometries, in turn, constrains the initial pressure. Substituting $K=0$ into Eq. (4) leads to

[TABLE]

which can be integrated to give

[TABLE]

where $f(p)$ is an arbitrary function. Within the unshocked region, where $e(\rho,p)=e_{1}$ and $p=p_{1}$, Eq. (24) can only be satisfied if either

[TABLE]

or

[TABLE]

Since we have $m=1$ or $2$ and $\rho_{1}$ is not constant (see Eq. (13)), Eq. (25) cannot be true so Eq. (26) must be applied to meet the conditions of conservation and the Noh problem conditions. Substituting the density solution from the unshocked region, Eq. (13), and $p_{1}=p_{0}$, into Eq. (26) provides

[TABLE]

implying that $p_{0}=0$ and $e_{1}=\text{constant}$ when $m\neq 0$. Thus, the curvilinear cases can only be studied provided that $p_{0}=0$ [3].

The primary focus of this paper is the planar geometry case when $m=0$, where the solution is piecewise constant. An example for the cylindrical and spherical geometries will be demonstrated with the Noble-Abel EOS for infinitely strong shocks.

2.2 Equations of State

In its original incarnation [1], and nearly every application since [2], the Noh problem used the ideal polytropic gas EOS [4, 14], initialized at zero pressure and specific internal energy, to model the thermodynamic behavior of the fluid,

[TABLE]

where $\gamma$ is the polytropic constant. Typically $\gamma$ is set equal to $5/3$, which corresponds to a monatomic ideal gas; however, it can be viewed simply as an adjustable parameter and empirically set to other values depending on the fluid and the thermodynamic path being modeled [4].

While Eq. (28) admits a particularly simple solution to the Noh problem, its behavior (as well as the conventional initial conditions) is far removed from most materials of practical interest. Thus there is significant value, from a verification and validation perspective, in extending the domain of applicability of the Noh problem to more realistic EOSs and initial conditions. We consider here three more complex EOSs: the stiff gas, useful for metals and liquids; the Noble-Abel, useful for explosive gases; and the Carnahan-Starling equation for a hypothetical atomistic model fluid composed of rigid, noninterpenetrating, hard spheres.

The stiff gas EOS [3, 5, 10, 11, 12, 13, 14], may be viewed as either a modification of the ideal polytropic gas EOS or a simplification (Taylor expansion in density) of the Mie-Grüneisen EOS when compressions are modest. In its most general form, it may be written

[TABLE]

where $p_{0}$, $\rho_{0}$, and $e_{0}$ denote the pressure, density, and specific internal energy at some convenient reference condition (e.g., the initial conditions of the problem of interest) and $\gamma$ is an adjustable constant such that $\Gamma=\gamma-1$ plays the role of an effective Grüneisen parameter [15].

If we choose our reference state such that $p_{0}$ and $e_{0}$ are zero, then Eq. (29) simplifies to the more common form

[TABLE]

This form of the stiff gas equation has been found to be qualitatively accurate for a variety of metals. Though similar in appearance to the ideal gas, Eq. (28), it enables a condensed-phase material to go into tension as well as compression, and allows use of a physically realistic bulk modulus though the parameter $c_{s}^{2}$. Weixin [16] has used this form of the stiff gas to model Al, Fe, Cu, and W and reports appropriate values for $c_{s}$ and $\gamma$. Harlow and Pracht [10] have used it to model Fe and Al using somewhat different values for $c_{s}$ and $\gamma$.

Alternatively, if one defines new parameters

[TABLE]

then Eq. (30) may be reexpressed as

[TABLE]

which is the form used by Courant to model water [14] and by Cole to model shock compression of salt water up to 8 GPa [13].

The significance of Eq. (29) and its simplifications, Eq. (30) and Eq. (33), is that the EOS remains simple enough algebraically to permit useful analytical manipulations, such as an exact solution to the Noh problem [3, 5], but is realistic enough to provide qualitatively accurate predictions for shock compression of liquids and metals, thus enabling code verification studies on physically realistic materials. In the remainder of this work, we focus on the Eq. (30) form, where the density, $\rho_{0}$, is set to the initial density [11].

The second EOS we consider is the Noble-Abel equation [17, 18, 19, 20, 21] (also sometimes referred to as the Clausius I [22, 23] or Hirn EOS [24, 25, 26]), which takes the form

[TABLE]

where the adjustable parameter $b$, referred to as the co-volume (units of volume per mass), is intended to account, approximately, for the finite size of the molecules (compare with Eq. (28)) [23, 25]. The original idea for this type of correction to the ideal gas traces its roots to Bernoulli [27, 28], and reemerges in Hirn’s, Clausius’s, and van der Waals’s subsequent EOS developments [42, 24, 25, 27] in the nineteenth century. For dilute gases, $b$ may be deduced from the second virial coefficient, and it can take on negative as well as positive values if the temperature is low enough that intermolecular attractive forces decrease the pressure below the ideal gas value. In hot or strongly compressed fluids, where it is generally more accurate, it takes on positive values (subject to the constraint $\rho b<1$) and functions to increase the pressure above the ideal gas value.

Eq. (34) has found extensive applicability in ballistics modeling [17, 18, 19], where it provides a simple and reasonably accurate EOS for propellant gases at the high densities and temperatures experienced in guns. It has also been used to model the release of hydrogen into air from high-pressure tanks [20], and in modeling the transfer and discharge of various other compressed gases [21].

The third EOS we consider is the Carnahan-Starling equation [29, 30, 31]. In contrast to the polytropic ideal gas, stiff gas, and Noble-Abel EOSs, which can be empirically fitted to a variety of fluids, the Carnahan-Starling EOS was developed for a very specific model of a prototypical atomisitic fluid, one composed of noninterpentrating, hard spheres. The hard-sphere fluid has long been of interest in the field of statistical thermodynamics, in part because the simple form of the interatomic potential greatly reduces calculations of properties such as virial coefficients and radial distribution functions; in part because the specific internal energy is exactly that of a monatomic ideal gas; and in part, because x-ray scattering experiments have shown that the local structure of dense real fluids, including liquid metals, is dominated by the harshly repulsive component of the interatomic potential, which can be well approximated by a hard sphere of an appropriately chosen diameter. Viewed this way, the hard-sphere fluid becomes a simple reference system about which the effects of interatomic attractions present in real fluids can be added perturbatively [30, 43].

Carnahan and Starling [29] observed an approximate, serendipitous recursion relationship between the second through the sixth virial coefficients of the hard-sphere fluid which, when postulated to apply equally well to the higher virial coefficients, allowed the infinite series to be summed exactly. The resulting EOS takes the form

[TABLE]

where $Z(\eta)$ is the compressibility factor

[TABLE]

and $\eta=b\rho$ is the packing fraction of the spheres. Here again, $b$ denotes the co-volume, now interpreted as the volume per mass occupied by the rigid spheres.

In extensive comparisons with Monte Carlo and molecular dynamic simulations of hard-sphere fluids, Eq. (35) has been found to be remarkably accurate over nearly the entire fluid range, from the dilute, ideal gas limit, $\eta\to 0$, up to packing fractions, $\eta\simeq 0.45$ [29, 30, 31].

2.3 Ideal Gas Solution

The ideal gas EOS, Eq. (28), can be rearranged in terms of energy as

[TABLE]

Ramsey et al. [3] discuss how the general solution to the Noh problem (Eqs. (13)-(21)) can be combined with the EOS to calculate the shock speed $D=1/a$. If we consider only the planar ($m=0$) case then the shock location is

[TABLE]

where

[TABLE]

The density for planar ideal gas solution for the unshocked region is $\rho_{1}=\rho_{0}$ and the density for the planar ideal gas solution for the shocked region is

[TABLE]

As shown in Section 2.1.1, for $m=1$ or $2$ solutions exists only in the infinitely strong shock case ($K_{0}=0$, $p_{0}=0$) [3]. The shock location is

[TABLE]

The density in the unshocked region is

[TABLE]

and the density in the shocked region is

[TABLE]

2.4 Stiff Gas Solution

The stiff gas EOS, Eq. (30), may be written in terms of energy as

[TABLE]

Combining the EOS with the general solution, the shock speed is

[TABLE]

where

[TABLE]

The density in the unshocked region is $\rho_{1}=\rho_{0}$ and the density in the shocked region is

[TABLE]

2.5 Noble-Abel Gas Solution

The Noble-Abel EOS is Eq. (34) may be written in terms of energy as

[TABLE]

Using the general Noh solution and the EOS, the shock speed is

[TABLE]

where

[TABLE]

The density for the planar solution for the unshocked region is $\rho_{1}=\rho_{0}$ and the density for the planar solution for the shocked region is

[TABLE]

As with the ideal gas solution, for $m=1$ or $2$ Noble-Abel gas solutions exists only in the infinitely strong shock case ($K_{0}=0$, $p_{0}=0$) [3]. The shock speed is

[TABLE]

The density for the Noble-Abel gas in the unshocked region is

[TABLE]

and the density for the Noble-Abel gas in the shocked region is

[TABLE]

2.6 Carnahan-Starling Gas Solution

The Carnahan-Starling gas EOS, Eq. (35), may be rearranged in terms of energy as

[TABLE]

Combining this with Eqs. (13)-(21) the planar Carnahan-Starling fluid solution for the unshocked region is

[TABLE]

Given $u_{2}=0$, the variables, $a$, $\rho_{2}$, $p_{2}$ and $e_{2}$ can be solved implicitly from the system,

[TABLE]

For the convergence studies on the Carnahan-Starling EOS, we employed the Powell Hybrid Method using the fsolve function via the optimize package in Python [44]. The accuracy of the implicit solver ($\sim 10^{-8}$) makes the numerical solutions viable as quasi-analytic solutions for verification studies.

3 Verification

Verification studies compare numerical implementations with analytically derived solutions. These studies are done in order to evaluate the accuracy and performance of numerical techniques and algorithms, in this case, for a hydrocode based on conservation laws. The FLAG code has been subject to rigorous verification and validation research. Some of the more recent efforts include the studies of the invariance with respect to the heat conduction in the Coggeshall problem [45], the Guderley shock problem [46], self-similar shock wave propagation [47], the impact of artificial viscosity on the Sedov problem [48], the Sedov shock wave blast [49], and the motivation behind the Godunov-like staggered grid hydrodynamic approach used in FLAG [50]. The common theme of all these investigations (which is exemplary of the larger body of work performed on the code to date) centers around the use of the ideal gas EOS coupled to the code’s hydrodynamic solver. Thus the motivation of our work: examples of verification studies performed using the same code and compressible-flow solver, but coupled to simple non-ideal EOS classes.

3.1 FLAG Hydrocode

The Los Alamos Lagrangian hydrocode FLAG uses a compressible-flow algorithm featuring artificial viscosity and ancillary grid stability methods. The solver implements explicit, geometry-compatible, finite volume spatial full-Lagrangian staggered-grid hydrodynamics. This method uses control volumes with zones that are assigned by drawing lines between nodes. The densities, pressures, and specific internal energies live at zone centers and the velocities live at the nodes. Time integration features a predictor-corrector method based on pressure or energy projection for momentum evolution, as discussed by Burton [32] and Caramana et al. [35]. For smooth (i.e., shock-free) gas dynamics problems, the hydrodynamics solver is second-order accurate in space and time (further assuming a zero linear artificial viscosity coefficient) [32, 33, 34, 35].

FLAG is a high-resolution numerical scheme with accuracy limited by the presence of the shock, which by Godunov’s order barrier theorem restricts the studies to first-order spatial convergence [38, 39]. Our hydrocode studies use the Barton artificial viscosity model. In compression, the quadratic viscosity coefficient is 2.0 and the linear viscosity coefficient is 0.3. The linear viscosity coefficient in expansion is 0.3 [51, 52]. (In the context of the Noh problem, the FLAG algorithms become equivalent to Reimann-solver based methods [6, 7, 8, 9]). These constructs are used to capture the zero-thickness of the jump discontinuity on a grid of finite-thickness and the nonzero linear coeffiencient are also responsible for expectation that the spatial convergence is first order [36, 37]. Another factor that alters the formation of the shock profile is the wall heating produced by the shock reflecting at the origin. This is caused by the pressure initially being computed and then energy/temperature added to the system, via the EOS, to compensate for a shortage in density [53]. The convergence order will remain first order despite refining the spatial mesh size because of the presence of the artificial viscosity and Gudunov’s theorem.

3.2 Verification Studies

FLAG results for the pressure, density, and specific internal energy as a function of time, position, and mesh resolution are compared with the exact solutions for the ideal and the non-ideal gases. Representative results for the density as a function of time are shown in Figs. 3-7 for the Noble-Abel fluid as a function of initial pressure, co-volume, and mesh size. The plots were produced using the verification tool ExactPack [54, 55]. All results use the same initial density $\rho_{0}=$1\text{,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}, the same adiabatic coefficient $\gamma=5/3$ and are measured at time $t=$0.6\text{\,}\mu\mathrm{s}.

In 3-6 the points are the FLAG results at indicated mesh resolutions ($\Delta r=0.01$, $0.005$, $0.0025$, and $0.001\,25\text{\,}\mathrm{c}\mathrm{m}$) and the lines represents the exact solution. In Figs. 3-6, results are shown for the Noble-Abel EOS in the planar and spherical cases. Figure 6 is the corresponding relative error between the FLAG results and the exact solution corresponding to the studies in Fig. 5. For the cylindrical and spherical cases, note the simulation remains a one-dimensional with modifications for the converging geometry taken into consideration.

In every case, the code behaves consistently with the exact solution except at wall and the shock front – deficiencies common to many hydrocodes. The representative plots shown here for the Noble-Abel EOS display the same trends as the other EOSs. The simplest trend is the effect of the increase in initial pressure on the shock speed and density. In all cases, as the unshocked pressure increases, the shocked density decreases and the shock speed increases. Increasing the initial pressure increases the sound speed at all densities, which then leads the shock to move faster. Also, as the co-volume, $b$, increases in the Carnahan-Starling gas and Noble-Abel gas (Fig. 4, 5) or as the reference sound speed, $c_{s}$, increases in the material of the stiff gas, the density in the shocked region decreases and the shock speed increases. When the reference sound speed of the stiff gas increases, the shock moves faster for the same reason as when the initial pressure is increased. By increasing the co-volume in the Noble-Abel gas and the Carnahan-Starling gas, the fluid becomes less compressible resulting in a smaller density jump across the shock front and faster moving wave. Conversely, when these EOSs have a negative co-volume, the molecules attract each other resulting in higher shocked densities and lower shock speeds than the ideal gas.

The simulation results also reveal additional correlations between the various EOS parameters and the most prominent computed solution errors. From the relative error plots shown in Fig. 6, the computed errors in the density at $r=0$ and at the shock front decrease as the co-volume, $b$, increases. Similar results were found when increasing the EOS parameter co-volume and initial pressure for the other EOSs. At first glance, it might appear that the drop in relative error at the discontinuity is connected to the shock’s position relative to the origin; the closer the shock is to the origin, the more wall heating influences the precise shape of the shock front. Upon closer investigation, by looking at a shock at earlier and later times, we found that the distribution of the relative error does not change. Therefore, we can conclude that the compressibility or sound speed in the fluid is what changes the distribution of the relative error at the shock front, not simply the distance from the wall. Higher co-volumes, sound speeds, and initial pressures all are observed to reduce the magnitude of the error at the wall and the jump discontinuity.

For infinitely strong shock cases, shown for the Noble-Abel gas in Figs. 6, the wall heating increases with increasing curvilinearity (planar to spherical). In cylindrical and spherical geometries, the change in the compression term, $b$, more dramatically affects the wall heating than it does in the planar case due to the additional dissipation associated with the flow converging on an axis or a single point.

First-order spatial convergence with respect to mesh size is observed with $\Delta r=0.01$, $0.005$, $0.0025$, and $0.001\,25\text{\,}\mathrm{c}\mathrm{m}$. The verification tool ExactPack [54, 55] analyzes the results by calculating the L1-norm [56, 57]

[TABLE]

where $y_{E}$ is the exact solution and $y_{F}$ is the FLAG data. The weights $w_{i}$ are determined by the dimensionality of the problem. For the 1D Cartesian problem, they are the zone lengths. The L1-norm with respect to the mesh size is fitted on a log-log scale to a line in the form

[TABLE]

Figure 7 shows a typical result obtained from an ExactPack convergence order study for the Noble-Abel EOS.

The slope, $p$, is determined from Eq. (65) and is known as the convergence order fitting this equation,

[TABLE]

The convergence orders for all the EOS mentioned in this paper are presented in Tables 1-4. Tables 1-2 and 4 are planar geometry while varying the EOS parameter and the initial pressure. Table 4 shows the results for the cuvilinear geometries ($m=0,1,$ and $2$) using the Noble-Abel EOS and $p_{0}=0$.

The results of Tables 1-4 fit our expectations of FLAG’s performance. The convergence studies of the code in the majority cases were close to being spatially first-order accurate.

4 Conclusion

Our results illustrate how non-ideal EOSs can provide a more physically-realistic complement to the ideal gas in verifying complex hydrocodes. Proper verification studies that feature shocks typically exhibit first-order spatial convergence [36, 37]. Our studies demonstrate that the LANL hydrocode FLAG models the Noh solution with first-order accuracy for both ideal and several non-ideal fluids when compared to the exact solutions. Additionally, the decrease in the amount of error as the shock speed increases is significant. We observe that the higher the initial pressure and/or the larger the initial bulk modulus, the smaller the magnitude of the error both at the wall and the shock front.

Further analysis could be done by examining pointwise verification metrics such as the shock location as a function of the non-ideal parameters, examining other EOSs, varying the adiabatic constant, examining shorter time scales, and varying the initial density. In addition to modifying parameters, it would be interesting to compare to a verification study utilizing an Eulerian hydrocode such as the LANL code, RAGE [58]. Research has already been published on the stiff gas EOS using the RAGE hydrocode[59]. The Eulerian code was shown to give better results for the Noh problem than the Lagrangian code as discussed by Rider [53]. In our results, the lowest convergence orders were seen in the curvilinear cases so this should be the focus of the comparison study. Additionally, one might undertake multidimensional studies and see how the results compare to the 1D curvilinear cases. Finally, exact solutions can be derived by applying Lie group methods to other problems such as the Noh-2 problem [1], the Sedov problem [60, 61], the Guderley problem [62, 63, 46, 64], and the Sod Problem [65]. Verification studies would be needed to show that the FLAG hydrocode performed as expected for these problems since they are considered more complex than the Noh problem.

{acknowledgment}

This work was performed under the auspices of the United States Department of Energy by Los Alamos National Security, LLC, at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. Extensive use of the LANL verification tool ExactPack was made in this study. The authors gratefully acknowledge the support of the U.S. DOE Advanced Strategic Computing Program, in particular the Physics Verification Project and Physics and Engineering Models Project.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1322106. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.