Properties of the Turbulent Mixing Layer in a Spherical Implosion

The implosion phase in inertial confinement fusion (ICF) [1] can be characterized by aggressive turbulent mixing between the hot fuel and the shell material. The turbulent mixing is driven by the Richtmyer–Meshkov (RM) [2–4] and Rayleigh–Taylor (RT) [5,6] instabilities, and can act to degrade the overall yield of the process. The RM instability is explained by the deposition of baroclinic vorticity as a result of shock impact of a corrugated interface that separates gases of different densities. The RT instability occurs when a light fluid is accelerated into a heavier fluid across a perturbed interface, where the acceleration can either be constant or time-dependent. To address these issues, there have been a vast array of studies of such canonical instabilities, although most of the efforts have been directed at the simpler, planar configuration. The effect of spherical convergence on the nature of such instabilities has been treated analytically, as well as through numerical simulations. We provide a brief review later. In this paper, we have investigated the behavior of turbulent mixing when the flow is constrained by a spherically confined geometry, as is the case with spherical implosions.

The effect of geometric convergence, due to either a cylindrical or spherical geometry, on the growth of the RT instability was illustrated by Bell and Plesset [7,8], and collectively referred to as Bell–Plesset (BP) effects. At high convergence ratios, BP effects that are also termed “undriven growth” can significantly enhance the development of the underlying RT instability. In a series of papers [9–11], Mikaelian investigated the linear growth rates of Rayleigh–Taylor and Richtmyer–Meshkov instabilities in stratified cylindrical and spherical shells, using a potential flow approach. For cylindrical convergence, the BP effects are separable from the growth of driven instabilities such as RM/RT [11], while more complex expressions for the growth rate were obtained for the corresponding spherical problem [10]. Lombardini et al. [12,13] developed a linear model for instability growth in both cylindrical and spherical geometries, based on the RM impulsive model, in the limit of RM growth driven by a self-similar shock [14,15]. The effect of viscosity on the linear growth of an RT-unstable spherical interface was examined by Terrones and Carrara [16] in the limit of infinite density ratio between the fluids. Finally, Ramshaw [17] developed a kinetic energy based model coupled with a wavelength renormalization approach to obtain expressions for both the linear and nonlinear growth rates of interfaces in spherical geometry.

Youngs and Williams [18] defined an implosion problem in spherical geometry, as an idealization of the ICF process. In particular, Youngs and Williams [18] used the Lagrange-remap code turmoil [19,20] to perform high resolution, three-dimensional (3D) numerical simulations of the implosion in a spherical wedge, with multimode perturbations prescribed at the interface between the gases. Joggerst et al. [21] performed a code comparison study to examine the effect of different meshes on instability growth at various stages of the implosion defined by Youngs and Williams [18]. They found that the growth rates of short wavelength modes in particular were affected by the choice of mesh, and a grid that is aligned with the interface was less susceptible to the rise of spurious modes at late times. Lombardini et al. [13,22] used large eddy simulations to investigate RM-driven turbulent mixing in a spherical geometry, where the drive provided by a self-similar shock wave as defined in Refs. [14] and [15]. In this paper, we report the results from high-resolution simulations of the implosion problem [18] using the flash [23] code, and describe in detail the behavior of several turbulent quantities at various stages of the implosion history.

The rest of the paper is organized into numerical details and problem definition (Sec. 2), results and discussion (Sec. 3), and summary and conclusions (Sec. 4).

where $Rbdt$ refers to the time-dependent trajectory of the boundary, $R0$ is its initial location, and $ubd$ is a constant (0.2 s⁻¹). Figure 1 is a cross section of the spherical wedge $(π/2−π/8≤(θ,ϕ)≤π/2+π/8)$ investigated here, which shows the arrangement involving the three fluids, along with the corresponding dimensions of the problem (in cgs units). Periodic boundary conditions were employed in the lateral directions, while the surface at r = 15 cm was treated as an outflow boundary to allow the reflected rarefaction (and other acoustic waves) to exit the computational domain without triggering spurious reflections. A reflecting boundary is placed at r = dr, so that the singularity at r = 0 is excluded from the computational domain. The fluid properties are listed in Table 1. We use $R0$ as a length scale, the initial interface jump velocity $ΔU$ as a velocity scale, and $R0/ΔU$ as the time scale for normalizing the variables.

