2. Numerical Methods

The governing equations of the 3T electron-ion radiation hydrodynamics conservation laws in Eulerian form [Reference Mihalas and Mihalas22] 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. [Reference Brown, Preston and Singletonjr23]. The ion and electron thermal conductivities are calculated using the approach of Lee and More [Reference Lee and More24]. 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 [Reference Colgan, Kilcrease and Magee25] are used for the simulations, as determined by the TOPS code [Reference Abdallah and Clark26]. 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 [Reference Lyon27].

The governing equations presented above are solved by using the technique of operating splitting, wherein two sets of simplified equations are derived [Reference Gittings, Weaver and Clover1]. 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 [Reference Fatenejad, Fryxell and Wohlbier28, Reference Fryxell, Olson and Ricker29] 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. [Reference Ramsey, Brown, Nelson and Alme30] 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 [Reference Chiravalle, Barlow and Morgan31] 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 [Reference Chiravalle, Barlow and Morgan31]. 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 [Reference Chiravalle and Morgan32]. 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 [Reference Weseloh, Clancy and Painter33].

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 [Reference Chiravalle, Barlow and Morgan31]. 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 [Reference Caramana, Burton, Shashkov and Whalen34].

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. [Reference Dodd, Benage and Kyrala35]. 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 [Reference McNally, Rothe and Sharp36]. 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. [Reference Haan, Lindl and Callahan37] 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. [Reference Haan, Lindl and Callahan37] 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 [Reference McNally, Rothe and Sharp36]. 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 [Reference Hunt, Wray and Moin38]. 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 [Reference Thomas and Kares39]. 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. [Reference Le Pape, Berzak Hopkins and Divol40], 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 [Reference Clark, Weber and Milovich41, Reference Haines, Clark, Weber, Edwards, Batha and Kline42]. 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 [Reference McNally, Rothe and Sharp36]. 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 [Reference Grinstein, Saenz, Dolence, Masser, Rauenzahn and Francois2], 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 [Reference Clark, Weber and Milovich41, Reference Haines, Clark, Weber, Edwards, Batha and Kline42]. 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 [Reference Le Pape, Berzak Hopkins and Divol40], 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 [Reference McNally, Rothe and Sharp36]. 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.