This paper discusses a strategy to initialize a two-dimensional (2D) Reynolds-averaged Navier–Stokes model [LANL's Besnard–Harlow–Rauenzahn (BHR) model] in order to describe an unsteady transitional Richtmyer–Meshkov (RM)-induced flow observed in on-going high-energy-density ensemble experiments performed on the OMEGA-EP facility. The experiments consist of a nominal single-mode perturbation (initial amplitude a 0 10 and wavelength λ = 100 μm) with target-to-target variations in the surface roughness subjected to the RM instability with delayed Rayleigh–Taylor in a heavy-to-light configuration. Our strategy leverages high-resolution three-dimensional (3D) implicit large eddy simulations (ILES) simulations to initialize BHR-relevant parameters and subsequently validate the 2D BHR results against the 3D ILES simulations. A suite of five 3D ILES simulations corresponding to five experimental target profiles is undertaken to generate an ensemble dataset. Using ensemble averages from the 3D simulations to initialize the turbulent kinetic energy in the BHR model (K0) demonstrates the ability of the model to predict the time evolution of the interface as well as the density-specific-volume covariance, b. To quantify the sensitivity of the BHR results to the choice of K0 and the initial turbulent length scale, S0, we execute a parameter sweep spanning four orders of magnitude for both S0 and K0, generating a parameter space consisting of 26 simulations. The Pearson's correlation coefficient is used as a measure of discrepancy between the 2D BHR and 3D ILES simulations and reveals that the ranges 8 S 0 20 μm and 109 K 0 10 10 cm2/s2 produce predictions that agree best with the 3D ILES results.

In the context of inertial confinement fusion (ICF), the mixing between initially separated different materials plays an important role in the performance of capsule implosion.1 For example, the thermonuclear yield may be reduced when the outer ablator material mixes with the central hotspot of the capsule, thus limiting the efficiency of fusion reactions. Material mixing may originate from the growth of hydrodynamic instabilities, such as the Rayleigh–Taylor (RT),2,3 Richtmyer–Meshkov (RM),4,5 and Kelvin–Helmholtz (KH)6,7 instabilities. These instabilities develop due to small imperfections (perturbations) present at the surface of the materials (e.g., due to manufacturing defects), which may grow under the influence of shear (KH), material accelerations (RT), and/or shock waves (RM). The turbulence (which drives mixing) ensuing from these instabilities and its effect on capsule performance are still an active area of research.

Laboratory experiments, especially under high-energy-density (HED) conditions, face many challenges in studying turbulence in part due to the low repetition rate of the experiments—resulting in only partially converged statistics—and the wide range of temporal and spatial scales the experiments are performed under. Furthermore, due to limitations in the diagnostics, experiments typically only provide integrated [i.e., two-dimensional (2D)] measures of three-dimensional (3D) complex flow dynamics, potentially ignoring important physics and mechanisms. Relying on numerical simulations is, therefore, critical and plays a complementary role to the experiments.

Simulating KH-, RT-, or RM-driven mixing using direct numerical simulations (DNS) can be prohibitively expensive because all the spatial and temporal scales have to be resolved, making approaches such as Reynolds-averaged Navier–Stokes (RANS) often the more favorable engineering option for applications like ICF. To this day, no DNS has been performed for ICF even on the largest supercomputers, as the resolution requirements are too stringent.8 However, RANS approaches also face their own challenges: RANS is based on the Reynolds decomposition of a flow where mean quantities are intended to represent an average over an ensemble of realizations, which is often replaced by a spatial average due to the scarcity of ensemble datasets. Replacing ensemble averages by space averages may be appropriate for flows that are in homogenous-, isotropic-, and fully developed turbulent states in which spatial, temporal, and ensemble averaging are often equivalent. However, most HED hydrodynamic experiments involve transitional periods in which the flow is neither homogeneous nor isotropic nor fully developed but may contain large-scale unsteady dynamics; thus, the equivalency of averaging can no longer be assumed. Yet, RANS models often still require to be initialized in such states of turbulence, and knowing how and when to initialize them in a transitional state is, therefore, challenging and is still poorly understood.

The goal of this paper is to develop a strategy allowing the initialization of a RANS model to describe an unsteady transitional RM-induced flow. We seek to examine how ensemble-averaged quantities evolve during the transition to turbulence based on some of the first ensemble experiments repeated under HED conditions. Our strategy involves using 3D high-resolution implicit large eddy simulations (ILES) to supplement the experiments and both initialize and validate the RANS model. We use the Besnard–Harlow–Rauenzahn (BHR) model,9–12 specifically designed to predict variable-density turbulent physics involved in flows like RM. Previous studies have considered different ways of initializing the BHR model. In early generations of the model (original model,9 BHR2,11,13 and BHR310), initialization is done in a heuristic manner, in which the turbulent kinetic energy is initially set to be a fraction of the laminar kinetic energy of the flow—although the meaning of the chosen fraction remains unclear. Later, Rollin and Andrews14 developed a modal model for the RT instability based on the analytic work of Goncharov,15 in which non-linear evolution of RT modes are tracked and used to initialize the BHR model at a time the flow is assumed to have transitioned to turbulence. To the author's knowledge, this approach has only been reported by Braun and Gore12,16 for the classical RM and RT instability16 and reversed-gravity RT,12 the latter introduced in the latest version of the model (BHR4). Other approaches that do not involve heuristic guesses or a modal model have also proven useful. For example, Haines et al.17 have leveraged the use of ILES simulations to validate the BHR model on an inverse chevron shock tube experiment18 and showed that reasonable agreement is obtained when the initial value of the turbulent length scale, S0, is of the same order of magnitude as that of the standard deviation of the surface roughness of perturbations. Haines et al.19,20 have applied this methodology to laser-driven shear experiments and showed reasonable agreement with the experiments. Flippo et al.21 have also demonstrated experimentally that S0 is related to the initial surface roughness. Although the above-mentioned studies have laid down a path for BHR initialization strategies, most of these studies focus primarily on RT or KH instabilities and are conducted in conditions relevant to traditional fluid mechanics. Thus, it is not apparent whether these strategies apply to the RM instability and under HED conditions.

In this work, we perform ensemble 3D high-resolution ILES simulations of on-going RM HED experiments performed on the OMEGA-EP laser facility22 (the ModCons experiments). One the unique features of these experiments is that they are repeated multiple times to generate an ensemble dataset. Although the repetition rate is low—this campaign has taken ∼20 shots over 2 years—these experiments are one of the first ensemble datasets generated under HED conditions.

This paper is organized as follows: first, we describe the ModCons experiments and show the differences between each target in initial conditions and late-time experimental radiographs in Sec. II. Next, we show results from the 3D simulations by first describing the flow evolution from one individual simulation and then considering the ensemble suite of simulations to investigate the order by which the 3D fields are reduced to 2D fields in Sec. III. Using the 3D simulations, we initialize 2D simulations using the BHR model, which are compared to the 3D results in Sec. IV. To evaluate the results from the 2D simulations, a sensitivity analysis with respect to the model parameters S0 and K0 is performed in Sec. IV C, and conclusions are drawn in Sec. V.

The ModCons experiments are one of the experimental campaigns designed by Los Alamos National Laboratory in an effort to expand our current understanding of turbulence and mixing in HED conditions. Originally designed to study the RM instability in single-mode23 and multi-mode24 configurations, it has also served as a modeling platform for combined instabilities25–27 and laser-induced time-dependent accelerations.23 The supporting experiments in this paper focus on single-mode RM.

A typical ModCons target is shown in Fig. 1(a) and consists of a cylindrical layer of polyamide-imide (PAI of density 1.45 g/cm3) adjacent to a layer of foam (density 0.1 g/cm3), both of which having a 1100 μm diameter. For visualization purpose, a 100 μm thick strip of iodine-doped plastic is inserted in the PAI adjacent to the interface between the PAI and the foam. The iodinated plastic is preferentially opaque to the x ray used to radiograph the system. Ripples at the interface are then machined into the PAI and serve as perturbation profile described by a single-mode sinusoidal wave with wavelength λ = 100 μm (wavenumber k = 2 π / λ) and initial amplitude a 0 = 10 μm such that y ( x ) = a 0 sin ( k x ), as shown in Fig. 1(b) from the top view. The system is then irradiated from the brim of the PAI (called the ablator) by three OMEGA-EP laser beams depositing 14.3 kJ of energy, launching a shock wave with a post-shock pressure 10 Mbar. The shock arrives at the interface at 13 ns after the lasers turn on, which causes the RM-instability to grow. At 10 ns, the lasers turn off causing an expansion wave to travel into the system, which reaches the interface at 19 ns. Therefore, the system is RM-unstable alone for 6 ns. When this expansion wave reaches the interface, it decompresses and decelerates the interface, making the system RT-unstable,23,26 such that both RM and RT contribute to perturbation growth past 19 ns.

FIG. 1.

(a) Pre-shot radiograph of a ModCons target and (b) milled ripples shown from the top view.

FIG. 1.

(a) Pre-shot radiograph of a ModCons target and (b) milled ripples shown from the top view.

Close modal

A unique feature of the on-going ModCons experimental campaign is the ability to repeat the experiments multiple times (of the order of 20) generating some of the first HED RM ensemble datasets used to calculate average quantities. Although each realization of the experiments is characterized by the same nominal perturbation profile (single-mode sinusoidal wave described in the previous paragraph), they differ between each other in the random character of the surface roughness imposed by the machining process. Figure 2 shows an example of a rendered perturbation profile over one wavelength [Fig. 2(a)] obtained from point coordinates of one experimental target (target 1) along with its associated spectral content [Fig. 2(b)]. For comparison purpose, Fig. 2(b) also includes the spectral content of four other targets. The spectral content of each target is obtained by taking a Fourier transform of the perturbation profile, y ̂, averaged over two different directions: the first is in the direction into the page [z direction, solid lines in Fig. 2(b)] and the second is in the spanwise direction [x direction, dashed lines in Fig. 2(b)]. As expected, the Fourier transform for the z-averaged profile shows the presence of a dominant mode (low-wavenumber high-amplitude) corresponding to the nominal perturbation, along with higher harmonics (high-wavenumber small-amplitude) corresponding to the surface roughness. By design of the experiments, all targets have the same dominant mode but slightly different higher harmonics due to the randomness of the surface roughness. The Fourier transform of the x-averaged profile also has a dominant mode corresponding to the mean location of the perturbation but its amplitude is different for each target due to slight variations in the location of the mean interface. The amplitudes of the higher harmonics are similar to those in the z-averaged direction.

FIG. 2.

(a) Rendered perturbation profile over one wavelength for target 1. (b) Comparison of the spectral content between different targets when the perturbation profile is averaged in the z direction (solid line) and in the x direction (dashed line).

FIG. 2.

(a) Rendered perturbation profile over one wavelength for target 1. (b) Comparison of the spectral content between different targets when the perturbation profile is averaged in the z direction (solid line) and in the x direction (dashed line).

Close modal

The primary diagnostic for the experiment is line-integrated absorption radiography, and only one high quality radiograph can be taken per shot.23 The images are timed the same on every shot, and an ensemble dataset of density maps is obtained at that time, where each image is a slightly different instantiation of the system evolution due to physical variations in target fabrication, like surface roughness. Figure 3 demonstrates this time-matching procedure over six different realizations of the same experiment. The resolution of the density map is limited by the line-integration through the target (CHI strip thickness), the resolution, and the motion blur. Density fluctuation statistics can still be obtained out the images above a certain scale,28 but these diagnostics limitations motivate us to use ILES simulations to examine other averaged variables or ways of averaging.

FIG. 3.

Experimental radiographs taken at 40 ns after lasers turn-on for six different realizations of the same experiment.

FIG. 3.

Experimental radiographs taken at 40 ns after lasers turn-on for six different realizations of the same experiment.

Close modal

The purpose of using RANS approaches in applications like ICF is to be able to predict results in an engineering-relevant timeframe that would otherwise be impractical to obtain with approaches like DNS. RANS models typically perform averages over a fully turbulent state, where spatial, temporal, and ensemble averages are often equivalent. During the transition to turbulence, however, this equivalency is not straightforward and we seek to examine how averaged quantities evolve with different types of averaging during the transition. To this end, ensemble 3D ILES simulations are performed, which allows us to access both ensemble and spatial averages (we do not consider the temporal average since the flow is in a transitional regime). Moreover, the 3D ILES simulations allow for BHR-relevant parameters (such as the initial turbulent kinetic energy) to be directly computed, which are then used to initialize 2D BHR simulations (Sec. IV). The ensemble set consists of five 3D ILES simulations which simulate the RM-only portion of perturbation growth, i.e., for 6 ns after the arrival of the shock at the interface.

Simulating an entire experimental target as shown in Fig. 1(b) is computationally too expensive, so that only a 100 × 100 μm2 cross section is selected, as depicted by the red square in Fig. 1(b). The computational domain is shown in Fig. 4, where an initial perturbation separates a layer of PAI from a layer of foam. The initial perturbation is generated from the point coordinates of the selected experimental cross section obtained from a stereolithography (STL) file and then reconstructed into a closed surface, as shown in Fig. 2(a). Table I indicates the number of point coordinates used to reconstruct the closed surface for each of the five targets. The amplitude of the initial perturbation, defined as the halfway distance between peak and valley, is a 0 10 μm and its wavelength λ = 100 μm. All simulations are performed with the radiation-hydrodynamics code xRAGE,29 which solves the multi-material compressible Euler equations of fluid motion as follows:
ρ t + x j ( ρ U j ) = 0 ( ρ U i ) t + x j ( ρ U i U j + p δ i j ) = 0 ( ρ E ) t + x j [ U j ( ρ E + p ) ] = 0 , ( ρ Y n ) t + x j ( ρ U j Y n ) = 0 ,
(1)
where ρ is the material density, Ui (i = 1, 2, 3) is the three-dimensional material velocity, p is the pressure, E 1 2 U i 2 + e is the total energy (with e being the specific internal energy), Yn is the mass fraction of species n, and δij is the delta function. These equations assume that body forces, viscous stresses, and heat flux are negligible (further justification is given in the next paragraph). Although the materials are initially in the solid phase, they vaporize upon shock passage and are characterized using a modeled equation of state corresponding to polystyrene (at density 1.45 g/cm3 for the PAI and 0.1 g/cm3 for the foam) obtained from the SESAME database.30 Although the shock in the experiments is driven by laser irradiation onto the PAI, simulating such physics (e.g., via a laser model) in three dimensions is computationally prohibitive. Instead, the shock is artificially driven by setting up inflow conditions obtained from a one-dimensional precursor simulation with a laser model, giving the post-shock density ρ s = 4.5 g/cm3, post-shock pressure p s = 13.5 Mbar, and post-shock material velocity u p = 24 μm/ns. This procedure has been described in detail previously by Di Stefano et al.24 and is, therefore, not repeated here; however,  Appendix A shows profiles of the density, pressure, and velocity obtained from this precursor simulation. In all our simulations (except for the 1D simulation of  Appendix A), time is taken relative to the arrival time of the shock on the interface and not relative to lasers turn-on.
FIG. 4.

Schematic of the computational domain with initial flow field for the 3D ILES simulations.

FIG. 4.

Schematic of the computational domain with initial flow field for the 3D ILES simulations.

Close modal
TABLE I.

Number of point coordinates (obtained from experimental STL files) used to reconstruct the closed surface for each target.

Target number Target 1 Target 2 Target 3 Target 7 Target 9
Number of point coordinates  146 675  141 746  136 639  130 320  133 631 
Target number Target 1 Target 2 Target 3 Target 7 Target 9
Number of point coordinates  146 675  141 746  136 639  130 320  133 631 

The assumptions made to arrive at Eq. (1)—viscous effects, heat flux, mass diffusion, and laser-related physics (such as radiation and three-temperature effects) are negligible—are justified by the fact that the Reynolds number (Re) and Péclet number (Pe) are high and that laser-related physics are important mainly when simulating the laser-energy deposition on the target. To further justify these assumptions, we performed a supporting 2D simulation of the entire target with a laser model31–33 (data not shown), and concluded that effects of the lasers on the evolution of the bulk hydrodynamics are minimal; including such effects would not change the results of this paper. Concerning viscous effects, we estimate the Reynolds number R e 10 5 based on the post-shock velocity (24 μm/ns), initial perturbation amplitude (10 μm), and a plasma kinematic viscosity (0.014 cm2/s). The latter was previously used in HED KH experiments34,35 with similar materials than the ones of the present study. Concerning mass diffusion, the Péclet number—the ratio of advective mass transfer to diffusive mass transfer—is estimated as P e 10 5 based on the study of Bender et al.36 who simulated an RM flow also with similar materials than the ones of the present study.

The system of Eq. (1) is discretized using a second-order Godunov finite-volume scheme on a cartesian grid generated by an adaptive-mesh-refinement procedure. The latter is chosen to have four levels of refinement, where level 0 corresponds to the cell size away from the interface (set to 1 μm) and level 3 (which is the minimum cell size in the domain) corresponds to the cell size at the interface (set to 1/8 μm). Inflow and outflow boundary conditions are imposed at the entrance and exit of the shock tube, respectively, and the four boundaries on the sides are reflecting. Each 3D simulation employs approximately 108 cells and runs approximately for 104 CPU-hours over 180 cores. A fourth-order Runge–Kutta scheme is used to march the equations in time with a varying time step of the order of 0.1 ps.

Figure 5 shows the time evolution of the perturbation for an individual target at three different times. Because the shock travels from a heavy to a light material, the initial interface undergoes a phase inversion, i.e., the troughs of the perturbation become peaks and vice versa. The type of phase inversion (direct,37 indirect,38 or marginal25,26) is determined by the relative speed of the post-shock interface with respect to the speed of the incident shock. In our case, since the latter is approximately the same as the speed of the post-shock interface, the initial interface undergoes a marginal phase inversion. As a result, the interface is relatively flat immediately after the passage of the shock, as shown in Fig. 5(a). Growth of the small scale roughness is evident at that time. As the nominal portion of the perturbation grows into a large bubble-and-spike structure—where the spike corresponds to the heavy material penetrating the light one and the bubble to the light material penetrating the heavy one—the perturbations associated with the small scale roughness also grow into their own (but smaller) bubble and spike, as depicted by the arrows in Fig. 5(b). Eventually, large KH roll-ups form along the part of the interface between the nominal bubble and spike, as seen in Fig. 5(c).

FIG. 5.

Three-dimensional isosurfaces of density of an individual target at (a) t = 1.1, (b) t = 3.7, and (c) t = 6.0 ns.

FIG. 5.

Three-dimensional isosurfaces of density of an individual target at (a) t = 1.1, (b) t = 3.7, and (c) t = 6.0 ns.

Close modal

Figure 6 compares the time evolution of the amplitude of the nominal perturbation between the 3D simulations and the experiments of Di Stefano et al.23 and shows good agreement between the two. Error bars of ± 3 μm in amplitude and ±0.5 ns in time are added to reflect the uncertainty made during the measurements. The agreement between the simulation and the experiments gives us confidence that the bulk hydrodynamic parameters (such as the shock strength and bubble-and-spike evolution) can be reproduced even with our artificial initialization of the shock.

FIG. 6.

Perturbation amplitude vs time obtained from an individual 3D ILES simulation and the experiments.23 

FIG. 6.

Perturbation amplitude vs time obtained from an individual 3D ILES simulation and the experiments.23 

Close modal
To evaluate the state of turbulence in the flow, a useful quantity to consider is the spectrum of the turbulent kinetic energy (TKE) defined as
E ( κ ) = 1 2 U ̂ i U ̂ i * ,
(2)
where the overhat represents the Fourier transform and the star the complex conjugate of the Fourier transform. Note that Eq. (2) uses the instantaneous velocity instead of the fluctuations since the mean velocity only affects the spectrum by adding a zero frequency, which is not visible in a log –log plot. The spectrum is represented as a function of the radial wavenumber κ = k x 2 + k z 2. Equation (2) leads to a 3D field of E, which is reduced to a 1D function of κ by averaging the values of E contained in the cylindrical strip defined by the radius κ, length Ny (where Ny is the number of points in the y direction), and thickness r i κ r i + 1, where r is a sample radial coordinate with i = 1 : N x (Nx being the number of points in the x direction). Figure 7 shows the TKE spectrum as a function of the radial wavenumber at different times. At all times, the lowest wavenumber (i.e., largest wavelength) mode dominates; this mode corresponds to the large nominal perturbation present at all times and growing in a bubble-and-spike structure. The high-wavenumber modes ( κ 10 μm−1) do not change the TKE levels significantly over time, as evidenced by the merging of all the curves into one, indicating that lower wavenumber modes are the ones responsible for TKE increase in the flow. Indeed, there is an amplification of the TKE levels by all modes contained in the range 1 κ 10 μm−1, although this amplification is significantly reduced past t 2 ns for modes in the range 5 κ 10 μm−1, indicating an increase in TKE levels due to lower and lower wavenumber modes. These modes are representative of the formation of the large KH roll-ups in the flow. Although the flow is neither homogeneous, isotropic, nor fully turbulent, the theoretical framework of homogeneous isotropic turbulence is often used as a reference (in particular the Kolmogorov power law κ 5 / 3). This reference is added to Fig. 7 (solid black line) and shows that part of the spectrum adopts the κ 5 / 3 trend over time, indicating the presence of an inertial subrange in which the large vortical structures associated with the KH roll-ups break down and transfer their kinetic energy to smaller eddies. It is interesting to note that the spectrum suggests that a dissipation range may be also present even though no physical viscosity is present (since only the Euler equations are solved); we attribute this dissipation to numerical viscosity introduced by the ILES numerical framework. In  Appendix B, we show that although numerical viscosity is present, the ILES framework provides a reasonable description of the flow. The trends shown in Fig. 7 reveal that the flow is indeed in a transitional state but that there is evidence of an evolution toward the Kolmogorov scales of developed turbulence, thus further motivating the need to understand the evolution of averaged quantities (which are used in BHR) in a transitional state.
FIG. 7.

Turbulent kinetic energy spectrum at different times. The solid black line represents the κ 5 / 3 slope.

FIG. 7.

Turbulent kinetic energy spectrum at different times. The solid black line represents the κ 5 / 3 slope.

Close modal
One of the implications of the flow not being in a fully developed turbulent state is that spatial and ensemble averages are not equivalent. In the experiments, although the flow is three-dimensional, the measurements are two-dimensional due to constraints in the diagnostics. Therefore, the 3D simulations must be averaged to obtain 2D quantities that can be compared to experimental data and used as validating reference for 2D BHR simulations. In this section, we examine whether the method of averaging to reduce the 3D fields to 2D fields is important in obtaining quantities of interest relevant to the BHR model. Specifically, we examine the order by which the 3D fields are reduced to 2D fields in the calculation of the density-specific-volume covariance, b, relevant to BHR. We consider two cases: one case consists of space averaging first then ensemble averaging (referred to as space-ensemble and denoted · s e) while the other case consists of ensemble averaging first then space averaging (referred to as ensemble-space and denoted · e s). The quantity b is a measure of density fluctuations in the flow and is a pertinent quantity to choose because the density field is the most accessible data from the experiments. It is defined as
b ρ v ,
(3)
where v 1 / ρ is the specific volume, the primes represent fluctuating quantities, i.e., f f f (for a general variable f), and the brackets represent either a space average (denoted by · z) or an ensemble average (denoted by · N), where N is the total number of realizations. By definition, b is non-negative and its values range between 0 and 1, where b = 0 signifies that the flow is fully mixed. The procedure for calculating b constitutes a series of steps that are dependent upon the type of averaging used.

1. The space-ensemble case: · s e

The space-ensemble case corresponds to the method of averaging closest to the averaging method performed in the experiments. Figure 8 shows the steps followed to obtain a 2D field of b from the 3D fields of the density. The first row shows the 3D density field for each target, which is first averaged in space in the z direction, leading to a 2D field in the (x, y) plane for each individual instantiation (step 0, second row), similar to line-absorption radiography. The ensemble average is then taken leading to a 2D ensemble-averaged density field, ρ e ( x , y ) (step 1). The large bubble and spike of the nominal profile are similar for each target but the small scale structures are different. Ensemble averaging conserves the shape of the nominal perturbation but results in smearing the small-scale structures. In the limit N , these small-scale structures would disappear and only the nominal perturbation would remain.

FIG. 8.

Steps followed to generate a 2D field of b at t = 6.0 ns in the space-ensemble case ( · s e), which corresponds to the steps performed in the experiments. Each step is encircled in a red rectangle.

FIG. 8.

Steps followed to generate a 2D field of b at t = 6.0 ns in the space-ensemble case ( · s e), which corresponds to the steps performed in the experiments. Each step is encircled in a red rectangle.

Close modal

The third row shows the 2D field of the density fluctuations, ρ , obtained by subtracting ρ e ( x , y ) from the density field for each instantiation: ρ ( x , y ) = ρ ( x , y ) ρ e ( x , y ) (step 2). The same is performed on the specific volume to obtain v .

The fourth row shows the 2D field of the product ρ v (step 3), which is only positive because regions of positive ρ correspond to regions of negative v and vice versa.

The resulting 2D field of b(x, y) is obtained by ensemble averaging the fourth row and is shown in the last column of the fourth row (step 4).

2. The ensemble-space case: · e s

In the ensemble-space case, all operations are performed on the 3D fields until the last step to obtain a 2D field of b. As in the space-ensemble case, Fig. 9 shows the steps followed to obtain a 2D field of b from the 3D fields of the density. The first row shows again the 3D density field for each target, which is now ensemble averaged in 3D first, leading to a 3D ensemble-averaged density field, ρ e ( x , y , z ) (step 1).

FIG. 9.

Steps followed to generate a 2D field of b at t = 6.0 ns in the ensemble-space case ( · e s). Each step is encircled in a red rectangle. The pictures in the dashed box correspond to the case when the 3D fields of the product ρ v are space averaged for each target.

FIG. 9.

Steps followed to generate a 2D field of b at t = 6.0 ns in the ensemble-space case ( · e s). Each step is encircled in a red rectangle. The pictures in the dashed box correspond to the case when the 3D fields of the product ρ v are space averaged for each target.

Close modal

The second row shows the 3D field of the density fluctuations obtained by subtracting ρ e ( x , y , z ) from the density field for each instantiation: ρ ( x , y , z ) = ρ ( x , y , z ) ρ e ( x , y , z ) (step 2). The same is performed on the specific volume to obtain v .

The third row shows the 3D field of the product ρ v (step 3). Ensemble averaging the third row leads to a 3D field of b (step 4), which is then space averaged in the z direction to obtain a 2D field of b (step 5, fourth row, last column). There is also another way to obtain this 2D field of b, which is to space average the 3D product ρ v for each target, which leads to 2D fields of b for each target, which are then ensemble averaged to obtain the final 2D field of b. This case is encircled in a dashed box in Fig. 9.

Each step is summarized in Table II for both methods of averaging. Only the steps for ρ are outlined but the same procedure follows for v. The steps are included in red in Figs. 8 and 9. The 2D fields of b obtained from both methods of averaging are compared in Fig. 10 and shows that both methods lead to similar results. To measure their differences, we calculate the root mean square (RMS) of the difference between the b obtained from the space-ensemble case and the one obtained from the ensemble-space case, which gives RMS 0.01. This low RMS value indicates that the method of averaging performed on the 3D ILES simulations is not important, i.e., it does not affect subsequent comparison with 2D-BHR simulations. This result shows two things: first, when compared to 2D-BHR simulations, one method is not better than the other but both are valid comparison, and second, that both methods connect the 3D ILES simulations to a 2D representation valid for the experiment. Note that this low RMS value is obtained only with five instantiations; it is interesting to conjecture that the RMS would tend to zero as the number of instantiations increases.

TABLE II.

Steps followed to calculate b when considering the space–ensemble vs ensemble–space cases.

Quantity calculated Space–ensemble · s e Field dimension Ensemble–space · e s Field dimension
Step 0  ρ  ρ [ x , y , z , N ] z ρ [ x , y , N ]  2D  –  – 
Step 1  ρe  ρ [ x , y , N ] N ρ e [ x , y ]  2D  ρ [ x , y , z , N ] N ρ e [ x , y , z ]  3D 
Step 2  ρ = ρ ρ e  ρ [ x , y , N ] = ρ [ x , y , N ] ρ e [ x , y ]  2D  ρ [ x , y , z , N ] = ρ [ x , y , z , N ] ρ e [ x , y , z ]  3D 
Step 3  ρ v   ρ [ x , y , N ] v [ x , y , N ]  2D  ρ [ x , y , z , N ] v [ x , y , z , N ]  3D 
Step 4  ρ v   ρ [ x , y , N ] v [ x , y , N ] N b [ x , y ]  2D  ρ [ x , y , z , N ] v [ x , y , z , N ] N b [ x , y , z ]  3D 
Step 5  –  –  –  b [ x , y , z ] z b [ x , y ]  2D 
Quantity calculated Space–ensemble · s e Field dimension Ensemble–space · e s Field dimension
Step 0  ρ  ρ [ x , y , z , N ] z ρ [ x , y , N ]  2D  –  – 
Step 1  ρe  ρ [ x , y , N ] N ρ e [ x , y ]  2D  ρ [ x , y , z , N ] N ρ e [ x , y , z ]  3D 
Step 2  ρ = ρ ρ e  ρ [ x , y , N ] = ρ [ x , y , N ] ρ e [ x , y ]  2D  ρ [ x , y , z , N ] = ρ [ x , y , z , N ] ρ e [ x , y , z ]  3D 
Step 3  ρ v   ρ [ x , y , N ] v [ x , y , N ]  2D  ρ [ x , y , z , N ] v [ x , y , z , N ]  3D 
Step 4  ρ v   ρ [ x , y , N ] v [ x , y , N ] N b [ x , y ]  2D  ρ [ x , y , z , N ] v [ x , y , z , N ] N b [ x , y , z ]  3D 
Step 5  –  –  –  b [ x , y , z ] z b [ x , y ]  2D 
FIG. 10.

Comparison of the 2D field of b between the space-ensemble case (left) and the ensemble-space case (right).

FIG. 10.

Comparison of the 2D field of b between the space-ensemble case (left) and the ensemble-space case (right).

Close modal

Having highlighted the differences between space-averaging and ensemble-averaging of the 3D ILES simulations, we now validate 2D BHR simulations against the 3D ILES simulations. The initial perturbation profile in the 2D simulations corresponds to the nominal perturbation of the 3D profile, i.e., λ = 100 and a 0 = 10 μm without the surface roughness. Such initial profile corresponds to the profile obtained if the ensemble average of the full 3D profile (i.e., with the surface roughness) was taken over an infinite amount of realizations.

The BHR model is based on the Favre decomposition of the flow, i.e., a variable f, as in a Reynolds decomposition, is decomposed in a mean and fluctuating parts, f = f ̃ + f , but where f ̃ is a density-weighted average, i.e., f ̃ = ρ f N / ρ N, where · N represents the ensemble average. Decomposing all variables in Eq. (1) this way and averaging these equations over an ensemble of realizations leads to a set of averaged transport equations that the BHR model evolves.10 In practice, three parameters are required to initialize the BHR model: the initial TKE, K0, the initial turbulent length scale, S0, and the time at which the model is “activated,” tBHR. These parameters are evaluated by directly computing their values from the 3D ILES simulations. The TKE is defined as
K = 1 2 ρ u i u i N ρ N ,
(4)
where u i are the velocity fluctuations defined as
u i = U i U ̃ i .
(5)

First, the ensemble average of the 3D density field is computed, giving ρ N. Then, U ̃ i is computed by taking the product of the 3D density field with the 3D velocity field for each target, ρ U i, taking the ensemble average of the product ( ρ U i N), and dividing by ρ N to give U ̃ i. The 3D velocity fluctuations are then computed for each target via Eq. (5); then, the triple product ρ u i u i is also computed for each target and ensemble averaged to finally give a 3D field of the TKE.

Since K0 is a global quantity (i.e., it only takes one value), the 3D TKE field must be reduced to a single value, which is done by computing its volumetric average in the entire field, K ¯ = ( 1 / d x dydz ) K dxdydz, and its time evolution is shown in Fig. 11. The TKE sharply increases between the initial time and t 1 ns due to the passage of the shock over the interface (increasing velocity fluctuations) and then decreases with some slight oscillations. The fact that the TKE has a peak value indicates that most of the TKE is produced during the interaction of the shock with the interface, suggesting that turbulent motion becomes significant only after t 1 ns. Therefore, the model “activation” time is chosen to be tBHR = 1 ns and the peak value K ¯ = 2.8 × 10 9 cm2/s2 is chosen as initial value for K0. Concerning S0, Haines et al.17 suggested that S0 is related to the standard deviation of the initial material interface. Using the point coordinates of the initial 3D surface from an STL file, the standard deviation is about 0.2 μm, which is the value chosen for S0. Note that if an ILES simulation is not available, strategies using potential flow models have been previously adopted,12,14,39 in which the evolution (amplitude and velocity) of initial modes present along the interface is tracked using linear perturbation analysis or vorticity-based equations. These potential flow models are run until an appropriate time typically determined by a critical Reynolds number, above which the flow may be turbulent.

FIG. 11.

Time evolution of the volumetric average of the 3D TKE field.

FIG. 11.

Time evolution of the volumetric average of the 3D TKE field.

Close modal

The values S 0 = 0.2 μm, K 0 = 2.8 × 10 9 cm2/s2, and tBHR = 1 ns constitute our baseline case and is discussed below. The model is “activated” in the entire domain with constant values for S0 and K0, and the other model parameters (such as the turbulent mass flux and the density-specific-volume covariance) are initialized at zero.

Our starting reference in the validation process is the case when K 0 = 2.8 × 10 9 cm2/s2, S 0 = 0.2 μm, and tBHR = 1 ns. Figure 12 compares the density field, the b field, and the TKE field between the 2D BHR simulations and the 3D ILES simulations at different times. As in Fig. 9, the 3D ILES simulations are first ensemble averaged over all realizations and then spatially averaged in the z direction to obtain a 2D field in the (x, y) plane. Initially, the interface is sharp in the 2D simulations but has some thickness in the 3D simulations, as seen on the density field, which is the result of ensemble averaging only over five realizations. The b value in the 2D simulation is initialized at zero but is close to one over the thickness of the interface in the 3D simulation (i.e., the materials are segregated), reflecting the fact that the location of the mean interface is slightly different for each realization. The TKE has the initial value of K0 in the 2D case and zero in the 3D case (as seen on the logarithmic colorbar). Slightly after shock passage over the interface (t = 1.1 ns), the interface has flattened (with growth of the small scale roughness in the 3D case) and an expansion wave is reflected into the PAI, as seen by the high-density values ahead of the wave. Both b and the TKE increase in the 3D simulation across the interface and the shock due to an increase in density and velocity fluctuations associated with the surface roughness. In the 2D simulation, b is still close to zero everywhere and the TKE increases but is approximately two orders of magnitude lower than in the 3D simulation. At t = 3.7 ns, the bubble and spike are clearly visible in both the 2D and 3D simulations with a non-uniform distribution of b along the interface: higher values near the spike than along the rest of the interface, as seen in the 3D simulation. These higher values are again the result of the low number of realizations in the ensemble average and are not visible in the 2D simulation. The TKE in the 3D simulation also shows spatial variation along the interface but with relatively similar magnitude on each half of the domain. Evidence of shock refraction patterns near the interface is also visible but not in the 2D simulation. At t = 6.0 ns, the roll-ups have now formed, causing an increase in both b and the TKE in the roll-up regions, as seen in the 3D simulation. Although this increase is slightly visible in the 2D simulation, the magnitude is lower by approximately a factor of two for b and an order of magnitude for the TKE.

FIG. 12.

Comparison between the 2D BHR simulations ( S 0 = 0.2 μm and K 0 = 2.8 × 10 9 cm2/s2) and the z-averaged 3D ILES simulations at different times for the baseline case for (a) the density, (b) the b parameter, and (c) the TKE.

FIG. 12.

Comparison between the 2D BHR simulations ( S 0 = 0.2 μm and K 0 = 2.8 × 10 9 cm2/s2) and the z-averaged 3D ILES simulations at different times for the baseline case for (a) the density, (b) the b parameter, and (c) the TKE.

Close modal

From this comparison between the 2D and the 3D simulations, the 2D-BHR simulation overall predicts well the behavior of the nominal perturbation at all times, both in terms of the shape of the interface and the growth of the bubble and spike. The interface in the 2D simulation stays smooth over time whereas growth of the small scale structures is shown in the 3D simulation (due to the relatively low number of realizations in the ensemble average). Concerning the parameter b and the TKE, the agreement between the 2D and 3D simulations is not as good: the 2D simulation predicts a thinner area over which b and the TKE are distributed, lower values in magnitude, and does not capture the high-b regions near the spike seen in the 3D simulation. To understand the sensitivity of these results with the initial values of K0 and S0, we report next on a parameter sweep with respect to K0 and S0.

Figure 13 shows the density, the b parameter, and the TKE at t = 6 ns in the parameter space (S0, K0) from the 2D simulations. For comparison purpose, the fifth column shows the solution from the 3D simulation (ensemble averaged first then space averaged in the z direction) and is repeated in each row. Increasing S0 or K0 has the general trend of diffusing the interface. Increasing K0 too much tends to diffuse the interface so much that the roll-ups are no longer visible (e.g., when S 0 = 10 μm and K 0 = 10 12 cm2/s2), while increasing S0 too much tends to underpredict the magnitude of b and the TKE (e.g., when S 0 = 100 μm and K 0 = 10 9 cm2/s2). The high-b regions near the spikes are best predicted when S 0 = 1 μm and K 0 10 9 cm2/s2 but the interface is too thin compared to the 3D case. Increasing S0 to S 0 = 10 μm thickens the interface to thicknesses closer to the 3D case but the peak values near the spike are then reduced. Concerning the TKE, its distribution along the interface as well as its order of magnitude for S 0 = 10 μm and K 0 = 10 10 cm2/s2 relatively corresponds to the ones for the 3D case. These qualitative observations suggest that the ranges 1 S 0 10 μm and 109 K 0 10 12 cm2/s2 may be the most appropriate for comparison with the 3D results. To further refine these ranges, we conducted an additional six simulations in a refined region of the parameter space corresponding to S 0 = 3 , 5 ,   and   7 μm and K 0 = 10 11   and   10 12 cm2/s2, and are shown in  Appendix C. For all three values of S0, the case K 0 = 10 12 cm2/s2 results in an interface that is too diffuse and the peak value of b near the spikes is not reproduced. For K 0 = 10 11 cm2/s2, all three values of S0 give similar results and would be appropriate values to use.

FIG. 13.

(a) Density field, (b) b field, and (c) TKE field in the parameter space (S0, K0) at t = 6.0 ns. The fifth column shows the corresponding 3D result and is repeated in each row. The blank field for S 0 = 100 μm and K 0 = 10 12 cm2/s2 corresponds to the case where the simulation could not reach the specified end time due to excessive diffusion introduced by the BHR model.

FIG. 13.

(a) Density field, (b) b field, and (c) TKE field in the parameter space (S0, K0) at t = 6.0 ns. The fifth column shows the corresponding 3D result and is repeated in each row. The blank field for S 0 = 100 μm and K 0 = 10 12 cm2/s2 corresponds to the case where the simulation could not reach the specified end time due to excessive diffusion introduced by the BHR model.

Close modal
To quantify the above observations, we consider the following Pearson's correlation coefficient (a measure of correlation between two datasets) between the 2D and 3D (collapsed to 2D) fields for the b variable:
r = cov ( b 2 D , b 3 D ) σ b 2 D σ b 3 D = i ( b i 2 D b ¯ i 2 D ) ( b i 3 D b ¯ i 3 D ) i ( b i 2 D b ¯ i 2 D ) 2 i ( b i 3 D b ¯ i 3 D ) 2 ,
(6)
where i = 1 , , M, where M is the total number of pixels in the b image, and b ¯ is the mean value of b taken over the entire image. Since the b-fields are 2D but the summation 1D, the 2D images are flattened into a 1D array. The Pearson's coefficient represents the ratio of the covariance between the 2D and 3D datasets, cov ( b 2 D , b 3 D ), to their respective standard deviations, σ b 2 D and σ b 3 D. Its values range from –1 to 1 with 0 meaning the two datasets are uncorrelated. When the coefficient is 1 or –1, the two datasets are linearly correlated. Figure 14 shows the correlation coefficient in the parameter space (S0, K0) at t = 6.0 ns and reveals that values closest to 1 lie in the region 8 μm S 0 20 μm and 109 cm2/s2 K 0 10 10 cm2/s2. Note that this range of values for S0 slightly differs from the values found from the qualitative analysis above, but the order of magnitude is the same.
FIG. 14.

Pearson's correlation coefficient in the parameter space (S0, K0) at t = 6.0 ns. The white circles show the values of S0 and K0 at which the 2D simulations are conducted.

FIG. 14.

Pearson's correlation coefficient in the parameter space (S0, K0) at t = 6.0 ns. The white circles show the values of S0 and K0 at which the 2D simulations are conducted.

Close modal
The Pearson's correlation coefficient is helpful in examining the differences between the 2D and 3D images. To further quantify these differences, we investigate the sensitivity of the 2D predictions with respect to S0 and K0 separately. We consider two quantities as follows:
Mean TKE : TKE ¯ = 1 d x d y k ( x , y ) dxdy , Mean b : b ¯ = 1 d x d y b ( x , y ) dxdy .
(7)
Figures 15 and 16 show the time evolution of both TKE ¯ and b ¯ for different values of S0 (with fixed K0) and K0 (with fixed S0), respectively.
FIG. 15.

Time evolution of the volumetric average TKE when varying (a) S0 with K 0 = 10 9 cm2/s2 and (b) K0 with S 0 = 10 μm.

FIG. 15.

Time evolution of the volumetric average TKE when varying (a) S0 with K 0 = 10 9 cm2/s2 and (b) K0 with S 0 = 10 μm.

Close modal
FIG. 16.

Time evolution of the volumetric average b when varying (a) S0 with K 0 = 10 9 cm2/s2 and (b) K0 with S 0 = 10 μm.

FIG. 16.

Time evolution of the volumetric average b when varying (a) S0 with K 0 = 10 9 cm2/s2 and (b) K0 with S 0 = 10 μm.

Close modal

Figure 15 shows that S 0 = 1 μm and K 0 = 10 9 cm2/s2 give the closest match to the 3D ILES simulations since the behavior of the mean TKE widely deviates from the 3D simulations for the other values of S0 and K0. Figure 16 shows that the evolution of b disagrees by similar factors for all values of S0 but the trend is relatively following the one from the 3D simulations. This disagreement in the magnitude of b comes from initializing b = 0 in the 2D-BHR simulations, which is not the case as demonstrated by the non-zero value of b at all times. When varying K0, the same observation of wrong initial value but right trend applies, with the best agreement again for K 0 = 10 9 cm2/s2.

This analysis helps us conclude that there exists a case of a 2D-BHR simulation ( 8 S 0 20 μm and 109 K 0 10 10 cm2/s2) that is a good comparison to the 3D ILES simulations and ultimately to the experiments.

In this work, we develop a strategy to initialize a 2D RANS model (LANL's BHR model) in the context of on-going HED ensemble experiments of laser-driven RM with delayed RT instabilities. Our strategy leverages 3D ILES simulations to initialize BHR-relevant parameters used in 2D BHR simulations, which are subsequently validated against the 3D ILES simulations.

By performing a suite of five high-resolution 3D ILES simulations corresponding to five experimental target profiles, we are able to provide an ensemble dataset used to evaluate the averaging method performed in the experiments, the latter taking an ensemble average of measured 2D mix quantities. We show that the order by which our 3D results are reduced to 2D fields (i.e., either by taking a space average of the 3D fields first followed by an ensemble average, or by taking an ensemble average of the 3D fields first followed by a space average) is not important when calculating mix quantities. This result is key because it shows that a direct comparison between experimental data and subsequent 2D BHR simulations can be achieved without additional constraints for the BHR model.

The availability of a simulated 3D ensemble dataset further allows us to initialize 2D BHR simulations by directly computing the initial parameters needed for BHR (initial TKE, K0, initial turbulent length scale, S0) from the 3D ILES simulations. By initializing K0 with the peak value of the TKE at the time immediately after shock passage and S0 with a value of the same order of magnitude than the standard deviation of the initial surface roughness, good agreement is obtained with the 3D ILES simulations for the shape of the interface but not for the spatial distribution of b nor the TKE. We show that the latter can be reproduced by BHR by increasing both S0 and K0 but within a reasonable range: increasing S0 too much underpredicts the magnitude of b while increasing K0 too much overdiffuses the interface to the point where roll-up structures are no longer visible. By investigating the distribution of the Pearson's correlation coefficient in the parameter space (S0, K0) generated by 26 2D BHR simulations, we show that the ranges 8 S 0 20 μm and 109 K 0 10 10 cm2/s2 produce predictions that agree best with the 3D ILES results. This range for S0 is an order of magnitude higher than the standard deviation of the initial surface roughness, which is in agreement with the previous findings of Haines et al.17 This parameter sweep with respect to K0 and S0 allows us to show that a case of a 2D BHR simulation corresponding to the experiments exist, which gives us confidence in future direct comparison with experimental data.

This work sets a precedent for practical BHR predictions of future HED experiments in which only a few realizations can usually be afforded. Further validation against the ModCons experimental dataset is warranted and will require the implementation of a laser model in 3D to capture the late-time delayed RT physics from laser turn-off. On-going research is tackling this task and will be used in a future study.

We wish to acknowledge Noah Braun and Brian Haines for helpful conversations about the BHR model and details on numerical methods, as well as Fernando Zigunov for fruitful conversations. Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by Triad National Security, LLC for the National Nuclear Security Administration of U.S. Department of Energy under Contract No. 89233218CNA000001.

The authors have no conflicts to disclose.

Sam Pellone: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Alexander Martin Rasmus: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (equal); Writing – review & editing (supporting). Carlos Alex Di Stefano: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Supervision (supporting); Writing – review & editing (supporting). Elizabeth C. Merritt: Conceptualization (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Supervision (supporting); Writing – review & editing (supporting). Forrest W. Doss: Funding acquisition (supporting); Supervision (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To give more details about the inflow conditions used in the 3D simulations, we perform a 1D precursor simulation with a laser model to track the evolution of the density, pressure, and velocity induced by the laser-generated shock. The energy deposition by the laser onto the target is modeled using a laser package,31–33 which aims to reproducing the 10 ns, square-shaped drive pulse present in the experiments. Figure 17 shows pressure values in the (t, y) plane from the precursor simulation. Time is taken relative to when the lasers turn on. Initially, material ablation from lasers turn-on launches a shock in the PAI material, as seen by the high-pressure region behind the shock (10 Mbar pressure 15 Mbar). After further examination, the pressure immediately behind the shock is 13.5 Mbar. This shock reaches the PAI-foam interface at ∼13 ns, which causes the interface to be displaced and a shock to be transmitted in the foam material, as seen by the pressure jump in the latter. The interface location is obtained by tracking the position of a passive tracer added at the location of the initial interface. At 10 ns, the lasers turn off, as seen by the sharp decrease in the pressure, which launches an expansion wave in the domain. To obtain further information on the evolution of the 1D system, Fig. 18 shows the evolution of the density [Fig. 18(a)] and the velocity [Fig. 18(b)] profiles at four different times: the initial time (t = 0 ns), a time after lasers turn-on (t = 5 ns), a time after lasers turn-off (t = 12 ns), and a time corresponding to the end of the RM-only phase (which lasts 6 ns after lasers turn-off, t = 16 ns). The location of the PAI-foam interface (referred to as y int) is shown by black dashed lines. Since the shock does not arrive at the interface until t = 13 ns, the interface location is the same for all times and corresponds to its initial location expect for t = 16 ns. Initially, the density profile corresponds to the three-material system used to initialize the simulation: the PAI, the foam, and a background vacuum modeled as a near-zero density material. This background region is 4300 μm thick behind the PAI and serves as a region in which the laser drive is initialized. At this time, the velocity is zero everywhere since the lasers have not turned on yet. After the lasers turn on, a shock travels into the PAI with a tail of decreasing density corresponding to the ablation region of the PAI material, as seen at t = 5 ns. This tail is also visible on the velocity profile with large negative values corresponding to the ablated material traveling in the direction opposite of the incoming laser beams. Note the magnitude of the density and velocity of the shocked PAI material: ρ s 4.5 g/cm3 and u s 24 μm/ns, respectively. At t = 12 ns, there is a change of slope in the density and velocity profiles in the tail region between 300 y 380 μm, corresponding to the expansion wave generated by the lasers turn-off at t = 10 ns. This expansion wave travels in the same direction as the shock, as shown by the positive values of the velocity. At t = 13 ns, the shock arrives at the interface, causing a shock to be transmitted into the foam and an expansion wave to be reflected into the PAI. At t = 16 ns, this traveling shock and expansion are evident.

FIG. 17.

Pressure contours in the (t, y) plane from the 1D simulation with the laser model. The white solid line represents the interface location and the black dashed line represents the time at which the lasers turn off. Time is relative to lasers turn-on.

FIG. 17.

Pressure contours in the (t, y) plane from the 1D simulation with the laser model. The white solid line represents the interface location and the black dashed line represents the time at which the lasers turn off. Time is relative to lasers turn-on.

Close modal
FIG. 18.

Density (a) and velocity (b) profiles from the 1D simulation at four different times. The location of the interface is shown by black dashed lines.

FIG. 18.

Density (a) and velocity (b) profiles from the 1D simulation at four different times. The location of the interface is shown by black dashed lines.

Close modal
The numerical viscosity in the ILES context represents the numerical dissipation introduced by the numerics on the solution, and is relevant when calculating the effective viscosity, ν eff = ν + ν n, where ν is the kinematic molecular viscosity and νn is the numerical viscosity. From νeff, an effective Reynolds number can then be calculated R e eff = U L / ν eff, where U and L are characteristic velocity and length scales, respectively. An estimation of the numerical viscosity can be performed40–42 based on the dissipation rate of the TKE, ε, and the strain-rate tensor, S i j 1 2 ( U i x j + U j x i ),
ν eff = ε 2 S i j S i j ,
(B1)
where the sum over repeated indices is implied (i.e., double-dot product), and the brackets represent a volumetric average. Using dimensional arguments, ε = D U 3 / L, where D is a dimensionless ratio approaching a constant value when the Reynolds number is large enough ( D 0.5). Following the methodology of Zhou et al.,42  L is taken to be the thickness of the mixing zone (defined in terms of mixedness), and U is of the order of the root mean square turbulence intensity, U = ( 2 TKE 3 ) 1 / 2. To compare νeff to the dissipation associated with the BHR model, we consider the turbulent viscosity based on the values of the initial TKE (K0) and the initial turbulent length scale (S0), ν t = 0.28 S 0 K 0. Figure 19 shows the time evolution of the effective viscosity and Reynolds number [Fig. 19(a)], along with a map of the turbulent viscosity in the parameter space (S0, K0) [Fig. 19(b)]. The effective Reynolds number toward the end of the simulation is R e 10 5, indicating that the effective value of the viscosity is not as significant and the ILES approximation is reasonable. Around 1 ns (when the shock impacts the interface), there is a rapid decrease in the Reynolds number to R e 10 2, suggesting that numerical viscosity may have more of an impact on the solution around that time, and reflecting the dissipative nature of shock-capturing schemes in the ILES framework. When comparing the numerical viscosity to the turbulent viscosity [Fig. 19(b)], the latter is multiple orders of magnitude higher than the former, reflecting the fact that RANS models are more dissipative than an ILES description.
FIG. 19.

(a) Effective viscosity and Reynolds number vs time. (b) Value of the turbulent viscosity in the parameter space (S0, K0). The solid white line represents the maximum effective viscosity present in the ILES simulations.

FIG. 19.

(a) Effective viscosity and Reynolds number vs time. (b) Value of the turbulent viscosity in the parameter space (S0, K0). The solid white line represents the maximum effective viscosity present in the ILES simulations.

Close modal

Figure 20 compares the 2D density field, b field, and TKE field to the 3D simulations for the refined values of S0 and K0.

FIG. 20.

(a) Density field, (b) b field, and (c) TKE field in the refined parameter space 3  S 0 7 μm and K 0 = 10 11 and 1012 cm2/s2 at t = 6.0 ns. The fourth column shows the corresponding 3D result and is repeated in each row.

FIG. 20.

(a) Density field, (b) b field, and (c) TKE field in the refined parameter space 3  S 0 7 μm and K 0 = 10 11 and 1012 cm2/s2 at t = 6.0 ns. The fourth column shows the corresponding 3D result and is repeated in each row.

Close modal
1.
D.
Hicks
,
N.
Meezan
,
E.
Dewald
,
A.
Mackinnon
,
R.
Olson
,
D.
Callahan
,
T.
Döppner
,
L.
Benedetti
,
D.
Bradley
,
P.
Celliers
et al, “
Implosion dynamics measurements at the national ignition facility
,”
Phys. Plasmas
19
,
122702
(
2012
).
2.
L.
Rayleigh
, “
Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density
,”
Proc. London Math. Soc.
1
,
170
177
(
1882
).
3.
G. I.
Taylor
, “
The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I
,”
Proc. R. Soc. London, Ser. A. Math. Phys. Sci.
201
,
192
196
(
1950
).
4.
R. D.
Richtmyer
, “
Taylor instability in shock acceleration of compressible fluids
,”
Commun. Pure Appl. Math.
13
,
297
319
(
1960
).
5.
E.
Meshkov
, “
Instability of the interface of two gases accelerated by a shock wave
,”
Fluid Dyn.
4
,
101
104
(
1969
).
6.
W.
Thomson
, “
XLVI. Hydrokinetic solutions and observations
,”
Lond. Edinburgh, Dublin Philos. Mag. J. Sci.
42
,
362
377
(
1871
).
7.
H. L.
von Helmholtz
, “
XLIII. On discontinuous movements of fluids
,”
Lond. Edinburgh, Dublin Philos. Mag. J. Sci.
36
,
337
346
(
1868
).
8.
B. M.
Haines
,
R. C.
Shah
,
J. M.
Smidt
,
B. J.
Albright
,
T.
Cardenas
,
M. R.
Douglas
,
C.
Forrest
,
V. Y.
Glebov
,
M. A.
Gunderson
,
C. E.
Hamilton
et al, “
Observation of persistent species temperature separation in inertial confinement fusion mixtures
,”
Nat. Commun.
11
,
544
(
2020
).
9.
D.
Besnard
,
F. H.
Harlow
,
R. M.
Rauenzahn
, and
C.
Zemach
, “
Turbulence transport equations for variable-density turbulence and their relationship to two-field models
,” Report No. LA-12303-MS (
Los Alamos National Laboratory
,
Los Alamos, NM
,
1992
).
10.
J. D.
Schwarzkopf
,
D.
Livescu
,
R. A.
Gore
,
R. M.
Rauenzahn
, and
J. R.
Ristorcelli
, “
Application of a second-moment closure model to mixing processes involving multicomponent miscible fluids
,”
J. Turbul.
12
,
N49
(
2011
).
11.
A.
Banerjee
,
R. A.
Gore
, and
M. J.
Andrews
, “
Development and validation of a turbulent-mix model for variable-density and compressible flows
,”
Phys. Rev. E
82
,
046309
(
2010
).
12.
N. O.
Braun
and
R. A.
Gore
, “
A multispecies turbulence model for the mixing and de-mixing of miscible fluids
,”
J. Turbul.
22
,
784
813
(
2021
).
13.
K.
Stalsberg-Zarling
and
R.
Gore
, “
The BHR2 turbulence model: Incompressible isotropic decay, Rayleigh-Taylor, Kelvin-Helmholtz and homogeneous variable density turbulence
,”
Report No. LA-UR-11-04773
(
Los Alamos National Laboratory
,
Los Alamos, NM
,
2011
).
14.
B.
Rollin
and
M.
Andrews
, “
On generating initial conditions for turbulence models: The case of Rayleigh–Taylor instability turbulent mixing
,”
J. Turbul.
14
,
77
106
(
2013
).
15.
V. N.
Goncharov
, “
Analytical model of nonlinear, single-mode, classical Rayleigh-Taylor instability at arbitrary Atwood numbers
,”
Phys. Rev. Lett.
88
,
134502
(
2002
).
16.
N. O.
Braun
and
R. A.
Gore
, “
A passive model for the evolution of subgrid-scale instabilities in turbulent flow regimes
,”
Physica D
404
,
132373
(
2020
).
17.
B. M.
Haines
,
F. F.
Grinstein
, and
J. D.
Schwarzkopf
, “
Reynolds-averaged Navier–Stokes initialization and benchmarking in shock-driven turbulent mixing
,”
J. Turbul.
14
,
46
70
(
2013
).
18.
D.
Holder
and
C.
Barton
, “
Shock tube Richtmyer–Meshkov experiments: Inverse chevron and half height
,” in
Proceedings of the Ninth International Workshop on the Physics of Compressible Turbulent Mixing
(
University of Cambridge
,
2004
).
19.
B. M.
Haines
,
F. F.
Grinstein
,
L.
Welser-Sherrill
, and
J. R.
Fincke
, “Simulations of material mixing in laser-driven Reshock experiments,”
Phys. Plasmas
20
,
022309
(
2013
).
20.
B. M.
Haines
,
F. F.
Grinstein
,
L.
Welser-Sherrill
,
J. R.
Fincke
, and
F. W.
Doss
, “
Simulation ensemble for a laser-driven shear experiment
,”
Phys. Plasmas
20
,
092301
(
2013
).
21.
K.
Flippo
,
F.
Doss
,
E.
Merritt
,
B.
DeVolder
,
C.
Di Stefano
,
P.
Bradley
,
D.
Capelli
,
T.
Cardenas
,
T.
Desjardins
,
F.
Fierro
et al, “
Late-time mixing and turbulent behavior in high-energy-density shear experiments at high Atwood numbers
,”
Phys. Plasmas
25
,
056315
(
2018
).
22.
F.
Doss
,
T.
Desjardins
,
A. M.
Rasmus
,
E. C.
Merritt
,
S.
Pellone
,
J.
Charonko
,
A.
Martinez
,
J.
Levesque
, and
C. A.
Di Stefano
, “
Density variance dynamics in disparate shock tubes
,”
Los Alamos National Laboratory Technical report No. LA-UR-23-20147
(
2023
).
23.
C. A.
Di Stefano
,
F. W.
Doss
,
A. M.
Rasmus
,
K. A.
Flippo
, and
B. M.
Haines
, “
The modeling of delayed-onset Rayleigh-Taylor and transition to mixing in laser-driven HED experiments
,”
Phys. Plasmas
26
,
052708
(
2019
).
24.
C.
Di Stefano
,
A.
Rasmus
,
F.
Doss
,
K.
Flippo
,
J.
Hager
,
J.
Kline
, and
P.
Bradley
, “
Multimode instability evolution driven by strong, high-energy-density shocks in a rarefaction-reflected geometry
,”
Phys. Plasmas
24
,
052101
(
2017
).
25.
A. M.
Rasmus
,
C. A.
Di Stefano
,
K. A.
Flippo
,
F. W.
Doss
,
J. L.
Kline
,
J. D.
Hager
,
E. C.
Merritt
,
T. R.
Desjardins
,
W. C.
Wan
,
T.
Cardenas
et al, “
Shock-driven discrete vortex evolution on a high-Atwood number oblique interface
,”
Phys. Plasmas
25
,
032119
(
2018
).
26.
A. M.
Rasmus
,
C. A.
Di Stefano
,
K. A.
Flippo
,
F. W.
Doss
,
C.
Kawaguchi
,
J. L.
Kline
,
E. C.
Merritt
,
T. R.
Desjardins
,
T.
Cardenas
,
D. W.
Schmidt
et al, “
Shock-driven hydrodynamic instability of a sinusoidally perturbed, high-Atwood number, oblique interface
,”
Phys. Plasmas
26
,
062103
(
2019
).
27.
S.
Pellone
,
C.
Di Stefano
,
A. M.
Rasmus
,
C. C.
Kuranz
, and
E.
Johnsen
, “
Vortex-sheet modeling of hydrodynamic instabilities produced by an oblique shock interacting with a perturbed interface in the HED regime
,”
Phys. Plasmas
28
,
022303
(
2021
).
28.
S.
Kurien
,
F. W.
Doss
,
D.
Livescu
, and
K.
Flippo
, “
Extracting a mixing parameter from 2D radiographic imaging of variable-density turbulent flow
,”
Physica D
405
,
132354
(
2020
).
29.
M.
Gittings
,
R.
Weaver
,
M.
Clover
,
T.
Betlach
,
N.
Byrne
,
R.
Coker
,
E.
Dendy
,
R.
Hueckstaedt
,
K.
New
,
W. R.
Oakes
et al, “
The rage radiation-hydrodynamic code
,”
Comput. Sci. Discov.
1
,
015005
(
2008
).
30.
S. P.
Lyon
, “
Sesame: The Los Alamos National Laboratory equation of state database
,”
Report No. LA-UR-92-3407
(
Los Alamos National Laboratory
,
Los Alamos, NM
,
1992
).
31.
J.
Marozas
,
M.
Hohenberger
,
M.
Rosenberg
,
D.
Turnbull
,
T.
Collins
,
P.
Radha
,
P.
McKenty
,
J.
Zuegel
,
F.
Marshall
,
S.
Regan
et al, “
First observation of cross-beam energy transfer mitigation for direct-drive inertial confinement fusion implosions using wavelength detuning at the national ignition facility
,”
Phys. Rev. Lett.
120
,
085001
(
2018
).
32.
I.
Igumenshchev
,
F.
Marshall
,
J.
Marozas
,
V.
Smalyuk
,
R.
Epstein
,
V.
Goncharov
,
T.
Collins
,
T.
Sangster
, and
S.
Skupsky
, “
The effects of target mounts in direct-drive implosions on omega
,”
Phys. Plasmas
16
,
082701
(
2009
).
33.
B. M.
Haines
,
G. P.
Grim
,
J. R.
Fincke
,
R. C.
Shah
,
C. J.
Forrest
,
K.
Silverstein
,
F. J.
Marshall
,
M.
Boswell
,
M. M.
Fowler
,
R. A.
Gore
et al, “
Detailed high-resolution three-dimensional simulations of omega separated reactants inertial confinement fusion experiments
,”
Phys. Plasmas
23
,
072709
(
2016
).
34.
E.
Harding
,
J.
Hansen
,
O.
Hurricane
,
R.
Drake
,
H.
Robey
,
C.
Kuranz
,
B.
Remington
,
M.
Bono
,
M.
Grosskopf
, and
R.
Gillespie
, “
Observation of a Kelvin-Helmholtz instability in a high-energy-density plasma on the omega laser
,”
Phys. Rev. Lett.
103
,
045005
(
2009
).
35.
O.
Hurricane
,
V.
Smalyuk
,
K.
Raman
,
O.
Schilling
,
J.
Hansen
,
G.
Langstaff
,
D.
Martinez
,
H.-S.
Park
,
B.
Remington
,
H.
Robey
et al, “
Validation of a turbulent Kelvin-Helmholtz shear layer model using a high-energy-density omega laser experiment
,”
Phys. Rev. Lett.
109
,
155004
(
2012
).
36.
J. D.
Bender
,
O.
Schilling
,
K. S.
Raman
,
R. A.
Managan
,
B. J.
Olson
,
S. R.
Copeland
,
C. L.
Ellison
,
D. J.
Erskine
,
C. M.
Huntington
,
B. E.
Morgan
et al, “
Simulation and flow physics of a shocked and reshocked high-energy-density mixing layer
,”
J. Fluid Mech.
915
,
A84
(
2021
).
37.
Y.
Yang
,
Q.
Zhang
, and
D. H.
Sharp
, “
Small amplitude theory of Richtmyer–Meshkov instability
,”
Phys. Fluids
6
,
1856
1873
(
1994
).
38.
R. L.
Holmes
,
G.
Dimonte
,
B.
Fryxell
,
M. L.
Gittings
,
J. W.
Grove
,
M.
Schneider
,
D. H.
Sharp
,
A. L.
Velikovich
,
R. P.
Weaver
, and
Q.
Zhang
, “
Richtmyer–Meshkov instability growth: Experiment, simulation and theory
,”
J. Fluid Mech.
389
,
55
79
(
1999
).
39.
J.
Canfield
,
N.
Denissen
,
M.
Francois
,
R.
Gore
,
R.
Rauenzahn
,
J.
Reisner
, and
S.
Shkoller
, “
A comparison of interface growth models applied to Rayleigh–Taylor and Richtmyer–Meshkov instabilities
,”
J. Fluids Eng.
142
,
121108
(
2020
).
40.
J. P.
Boris
,
F. F.
Grinstein
,
E. S.
Oran
, and
R. L.
Kolbe
, “
New insights into large eddy simulation
,”
Fluid Dyn. Res.
10
,
199
(
1992
).
41.
F. F.
Grinstein
,
L. G.
Margolin
, and
W. J.
Rider
,
Implicit Large Eddy Simulation
(
Cambridge University Press
,
Cambridge
,
2007
), Vol.
10
.
42.
Y.
Zhou
,
F. F.
Grinstein
,
A. J.
Wachtor
, and
B. M.
Haines
, “
Estimating the effective Reynolds number in implicit large-eddy simulation
,”
Phys. Rev. E
89
,
013303
(
2014
).