The interface between the inner and shell materials was initialized with a multimode perturbation suggested by Youngs and Williams [18] (consisting of a superposition of cosine waves), with wavenumbers confined to $kmin=π≤(kθ,kφ)≤kmax=π/2dr$⁠. The initial amplitudes and phases are randomized with a Gaussian distribution, thereby avoiding local pileups of peaks and valleys. Following Ref. [18], the power spectrum of the perturbation field satisfies $p(k)=Ck−2$ for $kmin<(kθ,kφ)<kmax$⁠, and $p(k)=0$ otherwise. The constant C is chosen so that $s.d.2=∫p(k)dk$⁠, where $s.d.=0.005$ is the standard deviation of the perturbation field. Surface perturbations within a cell volume are converted to corresponding perturbations of the mass fraction Y. When more than a single fluid is present in a cell, the fluids are assumed to be numerically mixed leading to a corresponding cell mass fraction. Thus, we do not explicitly track the interface between the fluids, an approach that has been shown to be successful in the handling of flows with sharp discontinuities and shocks.

The resulting perturbation field is shown in Fig. 2. The radial cell containing the perturbation is further subsampled to adequately resolve perturbation growth during the linear stages. In Section 3, we present the results from three simulations with grid resolutions in (θ, φ, r) of $128×128×384$⁠, $256×256×786$⁠, and $512×512×1536$ zones. The adaptive mesh refinement mesh is locally defined to be 64 × 64 × 192 in the region r < 1.5 cm to avoid the significant computational costs associated with tracking the convergent shock in that region at late times. From an examination of results from the two approaches, we found that the radial profiles of pressure, density, and other quantities are in excellent agreement, and the maximum discrepancy in local values does not exceed 10% at any radial location and at any time.

We verify our numerical problem setup by comparing results from a one-dimensional (1D), unperturbed simulation (768 zones in r-) with the corresponding data from Ref. [18]. Figure 3 is a plot of the trajectories of the inner/shell interface, as well as the external boundary between the shell and outer material which follows Eq. (1). We also track the location of the driving shock in each case through local maxima of the pressure gradients. Thus, for τ ≤ 0.24, we plot the trajectory of the incident shock, while we track the transmitted and reflected shocks from the origin for τ > 0.24. The unperturbed, interface location from flash is compared with the corresponding results from Ref. [18], and shows excellent agreement in Fig. 3. Following the incident shock impact at τ = 0.24, the interface is initially RM-unstable, while the large initial radius likely implies that the effects of convergence will be limited at this stage. For the initial heavy → light interaction, the Atwood number associated with the inner/shell interface is −0.905. The shock impact results in a transmitted shock and a reflected rarefaction, while the interface appears to coast with a nearly constant radial velocity. However, the transmitted shock is reflected from the origin and shocks the collapsing interface several times in a light → heavy configuration. These interactions cause the interface to stagnate at τ ∼ 1.2, followed by a slight rebound. The flow development at this stage is complex and involves stages of RM instability immediately following each shock, interspersed with durations of RT instability subjected to variable acceleration. For the flow conditions investigated here, a radial convergence ratio of ∼4 is observed in the flash simulations. As shown in Ref. [21], the use of an aligned coordinate mesh in flash here eliminates the formation of spurious, grid-seeded modes even at high convergences.

The simulations were repeated in 3D with multimode perturbations described in Sec. 2, and at the three different mesh resolutions to establish grid convergence. In Fig. 4, we plot the time history of the extents of the mixing layer defined as the r-locations where the planar-averaged mass fractions achieve values of 1% and 99%. The flash results shown are at the maximum mesh resolution employed in this study (⁠$512×512×1536$⁠) and compared with the corresponding data from Ref. [18]. From Fig. 4, it is clear that the initial RM phase of evolution produces scant growth, since the small amplitudes of the imposed perturbations ensure that modes remain in linear growth at this stage. Significant amplitude growth is observed beginning at around τ ∼ 0.96 following the first reshock event, which re-energizes the nonlinear mixing layer. This is followed by aggressive growth of the mixing layer driven by a combination of reshock-driven RM growth, variable acceleration RT growth as well as BP effects due to the high convergences at late times.

