Design optimization for Richtmyer–Meshkov instability suppression at shock-compressed material interfaces

The interface between two distinct fluids of different densities is generally conducive to the generation of unstable fluid motion. For the specific case in which this interface is impulsively accelerated (e.g., from the passage of a shock wave), this unstable fluid motion is referred to as the Richtmyer–Meshkov instability (RMI). Under these conditions, perturbations in the interface begin to grow due to the RMI, which leads to the propagation of jets or spikes of the denser substance into the fluid of lower density. RMI is prevalent both in natural processes, such as supernova explosions,^1,2 and in technological settings, such as inertial confinement fusion (ICF) applications.^3–7 The burning plasma that was recently produced by Zylstra et al.⁷ in their landmark study represents a major step forward in the field of ICF. However, there are still several major challenges involved in attaining self-sustaining fusion energy. One major challenge encountered in ICF experiments, which generally use high-pressure shock waves to compress a fuel target, is asymmetry from mixing of the fuel with the surrounding capsule material during the implosion process. This mixing is often caused by instabilities, such as RMI, which can be seeded by imperfections in the fuel capsule. The mixing of fuel with the surrounding capsule material can significantly reduce the energy yield of these experiments. Thus, the suppression of perturbation growth caused by RMI is an important, yet challenging, goal.

In many applications, the denser substance at the interface of interest may be an elastic–plastic material (i.e., a solid). Under high stresses and strain rates, a solid can significantly deform and exhibit fluid-like behavior, which includes the development of instabilities such as RMI. For RMI at interfaces that initially involve a solid, the elasticity and yield strength of the solid can play a major role in the propagation of RMI spike growth. At a certain point, material elasticity and strength can stabilize the RMI and arrest further spike growth.⁸ Under certain circumstances, RMI spike growth may continue until these spikes rupture and form ejecta. The extreme pressure and temperature changes that occur during shock-wave compression can also lead to phase transitions (e.g., melting), which may further complicate unstable RMI behavior. Several experimental studies have investigated RMI involving a solid material to study material properties at extreme pressures and strain rates.^9–12 Theoretical and computational studies have also been conducted to provide a mechanistic understanding of RMI involving elastic–plastic materials. The landmark work of Piriz et al.⁸ describes an analytical model for characterizing RMI growth at solid–vacuum interfaces. Their model provides scaling laws for predicting the maximum amplitude growth of the RMI spikes and showed that this amplitude is inversely proportional to the yield stress of the solid material. Dimonte et al.¹³ built on the work of Piriz et al. by developing a scaling law for predicting RMI spike amplitude growth. This scaling law can also be used to infer the yield stress of a material by measuring the RMI spike amplitude growth. The more recent work of Chen et al.¹⁴ provides an analytical description of the dependency of elastic–plastic RMI perturbation growth on the Atwood number, a dimensionless number that describes the relative densities of the materials at the interface. (We describe the Atwood number in more detail in Sec. II of this paper.)

These experimental and theoretical studies help to provide a mechanistic understanding of RMI involving elastic–plastic materials occurring under various conditions. However, an effective method for suppressing RMI perturbation growth remains difficult to achieve. Several methods have previously been proposed for suppressing RMI perturbation growth. These include using magnetic fields to stabilize the RMI and prevent perturbation growth^15–17 when a plasma is present at the interface. Alternatively, Sano et al.¹⁸ propose using a density transition layer at the interface to suppress RMI. This density transition layer provides a continuous transition between the densities of the materials at the interface and appears to be effective at suppressing perturbation growth if a sufficiently thick transition layer is used. These methods for suppressing RMI perturbation growth tend to require specific conditions to be present and are not necessarily applicable to all applications involving RMI. Several recent studies, such as those of Henry de Frahan et al.,¹⁹ Feng et al.,²⁰ and Liang and Luo^21–23 have investigated RMI in scenarios involving several layers of different density fluids with coupled interfaces. These studies have focused on analyzing the effect of upstream interfaces with perturbations that create reverberating shock waves inside the fluid layers and have found that RMI and other hydrodynamic instabilities can be enhanced or suppressed in this manner.

For this study, we propose using a shape optimization approach²⁴ to reduce or suppress the RMI spike growth of a solid target material (with an initially perturbed surface in contact with atmospheric air) undergoing shock-wave compression. Our design optimization focuses on modifying the geometry of the target that is upstream of the initially perturbed interface where the primary RMI spike growth occurs. Note that we do not simply alter the geometry of the initially perturbed interface, as it may not be possible or practical to alter this interface in various engineering applications. For instance, this perturbed interface may arise in an engineering scenario from a limitation in manufacturing methods or may have desirable properties for a specific application. Our computational study demonstrates that reduction and partial suppression of RMI spike growth is achievable under specific conditions. The optimized designs that we have developed suppress RMI perturbation growth by intentionally generating another RMI that occurs upstream of the initially perturbed interface, which appears to counteract unstable spike growth at the interface of interest. We refer to this intentionally produced RMI as the rectifying RMI. We also investigate the underlying dynamics of why our optimized designs are effective at suppressing RMI spike growth by considering shock-wave behavior (such as reverberating shock waves) and vorticity deposition in the target material.

