Verification and Validation of Two Hydrodynamic Methods for Simulations of High Energy Density Physics Problems

Abstract

A 3D verification and validation suite of test problems is presented and used to evaluate hydrodynamic methods within a radiation hydrodynamics code, xRAGE. These test problems exercise different levels of complexity, building towards ICF problems which in addition to hydrodynamics also include three temperature plasma physics, thermal conduction, and radiation diffusion. Among the problems in the test suite are the Kidder ball problem, the Verney shell problem, and a 5-material compression problem, which exercise different purely hydrodynamic methods implemented within xRAGE. There is excellent agreement between 2D and 3D XRAGE simulation results and between the xRAGE results and the benchmark solutions. Two 3D ICF test problems are also presented, based on an OMEGA direct drive capsule experiment and on a NIF indirect drive capsule experiment. It is demonstrated that the newer unsplit hydrodynamic method in xRAGE produces more vorticity relative to the older default method. For the indirect drive capsule, the 3D simulations are in reasonable agreement with the experimental values of ion temperature and neutron production.

1. Introduction

A plethora of high energy density physics (HEDP) experiments including inertial confinement fusion (ICF) capsule experiments have been undertaken over the past decade which have included 3D geometric features and defects, creating a need for verification and validation of 3D capabilities within the Eulerian arbitrary mesh refinement (AMR) radiation hydrodynamics code xRAGE [1]. Many HEDP experiments involving laser platforms such as OMEGA and NIF require 3D simulation capabilities in order to fully understand the complex interaction of various physical processes and accurately predict performance metrics. Implicit large eddy simulations (ILES) in 3D with xRAGE have provided insight into instability growth and the transition to turbulence for problems such the Taylor-Green vortex [2]. Furthermore, certain kinds of ICF capsule experiments have prominent 3D features such as DIME capsules with equatorial trenches [3], MARBLE capsules with 3D porous foam [4], and NIF high foot capsules that include a tent supporting structure and fill tubes [5]. These 3D features are the focus of ongoing studies to better quantify their effect on capsule performance.

Numerous works have been completed for verification and validation of the numerical methods for hydrodynamics, radiation diffusion, three temperature (3T) plasma physics, and thermal conduction within xRAGE. A few prior examples of verification and validation are provided here. The default xRAGE hydrodynamic method has been tested using a diverse set of analytic test problems as well as detailed experiments involving Richtmyer-Meshkov instabilities in planar shock tubes [1, 6, 7]. In addition, when coupled with other appropriate HEDP models such as thermal conduction, the default hydrodynamic method has also been validated against several laser-driven experiments including 3D supersonic jets, Mach reflection experiments, shear layers, and collapsing cylinders with single-mode Rayleigh-Taylor instabilities [8–13]. Augmenting the large amount of validation data available from laser-driven experiments, several researchers have developed analytic problems to further test HEDP methods in radiation hydrodynamics codes. For instance, McClarren and Wohlbier developed a 3T variant of the Su-Olson radiative shock wave problem with a spherical energy source, which they then used to test the 3T plasma physics method in xRAGE in three dimensions [14]. Furthermore, xRAGE has been used extensively for 3D simulations of OMEGA direct drive ICF capsules [15, 16] and NIF indirect drive ICF capsules including MARBLE capsules with 3D porous foam [4].

Arguably, the most important algorithm within xRAGE is the hydrodynamics method, as it is responsible for the evolution of material position which ultimately sets the environment for all other physical processes in the simulation. This work explores three distinct hydrodynamic methods that have been implemented within xRAGE for 3D simulations: (1) the default Godunov-type hydrodynamic method that is directionally split [1], (2) the default method using a volume of fluid (VOF) interface treatment, and (3) the newer directionally unsplit hydrodynamic method which is a higher-order Godunov-type approach incorporating quadratic polynomial reconstructions for velocity [17]. The default directionally split algorithm is a variation of the second-order MUSCL and PPM formulations [18, 19] for solving the Euler equations, where velocity is linearly reconstructed and time-centered within a zone. Left and right velocity values at a zone boundary from the linear reconstructions are used to solve an acoustic Harten-Lax-van Leer (HLL) Riemann problem [20], which is then used to find mass, momentum, and energy fluxes at that zone boundary. The unsplit method dispenses with the directional splitting of the default method, uses a method of lines approach for time discretization, and solves a different acoustic Riemann problem based on the Harten-Lax-van Leer-contact (HLLC) approach [21]. The unsplit method has proven to be superior to the default method for simulating certain problems such as the Sedov blast wave problem, but the numerical implementation of this method within xRAGE is not compatible with the VOF algorithm developed for the default method and is thus limited when simulating sharp material interfaces.

The governing equations of radiation hydrodynamics are presented in Section 2, and a few general details about the numerical implementation within xRAGE are discussed including choices for opacity and thermal diffusivity models and the method of operator splitting for separating hydrodynamic and diffusive contributions to the governing equations, as well as the numerical technique for updating electron and ion energies due to the electron-ion energy exchange term. This is followed by the development of a comprehensive test suite including new test problems to evaluate these methods. In Section 3, all three hydrodynamic methods are used to simulate the isentropic spherical compression problem, which is a variant of the Kidder ball problem. The default method with VOF was coupled with a simple elastic-plastic material strength model in Section 4 in order to simulate the Verney shell collapse problem. A more sophisticated 5-material compression problem that uses the same simulation approach is presented in Section 5. Finally, in Sections 6 and 7, the default and unsplit methods without VOF are coupled with 3T physics, thermal conduction, and radiation diffusion, in order to simulate direct drive and indirect drive ICF capsule test problems.

2. Numerical Methods

