Atomic scale mixing for inertial confinement fusion associated hydro instabilities

Abstract

Hydro instabilities have been identified as a potential cause of performance degradation in inertial confinement fusion (ICF) experiments. We study instabilities associated with a single Richtmyer–Meshkov (RM) interface in a circular geometry, idealized from an ICF geometry. In an ICF application, atomic level mix, as an input to nuclear burn, is an important, but difficult to compute, variable. We find numerical convergence for this important quantity, in a purely hydro study, with only a mild dependence on the Reynolds number of the flow, in the high Reynolds number limit. We also find that mixing properties show a strong sensitivity to turbulent transport parameters; this sensitivity translates into an algorithmic dependence and a nonuniqueness of solutions for nominally converged solutions. It is thus a complication to any verification and validation program. To resolve the nonuniqueness of the solution, we propose a validation program with an extrapolation component, linking turbulent transport quantities in experimental regimes to mildly perturbed turbulent transport values in ICF Reynolds number regimes. In view of the observed solution nonuniqueness, the validation program and its justification from the results presented here, has a fundamental significance.

Introduction

Hydro instabilities are known to limit inertial confinement fusion (ICF) neutron production [1]. See also the more recent discussions in Ref. [2]. For predictions of neutron production, the atomic level description of temperature and mixture concentrations provide the main input to a thermonuclear (TN) burn computation [3]. The interest in mix is brought about by current experiments of ICF processes. For instance, there is an observed degradation of 40–75% in neutron production, attributed to atomic mixing, estimated [4] based on experiments conducted on the Omega laser. In addition, it is proposed that low mixing is required to achieve ignition, specifically, the required maximum level of contaminated fuel from the ablator material must be less than an estimated 25%–40% [1]. Modal growth factors of 10–100 or more are reported in the linearized analysis of [1], suggesting a role for the nonlinear mode coupling and even perhaps the fully developed turbulent mixing regimes, beyond the ablation stabilized weakly nonlinear theory reviewed in Ref. [5]. ICF simulations do not predict correctly the observed dynamics of ICF capsules, which exhibit significant yield degradation relative to design. Thus a re-examination of the validation process would appear to be appropriate. Experimental analysis [6] of the ICF hot spot mix is based on Ge doped ablator material found within the hot spot at ignition time. Only mix at a level as predicted by design analysis was observed.

Known ICF instabilities or asymmetries occur primarily in four distinct stages during an ICF process. They begin with (1) the laser drive and are followed by (2) instability growth at the ablation surface. In addition, (3) Richtmyer–Meshkov (RM) instabilities occur at the boundary of ablator to deuterium-tritium (DT) ice and the boundary of DT ice to DT gas. These propagate inward to define initial conditions for the (4) deceleration phase, which is known to be strongly RT unstable. The second of these instabilities is an ablation modified RT instability, the third is RM and the fourth is RT in nature. In Ref. [7], velocity perturbations of 1% at the beginning of the deceleration phase (4) are considered to be dangerous for the performance of the capsule. The Verification and Validation (V&V) methods of this paper are applicable to the second, third and fourth instabilities mentioned above. Inaccuracies in prediction of the yield cliff suggests a possible need for improved modeling of the ablation phase, such as the inclusion of subgrid scale (SGS) models as proposed in Ref. [2], and the value in reexamination of any of the four instability stages, as well as other issues not discussed here.

Here we address V&V aspects of a hydro instability study related to these broad ICF concerns. We study the RM instability, which is a primary ICF hydro instability, in its relation to atomic scale mix; the conclusions have implications for RT V&V as well.