In this section, we describe several characteristics of the RMI that are important for understanding the methods that we use for suppressing this instability. We consider the case where the impulsive acceleration that generates the RMI is caused by an incident shock wave moving in a normal direction to the initially perturbed interface (i.e., the x-direction in Fig. 1). One common description of the mechanism that leads to RMI perturbation growth comes from considering vorticity deposition $ω→$ as the incident shock wave refracts at the perturbed interface. Note that vorticity represents the curl of the velocity vector (i.e., $ω→=∇×u→$⁠). This vorticity deposition generates counterclockwise (positive) vorticity on the one side of an interface perturbation and clockwise (negative) vorticity on the other side of the perturbation (see diagrams in Fig. 1).^{25,26,45–47} This vorticity imbalance causes the interface perturbations to deform, leading to RMI perturbation growth. For single-mode perturbations, this perturbation growth generally occurs as long narrow jets or spikes with a rounded bubble structure in between these spikes as shown in Fig. 1.

The general behavior of RMI perturbation growth at an interface is strongly dependent on a dimensionless ratio known as the Atwood number. The Atwood number quantifies the difference in density between the two interface substances prior to shock loading and is defined as^14,26

where ρ₁ is the density of the upstream material through which the shock wave is initially moving and ρ₂ is the density of the other substance or material at the interface. If the magnitude of the Atwood number is close to unity (i.e., $|A|≈1$⁠, which is the case for solid–vacuum interfaces), this indicates a large difference in the density of the two substances or materials. If the magnitude of the Atwood number is close to zero (i.e., $|A|≈0$⁠), the densities of the materials are approximately equal and this generally represents a more stable configuration. The developmental behavior of RMI spike growth is dependent on the Atwood number. For instance, if A > 0, then the interface perturbation spikes generally grow in the x-direction toward the lower-density material as shown in Fig. 1(a). However, if A < 0, a phase inversion is likely to occur as the perturbations begin to deform as shown in Fig. 1(b) and it is these phase-inverted spikes that grow into the lower-density material.²⁶ For our analysis, we focus on suppressing RMI spike growth in scenarios where the density of the upstream material ρ₁ is significantly greater than the density of the other substance at the interface ρ₂ (i.e., scenarios where $A≈−1$⁠).

Additionally, the shock impedance of the materials at the interface plays a crucial role in the shock-wave behavior and corresponding vorticity deposition. The shock impedance of a given material is defined as $Z=ρius$⁠, where ρ_i is the initial density (prior to shock loading) of the material and u_s is the shock velocity magnitude within the material. According to Nellis et al.,²⁷ if we define Z₁ as the impedance of the upstream material (through which the shock wave is initially moving) and Z₂ as the impedance of the downstream material, then the following is true. If Z₁ < Z₂, then the shock wave is at least partially reflected when it reaches the interface between the two materials and reverberates back through material Z₁ at higher pressure. Otherwise, if Z₁ > Z₂, then a rarefaction wave is produced in material Z₁ when the shock wave reaches the interface and the pressure is isentropically released. This concept will be important in understanding the shock-wave behavior and corresponding vorticity deposition for our results in Sec. III.

The simulation setup that we consider in our analysis involves a high-velocity copper impactor striking a stationary copper target to generate a shock wave that passes through the target. We generally use an impactor velocity of $Vi=$ 0.2 cm/μs in our simulations. The copper target is in contact with air on the side of the target opposite to the side where the impact occurs [see diagram in Fig. 2(a)] and has an Atwood number of $A≈−1$⁠. For convenience, we will henceforth refer to the side of the target struck by the impactor as the left-hand side and to the opposite side of the target as the right-hand side. The right-hand side of the target includes non-smooth perturbations along the interface with the air. For these initial perturbations, we have chosen to use a periodic single-mode sinusoidal profile. During compression, these perturbations first invert phase and then begin to grow after the shock passes through the copper–air interface on the right-hand side of the target. The conditions chosen for these simulations are similar to gas-gun dynamic-compression experiments being conducted at the High Explosives Application Facility (HEAF) at Lawrence Livermore National Laboratory (LLNL) and the Dynamic Compression Sector user facility at Argonne National Laboratory.^28,29 These experiments are currently ongoing and will be published in a future article.

The dimensions of the copper impactor and target in our simulations are listed in Table I. Note that in our simulations, we use N_r = 10 perturbation periods on the right-side of the target and zero-velocity (rigid wall) boundary conditions in the y-direction on the top and bottom sides of our simulation domain. While we could also have used periodic boundary conditions with N_r = 1 in our simulations, we wanted to investigate potential edge effects. However, we found the edge effects to be largely negligible in our simulations (see simulation density plot in Fig. 3). For this reason, we only show a reduced number of perturbation periods (e.g., N_r = 2) for our simulation results in Sec. III. For the left and right sides of our simulation boundaries, we use zero-velocity boundary conditions in both the x- and y-directions.