The governing equations of the 3T electron-ion radiation hydrodynamics conservation laws in Eulerian form [22] which are implemented within xRAGE are presented below, in a simplified form assuming the ions are composed of a single species:

(1) $\begin{array}{c}\frac{\partial \rho }{\partial t}+\nabla \cdot \rho u=0,\frac{\partial \rho u}{\partial t}+\nabla \cdot \left(\rho u\otimes u+P_e+P_i\right)=0,\\ \frac{\partial \rho e_i}{\partial t}+\nabla \cdot \left(\rho e_{i}u\right)-\nabla \cdot \left({\kappa }_i\nabla T_i\right)+P_i:\nabla u={\gamma }_{ei}\left(T_e-T_i\right)+S_i,\\ \frac{\partial \rho e_e}{\partial t}+\nabla \cdot \left(\rho e_{e}u\right)-\nabla \cdot \left({\kappa }_e\nabla T_e\right)+P_e:\nabla u={\gamma }_{ei}\left(T_i-T_e\right)+c{\sigma }_a\left(E_{r}a{T}_e^4\right)+S_e,\\ \frac{\partial \rho E}{\partial t}+\nabla \cdot \left(\left(\rho EI+P\right)\cdot u\right)-\nabla \cdot \left({\kappa }_i\nabla T_i+{\kappa }_e\nabla T_e\right)=c{\sigma }_a\left(E_{r}a{T}_e^4\right)+S_i+S_e,\\ \frac{\partial E_r}{\partial t}+\left(\frac{4}{3}\right)\nabla \cdot \left(u{E}_r\right)-\nabla \cdot \left(\kappa \nabla E_r\right)=-c{\sigma }_a\left(E_r-a{T}_e^4\right)+\left(\frac{1}{3}\right)u\cdot \nabla E_r,\end{array}$

where ρ is the mass density, u is the fluid velocity, $P=P_e+P_i$ is the pressure tensor (equal to minus the Cauchy stress tensor); P_e is for the electrons and P_i for the ions, e_i and e_e are the specific internal energies of the ions and electrons (both per unit ion mass), $E=e_i+e_e+v^2/2$ is the total specific energy, T_i and T_e are the ion and electron temperatures, ${\gamma }_{ei}$ is the electron-ion coupling term, S_i and S_e are ion and electron internal energy density sources, κ_i and κ_e are the ion and electron thermal conductivities, c is the speed of light, σ_a is the absorption opacity, E_r is the frequency averaged energy in the radiation field, and κ is the radiation diffusion coefficient. There are equations for mass, momentum, and total energy conservation, as well as an equation for conservation of radiation energy density, E_r. In addition, there are separate energy equations for electrons and ions. For simplicity, the equations mentioned above describe a single ion species; however, the implementation within xRAGE has been generalized to accommodate an arbitrary number of ion species, and for this more general case, the species volume fractions within a zone are evolved so that pressure and temperature equilibrium is enforced among the species within that zone. xRAGE also has the option not to enforce pressure-temperature equilibrium within a zone, and this capability is exercised for the 5-material compression problem in Section 5. The coupling of energy between electrons and ions is an important part of the 3T plasma physics model. ${\gamma }_{ei}$ is determined using the method of Brown et al. [23]. The ion and electron thermal conductivities are calculated using the approach of Lee and More [24]. In calculating the ion-electron coupling term and the thermal conductivities, the effects of partial ionization are neglected for the test problems considered in this work. For the ICF test problems in Sections 6 and 7, tabulated values of absorption opacity from the OPLIB database [25] are used for the simulations, as determined by the TOPS code [26]. The governing equations mentioned above are written for the more general case where the off diagonal terms of the pressure tensor are not zero, which is required for problems involving material strength which have deviatoric stress such as the Verney problem in Section 4 and the 5-material compression problem in Section 5. For the ICF test problems considered in Sections 6 and 7, it is assumed that the pressure tensor takes the simplified form of a scalar pressure, $p=p_i+p_e$ , times the identify matrix, P = pI. It is further assumed that p is a function of density and temperature, represented by an equation of state (EOS) model. xRAGE can use several different EOS models including an ideal gas model and a SESAME tabular EOS model [27].

The governing equations presented above are solved by using the technique of operating splitting, wherein two sets of simplified equations are derived [1]. The first set of simplified equations includes the advection terms in the governing equations, but excludes source terms, ion-electron coupling, thermal conductivity, and radiation diffusivity. This set is solved first, using either the default directionally split method with the HLL Riemann solver or the unsplit method that uses the HLLC Riemann solver. Then, a second set of simplified equations is solved, and these simplified equations exclude the advection terms but include the 3T physics such as source terms, ion-electron coupling, thermal conductivity, and radiation diffusivity. The first set of simplified equations in the operator splitting scheme is given below, for the case of the simplified pressure tensor, P, containing only diagonal terms.

(2) $\begin{array}{c}\frac{\partial \rho u}{\partial t}+\nabla \cdot \left(\rho u\otimes u\right)+\nabla P=0,\frac{\partial \rho e_i}{\partial t}+\nabla \cdot \left(\rho e_{i}u\right)+p_i\nabla \cdot u=0,\\ \frac{\partial \rho e_e}{\partial t}+\nabla \cdot \left(\rho e_{e}u\right)+p_e\nabla \cdot u=0,\\ \frac{\partial \rho E}{\partial t}+\nabla \cdot \left(\left(\rho E+p\right)\cdot u\right)=0,\\ \frac{\partial E_r}{\partial t}+\left(\frac{4}{3}\right)\nabla \cdot \left(u{E}_r\right)=0.\\ \frac{\partial \rho }{\partial t}+\nabla \cdot \rho u=0,\end{array}$