The main thrust of the paper is to offer scientific data in support of a proposal for validation with extrapolation. We use a single code with no free parameters, validated in common for RT and RM data. This goal, infeasible for Reynolds Averaged Navier–Stokes (RANS), is a natural objective for Large Eddy Simulations (LES), which use dynamic SGS models. Validation for RT was previously considered and is summarized in Sec. 1.4. The RT case is more sensitive to transport (molecular, turbulent or numerical) than are RM instabilities. RT is thus more demanding for validation. This difficulty can be inferred from the fact that (in the authors' opinion) the FT/LES/SGS code discussed here appears to be unique as a multimode RT validated compressible LES code. For RM, we rely on an extensive code comparison study [8] with RAGE [9]. Generally excellent agreement was found for most measures of comparison, however, exceptions were noted for the temperature field (resolved in favor of our code FronTier (FT)). In the case of high density ratios, such as those that occur at the hot/cold DT interface in an ICF capsule, the edges of the mixing zone showed slow numerical convergence, especially for RAGE, indicating potential problems for commonly used resolution levels. Translated into the ICF context, this raises the question of simulation velocity perturbations at the beginning of the deceleration phase and perhaps of penetration of cold DT gas from the original DT ice layer of an ICF capsule into the ICF hot spot, or of ablator material (as has been observed experimentally), any of which can potentially degrade performance. Here we are mainly concerned with the extrapolation from the experimental regime to an ICF regime or to still higher Reynolds numbers. For the purpose of this paper, we understand the experimental regime to be Re ≈ 3.5 × 10⁴, typical of laboratory RT experiments and the ICF regime to be Re ≈ 6 × 10⁵, with higher values of Re also considered. We observe that the obvious step of RT validation using NIF or Omega laser data does not address the multimode, mode coupling RT growth stage, as the laser experiments performed to date have not reached this time regime.

The two main steps suggested here (RT validation and assessment of high Re extrapolation, with a zero parameter, theoretically constrained model) are not part of common V&V methodology. Most multiphysics codes fail for RT validation in the LES regime, and generally ICF codes do not employ dynamic subgrid scale models, which allow for smooth extrapolation from experimental to ICF Reynolds numbers. The validation proposed here is suggested as a complement to, or preliminary to, and not a replacement for, validation as presently performed. The most fundamental validation (as yet apparently still an open question) is comparison of simulations against ICF capsule performance directly.

The issue we raise is that the atomic mix properties are sensitive to the turbulent transport coefficients, and through them to the numerical algorithm. Stated more directly, apparently converged solutions are nonunique. They depend in an essential manner on details of the numerical algorithm. This unfortunate fact is documented here. However, we find only a mild variation in these coefficients in passing from an experimental regime to the ICF regime, as far as purely hydro issues are concerned. Thus a mild extrapolation from experimentally validated transport coefficients is feasible, yielding a significant constraint on the turbulent transport terms. We bridge this gap using algorithmic choices supported and constrained by theory and mathematical verification. We propose a program to carry out this objective, based on (a) a numerical algorithm with front tracking, to limit numerical thermal conduction and concentration diffusion and on (b) dynamic subgrid models to assure a theoretically consistent parameter free setting of the turbulent transport coefficients.

The context of the present study is a circular Richtmyer–Meshkov fluid instability problem; in view of the twin computational requirements of extreme levels of mesh refinement and of a parametric variation study involving multiple simulations, we consider this problem in 2D. At high Re, with ionized molecules, we consider a fluid with plasma like properties. These particles may be strongly coupled to a radiation field, while the thermal conduction between ions may be mediated by the ion-electron thermal coupling and the faster electron thermal diffusion. To model at a purely hydro level, we assume a highly heat conductive Prandtl number Pr = 10⁻⁴ and a Schmidt number Sc = 1. We define the high Re limit to be taken as viscosity ν → 0, with fixed values for Sc and Pr. The parameter study analyzes a series of mesh levels and Re values, up to Re = 6 × 10⁷ and Re = ∞, extending a study [10] of the same problem. In future work, we will include plasma and radiation effects.

Here, we explain a fundamental V&V challenge and present scientific data which supports a V&V route to bypass these obstacles. The challenge is the nonuniqueness and algorithmic dependence of nominally converged solutions for turbulent mixing problems. The route to bypass this obstacle is the combined use of parameter free dynamic SGS models and control over numerical diffusion through front tracking.

To validate this algorithmic choice, we use a two step V&V method. The main step is RT validation against the wealth of RT experimental data at Re ≈ 3.5 × 10⁴. This data is highly sensitive to transport (molecular, turbulent, numerical) to the point that the FT/LES/SGS code used here appears to be the unique compressible multiphysics LES code to have been validated in a zero parameter manner against this data. This first step was reported in earlier publications [11], [12] and is summarized here.

The second step in our V&V method is an extrapolation from Re ≈ 3.5 × 10⁴ to the ICF regime of Re ≈ 6 × 10⁵ or higher up to Re = ∞. This step leads to the main technical result of this paper. In essence, we show (a) that the consequences of the extrapolation are not large in their influence on numerical simulations, and (b) that there is an available theory to constrain the extrapolation. The theory [13], which is the basis of dynamic SGS models, is not new to this paper.

Thus our main results can be stated as (a) numerical convergence of the turbulent mixing statistics in question and (b) a small magnitude of the perturbation in these statistics due to a change of Re.

The numerical study is based on the front tracking code FronTier [14], [15], enhanced with an LES turbulence model having dynamic SGS terms [13]. We refer to this algorithm as FT/LES/SGS. These two algorithms are built upon a hydro package (compressible, with general equation of state) using the MUSCL [16] algorithm. To facilitate coupling of these features with other hydro packages, we have developed an API. To avoid complications with multispecies diffusion [17], we assume mixing of two fluids only.

With front tracking, our aim is to achieve enhanced resolution of interfaces between fluids and maintain this resolution throughout the simulation. This allows us to minimize the numerical mass diffusion associated with steep concentration or thermal gradients. We see improved results even for coarser grids, since the discontinuity that normally leads to the smearing in coarser grids is dynamically tracked, preventing one fluid from numerically interacting across the interface with the other.

We accomplish this goal by creating a separate, lower dimensional mesh along the front (interface) which separates the fluids, for example following an iso-temperature contour within a sharp thermal gradient. The front then separates fluids with distinct physical characteristics. During each time step we propagate the front's mesh to its new location based on the surrounding interior states. The front points are propagated in a normal direction, utilizing a Riemann problem formulation to resolve the different velocities obtained via one sided extrapolation from interior states to the front location.

When interfacial tangling occurs, we utilize algorithms to untangle and redistribute points on the front mesh, allowing for a logical reconstruction of a tangled interface. For more information on this process and the FronTier front tracking algorithm, we refer the reader to [15], [18], [19].

The front mesh is overlayed on a normal Eulerian grid. Where it intersects cells, we cut those cells into multiple pieces to ensure separation of the distinct fluids across an interface, illustrated in Fig. 1.

When updating interior states on a cut cell, ghost cells constructed from states defined in neighboring cells on the same side of the front allow for a 1-sided update, consisting only of states from those cells on the same side of the discontinuity interface.

This construction allows us to control numerical diffusion, especially with complicated interfaces, as arise in hydro instabilities in ICF processes.

The purpose of subgrid models is to capture the effect of the unresolved scales (below the grid level) on those that are resolved. We start with the interpretation that the numerical solution values represent grid cell averaged quantities. Conventionally, in the derivation of filtered Navier–Stokes equations for use in Large Eddy Simulations, there is a convolution average by a positive weight function (the filter), and the equations are derived for the filtered quantities. We omit this filter step, and use the cell averaged quantities directly. The average of the nonlinear terms introduces new unknowns into the equations for which new equations, called closures, are required. These show up as flux type terms, and are the origin of the turbulent transport. The closure term must be modeled; often a solution gradient is used for a turbulent flux, with an undetermined coefficient. With dynamic subgrid models [13], [20], [21], the otherwise missing coefficient is determined by the numerical solution itself, and varies with space-time according to local flow conditions. The determination is achieved by consideration of the closure on the current grid and on a once coarser grid, and by an asymptotic assumption that the model coefficient is a mesh convergent or asymptotically mesh independent quantity, so that a common value can be used between the two grid levels considered. The advantage of cell averages as opposed to a filter is the absence of loss of spatial resolution and of solution fluctuations averaged over by the filter.

The API is presently under development http://www.ams.sunysb.edu/∼rkaufman/api.

Presently, it links a client code to the FT/LES/SGS front tracking library. All functions in the API are designated as belonging to the API or to the client. Those belonging to the client depend on client specific information, such as the data layout of solution state values. For these functions, the API includes a reference implementation in a few regular cases, such as for regular rectangular grid based state data. The main point of the API functions is to communicate data between the client and the front library. This communication is, of course, two directional. From the client data, we extract velocities, used to advance the front in time. From the front library, we extract front location information, used in a client function to label any given finite difference stencil as regular or irregular, according to whether the stencil does not or does cross a tracked front. Finally, a client function will offer modified (tracked) stencil state values, to affect ghost cells in a neighborhood of the front. The modified stencil of state values will be differenced using the usual client difference algorithm.

FT/LES/SGS has achieved systematic agreement with a wide range of Rayleigh–Taylor experiments, including mixing of both concentration and thermal fields. Uncertainty over initial conditions (not recorded) has been resolved, by backward solution of the fluid equations from an early time (recorded) to an initial time (not recorded), and with uncertainty quantification for the possible errors in the reconstruction [22]. Agreement is precise enough to distinguish between distinct experiments and their distinct values of α. See Table 1.

FronTier was used in a three way code comparison (with RAGE and PROMETHEUS), with experiment and with theory [26] for single mode Richtmyer–Meshkov instabilities. Extensive comparison of FT/LES/SGS with the already validated RAGE code for multimode Richtmyer–Meshkov instabilities, before and after reshock, gave generally excellent results [8], with exceptions related to the thermal field as noted above. In view of the more difficult RT thermal validation and the absence of thermal RM experimental data, we plan to proceed to the next stage of our studies, which is to bring the simulations closer to an ICF experimental plane.

In Sec. 2, we expand on a method for stochastic convergence of turbulent solutions we previously introduced and lay the framework for the convergence strategy of our algorithm. In Sec. 3, we study the convergence of the solution PDFs and their variation with Re. In Sec. 4, we show that the PDFs are sensitive to the definition of turbulent transport. In other words, we establish numerical nonuniqueness for turbulent mixing in the high Re limit. In Sec. 5 we outline our ideas for the validation of numerical studies of hydro instabilities related to ICF performance.