Our two-dimensional simulations are constructed using the hydrodynamics code MARBL,^30,31 developed at LLNL, which numerically solves the conservation of mass, momentum, and energy equations. MARBL is an arbitrary Lagrangian–Eulerian code that uses a high-order finite element approach to increase accuracy of the numerical calculations.³² For the materials in our simulations, we use the Livermore Equation of State (LEOS) library, which was also developed at LLNL. The material equations of state for copper, air, and lucite (a low-density plastic resin that we use in some of our simulations) come from LEOS tables 290,³³ 2260, and 5070, respectively. They are based on the quotidian equation of state with extension formalism.^34,35 For copper, we use the Steinberg–Guinan strength model³⁶ with a shear modulus of 0.477 Mbar, an initial yield stress of $1.2×10−3$ Mbar, and a maximum work hardening stress of $6.4×10−3$ Mbar. We neglect the strength of lucite in our simulations.

For the hydrodynamics simulations used in our study, we use a largely uniform mesh of approximately 83 000 elements with a grid spacing of $Δx=0.0208$ cm in both the x- and y-directions. We found a maximum difference of only approximately 2% between the RMI spike growth at this mesh resolution and at higher resolutions, which indicates that our simulations are largely converged. In  Appendix A, we also compare simulation results using this setup to the scaling law developed by Dimonte et al.¹³ to demonstrate that the material models we use (e.g., the strength and failure models for the copper impactor and target) are accurately capturing the RMI behavior in the high-pressure regimes that we consider.

Now that we have described our simulation setup, we outline the methodology that we use for reducing and suppressing RMI spike growth. For our study, we focus on investigating designs that involve modifying the geometry of the left-hand side of the target. This allows us to tune the dynamics of the shock wave and response of the target material to the shock wave. The geometry of the right-hand side of the target, which includes the initial perturbations that seed the RMI, remains unchanged in our simulations. As previously discussed, it may be necessary or unavoidable for a specific engineering application to have these perturbations, but highly undesirable for RMI spike growth to occur because of these perturbations.

One promising parameterization that we have developed is to use piecewise cubic Hermite interpolating polynomials (PCHIP) to shape the geometry of the left-hand side of the target. For this parameterization, we select several equidistant adjustment points along the axis parallel to the left-hand side of the target (i.e., the y-direction) $yk, k=0,1,2,…,N$ as well as corresponding x-direction values at these points $f(yk)=fk$⁠. Between these adjustment points, we interpolate using a polynomial equation of the following form:

where $t=(y−yk)/(yk+1−yk)$ and $Hij(t)$ are Hermite basis functions and are defined as

In Eq. (2), the $fk′$ term represents the slope at each of the points y_k. For our application, we have chosen to preserve monotonicity between the y_k points using the method developed by Fritsch and Carlson³⁷ and set this slope equal to

where $hk=yk+1−yk$ and $dk=(xk+1−xk)/hk$ represent the internal slopes at each point y_k. The variables w₁ and w₂ are weighting factors and are set equal to $w1=2hk+hk−1$ and $w2=hk+2hk−1$⁠. If the slopes d_k and $dk−1$ are opposite in sign or one of them is equal to zero, then $fk′$ is simply set equal to zero.

Using the selected adjustment points y_k, the corresponding values f_k, and interpolating between these points using Eqs. (2) and (4), we can construct a half wavelength profile, which can then be mirrored to produce a full wavelength profile. This profile is then repeated to produce a periodic interface on the left-hand side of the target, which we refer to as a PCHIP profile. Figure 4 shows what a single period of a potential PCHIP profile may look like for four y_k points. This allows us to create a periodic wave profile that can achieve many different wave profile shapes that depend on the values chosen for f_k. The PCHIP profile also preserves monotonicity between the y_k points and is first derivative continuous. An alternative parameterization method that we have also tried is to interpolate between the adjustment points y_k using simple linear interpolation rather than PCHIP polynomials. Piecewise linear interpolation produces a less complex profile on the left-hand side of the target that is simpler to manufacture using standard machine tools (e.g., for an experiment) because it only uses straight lines to produce the profile. We show an optimized result in Sec. III that uses piecewise linear interpolation in addition to results that use PCHIP interpolation.

Note that these profiles create gaps or cavities between the target and the impactor (see Fig. 4) that can be filled with either air or another material. We typically fill these cavities in our simulations using lucite, a low-density plastic resin, that helps to maintain numerical stability. However, the difference in the final results between using air or lucite to fill these cavities is generally small.

Using the PCHIP profile (or alternatively the piecewise linear profile), we can select several parameterizations to generate our design space. For our study, we adjust four values [i.e., $f(y0), f(y1), f(y2)$⁠, and $f(y3)$] and the number of periods on the left-hand side of the target N_l (or equivalently the wavelength of the repeated profile). Note that we equidistantly space the four points y_k such that they span an entire half period. This creates a total of five parameters for our design space that can be optimized to reduce or suppress RMI spike growth.