The second set of simplified equations in the operator splitting scheme is

(3) $\begin{array}{c}\frac{\partial \rho e_i}{\partial t}-\nabla \cdot \left({\kappa }_i\nabla T_i\right)={\gamma }_{ei}\left(T_e-T_i\right)+S_i,\frac{\partial \rho e_e}{\partial t}-\nabla \cdot \left({\kappa }_e\nabla T_e\right)={\gamma }_{ei}\left(T_i-T_e\right)+c{\sigma }_a\left(E_r-a{T}_e^4\right)+S_e,\\ \frac{\partial \rho E}{\partial t}-\nabla \cdot \left({\kappa }_i\nabla T_i+{\kappa }_e\nabla T_e\right)=c{\sigma }_a\left(E_r-a{T}_e^4\right)+S_i+S_e,\\ \frac{\partial E_r}{\partial t}-\nabla \cdot \left(\kappa \nabla E_r\right)=-c{\sigma }_a\left(E_r-a{T}_e^4\right)+\left(\frac{1}{3}\right)u\cdot \nabla E_r.\end{array}$

This second set of equations is further split so that the contribution to the change in electron and ion specific energy densities due to the electron-ion coupling term is treated separately from the other contributions. An exponential electron-ion relaxation scheme is used to account for the contribution from the electron-ion coupling term in xRAGE simulations of the ICF test problems in this work. This exponential electron-ion relaxation scheme was developed for the FLASH radiation hydrodynamics code [28, 29] and is now briefly described. The scheme assumes that during a given time step from n to n + 1, the electron and ion-specific heats, C_Ve and C_Vi, remain constant. In terms of the specific heats and temperatures, the ion and electron energy equations can be written as

(4) $\begin{array}{c}\rho C_{Vi}\frac{\partial T_i}{\partial t}={\gamma }_{ei}\left(T_e-T_i\right),\rho C_{Ve}\frac{\partial T_e}{\partial t}={\gamma }_{ei}\left(T_i-T_e\right).\end{array}$

The solution of the discrete form of these equations in the time interval from n to n + 1 is facilitated by defining two new variables:

(5) $\begin{array}{c}\Delta =T_e-T_i,\sum =\frac{\left(C_{Ve}{T}_e+C_{Vi}{T}_i\right)}{\left(C_{Ve}+C_{Vi}\right)}.\end{array}$

Using these variables, the electron and ion temperatures are advanced from n to n + 1 using the following difference equations with $\delta t=t^{n+1}-t^n$ :

(6) $\begin{array}{c}{\Delta }^{n+1}={\Delta }^n\exp \left(-\frac{{\gamma }_{ei}^n}{{\rho }^n}\left(\frac{1}{C_{Ve}^n}-\frac{1}{C_{Vi}^n}\right)\delta t\right),{\sum }^{n+1}={\sum }^n.\end{array}$

The electron temperature at n + 1 is determined from equation 6 as follows, and the corresponding ion temperature is then found using the definition of Σ:

(7) $\begin{array}{}{T}_e^{n+1}=\frac{C_{Vi}^n}{C_{Ve}^n+C_{Vi}^n}\left({\Delta }^{n+1}+{\Sigma }^n\right).\end{array}$

3. Isentropic Spherical Compression

The isentropic spherical compression problem is the simplest case in the present verification and validation test suite, and this problem involves pure hydrodynamics with an ideal gas EOS model. The isentropic spherical compression problem is a variant of the Kidder ball problem with an initial Gaussian density profile. The isentropic spherical compression problem involves an ideal gas with γ = 5/3 and specific heat of $3.6\times {10}^{11}erg/\left(geV\right)$ , having the following initial conditions: a linear velocity field, a Gaussian density profile, and a constant specific internal energy. Ramsey et al. [30] derived a self-similar solution for this problem that has a spatially-uniform but time-varying specific internal energy. The boundary conditions include three planes of symmetry encompassing the origin and an outer spherical boundary that moves inward with the prescribed initial velocity above, directed towards the origin. This problem is simulated with the outer boundary internal to the mesh and with a separate material region outside that boundary. The inner material is an ideal gas consisting of a 3 cm radius sphere where the initial conditions are prescribed as discussed above. The outer material or background region has the same ideal gas equation of state, but with zero velocity, a temperature of 1 eV and a pressure of $1.96\times {10}^{9}erg/c{m}^3$ . The inner gas region collapses onto the origin without shocks.

The goal is to assess the accuracy of the xRAGE solution over the inner material region to a radius of 2 cm. Density values from a series of fixed tracer particles at different solid angles within a radius of 2 cm from the origin and at a simulation time of 0.5 μ s are shown in Figure 1. Both 2D and 3D xRAGE simulation results using the default method with VOF are presented, revealing excellent agreement between the analytic solution and the xRAGE results. A uniform mesh with spatial resolution of 0.01 cm was used for these simulations. The simulations were initialized using a link file, containing 288,000 triangular elements that represent the initial state of the collapsing gas region. 2D simulation results using the default method without VOF and the unsplit method are compared in Figure 2. There is good agreement between these two hydrodynamic methods, and even without any kind of interface treatment, these methods show remarkably good symmetry with relatively little scatter in density at any given radial distance, across seven different solid angles. These simulations have a uniform spatial resolution of 0.00125 cm.

Figure 1: 2D and 3D density results from xRAGE for the isentropic spherical compression problem at t = 0.5 μ s, compared against the analytic solution. These simulations used the default hydrodynamic method with VOF at 0.01 cm spatial resolution.