Figures 5(a)–5(e) represent isosurfaces of mass fraction of the inner fluid satisfying the condition <Y_inner> = 0.5, where the <•> denotes planar-averaging along (θ, φ). At τ = 0.87, the interface growth favors shorter wavelengths, which have the highest linear RM growth rates. The shape of the mass fraction isosurface immediately following reshock is illustrated in Fig. 5(b), and shows that once individual modes have reached nonlinear saturation (h ∼ λ), larger wavelengths are formed through mode-coupling. By τ = 1.08 (Fig. 5(c)), we observe the onset of cross-modal interactions and mode coupling which will eventually lead to mixing. Figures 5(d) and 5(e) show the turbulent mixing layer at late times, when multiple reflected shocks have passed through the interface. At these late times, the isosurface is characterized by significant levels of mixing, and asymmetry between the bubble and spike fronts. The bubble-spike asymmetry is a result of the large density contrasts between the fluids, but also due to convergence effects which favor inward pointing spikes.

We plot the angular-averaged mass fraction (inner fluid) profiles <Y_inner> across the mixing layer in Fig. 6(a), corresponding to the times shown in Fig. 5. When the radial coordinate is scaled as shown in Fig. 6(a), the profiles are anchored at (r − r_c)/W ∼ 0.5. At late times, the flow appears to settle down to a self-similar state, evidenced by the collapse of the mass fraction profiles for τ > 0.96. As the reflected shocks repeatedly impact the flow at late times, the underlying self-similarity of the flow does not appear to be affected. The corresponding product <Y_innerY_shell> is plotted as a radial function in Fig. 6(b). At late times, <Y_innerY_shell> peaks at the centerline at ∼0.18, where a value of 0.25 would indicate perfect mixing. Once again, the profiles collapse at late times (τ > 0.96) suggesting that the flow reaches a self-similar state when subjected to mixing from repeated shocks and RT.

For RM flows, the amplitude of the mixing layer is typically defined as $r99%‐r1%$⁠, where $r99%$ (⁠$r1%$⁠) denotes the radial location of the surface <Y_inner> = 99% (1%). In Fig. 7, we plot the time evolution of the mixing layer amplitude from flash simulations with mesh resolutions of $256×256×768$ and $512×512×1536$⁠. The results are compared with the corresponding data from Ref. [18] and show good agreement. Consistent with the images of the isosurfaces discussed earlier, the growth rate of the mixing layer falls in the linear regime until the first reshock event τ ∼ 0.96. Significant nonlinear growth and mixing is observed at late times, due to the combination of RM/RT/BP effects. In Fig. 7, we also show results from the simulations at lower resolutions (128 zones and 256 zones). At late time, following reshock, the simulations show good mesh convergence and are in agreement with the results from Ref. [18]. At early times, the simulation with 128 zones (and to a lesser extent the 256 zone calculation) shows a larger mixing width compared with the highest-resolution case. At these early times, the growth of the mixing width is due to the linear growth of individual modes, so that the larger width observed in the poorly resolved (128 zone) case is attributable to the numerical diffusion dominating over the linear growth.

where ℘(t) and δ(t) represent the actual chemical product thickness and the maximum thickness product formed with complete mixing, respectively [22]. Thus, the measure in Eq. (2) compares large-scale stirring with small-scale mixing that is associated with a numerical diffusivity. Note that Θ(t) is related to the intensity of segregation I(t) through Θ(t) = 1 − I(t), and can change from 0 → 1 to indicate variation from segregated to completely mixed states. In Fig. 8, we plot the evolution of the atomic mix parameter Θ(t), as well as ℘(t) and δ(t) from our highest resolution simulations. The high values of Θ(t) observed at early times (τ < 0.96) during which the linear mixing layer exhibits little mixing are an artifact of the discrete representation of the initial interface in our simulations [$512×512×1536$]. At early times, when the interface is confined within a cell, the only contribution to the integral in Eq. (2) is from the presence of both fluids in a given cell volume [25]. After reshock (τ > 0.96), the nonlinear growth spreads the interface across multiple mesh zones, so that the atomic mix becomes a more reliable indicator of the fraction of molecularly mixed fluid within the mixing layer. At late times, the atomic mix parameter approaches a value of ∼0.7, consistent with earlier findings of mixing in RT and RM flows.