For the design optimization analysis, we must first define a performance metric or objective function that we intend to minimize such that RMI spike propagation is reduced or suppressed. One objective function that we have selected for this purpose is

where N_r is the number of periods of the initial perturbations on the right-hand side of the target, $hsp,i$ is the amplitude of the ith RMI spike, and $hbu,i$ is the amplitude of the ith RMI bubble. The initial locations of the tracer particles that we use to determine $hsp,i$ and $hbu,i$ are shown in Fig. 4. The objective function in Eq. (5) represents an average amplitude difference between RMI spikes and bubbles across all perturbation periods and is useful for ensuring that the right-hand side of the target is as flat as possible at a specified time.

For certain applications, it may be desirable to ensure not only that RMI spike propagation is suppressed but also that the target remains as compact as possible. For this purpose, we have also investigated another objective function that is useful for maintaining target compactness by quantifying the deviation from the target's center-of-mass, which is computed via

where ρ denotes the material density at point x and $x*$ denotes the center-of-mass. The integral is computed only over the target material Ω_T. We use a value of p = 6, which can be considered a soft-max function (i.e., a smooth estimate of the maximum deviation of material about the center-of-mass). A diagram of the properties that we use in the objective function of Eq. (6) is illustrated in Fig. 5. This new metric does not utilize tracers as in the objective function of Eq. (5), but rather it requires images of material density to be saved to disk and post-processed. This post-processing involves determining the x location for some finite-sized volume of the target and the corresponding density $ρ(x)$ of this volume. This is repeated to cover the entire target volume Ω_T, which allows us to calculate the integral of Eq. (6).

The objective of our optimization is to find the design that minimizes the objective function [i.e., either Eq. (5) or Eq. (6)] at a specific time $Δt$ after the impactor strikes the target or minimizes the maximum objective function value over a specified time duration (e.g., between 0 and 7 μs after impact). We have somewhat arbitrary selected $Δt=7$ μs as the final time in our analysis. This allows sufficient time for the shock wave to pass through the perturbed interface on the right-hand side of the target and additional time for RMI spike growth to occur. Another reason that we have chosen to set the final simulation time equal to $Δt=7$ μs is because the in-progress independent experiments on which our simulation setup is based generally take measurements up until approximately 7 μs after the impactor strikes the target.

Our design optimization problem is largely defined by a mapping between the parameterization of Eq. (2) and the objective function of either Eq. (5) or Eq. (6). We presume that this is a smooth, but not necessarily convex, mapping. The design optimization problem is unconstrained other than simple bounds on the design parameters. We chose the search bounds for the PCHIP parameters to be within the range $f(yn)∈[−0.25 cm,0.25 cm]$⁠. These bounds correspond to a distance from the left-hand side of the target and were chosen to be sufficiently small to prevent the target from becoming too thin and rupturing during the RMI spike growth phase of the simulation. The parameter bounds for the number of periods of the profile on the left-hand side of the target is within the range $Nl∈[0.25,2]Nr$⁠, where N_r is the number of periods on the initially perturbed right-hand side of the target. In summary, our optimization problem can be written in the standard optimization formula as

We apply a gradient-free design optimization algorithm known as differential evolution optimization (DEO) to the optimization problem defined in Eq. (7), which is described further in  Appendix B.

We use a hierarchical optimization process that involves two primary stages where a set of parameters are optimized using the DEO algorithm. The first DEO stage involves identifying which of the five parameters [i.e., N_l, $f(y0), f(y1), f(y2), f(y3)$] play a dominant role in reducing the objective function and which parameters can simply be set equal to a boundary value. From the first DEO stage, we found that using $Nl=Nr/2$ is generally the most effective for suppressing RMI spike growth. For most of our design optimization processes, we also found that setting $f(y0)$ and $f(y3)$ equal to the minimum and maximum bounds, respectively [i.e., $f(y0)=−0.25 cm$ and $f(y3)=0.25 cm$], appears to be effective for reducing RMI spike growth. This leaves us with only two parameters for the second optimization stage [i.e., $f(y1)$ and $f(y2)$] that we again optimize using the DEO algorithm. The number of parameters D, approximate number of DEO generations G, and approximate number of objective function $ϕ$ evaluations for each DEO stage are summarized in Table II. Note that the number of objective function evaluations corresponds to the number of times the hydrodynamics simulation must be run.

We compare the optimized designs developed using this optimization hierarchy to the baseline case where the left-hand side of the target is completely flat. For all density plots that we show in this section, the density values map to the color bar shown in Fig. 6. Figure 7 shows density plots for the baseline case at several time instants after the impactor strikes the target. Figure 7(a) displays the initial target configuration at the instant that it is struck by the impactor at $Δt=$ 0 μs (just prior to deformation). After $Δt=$ 1.75 μs, a phase inversion of the interface perturbations has clearly occurred as shown in Fig. 7(b). At later times [see Figs. 7(c)–7(e)], the inverted perturbations continue to grow into long narrow spikes due to the RMI.