Figure 2: 2D density results from xRAGE for the isentropic spherical compression problem at t = 0.5 μ s, using the default method without VOF (red) and the unsplit method (blue). These simulations have a uniform spatial resolution of 0.00125 cm. For each value of radial distance, simulation results at seven distinct solid angles are presented.

4. Verney Shell Collapse

The Verney test problem [31] exercises material strength and involves an isentropic imploding steel shell with an initial inner radius of 8 cm and a thickness of 0.5 cm. The initial velocity spatial profile is chosen to create a constant density implosion that converts the initial kinetic energy into internal energy. A simple strength model consisting of a constant shear modulus of 0.895 Mbar and a constant yield stress of 0.05 Mbar is used together with the SESAME 4272 EOS for steel. The background air is represented by the SESAME 5030 EOS. Similar to the isentropic compression problem, a link file with triangular elements, representing the initial state of the steel shell, is used to setup the xRAGE simulations.

As the shell collapses kinetic energy is converted into internal energy through the plastic deformation of the shell. The internal energy of the steel shell as a function of time was derived using an analytical solution approach [31]. Both 2D and 3D xRAGE results using the default method with VOF at a spatial resolution of 0.01 cm are shown in Figure 3. There is excellent agreement between the 2D and 3D xRAGE simulations and between the xRAGE simulations and a 3D pure Lagrangian simulation using FUEL for the total internal energy of the shell as a function of time. FUEL implements a finite element staggered grid method for Lagrangian hydrodynamics [32]. At stagnation time, which occurs near 55 μ s, the calculated internal energy using xRAGE is virtually identical to the analytical solution. The 3D xRAGE simulation contained a uniform mesh with approximately 729 million zones. The VOF method for 3D simulations in this work incorporates the “onionskin” approach for ordering materials in mixed cells. This method was implemented within xRAGE for 3D simulations, following the approach used in Pagosa, a separate 3D Eulerian hydrodynamics code [33].

Figure 3: Calculated total internal energy of the Verney shell versus time with 2D xRAGE results is shown in red, 3D xRAGE results in green, and 3D Lagrangian results from FUEL in blue. The analytic solution is shown as the solid black line. The xRAGE simulations used the default hydrodynamic method with VOF at 0.01 cm spatial resolution.

Contours of velocity magnitude and density from the 2D xRAGE simulation are presented in Figure 4, showing a high degree of spatial symmetry, both at the initial time and at the final time of 55 μ s. Using a link file to setup the initial geometry is essential in order to obtain this level of symmetry. The link file represents the shell with triangular elements, generated from a conformal mesh of the shell with 1/64 degree angular resolution. This amount of angular resolution is particularly important in order to represent the initial velocity profile within the shell.

Figure 4: Density and velocity profiles from the 2D xRAGE simulations, using the default method with VOF at a spatial resolution of 0.01 cm.

5. 5-Material Compression Problem

The Verney shell problem involves an isentropic, smooth collapse, whereas the 5-material compression problem has a more complex collapse, involving contact discontinuities, release waves, and shocks. The 5-material problem has two steel cylinders that are separated by a high pressure gas region as shown in Figure 5. The steel cylinders are driven by the high pressure gas [31]. The outer steel cylinder has a thickness of 0.25 cm and the inner shell is 0.5 cm thick, with both cylinders having a length of 5 cm long. The high pressure gas region between the cylinders provides the energy to compress the inner steel cylinder and to expand the outer steel cylinder. Contained within the inner cylinder is a low pressure gas region. Both gas regions were represented using the SESAME 5030 EOS in xRAGE simulations and were initialized according to the conditions specified in Figure 5. Both steel cylinders were represented by the SESAME 4272 EOS with the elastic, perfectly-plastic constitutive model from Section 4. This test problem is simulated to a time of 5 μ s.

Figure 5: The geometry of the 5-material compression problem.

2D and 3D xRAGE simulations of this problem were performed using the default method with VOF on an AMR mesh with finest resolution of 0.00625 cm within the steel shells. Unlike other simulations presented in this work, pressure-temperature equilibrium was not enforced in mixed cells, such that each material has a distinct temperature. Both 2D and 3D simulations give a total energy (internal plus kinetic) of about 2.7 MJ for the inner steel cylinder at 5 μ s. The internal energy calculated during the compression of the inner steel cylinder is presented in Figure 6. There is excellent agreement between 2D and 3D xRAGE simulations for the internal energy as a function of time. Both xRAGE simulations agree fairly well with a 2D axially symmetric Lagrangian calculation using FLAG, a Lagrangrian hydrodynamics code that uses a compatible hydrodynamic algorithm which conserves total energy [34].

Figure 6: 2D and 3D xRAGE simulation results for the internal energy of the inner steel shell, compared with 2D Lagrangian simulation results from FLAG.

6. Direct Drive ICF Problem

VOF is essential for problems with material strength, at low temperatures, in order to control artificial diffusion of materials in mixed cells; however, this is much less of an issue at high temperatures such as those encountered in ICF. The default and unsplit methods without VOF have been applied to simulate two ICF test problems. The first ICF test problem discussed here is based on a deuterium-filled capsule experiment, OMEGA 50997, as described by Dodd et al. [35]. This test problem involves a SiO₂ glass shell that is filled with a mixture of D₂ and ³He, where 13.5 kJ of energy is deposited within the SiO₂ shell during 1 ns, approximating the direct drive laser absorption process in the 50997 capsule experiment. In simulations of this test problem, the SiO₂ glass shell is represented by the SESAME 7383 EOS.