The behavior of the widths ℘(t) and δ(t) from the lower resolution simulation (128 zones) is consistent with the trend observed in Fig. 7 for the mixing layer amplitude—discrepancies at early time due to the dominance of numerical diffusion over linear growth of modes, and convergence to the high-resolution results at late times when the flow is dominated by nonlinear growth and turbulence. Simulations at the lower mesh resolution show the same qualitative trend in the behavior of the atomic mix parameter, while the transition from the initial (artificial) Θ value occurs at different times in each case. Specifically, the higher resolution simulations show this transition earlier, since these calculations had much smaller mesh sizes that the nascent flow can outgrow at an earlier time. When this occurs, the value of Θ transitions from the artificial initial value to a measure that is reflective of the fraction of atomically mixed fluid within the flow.

where $vivj$ are components of the Reynolds stress tensor. For isotropic turbulence, B_ij = 0, while strong anisotropy suggests B_ij → 2/3 in the dominant direction (and B_ij → −1/3 in the other directions) [26]. The diagonal components of the anisotropy tensor are plotted as a function of the nondimensional radial coordinate at stagnation (Fig. 10(a)) and at the end of the simulation (Fig. 10(b)). Both plots suggest strong anisotropy in the radial direction with B_rr ∼ 0.27±0.01 and 0.38±0.01 at τ ∼ 1.22 and τ ∼ 1.44, respectively. In contrast, the angular components are restricted to 18.5% of the total energy at those times. However, the levels of anisotropy observed here are lower than in other directed flows such as RT [25] and could be due to the repeated realignment (scrambling) of the mixing layer from multiple reshock events.

In variable density flows such as RT and RM turbulent flows, quantities such as the normalized mass flux $ai$ and the density specific volume correlation parameter $b$ are important in understanding turbulent transport within the mixing layer [26,27]. In particular, in the non-Boussinesq limit, these quantities can develop asymmetries across the mixing layer, in proportion to the density contrast between the fluids. The normalized mass flux in the ith direction is defined as $ai=<ρ′vi′>/<ρ>$ and has been recognized to be a critical quantity in the conversion of potential energy to kinetic energy in buoyancy-driven flows [26,27]. In Fig. 11(a), we plot the cross-stream profiles of the normalized mass flux in the radial direction (⁠$ar$⁠), for the two highest resolved cases at the end of the simulation (τ ∼ 1.44). The obtained $ar$ profiles are asymmetric and skew toward the spikes, suggesting greater mass flux associated with the spikes due to the large density differences. This velocity increase on the spike side is also compounded by greater geometric convergence experienced by the spikes, which are directed inward toward the origin.

In Refs. [22] and [27], the authors show that the density specific volume correlation parameter $b=‐ρ′1/ρ′$ can modify the production term in the radial mass flux equation. In Fig. 11(b), we plot the radial profile of the b parameter at τ ∼ 1.44, and at the two highest mesh resolutions. The peak values obtained here (∼0.28) are similar to values reported by Lombardini et al. [22] for their heavy → light simulation in a spherical geometry. Once again, the profiles are skewed toward the light side of the mixing layer, suggesting that the flow is highly non-Boussinesq and requiring the solution of variable density equations in modeling transport [27,28].

where “ $̂$ ” indicates a two-dimensional Fourier transform operation performed on a (θ, φ) plane of statistically homogeneous data, while “*” represents the corresponding complex conjugate. Thus, the two-dimensional fast Fourier transform operations were performed on (θ, φ) planes corresponding to the radial location $r|Yinner=0.5$⁠, while the double integral in the above equation indicates averaging over a circular region in wavenumber space, thus resulting in a 1D kinetic energy spectrum. In Fig. 12, we plot energy spectra from three critical moments in the flow, viz., following the initial shock as modes achieve nonlinearity (Fig. 12(a)), immediately following reshock (Fig. 12(b)) and at late time (Fig. 12(c)). In each case, the kinetic energy spectra are compensated by the coefficient k^5/3 to highlight the inertial range, while the wavenumber axis is nondimensionalized with R(t), the radius of the corresponding unperturbed interface at that time.