Our optimized designs aim to halt and suppress the significant RMI spike growth seen in Fig. 7. Figure 8 shows density plots (at the same time instances after impact as in Fig. 7) of one of our designs using the DEO algorithm to optimize the PCHIP parameters $f(y1)$ and $f(y2)$⁠. Equation (5) is used as the objective function, where this function is minimized at $Δt$ = 7 μs. The cavities created by the PCHIP profile are filled with lucite as previously mentioned. The PCHIP profile on the left-hand side of the target [see Fig. 8(a)] creates an entirely new RMI, which we refer to as the rectifying RMI. This rectifying RMI is equivalent to a positive Atwood number (i.e., A > 0) RMI where no perturbation inversion occurs and the spikes grow left toward the impactor. The vorticity created by this new instability largely counteracts/overtakes the vorticity on the right-hand side of the target and halts or redirects RMI spike growth on this side of the target. We further investigate shock-wave behavior and vorticity deposition for our optimized designs in Sec. III E.

An alternative design is shown in Fig. 9. The cavities between the RMI spikes on the right-hand side of the target [see Fig. 9(c)] are “pinched” together by $Δt=$ 7 μs. This leads to the largely flat interface seen in Fig. 9(e). Similar to the design shown in Fig. 8, the design in Fig. 9 also uses the PCHIP parameterization on the left-hand side of the target and uses a form of the objective function of Eq. (5). However, we make two modifications to the way that the objective function is evaluated. Instead of minimizing the objective function only at $Δt$ = 7 μs, we minimize the maximum value of the objective function over the entire time duration between $Δt$ = 0 and 7 μs. This helps to maintain a relatively flat interface on the right-hand side of the target throughout the simulation. In evaluating the objective function of Eq. (5), we also found it useful to use the average of every other tracer pair $i=1,3,5,…,Nr$ rather than the average of all tracer pairs $i=1,2,3,…,Nr$ as we did for the design in Fig. 8. Using every other tracer pair appears to be more conducive to producing the “pinching” effect that we observe for the optimized design in Fig. 9.

Another alternative optimized design that we have found to work well for suppressing RMI spike growth on the right-hand side of the target is shown in Fig. 10. This design uses piecewise linear segments (instead of PCHIP) to produce the profile on the left-hand side of the target and involves using a layer of lucite between the impactor and the target [see Fig. 10(a)]. We again use the DEO algorithm to optimize two parameters, which are the thickness of the lucite layer and the angle of the grooves on the left side of the target. The effect is similar to the designs in Figs. 8 and 9 where a rectifying RMI is created with a positive Atwood number. As seen in Fig. 10(e), this design is also fairly effective at suppressing RMI growth on the right-hand side of the target. Additionally, this design is relatively easy to manufacture using standard machine tools as previously discussed and is a good initial candidate for an experimental study.

The above results illustrate that RMI, as defined by the objective function in Eq. (5), can be reduced to near zero by optimizing the profile of the left-hand side of the target. Note, however, that the results in Figs. 8–10 show that the left-hand side of the target becomes quite distorted as large bubbles form on this side of the target. This may be inappropriate for some applications. As a second design experiment, we investigate the target compactness objective function of Eq. (6) in an attempt to both reduce RMI and also keep the target relatively compact.

In this new study, the target and impactor materials, dimensions, and initial velocities are identical to those in the previous study. Likewise, the same PCHIP parameterization is used on the left-hand side of the target. Note that in this study the area (and mass) of the target is kept constant, meaning that the PCHIP polynomial coefficients are not independent because of a constraint on the area of the target. The DEO algorithm is again used to find the optimal design. The final design for this new study displays some interesting characteristics and is shown in Fig. 11(a). The design has a low-amplitude high-frequency nature to it, which is quite different than the previous designs in Figs. 8–10. The results at several other times are shown in Figs. 11(b)–11(e). The results in Fig. 11 show that this design has reduced RMI spike formation compared to the baseline case in Fig. 7, although the reduction is not as significant as that in Figs. 8–10. The design in Fig. 11 may also be difficult to manufacture due to its high-frequency nature. However, the left-hand side of the target no longer has large bubbles and the target is clearly more compact, which was the primary design objective.

Figure 12 provides a more quantitative comparison of the spike growth on the right-hand side of the target for the different designs in Figs. 7–11 by plotting the difference between the spike amplitude and the bubble amplitude (i.e., $|hsp−hbu|$⁠) of each design. Tables III and IV show $|hsp−hbu|$ and the percent reduction for the designs in Figs. 7–11 at 7 μs after impact occurs as well as at 4.5 μs after the phase inversion of the perturbations occurs. As is clear from Fig. 12 and Tables III and IV, the optimized designs all significantly reduce the RMI spike growth below that of the baseline case in Fig. 7 (i.e., design 0 in Fig. 12).