When coupled with radiation diffusion, 3T physics, and thermal conduction, the default and unsplit hydrodynamic methods give similar burn-averaged ion temperatures, both when the hot spot forms and when the fusion reaction rate is maximum (“bang time”), as shown in Figure 7. The burn-averaged ion temperature is calculated by taking a weighted average of the ion temperature across the entire mesh, using available deuterium fusion reaction rate data to find the weight for each zone [36]. Highly resolved 1D simulation results are presented at a uniform spatial resolution of 0.125 μ m. After 0.8 ns, there is a noticeable difference in burn-averaged ion temperature between the two methods, with the unsplit method producing a temperature that is 0.8 keV higher at 1.2 ns. This difference in temperature translates into a difference of about 13% in the total number of neutrons produced from DD fusion reactions.

Figure 7: An ICF test problem based on an OMEGA capsule experiment: (a) the initial geometry of the capsule and (b) the calculated burn-averaged ion temperature from separate highly resolved 1D simulations using the default and unsplit hydrodynamic methods in xRAGE.

In addition to 1D simulations, 3D simulations of this problem have been completed using a 2D-linked-into-3D approach where 2D axially symmetric xRAGE simulation results at 0.5 ns are used to initialize a 3D simulation on a truncated spatial domain consisting of a 600 μ m by 1200 μ m by 600 m sized box region, including both poles of the capsule which are on the y-axis and the-y-axis, respectively. The x-y and y-z planes are symmetry planes in the 3D simulation. The mesh for the 2D simulation consists of 4000 by 8000 square zones, with each zone having a uniform edge size of 0.25 μ m. The corresponding 3D mesh includes a base mesh of 300 by 600 by 300 cubic zones, with each zone having a uniform edge size of 2 μ m. In addition to the base mesh, there are 2 spherical regions where the mesh is further refined. The first region of mesh refinement is a sphere with a radius of 256 μ m, centered at the origin. Within this region, the mesh resolution is refined so that zones have a uniform edge size of 1 μ m. Nested within the first mesh refinement region is an additional region of mesh refinement, consisting of a sphere with a radius of 128 μ m, centered at the origin. Within this region, zones have a uniform edge size of 0.5 μ m. There are approximately 90 million zones on the 3D mesh. 2D simulation results are coarsened from 0.25 μ m to 2 μ m in order to enable mapping onto the 3D mesh. The link time occurs before the formation of a hot spot at the center of the capsule at about 0.6 ns. For 3D simulations, the inner and outer surfaces of the shell include a complete spectrum of roughness, determined from experimental surface measurements of CH cryogenic capsules reported by Haan et al. [37] and presented here in Figure 8. The CH inner and CH outer radial perturbations, presented in Figure 8(b) are used here for 3D simulations of the SiO₂ shell in the direct drive ICF test problem.

Figure 8: Surface roughness for the direct drive ICF test problem: (a) the 1D power spectrum of surface roughness as a function of Fourier mode number for generic cryogenic ICF capsules as reported by Haan et al. [37] and (b) the corresponding radial surface perturbations, used in xRAGE simulations, as functions of polar angle. The radial surface perturbations in (b) were determined by summing the power spectrum from (a) over all mode numbers.

In general, the calculated vorticity contours are qualitatively similar for both the default and unsplit methods as shown in Figure 9, but the unsplit method generates significantly more vorticity from the hot spot and also more vorticity at the interface between the shell and the gas, especially at 0.8 ns. For both methods, a region of vorticity created by the hot spot is evident, consisting of a series of compact azimuthal vortex rings, and this region moves along the-y-axis during the course of the simulation, in a direction corresponding to the polar angle of zero degrees in Figure 8(b), where the initial surface perturbations on the capsule surface are the largest. The burn-averaged ion temperatures for both the default and unsplit methods are presented in Figure 10 and unlike the corresponding 1D simulation results in Figure 7(b), where the two temperature profiles are virtually identical at bang time, for the 3D simulations, the unsplit method produces a temperature that is about 0.6 keV less than the default method at bang time which occurs at 0.72 ns for the 3D simulations. A possible explanation for the lower temperature observed with the unsplit method is that the greater vorticity produced by the unsplit method leads to more mixing of the shell material into the gas region, and this additional shell material lowers the temperature of the gas.

Figure 9: Contours of vorticity magnitude for the direct drive ICF test problem using both the default and unsplit methods at four different times.

Figure 10: Burn-averaged ion temperature profiles from 3D simulations of the direct drive ICF test problem.

Table 1 presents calculated values for peak burn-average ion temperature at bang time and total number of DD neutrons produced. The DD neutron production is determined using species number densities and the ion temperature from the xRAGE simulations, together with available DD fusion reaction rate data [36]. The calculated neutron production value from 1D or 3D simulations using either the default or unsplit method is less than half of the value observed in the 50997 experiment. This reflects the limitations of the simple energy source term used in the simulations; a more accurate model of the laser energy absorption is needed for better agreement. A 1D simulation using the unsplit method produces a burn-averaged ion temperature of about 7.6 keV at bang time and about $1.8\times {10}^{10}$ neutrons from DD fusion reactions. In 3D, the unsplit method produces a burn-averaged ion temperature at bang time of 7.2 keV and about $1.4\times {10}^{10}$ neutrons. Additional sources of mixing at the shell and gas interface that are present in the 3D simulation produce a smaller value of burn-averaged ion temperature and less neutrons from DD fusion reactions.

Table 1: Comparison of ion temperature and total DD neutron production for different simulations of the direct drive ICF test problem.