The energy spectrum shown in Fig. 12(a) was obtained at τ = 0.77, following the incident shock and still shows the imprint of the initial conditions. At this early time, a narrow band of modes appears to have evolved to nonlinearity (h ∼ λ), while the long wavelengths have not filled in. At τ = 0.95 (Fig. 12(b)), the appearance of an inertial range is evident as the flow has clearly transitioned to turbulence following the reshock event. However, the inertial range at this time is narrow and represents a flow developing into an eventual self-similar state. Figure 12(c) is obtained from the $r|Yinner=0.5$ surface at late time (τ = 1.44) and shows an inertial range spanning a decade with a slope approaching k^−5/3. The broadband nature of the spectrum and the existence of a substantial inertial range suggest that the flow exhibits a self-similar state at these late times. From the kinetic energy spectra, we determine the inner viscous length scale λ_ν as the intersection point between the fiducials corresponding to the inertial and the dissipation ranges. For τ = 1.44, we find k_ν R(t) ∼ 220, where k_ν = 2π/λ_ν and R(t) is the location of the corresponding unperturbed interface. The outer length scale is given by δ ∼ W(t), where W(t) is the amplitude of the turbulent mixing layer, so that the effective Reynolds number is obtained as follows [30–33]: $Re=50δλν4/3∼1.1×10$⁴. A similar analysis conducted at τ = 0.95 yields k_ν R(t) ∼ 410 and Re ∼ 1.25 $×10$³.

We have discussed results from numerical simulations of a spherical implosion of relevance to the ICF application. The problem statement was originally defined by Youngs and Williams [18] and also serves as a benchmark for evaluating different numerical algorithms. The simulations were performed here with the widely used flash code, for two different conditions: unperturbed interface and an interface initialized with multimode perturbations. The flow conditions and initialization were chosen to match parameters in Refs. [18] and [21], while the problem setup had to be modified to introduce an outer layer of fluid that could serve to support the driving shock. Thus, the interface between the outer fluid and shell material in the Eulerian simulations discussed here serve the purpose of the moving boundary in the ALE code turmoil [19,20]. Nevertheless, we find several quantities of interest including the radial trajectory of the unperturbed interface and the mixing width of the turbulent flow to be in very good agreement with the results from Ref. [18].

The implosion of a spherical interface is susceptible to multiple instabilities, resulting in a complex evolution over time. The initial growth is driven by the convergent shock and can be attributed to RM growth of linear modes. At late times, the nonlinear mixing layer is repeatedly impacted by shocks reflected at the origin, interspersed with stages of RT-driven growth. This is also reflected in the atomic mix parameter, which saturates to ∼0.7 at late times following the interactions of the interface with the reflected shocks. In addition, at the convergence ratios achieved here, BP effects are expected to play a role, thus contributing preferentially to the growth of the spike front. At late times, the flow evolves to higher levels of anisotropy, with much of the turbulent kinetic energy residing in the radial component.

In the future, we plan to extend this work by comparing results from the flash simulations with turbulent mix models [9,27,28]. Future studies will also include detailed investigations on the effects of the convergence ratio, the initial conditions, and the shock strength on the development of the turbulent mixing layer properties discussed in this paper.

This work was supported in part by the (U.S.) Department of Energy (DOE). FLASH was developed by the DOE-sponsored ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. The authors are grateful to Mr. Karkhanis for assistance in preparing the figures.

• •

  Los Alamos National Laboratory (Grant No. DE-AC52-06NA2-5396).

B_ij =

component of anisotropy tensor

e =

internal energy

h =

amplitude of perturbation

$I(t)$ =

intensity of segregation

$k$ =

wavenumber of the perturbation

p =

pressure

$P(k)$ =

power spectrum of perturbation field

r =

radial coordinate

$Rbd$ =

radial location of interface between shell and outer fluids

$R0$ =

initial radial location of interface between inner and shell material

s.d. =

standard deviation of random perturbation field

$ubd$ =

radial velocity of boundary

W =

mixing width

Y =

mass fraction

γ =

adiabatic index

$δ(t)$ =

maximum mixing product width

θ =

azimuthal angular coordinate

$Θt$ =

atomic mix

ρ =

fluid density

$τ$ =

nondimensional time scale

φ =

zenithal angular coordinate

$℘t$ =

total product width