The designs that result from our optimization process have several key features. The first key feature is the lucite that fills the cavities created by the optimized profile on the left-hand side of the target. This low-density lucite creates conditions that are similar to a RMI scenario where the Atwood number is positive (i.e., $ρ1<ρ2$ and A > 0) and produces the rectifying RMI. For the rectifying RMI that we generate in our simulations, the Atwood number is equal to A = 0.7657 (where $ρ1=1.186$ g/cm³ for lucite and $ρ2=8.938$ g/cm³ for copper). In cases where A > 0, a perturbation inversion does not occur as it does for the case where A < 0 (see diagrams in Fig. 1).²⁶ Instead, the protrusions created by the optimized profile on the left-hand side of the target begin to grow into the lower-density lucite. However, these protrusions have very little or no room to grow horizontally and almost immediately press up against the impactor. This causes the growing protrusions to immediately start rolling out or “mushrooming” [see Fig. 13(b)]. Nevertheless, these protrusions continue to grow (possibly penetrating into the impactor) and this unstable growth puts stress on the target material. The lucite cavity region also has lower shock impedance Z than the copper target and creates the conditions necessary for compression waves to reverberate through the target and contributes to vorticity deposition that produces the rectifying RMI. (We discuss this further in Sec. III E.)

Another key feature of our optimized profiles is a sharp corner or indentation that protrudes into the target material on the left-hand side of the target. This sharp point acts as a stress concentration point (i.e., the stress in this region is higher than in the surrounding region in the target). In a standard RMI scenario, the stress created by the unstable spike growth leads to the formation of bubbles between the spikes. A stress concentration point between two rectifying RMI spikes (on the left-hand side of the target) appears to exacerbate the bubble effect, which creates a relatively large bubble that penetrates deep into the target as it develops [see Fig. 13(b)]. This exacerbated RMI bubble has the effect of putting vertical stress on the perturbations on the right-hand side of the target. This vertical stress alters the direction of vorticity along the right-hand side of the target and redirects the velocity of the growing perturbations on the right-hand side of the target in a more vertical direction. This counteracts much of the RMI spike growth that would normally arise on the right-hand side of the target.

A third key feature is the difference in the number of periods between the optimized profile on the left-hand side of the target and the perturbations on the right-hand side of the target. As previously mentioned, we found that using half the number of periods on the left-hand side of the target (i.e., $Nl=Nr/2$⁠) to be the most effective for suppressing spike growth on the right-hand side of the target. This is likely due to the pinching and flattening effect that the exacerbated bubbles from the rectifying RMI have on the right-hand side perturbations.

The shock-wave behavior plays an important role in the vorticity deposition throughout the target and is necessary for fully understanding why our optimized profiles are effective at suppressing RMI spike growth. For the baseline case shown in Fig. 7, the shock-wave behavior is relatively simple. When the impactor strikes the target, a normal shock wave is produced that moves to the right through the target. Once this shock wave reaches the air at the perturbed interface on the right-hand side of the target, the shock wave refracts and a rarefaction wave is produced that releases the pressure in the target. This rarefaction wave is produced because the low-density air has a much lower shock impedance Z than the copper target (see shock impedance discussion in Sec. II A). This results in vorticity deposition on only the right-hand side of the target, as seen in Fig. 14, which shows the vorticity at several time instants for the baseline case. Note that red contours imply positive vorticity (counterclockwise) and blue contours imply negative vorticity (clockwise). As is clear from Figs. 14(b) and 14(c), vorticity is high at the base of the growing RMI spikes and acts in opposite directions (clockwise and counterclockwise) on each half of these spikes. The balance of opposite direction vorticity for each spike is correlated with the growth of the spike in the horizontal direction. In Fig. 14, there is essentially zero vorticity deposition on the left-hand side of the target, which is flat for the baseline case.

For the optimized design shown in Fig. 9, the shock-wave behavior is slightly more complicated and the vorticity deposition on the left-hand side of the target is significant. Figure 15 shows a plot of pressure at several early time instants between $Δt=0.75 μs$ and $2.25 μs$ for the design in Fig. 9. When the impactor strikes the target, a normal shock wave is initially produced. However, the shock wave travels at higher velocity and at higher pressure through the copper (which has higher shock impedance Z than the lucite cavity region) [see Fig. 15(a)]. The higher velocity portions of the shock wave move through the copper and lucite and eventually meet to produce a high-pressure region near the interface between the lucite cavity region and the copper target [see Fig. 15(b)]. This high-pressure region reflects off this interface because lucite has a lower shock impedance Z than the copper target and reverberates radially outward [see Fig. 15(c)]. This radial compression wave then interacts with the corresponding wave from the adjacent periods to produce a high-pressure region within the copper protrusions [see Fig. 15(d)]. Eventually, after significant vorticity deposition due to these reverberating shock waves, rarefaction waves release the pressure in the target.

This vorticity deposition produces the rectifying RMI on the left-hand side of the target. Figure 16 shows the vorticity at several time instants for the optimized design in Fig. 9. As the rectifying RMI is generated on the left-hand side of the target, the corresponding bubble penetrates deeply into the target and affects the vorticity on the right-hand side of the target. The result is that the balanced vorticity for each spike (which occurs in the baseline case) becomes unbalanced, causing the spikes on the right-hand side of the target to move more vertically and reducing growth in the horizontal direction. This creates the pinching effect seen in Figs. 9(d) and 9(e) that contributes to the flattening of the right-hand side of the target.