While mixing between the shell material and the gas occurs in both 1D and 3D simulations, two specific mechanisms for this mixing are clearly identified in the 3D simulations: (1) an intense region or sheet of vorticity along the entire gas/shell interface and (2) the development of several distinct azimuthal vortex rings at discrete locations on the gas/shell interface. The later mechanism is clearly illustrated in Figure 11, which compares isosurfaces where the Q-criterion is ${10}^{18}{\text{s}}^{-2}$ . Q-criterion is a useful way to visualize regions of vorticity [38]. In order to create isosurfaces of Q-criterion, the simulated velocity values within a 256 μ m by 512 μ m by 256 μ m sized box region on the computational mesh, centered at the origin, are mapped from the block AMR computational mesh onto a uniform mesh with spatial resolution of 0.5 μ m. Velocity gradients and corresponding values of Q-criterion are calculated on this uniform mesh. The axis of cylindrical symmetry is along the y axis, and the x-y plane and the y-z plane are symmetry planes. At 0.7 ns, there is a region of hot spot vorticity close to the origin, and as one moves away from the origin, there is an outgoing shock, which is contained entirely within the gas. Beyond, the outgoing shock is the SiO₂ shell, and several vortical structures are present within the SiO₂ shell. The unsplit method produces a series of azimuthal vortex rings within the shell, centered around the-y-axis, whereas the default method has single vortex ring in the same location. In addition, with both the unsplit and default methods, additional vortical structures appear at 0.7 ns, roughly halfway between the x-y and y-z symmetry planes, as indicated in Figure 11 by the arrows pointing toward these regions of vorticity within the shell. These vortical structures are three-dimensional in nature, and are no longer present at 0.8 ns. At later times, the outgoing shock expands into the SiO₂ shell, and a series of azimuthal vortex rings is present at the interface between the gas and the SiO₂ shell. The unsplit method develops azimuthal vortex rings at several locations on the gas/shell interface, whereas the default method produces a single azimuthal vortex ring, which is most evident at 0.9 ns as illustrated in Figure 11. Azimuthal vortex rings of this kind have been observed in previous simulations of ICF capsules [39]. In addition to azimuthal vortex rings, the unsplit simulation also produces polar vortex rings, which are clearly visible at 0.8 ns and 0.9 ns as indicated in Figure 11. These polar vortex rings are unique to the unsplit simulation. The default simulation produces vortex tubes at a few discrete polar angles which do not develop into coherent rings.

Figure 11: Isosurfaces where the Q-criterion is equal to ${10}^{18}{\text{s}}^{-2}$ at four different times within a 256 μm by 512 μm by 256 μm sized box region for both the default and unsplit methods.

7. Indirect Drive ICF Problem

Both the default and unsplit methods have also been applied to simulate a more challenging indirect drive ICF test problem requiring additional physics including multigroup radiation diffusion, to transport X-ray energy from the cylindrical hohlraum to the target capsule. The capsule in this test problem is based on a NIF cryogenic capsule experiment, N170601, as described by Le Pape et al. [40], involving a high density carbon (HDC) shell that is 70 μ m thick with an outer radius of 980 μ m. The HDC shell surrounds a cryogenic deuterium-tritium (DT) layer with mass of 0.13 mg, and DT gas fills the center of the capsule. There are 67 energy fluxes which are specified on the spatial boundary as functions of time, and these were formulated in previous studies using a separate radiation hydrodynamics code, HYDRA, in order to match certain tuning data collected from dedicated NIF shots related to N170601 such as early-time shock propagation through surrogate targets, inflight shell implosion velocity, low-mode shell shape, and bang time [41, 42]. These boundary energy fluxes are applied here for xRAGE simulations without modification.

A hot spot forms at the center of the capsule at about 7.7 ns, and the peak DT fusion reaction rate occurs at 8.33 ns, as indicated in Figure 12, which shows 1D simulation results at a uniform spatial resolution of 0.25 μ m. The burn-averaged ion temperature is calculated by taking a weighted average of the ion temperature across the entire mesh after the simulation is completed, using available DT fusion reaction rate data to find the weight for each zone [36]. Fusion energy is not deposited in the DT gas or cryogenic layer during the simulation. When coupled with multigroup radiation diffusion, 3T physics, and thermal conduction, the default and unsplit hydrodynamic methods give virtually identical burn-averaged ion temperatures as a function of time during hot spot formation, and give similar spatial profiles of ion temperature, density, and mass concentration of the cryogenic layer at 8.3 ns.

Figure 12: Results from highly resolved 1D simulations of the indirect drive ICF problem at 0.25 μm uniform spatial resolution: (a) burn-averaged ion temperature during formation of the hot spot and (b) spatial profiles at 8.3 ns of density normalized by 380 g/cm³, ion temperature normalized by 8 keV, and cryogenic layer mass concentration.

While 1D simulations are useful for evaluating the default and unsplit methods together with multigroup radiation diffusion and 3T physics, such simulations are not adequate to study the development of vortex rings within the capsule and the transition to turbulence. At a given spatial resolution, the unsplit method is expected to generate more vorticity, compared to the default method, based on prior 3D simulations of the Taylor-Green vortex in the literature [2], as well as the 3D simulation results for the direct drive capsule problem from the previous section. 3D simulations of this indirect drive ICF problem have been completed using the same 2D-linked-into-3D approach as the direct drive problem, where a 2D simulation is mapped onto a spatial domain consisting a 800 μ m by 1600 μ m by 800 μ m sized box region, including both poles of the capsule which are on the y-axis and the-y-axis, respectively. The x-y and y-z planes are symmetry planes in the 3D simulation. The 3D computational methodology begins with a 2D axially symmetric simulation, including multigroup radiation diffusion with 67 energy groups, representing the energy flux from the hohlraum in the experiment. 1D and 2D simulations use multigroup diffusion to approximate the radiation transport between the hohlraum and the target capsule, utilizing 67 energy flux terms on the computational boundary. Low-density, optically thin, ablated plasma from the capsule is present between the hohlraum and the capsule, and these conditions require corrections to multigroup diffusion. The separate HYDRA simulations which determined the 67 energy flux terms also used multigroup diffusion but implemented two additional corrections together with multigroup diffusion, a mean free path to boundary correction and disabling the flux limiter, in order to achieve reasonable agreement with HYDRA simulations of the shell trajectory and shock timing experiments using radiation transport [41, 42]. The boundary energy flux functions are identical for both the 1D and 2D xRAGE multigroup radiation diffusion simulations. However, the size and shape of the spatial boundary are different for the 1D and 2D simulations. In 1D, the spatial boundary occurs at a spherical radius of 5 mm, but in 2D, the spatial boundary in the x-y plane consists of a 5 mm by 10 mm box, with a similar boundary for the y-z plane. The spatial boundary of the 2D simulation here differs from the corresponding shape of the hohlraum in N170601 [40], which is a cylinder that has a radius of 3.1 mm, as compared with 5 mm in the simulation, and a length of 11.3 mm, as compared with 10 mm in the simulation. A larger spatial domain was chosen for the indirect drive test problem, relative to the experiment, in order to allow for more expansion of the ablated plasma from the capsule, before this material impacts the computational boundary. Also, for the 2D simulation, which has a uniform spatial resolution of 0.5 μ m, both the inner and outer surfaces of the ablator, as well as the inner surface of the cryogenic DT ice layer, include a complete spectrum of roughness, represented by the radial perturbations presented in Figure 8(b). At 8 ns, the 2D simulation results are mapped onto a 3D block AMR mesh with 0.5 μ m finest resolution within the cryogenic DT layer and the central DT gas region. There are about 70 million zones on the 3D block AMR mesh. The 3D simulation uses grey radiation diffusion and runs until 9 ns, with bang time occurring at 8.36 ns.

The indirect drive problem uses a 2D mesh with 10000 by 20000 square zones, with each zone having a uniform edge size of 0.5 μ m. The 3D mesh has a block AMR structure similar to the 3D mesh for the direct drive problem. There is a base mesh consisting of 200 by 400 by 200 cubic zones, with each zone having an edge size of 4 μ m. There are 3 spherical regions where the mesh is refined further. The first mesh refinement region is a sphere with a radius of 512 μ m, centered at the origin. Zones within this region have an edge size of 2 μ m. Nested within this region is a second mesh refinement region, which is a sphere with a radius of 256 μ m, centered at the origin. Zones within this region have an edge size of 1 μ m. Nested within the second mesh refinement region is a third mesh refinement region, consisting of a sphere with a radius of 128 μ m, centered at the origin. Zones within this region have an edge size of 0.5 μ m.

In 3D simulations at 8.3 ns, one can see in Figure 13 that both the default and unsplit methods produce compact vortex rings from the hot spot, which move along the-y-axis, in a direction corresponding to the polar angle of zero degrees in Figure 8(b), from the origin to the ice/carbon interface as the simulation progresses. In addition to these regions of vorticity on the y-axis, the default method produces additional vortex rings from the hot spot, which are located at polar angles of 45 degree and 135 degrees, near the ice/carbon interface. Arguably, these additional vortex rings arise in part from numerical artifacts due to the directionally split nature of the default hydrodynamic method. The unsplit method does not have these features. Later in time at 8.5 ns, the unsplit method generates significantly more vorticity at the ice/carbon interface. In addition, the unsplit method produces more vorticity from the hot spot which moves along the-y-axis, penetrates the ice/carbon interface, and then overtakes the outgoing shock in the carbon shell. The behavior of the hot spot vorticity together with the larger amount of vorticity at other locations along the ice/carbon interface leads to greater mixing between the ice and the carbon regions, relative to the default method. This behavior is similar to what was observed for the direct drive ICF test problem in the previous section.

Figure 13: Contours of vorticity magnitude for the indirect drive ICF test problem using both the default and unsplit methods at four different times.

The burn-averaged ion temperature profiles from the 3D simulations are presented in Figure 14. The peak temperature near bang time for the default method is about 5.7 keV, whereas it is about 5.9 keV for the unsplit method. These peak temperatures are slightly larger than what is observed in corresponding 1D simulations at higher resolution, where the peak temperature is about 5.5 keV for both the default and unsplit methods. In 1D, the temporal profiles of temperature from the two methods closely follow each other as shown in Figure 12(a). For the 3D simulations, the unsplit method produces slightly higher temperatures leading up to bang time, but then the two temperature profiles come into close agreement after 8.4 ns.

Figure 14: Burn-averaged ion temperature profiles from 3D simulations of the indirect drive ICF test problem.

Table 2 compares 1D and 3D simulated peak burn-averaged ion temperature near bang time, and the corresponding total number of neutrons produced from DT fusion reactions, which is calculated using values of species number densities and ion temperature from the xRAGE simulation together with available DT fusion reaction rate data [36]. The 1D results presented in Figure 12 used multigroup diffusion for the entire duration of the simulation. In order to facilitate comparison with the 3D simulations, the 1D results presented in Table 2 follow a similar computational approach as the 3D simulations, where multigroup radiation diffusion is turned off at 8 ns and grey diffusion is then used. Using this computational approach, both 1D simulations produce about $2.2\times {10}^{16}$ neutrons from DT fusion reactions, which is considerably more than what was observed in the N170601 experiment. The neutron production is only $1.3\times {10}^{16}$ in 1D simulations where multigroup diffusion is active after 8 ns. The peak burn-averaged ion temperature near bang time in 1D simulations where grey diffusion is used after 8 ns is about 5.9 keV, which is 5% higher than in corresponding simulations with multigroup diffusion. In 3D, the default and unsplit methods are in fairly good agreement with each other, with the unsplit method producing larger values of peak burn-averaged ion temperature and neutron production. Both the default and unsplit methods produce burn-averaged ion temperatures near bang time that are about 20% larger than the experimental value. The calculated value of burn-averaged ion temperature from a 1D simulation agrees well with the corresponding value from a 3D simulation when using the unsplit method for both simulations, as indicated in Table 2. However, the 3D simulation produces a smaller number of neutrons due to additional mixing at the ice/carbon interface. In 3D vortex, rings at the ice/carbon interface mix material from the carbon shell into the DT ice region, decreasing the overall number of DT fusion reactions. The test problem here is a simplified variant of N170601 in which the energy of alpha particles generated by DT fusion reactions is not deposited within the computational domain. In order to better understand how heating from alpha particles would affect the results in Table 2. a separate 1D simulation was completed that deposited the energy of alpha particle generated by fusion reactions uniformly within the DT ice layer. This simulation produced about 30% more neutrons than corresponding simulations without energy deposition from alpha particles.