Our results show that it is possible to reduce and suppress RMI growth using gradient-free design optimization methods. Our optimized designs achieve suppression of RMI spike growth by using an optimized profile on the impactor side of the target (the side opposite to the initially perturbed interface) that generates a rectifying RMI. This rectifying RMI nearly completely counteracts RMI spike growth of the initially perturbed interface. In our analysis, we have identified several key features of our optimized designs that are necessary for suppressing RMI spike growth at the perturbed interface. These include several geometric and other features of the optimized profile that are particularly favorable for generating this rectifying RMI. We have also investigated shock-wave behavior and vorticity deposition for our optimized designs and found that the optimized profiles significantly increase vorticity deposition on the impactor side of the target by producing reverberating shock waves in the target. These reverberating shock waves play a significant role in instigating the rectifying RMI.

Overall, the methods that we have developed are very successful at suppressing RMI spike growth at an initially perturbed interface that has been impulsively accelerated using a shock wave. Future work will involve applying the design optimization methodologies that we have developed in this work to other RMI scenarios involving different materials, target geometries, and compression conditions. This may include a multi-objective optimization where the objective functions in Eqs. (5) and (6) are simultaneously minimized to reduce RMI spike growth while maintaining target compactness. Another area for future research may be to use design optimization to exacerbate RMI spike growth (rather than suppressing it as we have done in this work). Additionally, shock-compression experiments will be needed to further validate the results that we have presented in this study and possibly help guide future computational work.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and was supported by the LLNL-LDRD Program under Project No. 21-SI-006. We also thank Y. Choi, J. R. Gambino, A. B. Kostinski, R. N. Rieben, and W. J. Schill for helpful discussions and guidance related to the research presented in this paper.

The authors have no conflicts to disclose.