Table 2: Comparison of ion temperature and total DT neutron production from different simulations of the indirect drive ICF problem. 1D simulations use grey diffusion after 8 ns.

While the vorticity production may explain why 3D simulations produce less neutrons than corresponding 1D simulations, there is no evidence that transition to turbulence occurs within the DT ice and DT gas regions. Vortex rings are generated at the ice/carbon interface, but these remain largely coherent from 8 ns to 9 ns as shown in Figure 15. The density of the DT ice is shown in Figure 16 at both 8.3 ns and 9 ns. Rayleigh-Taylor instabilities are clearly visible at the interface of the DT ice and the carbon shell. But these features are coherent and not turbulent. If transition to turbulence does take place in these simulations, it occurs within the carbon shell of the capsule and not within the DT ice or DT gas regions.

Figure 15: Isosurfaces where the Q-criterion is equal to ${10}^{18}{\text{s}}^{-2}$ at four different times within a 256 μm by 512 μm by 256 μm sized box region for the indirect drive ICF test problem using both the default and unsplit methods. The spatial resolution within the box region is 0.5 μm.

Figure 16: Contours of DT ice density for the indirect drive ICF test problem using both the default and unsplit methods. The top row shows results at 8.3 ns and the bottom row shows results at 9 ns.

8. Conclusions

Recent 3D xRAGE simulations of a set of verification and validation test problems have produced promising results. Among these problems are the Kidder ball problem, the Verney shell problem, and a 5-material compression problem. There is excellent agreement between 2D and 3D xRAGE simulation results and between the xRAGE results and the benchmark solutions for all three of these problems. In addition to using the default method with VOF, several other less common simulation options were explored which enabled the xRAGE simulations to closely match the benchmark solutions. For both the Kidder and Verney problems, it was advantageous to start the simulation using a link file with triangular elements that are mapped onto a uniform mesh in xRAGE. Furthermore, it was essential not to enforce temperature-pressure equilibrium in mixed cells, so that each material has a distinct temperature, when simulating the 5-material problem in order to closely match the corresponding Lagrangian simulation result.

ICF experiments are uniquely well suited for a verification and validation suite because they involve surface perturbations, convergence, and instability growth, all in a high temperature regime. A strong focus of this work is to compare the default and unsplit hydrodynamic methods for ICF test problems, coupling hydrodynamics with 3T physics, thermal conduction, and radiation diffusion. A 3D ICF verification problem was created, based on an OMEGA direct drive capsule experiment. In 1D simulations, the default and unsplit hydrodynamic methods produce nearly identical burn-averaged ion temperatures when the hot spot forms and up until bang time. The peak ion temperature near bang time in the simulations is similar to what is observed in the 50997 capsule experiment. In 3D simulations, the unsplit method generates more vorticity than the default method, causing additional mixing between the gas and the shell, resulting in a lower peak burn-averaged ion temperature near bang time that is in better agreement with the experiment. Vortex rings are generated within both the shell and the gas regions, and these rings remain coherent during the simulations.

A more challenging indirect drive test problem, using multigroup radiation diffusion, was created based on a cryogenic DT capsule experiment on NIF. In highly resolved 1D simulations, the default and unsplit hydrodynamic methods give virtually identical burn-averaged ion temperatures as a function of time during hot spot formation and give similar spatial profiles of ion temperature, density, and mass concentration of the cryogenic layer near bang time. In 3D simulations, both the default and unsplit methods produce compact vortex rings from the hot spot, which move along the-y-axis, from the origin to the ice/carbon interface as the simulation progresses. In addition, a sheet of vorticity is generated at the ice/carbon interface, which is more intense using the unsplit method. Vorticity production mixes material from the carbon shell into the DT ice and DT gas layers. Vorticity production may explain why the 3D simulations produce less neutrons than corresponding 1D simulations. When initialized from a 2D multigroup simulation with uniform 0.5 μ m spatial resolution, the 3D simulations are in reasonable agreement with the experimental values of ion temperature and DT fusion neutron production. There is no evidence of transition to turbulence within the DT ice and DT gas regions in these simulations, as the vortex rings at the ice/carbon interface remain largely coherent.

Future work should explore the sensitivity of the 3D simulation results for both the direct drive and indirect drive problems to additional features such as the fill tube. In addition, future simulations should also attempt to induce transition to turbulence within the central gas region in these two problems by introducing 3D spatial perturbations when mapping from 2D to 3D, following the approach used by Haines et al. [15]. The direct drive problem can be improved further by including laser energy deposition, which would be more realistic than the present energy source term. The two ICF problems presented here would be useful test cases for evaluating plasma viscosity and transport models that have been discussed in the literature.