Dane M. Sterbentz: Conceptualization (equal); Formal analysis (lead); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (lead); Writing – review and editing (equal). Charles F. Jekel: Conceptualization (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Writing – review and editing (equal). Dan White: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (supporting); Writing – review and editing (equal). Sylvie Aubry: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – review and editing (equal). Hector Lorenzana: Conceptualization (equal); Supervision (supporting); Validation (supporting); Writing – review and editing (equal). Jonathan Belof: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review and editing (equal).

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

To demonstrate that our material equations of state and strength models are appropriate for the high-pressure regimes that we consider in this study, we compare simulation results obtained from the hydrodynamics code MARBL to an analytical scaling law developed by Dimonte et al.¹³ The Dimonte et al.¹³ relation provides a correlation between the maximum amplitude of the RMI spike and the peak growth rate of the RMI spike. Under certain compression conditions, the RMI perturbation growth is halted due to the strength of the solid material at the interface (e.g., copper). This implies that there is a maximum spike amplitude $max(hsp)$ that can be obtained for a given set of compression conditions. The plot in Fig. 17(a) shows the spike amplitude $hsp$ vs time for our simulation setup using several different impactor velocities. As seen in this plot, the RMI spike growth is largely arrested after a certain period of time and the spike amplitude levels off to an approximately constant value.^8,13 The Dimonte et al.¹³ analytical relation provides a means of predicting this maximum spike amplitude $max(hsp)$ for a given set of material and compression conditions.

The Dimonte et al.¹³ relation draws from the Piriz et al.⁸ relation, which is derived by integrating the equation of motion for the linear evolution of the perturbed material interface. Piriz et al. assume that the material behaves as a Hookean solid (i.e., as a linearly elastic medium) below the yield criterion. These relations are also developed under the assumption that the upstream material through which the shock wave is initially moving can be approximated as behaving in an elastic–perfectly plastic manner (although we use the Steinberg–Guinan strength model³⁶ in our simulations) that follows the Prandtl–Reuss flow rule with the von Mises yield criterion.⁸ According to Dimonte et al.,¹³ a RMI strength parameter h_Y can be written as a combination of other properties to obtain the following dimensionless form:¹³ 

where k is the wave number of the interface perturbations, ρ₁ is the density of the upstream material (prior to shock loading), Y is the maximum yield stress of this material (we use $Y=6.4×10−3$ Mbar for copper), and $ḣsp0$ is the peak growth rate of the RMI spike, such that $ḣsp0=max(dhsp/dt)$⁠. The wave number k is defined as $k=2π/λ$⁠, where λ is the wavelength of the interface perturbations. The maximum spike amplitude h_sp can be written as a linear function of $khY$ as

where β₀ and β₁ are constants that are dependent on the system under consideration. Dimonte et al. found that $β0≈0.08$ and $β1≈0.24$ work well for scenarios where $A≈−1$⁠, such as the scenario that we consider in this work.

We use Eqs. (A1) and (A2) (setting $β0=0.08$ and $β1=0.24$⁠) and compare this analytical relation with our numerical simulation. Figure 17(b) shows a comparison between the analytical relation of Eq. (A2) and the maximum spike amplitude $khsp$ for several simulation results. As is clear from the results in Fig. 17(b), our simulations appear to be well approximated by the analytical relation of Eq. (A2). Note that one caveat to the analytical relation of Eq. (A2) is that it is only valid for $khY$ below a certain point. For values of $khY$ above this point, the RMI spikes grow quite rapidly and nonlinearly and the assumptions used to develop the relation in Eq. (A2) tend to break down. This relation is also only strictly applicable when the strength model is approximated as being elastic–perfectly plastic. However, our simulations, which use the Steinberg–Guinan strength model for copper, appear to be well approximated by the Dimonte et al. relation.

Gradient-free optimization methods provide a relatively straightforward means of exploring the parameter space without requiring the gradient of the objective function. They are also capable of searching for a global minimum, whereas gradient-based methods are often restricted to finding local minima. There exists a vast variety of gradient-free methods, and most employ a directed random search that attempts to balance exploration vs exploitation (also called diversification vs intensification). It is important to note that “No Free Lunch” theorems³⁸ show that most gradient-free methods perform as well on average as any other method. Therefore, we have not exhaustively investigated gradient-free methods, but have instead simply selected differential evolution optimization (DEO),^39,40 an evolutionary algorithm that shares some similarities with a genetic algorithm. Several studies in the prior literature have found DEO to be useful for solving design optimization problems that also involve hydrodynamics simulations.^41–43 We use a variation of the DEO algorithm in the SciPy.optimize software package in our analysis.⁴⁴ A brief description of this algorithm is provided in the following paragraphs, but additional details on DEO can be found in the work of Storn and Price,³⁹ the developers of this algorithm.

The first step of the DEO algorithm is to generate an initial set or population of trial parameter vectors using Latin hypercube sampling (LHS) to ensure that these values sufficiently span the parameter search space. We designate these trial parameter vectors as $x→i,G$⁠, where the subscript i is the index of a given trial vector within the population and G is the iteration or “generation” number. The elements of these parameter vectors correspond to the individual parameters of interest. The developers of the DEO algorithm recommend using a population size containing 5D to 10D parameter vectors, where D is the number of parameters used.³⁹ These trial parameter vectors $x→i,G$ are then individually substituted into hydrodynamics simulations as described in Sec. II B, which can all be run in parallel for the entire population of trial parameter vectors for a given generation of the DEO algorithm. These hydrodynamics simulations output data for computing the objective function (e.g., tracer data for the amplitude of the spike $hsp$ and bubble $hbu$ for each period of the perturbed interface). These data are then substituted into the objective function of Eq. (5) or (6) to determine the relative fitness of each trial parameter vector for reducing RMI spike growth.

The trial parameter vector that minimizes the objective function the most is selected from among the population of vectors. We refer to this parameter vector as the “best” parameter vector $x→best,G$⁠, and this vector is used in the next generation of the DEO algorithm to produce a set of “mutant” vectors in the following way. Two parameter vectors are individually selected at random from the population of the previous generation. These two parameter vectors (which we designate as $x→r1,G$ and $x→r2,G$⁠) are subtracted to produce a differential variation,³⁹ which is scaled by a weighting factor F. The scaled differential variation is then added to the best parameter vector to produce a mutant vector as shown in the following equation:

This is repeated until a full population of mutant vectors for the next generation $v→i,G+1$ is produced.

After generating the mutant vectors, each trial parameter vector in the previous generation population undergoes a round of “crossover” with the mutant vectors. For the crossover process, we refer to the trial parameter vectors of the previous generation $x→i,G$ as the target vectors. For each element in a target vector, a random number $rand(j)$ between 0 and 1 is chosen. If this random number is less than or equal to the crossover probability CR, then the parameter value from the mutant vector is used in the trial vector for the next generation $x→i,G+1$⁠. Otherwise, the parameter value for the trial vector remains the same as the value from the previous generation. The crossover process for a single vector element $xji,G+1$ can be represented as

where the subscript j is the index of a given element within the parameter vector. An illustration of the crossover process is also provided in Fig. 18. Note that the choice of the weighting factor F and the crossover probability CR are generally application specific. For our application, we set F = 0.5 and CR = 0.75.

This new set of trial parameter vectors $x→i,G+1$ is then used as the trial population for the next generation of the DEO algorithm and are individually substituted into hydrodynamics simulations. The hydrodynamics data are again used to find the best parameter vector $x→best,G$ for the current generation using the objective function of Eq. (5). This process is iteratively repeated until a stopping criterion is reached. This stopping criterion can be specified as a specific number of generations or as a convergence criterion such as

where $σϕi$ is the standard deviation of the objective function across all parameter vectors of the current generation, $⟨ϕi⟩$ is the average value of the objective function across all parameter vectors, and ϵ is a specified convergence tolerance that we set equal to $ϵ=0.01$⁠. Once the stopping criterion is reached, the best parameter vector for the current generation is output and ideally contains the optimal set of parameters for minimizing the objective function. This constitutes the DEO algorithm that we apply to the RMI reduction design problem.