Ejecta are typically formed at the free surface of a target metal when it is loaded by explosives, ballistic plate impact, or a laser source. The dynamics of material transport resulting from particle ejections at the metallic free surface are of great interest in high-energy-density applications.^1–4 The problem of mass ejections from shocked free surfaces is complex and involves multiple physical phenomena, yet it has been shown that significant progress can be made by treating the ejection process as a limiting case of the Richtmyer–Meshkov (RM) instability.^5,6 In recent years, this approach has been successfully employed in developing models for the mass source,^7,8 extension to nonstandard shapes,^7,9 predicting the velocities of the bubble cavities,¹⁰ developing self-similar models of ejecta growth,^11,12 and even using ejecta hydrodynamics to diagnose the yield strength of a material.^13,14 The objective of this paper is to build on this progress to (i) propose an empirical model for ejecta velocities, (ii) suggest an approach to apply the model to doubly shocked metals, and (iii) validate these ideas using high-fidelity, continuum numerical simulations. By focusing solely on the hydrodynamic aspects, we have avoided complicating factors such as material defects, spallation, and uncertainties in the equations of state.

There has been significant recent interest in the ejecta problem, encompassing theory,^7–13,15,16 simulations,^7,13,16–22 and experiments.^{13–15,23–27} For recent detailed reviews, we refer the reader to Refs. 28 and 29. In contrast, there have been few detailed studies of the corresponding double-shock problem. Charakhch’an³⁰ performed numerical simulations of doubly shocked ejecta in Al, which had been rendered in a liquid state following the strong shocks. Continuum simulations of doubly shocked Pb samples were investigated in Ref. 19 for initial perturbation properties matching machined grooves on coupons used in experiments. In Ref. 31, we extended a recently proposed ejecta source model to the multiply shocked situation, and validated the approach using continuum simulations. Researchers from Los Alamos National Laboratory (LANL) reported measurements using a novel two-shockwave tool^24,25 that was used to study ejecta from multiple shocks emanating from an Sn target. Measurements from the experimental campaign included cumulative ejecta mass and time-resolved velocities from both shock events.

The RM instability occurs when a material interface separating different fluids is impulsively accelerated, usually by a shockwave that has just crossed the interface boundary. Thus, perturbations imposed at the interface are amplified as a result of the baroclinic vorticity deposited by the shock–interface interaction. If the fluids straddling the boundary have densities ρ_A and ρ_B, an important parameter governing the growth of the interface is the Atwood numberThen, hydrodynamic aspects of the ejecta phenomenon may be understood by treating it as a limiting case of the RM instability, where A → −1 (corresponding to the RM case where the tenuous material has a vanishing density). The resulting instability growth occurs in stages and depends on the initial interface perturbation amplitude h₀ relative to its wavenumber k ≡ 2π/λ. When kh₀ ≪ 1 [or kh(t) ≪ 1], the growth rates are obtained by linearizing the perturbation equations, and this situation is therefore referred to as linear growth. When either condition is violated, higher-order harmonics become significant, resulting in nonlinear growth of the interface. These harmonics shape the interface into distinct bubble cavities penetrating the heavy fluid and spike channels of heavy fluid advancing into the light fluid. In the ejecta limit, bubbles are rounded wells advancing into the donor material, while spikes are long and narrow jets, which are the ejecta.

In this paper, we address the issue of ejecta velocities in materials that are subjected to two shocks. The problem has gained recent relevance following the successful double-shock experimental campaign conducted at LANL^24,25 through the development of a unique two-shockwave physics package. The remainder of the paper is organized as follows. In Sec. II A, we review previous models for bubbles and spikes, and suggest an empirical model for ejecta velocities that addresses some shortcomings of earlier efforts. In Sec. II B, we discuss shape effects, and then, in Sec. II C, we suggest a protocol for applying the model to the double-shock problem by properly accounting for these effects. In Sec. III we describe the numerical methods and problem setup. The results, which include a detailed comparison of our simulation results with the model over a wide range of conditions, are presented in Sec. IV A, and we present a preliminary comparison with the experimental results of Refs. 24 and 25 in Sec. IV B. We conclude with a summary in Sec. V.

We briefly review previous efforts to describe linear and nonlinear growth rates associated with the RM instability and the ejecta phenomenon. A more detailed summary can be found in recent papers on this subject (see, e.g., Ref. 32). Throughout this paper, we adopt the following nomenclature to clarify the different quantities under discussion: the symbols h and V represent interface amplitudes and velocities, and the ejecta mass is denoted by m. The states of these quantities are indicated by subscripts and superscripts. Specifically, initial states are indicated by the subscript 0, while pre- and post-shock states are indicated by superscripts − and +, respectively. However, for the states immediately preceding and following the second shock, we use the superscripts −− and ++, respectively. Finally, when discussing quantities derived from theoretical models, we use the authors’ initials as superscripts to identify the models.

In the ejecta limit (A → −1), the linear growth rate of interfacial perturbations amplified by an incident shock with a velocity W_i is given by^5,33V0MB=FCkh0−ΔU,(1)where ΔU is the RM interface velocity (or free-surface velocity in the case of ejecta), and F^C = 1 − ΔU/2W_i accounts for shock compression of the interface. When the interface transitions to nonlinear growth followed by differentiation into bubbles and spikes, Eq. (1) is no longer adequate for describing the growth rates. Mikaelian¹⁰ used the potential flow approach (first introduced by Layzer³⁴ for such flows) to obtain an expression for the time-dependent bubble velocity:VbuKM(t)=Vbu,01+32Vbu,0kt,(2)where V_bu,0 represents an initial bubble velocity, which could itself be compromised by nonlinear effects when kh0−>1. To address this effect, Buttler et al.²³ suggested using Vbu,0WB=FbuNLV0MB, where FbuNL is a nonlinear correction factor that was obtained as a fit to higher-order expressions based on Padé approximants³⁵ and is given byFbuNL=11+16kh0−.(3)The corresponding spike model was developed by Zhang,³⁶ starting from potential flow theory, but taking the spike curvature as being opposite in sign to the bubble curvature: ξ_spike = −ξ_bubble = kh₀/2. Asymptotically, this leads to a terminal velocity for the spike tip, which is obtained from the model asVspQZ=V03kh0−+13kh0−+1,(4)where nonlinear effects arising from finite-amplitude initial conditions are encapsulated through the term under the square root. Once again, Zhang³⁶ did not specify a choice for the initial velocity V₀ in Eq. (4), but for kh0−≥1 it is expected the initial growth rate will be decreased from the ideal. Dimonte and Ramaprabhu³⁷ and Buttler et al.²³ suggest a nonlinear correction factor for the initial spike velocities (obtained as a curve fit to the more detailed expression consisting of Padé approximants given in Ref. 35),FspNL=11+(12kh0−)2.(5)Thus, the initial bubble and spike velocities in Eqs. (2) and (4) must first be computed as suggested by Buttler et al.²³ according toVbu0/sp0WB=Fbu/spNLV0MB.(6)Spike growth can also be inhibited at high Mach numbers of the incident shock, M_i, as suggested by Dimonte et al.¹⁶ and observed in numerical simulations.³² Pending a fully compressible potential flow model for spikes, Dimonte et al.¹⁶ proposed an empirical model for the asymptotic spike velocity that depends on M_i, kh₀, and the initial bubble velocity according toVbu0GD=V0MB1+ϕbukh0−,VspGD=|Vbu0GD|ϕsp2kh0−+13Mikh0−+1.(7)The suggested values for the empirical coefficients are ϕ_bu ∼ 0.45, and ϕ_sp = 3.5γ_A/(1 + γ_A), where γ_A is the adiabatic index of the donor metal. While Eq. (7) predicts spike velocities at finite values of M_ikh₀, it does not approach the exact incompressible model (4) in the limit of M_ikh₀ → 0.

To address the shortcomings of existing models, and to explain observed ejecta velocities across a wide range of perturbation groove amplitudes and incident shock Mach numbers, we recently proposed³² an empirical model for the asymptotic spike velocity:Vsp=VMB3kh0−+13Mi(kh0−)2+1.(8)When expressed in the form of Eq. (8), our model has the advantage that it avoids the ambiguities associated with the calculation of the initial velocity V_sp0, and can instead be computed entirely in terms of well-defined quantities such as M_i and kh0−. When recast in the form of Eq. (4), the proposed empirical model for ejecta is expressed asVsp=FspNLV0MB3kh0−+13kh0−+1,(9)where FspNL can be interpreted as a prefactor that accounts for both nonlinearities in the initial conditions and compressibility, and is now given byFspNL=3kh0−+13Mi(kh0−)2+1.(10)As shown in Ref. 32 and in Sec. III of the present paper, Eq. (8) accurately captures ejecta behavior due to nonlinearities and compressibility effects without the need for tunable coefficients, but reverts to the same forms as the exact incompressible, potential flow model result of Ref. 36 in the limit of kh₀ → 0 or M_ikh₀ → 0. For M_i → 1, the models diverge for kh₀ ≫ 1, with significant evidence from experiments and simulations that Eq. (8) more accurately explains the data reported in Refs. 15, 23, 26, and 27. Note that in Eqs. (8) and (10), compressibility effects are coupled to the nonlinear term through the incident shock Mach number, so that even at low values of M_i, these effects may be significant at large enough values of kh₀. Indeed, this is the situation encountered in several experiments that use target coupons with groove amplitudes kh₀ ≫ 1, or in doubly shocked flows where the interface perturbation from the first shock has matured to a nonlinear amplitude prior to arrival of the second shock.

For targets with nonsinusoidal perturbation grooves, Cherne et al.⁷ argued that the longest wavelength associated with the perturbation is responsible for sourcing most of the mass. In practice, this can be accounted for by defining an effective wavelength λ_eff = A_sh/h₀, where A_sh is the missing area under the interface profile (i.e., negative space) and h₀ is the half-amplitude. Thus, λ_eff can be considered the wavelength of an equivalent sinusoid that sources the same amount of mass as the profile under consideration, and when used in source models was found to accurately predict the ejecta areal mass for different shapes. The same approach can be extended to velocities, and was found to be successful in Ref. 32. Thus, the bubble velocity model of Ref. 10 can be written asVbu(t)=(Vbu0WB)eff1+t/τeff,(11)where (Vbu0WB)eff is the initial bubble velocity with the nonlinear correction given in Eq. (3) and evaluated using λ_eff. In Eq. (11), τ_eff represents a nonlinear time scale over which the bubble saturates, and is given by τeff=λeff/3π(Vbu0WB)eff. The corresponding spike velocity from our proposed model [or alternatively Eq. (4)] can be extended to nonstandard shapes usingVsp=(FspNL)eff(V0MB)eff3(kh0−)eff+13(kh0−)eff+1,(12)where terms with the subscript eff are computed using the effective wavelengths and wavenumbers, respectively. Equations (11) and (12) were validated using continuum hydrodynamics and molecular dynamics simulations in Ref. 32, for interfaces with sinusoidal, Gaussian, flycut, and chevron profiles.

When an interface is subjected to multiple shocks, our hypothesis is that it is the bubble profile (amplitude and effective wavelength) at the time of arrival of each shock that explains the ejecta mass and velocity. In Ref. 31, we validated this hypothesis for ejecta mass using data from simulations and experiments, and we address the issue of velocities in the current work. The bubble surface can acquire a nonsinusoidal shape by the time of arrival of the second shock, owing to the influence of higher harmonics in the perturbation equations. Specifically, we find the nonlinear bubble most closely resembles the shape carved out by a flycutter machining tool, which we refer to as a “flycut.” The equation for a flycut graph is parametrized by three independent parameters ^{[Fig. 1(a)]}, namely, the bubble amplitude h_bu [the distance between the free surface—the dashed line in Fig. 1(a)—and the bubble apex], the wavelength of the original sinusoidal perturbation λ, and R, which is the radius of an off-center circle centered at C and osculates with the bubble surface. The flycut profile is then given byh(x)={hbu,0≤x≤12λ−b,R−hbu−R2−(x−12λ)2,12λ−b<x<12λ+b,hbu,12λ+b≤x<λ.(13)where x is the coordinate aligned with the shock direction and b=R2−(R−2hbu)2 is the chord length of the circle at a distance of 2h_bu from the bubble tip ^{[Fig. 1(a)]}. We find that representing the bubble profile with a sinusoidal or parabolic waveform³¹ under-counts the missing area of the perturbation and hence the subsequent ejecta areal mass. We show this in Fig. 1(b), where we compare the nonlinear bubble surface from one of our simulations with various potential waveforms. The solid vertical line represents the location of the free surface also obtained from the simulation. As shown in Fig. 1(b), a sine perturbation undercounts the missing area of the bubble and would therefore under-predict the amount of ejecta that this bubble would produce when shocked.³¹

When working with data from simulations and experiments, it is more convenient to use a lateral bubble length scale rather than R to define the flycut profile. Thus, we choose the bubble half-width L_bu (Fig. 1) as the positive y intercept of the bubble profile, with the free surface as the third independent parameter in defining the flycut profile. If the corresponding spike half-width is L_sp, then 2(L_bu + L_sp) = λ. From Fig. 1(a), it is clear thatR=12(hbu+Lbu2hbu),(14)so the flycut profile may be specified instead in terms of (h_bu, L_bu, λ). The bubble half-width is directly available from simulation data or by computing the streamlines from recently proposed similarity models¹¹ of ejecta (see Sec. IV B). The missing area is the area bound by the bubble tip and the flat surface [at a distance 2h_bu from the bubble tip in Figs. 1(a) and 1(b)] and is given byAsh=θR2−b(R−2hbu),(15)where θ = arcsin(b/R). Thus, our prescription for analyzing doubly shocked ejecta is as follows: From the nonlinear bubble surface just prior to the second shock, obtain the bubble half-width Lbu−−. The radius of the osculating circle R and the corresponding chord length b are computed from Eq. (14) using Lbu−− and the pre-second-shock bubble amplitude hbu−−. The effective perturbation wavelength λ_eff for use in Eqs. (11) and (12) is then computed as λ_eff = A_sh/h₀ using the missing area given in Eq. (15) and hbu−−. Finally, we note that the above analysis is only valid for bubbles that satisfy R>2hbu−−.

In the next section, we compare bubble and spike velocities from our numerical simulations with the models given in Eqs. (2), (4), (7), and (8), but evaluated for the second-shock parameters using the effective wavelength λ_eff as described above. For instance, Eqs. (4) and (8) were evaluated by recasting in the form of Eq. (12) in terms of the corresponding eff parameters.

We validate the above ideas through detailed numerical simulations using the continuum hydrodynamics code FLASH.^38,39 Details of the numerical methods used in FLASH are provided elsewhere,^38,39 and only a brief summary is given here. FLASH solves the governing Euler equations in conservative form, while we use an ideal gas equation of state in our simulations. The hydrodynamics are handled through the piecewise-parabolic method (PPM).⁴⁰ Adaptive mesh refinement (AMR) strategies are used, where the allocation of mesh resources is optimized based on the gradients of pressure and density variables.

The problem setup employed in our simulations in shown in Fig. 2, with terminology adapted from Refs. 37 and 41. The numerical domain represents a computational shocktube, in which a planar incident shock initialized in fluid A traverses along the positive x-coordinate direction. A material interface separates the two fluids A and B, and supports an initial perturbation with an amplitude h₀ and wavenumber k. As the shock traverses the interface from heavy (A) to light (B) fluid, a transmitted shock and a reflected rarefaction are generated. The FLASH simulations were performed with an adaptive mesh resolution equivalent to 256 zones/λ, where λ was chosen to be 1 cm for all cases. Periodic boundary conditions were used on the domain boundaries at y = 0 and y = 2λ, while outflow conditions implemented on the boundaries in the shock direction ensured the egress of acoustic waves and shocks without leading to spurious oscillations within the flow region. In our simulations, we approximated the hydrodynamic response of shocked Sn by an equivalent γ-law fluid with ρ_A = 7.3 g/cm³ and γ_A = 3. It has previously been shown¹⁶ that for these conditions, the modeled γ-law fluid reproduced the U_s–U_p relationship for Sn over particle velocities ranging from 0 to 5 km/s. The light fluid (B) was chosen to have ρ_B = 1.22 × 10⁻³ g/cm³ and γ_A = 5/3, resulting in an Atwood number A → −1.

A summary of all the simulations reported in this work is given in Table I. Cases 1–5 were initialized with a sinusoidal perturbation with scaled initial amplitudes of kh0−=0.12, while the arrival time for the second shock was varied to generate interfaces with different initial amplitudes kh0−− for the second shock. Case 6 was initialized with a chevron interface with a waveform given byand is representative of grooved perturbations on target coupons used in dynamic experiments. Table I also includes the flycut parameters associated with the pre-second-shock bubble surface, namely, the bubble half-width Lbu−−, chord length b, and effective wavelength λ_eff.

The x–t diagram in Fig. 3 shows the sequence of significant events captured by a typical simulation. An incident shock IS₁ with M_i = 1.6 is initialized in fluid A with a piston velocity, and its interaction with the material interface is labeled SI₁ in the figure. The shock–interface interaction leads to a transmitted shock TS₁ and reflected rarefaction RR₁, while the interface acquires a jump velocity ΔU⁺. As the RR₁ wave leaves the domain through the x_l boundary, the guard cells adjacent to that boundary are filled with material corresponding to the conditions of a second shock IS₂. In this manner, a clean second shock can be generated, free of perturbations from the passage of exiting acoustic waves. Furthermore, the timing of the second shock (and thereby the interface amplitude prior to second shock) can be precisely controlled, enabling a systematic study of the effects of nonlinearities on the ejecta velocities. The products of the second shock interaction with the interface are labeled TS₂ and RR₂, and the jump velocity is ΔU⁺⁺.

In Figs. 4–6, we present qualitative results from our simulations by plotting the time evolution of density contours at key stages of the flow development. The contours in Fig. 4 correspond to case 1, for a sinusoidal initial perturbation with an initial amplitude kh0−=0.12 and a second-shock amplitude khbu−−=0.45. The images are plotted in units of scaled time (t − t_SI2)/τ⁺⁺, where t_SI2 is the time instance of second-shock arrival and τ⁺⁺ is the bubble nondimensional time scale as defined earlier, but evaluated here for second-shock conditions. The initial configuration corresponding to a scaled time of −2.2 units is displayed in Fig. 4(a) and also shows the planar shock positioned immediately next to the sinusoidally perturbed interface. Note that the dashed white line in each figure indicates the location of the unperturbed interface obtained from companion simulations under the same shock conditions.

Immediately following the shock interaction, the perturbed interface undergoes a phase inversion due to the negative growth rate resulting from the A < 0 condition.⁴² The growth from the first-shock interaction is shown in Figs. 4(b) and 4(c), while Fig. 4(c) also shows the arrival of the second shock at a scaled time (t − t_SI2)/τ⁺⁺ = 0. Following the second-shock interaction, a second phase inversion is observed in Fig. 4(d) corresponding to a nondimensional time of ∼2.2. The subsequent nonlinear growth and saturation of bubbles are summarized in Figs. 4(e)–4(g). At late times, bubbles converge to a rounded tip, while spikes are narrow and pointed and appear to grow indefinitely with a terminal velocity. In our simulations, the bubble velocities decay asymptotically as ∼1/t, in agreement with hydrodynamic models, while experiments suggest bubble growth is arrested at a finite amplitude, possibly owing to other physics not included here. Relatedly, in the simulations that correspond to idealized conditions, the spike jets do not suffer a break up at late times, since surface tension effects were not included. In spite of these late-time differences with experimental observations, such simulations are valuable since they provide a reliable description of the underlying hydrodynamic aspects that determine gross aspects of bubble and spike growth, including their velocities. Furthermore, as reviewed in Sec. II (and elsewhere in Refs. 7 and 16), models for bubble/spike velocities and source models for ejecta mass are based on purely hydrodynamic potential flow theory.¹⁰ The simulations summarized in this paper and in Refs. 7 and 16 can thus be directly compared with the results of such models.

The density contours in Fig. 5 correspond to case 5, which had a scaled initial amplitude kh0−=0.5 ^{[Fig. 5(a)]}, while the bubbles were allowed to mature to a scaled amplitude khbu−−=1.3 before second-shock impact ^{[Fig. 5(c)]}. Since the bubbles had evolved over ∼6.8τ⁺⁺ before reshock, additional features due to secondary instabilities are visible in Fig. 5(c). Following the second-shock impact, the bubble shape inverts and the secondary features are amplified from the baroclinic vorticity deposited at the interface into the smaller spikes located at y = λ/2 and 3λ/2. In spite of the appearance of secondary instabilities, the dominant features eventually establish themselves, as observed in Figs. 5(a)–5(c). Similar to case 1, these features include a concave bubble tip and a narrow spike jet whose shape is determined by the influence of higher harmonics. Finally, qualitative results from the chevron case are presented in Figs. 6(a)–6(e). The interface properties at first and second shock were kh0−=0.5 and khbu−−=1.3, since the bubble evolved to τ⁺⁺ ∼ 7, at which time they were reshocked. While the interface shapes immediately following the first and second shocks are different owing to the presence of higher harmonics in the Fourier series representation, these effects subside quickly, leaving the primary mode to dominate at late times in a manner similar to cases 1 and 5 [the amplitude of the nth term in the Fourier series decays as ∼(2n + 1)⁻², as shown in Ref. 7]. Note that the chevron perturbation has the same missing area as a sinusoid of the same wavelength, given by A_sh = h₀λ, so that λ_eff = λ.

From images such as those in Figs. 4–6, we infer the amplitudes of bubbles and spikes using the following approach. The bubble tip position x_bu is defined as the x location where the planar-averaged (y) density dropped below 90% of the release density of the liquid metal. Similarly, the spike position x_sp is identified as the x location where the planar-averaged (y) mass fraction of the heavy fluid dropped below 1%. The location of the free surface x_fs is obtained by requiring thatm(t)=∫xbu(t)xfs(t)[ρrel−〈ρ(x)〉]dx=∫xfs(t)∞〈ρ(x)〉dx,(17)Eq. (17) is solved iteratively for x_fs(t), which is then used to calculate the bubble and spike amplitudes according to hbu/sp=|xbu/sp−xfs|. We plot these quantities and corresponding velocities for cases 1, 5, and 6 in Figs. 7–9, respectively.

Figures 7(a) and 7(b) show the evolution of the scaled bubble and spike amplitudes kh_bu and kh_sp as functions of the nondimensional time. We take the initial bubble (spike) amplitudes as positive (negative) in sign, so that the second-shock interaction leads to sign reversals in these quantities. The bubble amplitude khbu++ in Fig. 7(a) shows a saturation-like behavior for large values of (t − t_SI2)/τ⁺⁺. As can be seen in Fig. 7(b), spike amplitudes grow linearly and without bound in the absence of mitigating factors such as material strength. The bubble and spike velocities plotted in Figs. 7(c) and 7(d) are scaled with the ideal linear growth rate from Eq. (1), which has been evaluated here for second-shock conditions and using λ_eff obtained from fitting the second-shock bubble surface to a flycut using the procedure described in Sec. II B. When the potential flow model of Ref. 10 given in Eq. (2) is evaluated using λ_eff , it is in excellent agreement with the scaled bubble velocities from our FLASH simulations.

The corresponding scaled spike velocities are shown in Fig. 7(d), and it can be seen that an asymptotic velocity is achieved following a brief rise time. We compare our simulation results for the spike velocity with the potential flow model of Ref. 36 as well as the empirical models of Ref. 16 and Eq. (8). Once again, all models were evaluated using λ_eff computed from a flycut fit to the bubble surface at the instant of the second shock. The empirical model of Ref. 16 was evaluated using ϕ_sp corresponding to γ_A = 3, to match the ideal gas representation of the shocked liquid metal in our simulations. The empirical models are in excellent agreement with the simulations, while the incompressible potential flow model of Ref. 36 slightly under-predicts the late-time velocities owing to the finite values of Mikh0−− at the second shock.

Figure 8 plots amplitudes and velocities from case 5, where the bubble had evolved to khbu−−=1.3 at second-shock impact. The results from this case follow the same trends as the lower-amplitude data presented in Fig. 7, although the presence of secondary instabilities shown in Fig. 5 introduces some transient aspects in the spike behavior. As the bubble initial amplitude khbu−− was larger in this case, the time to nonlinear saturation, t_nl/τ, increases, as observed in Fig. 8(a). When evaluated with the effective wavelength λ_eff for a fitted flycut surface, the potential flow model of Ref. 10 gives excellent agreement with the simulation data ^{[Fig. 8(c)]}. Spike velocities from this case are shown in Fig. 8(d), and exhibit prolonged transients following the second-shock event, owing to the presence of secondary instabilities discussed earlier. Since the secondary spikes ^{[Figs. 5(e)–5(g)]} occur at a smaller scale (higher k_eff), a period of higher growth rates is observed at first [for 0≤(t−tSI2)/τ++≤10], as expected from Eq. (1). Eventually [at (t−tSI2)/τ++>10], the secondary structures are dominated by the fundamental mode, and the observed nonlinear velocity saturates to a constant value that agrees with the empirical models [Eqs. (7) and (8)] when evaluated using the procedure discussed in Sec. II. Once again, the incompressible potential flow model of Ref. 36 slightly under-predicts the observed spike velocities, owing to the finite values of Mikh0−− for this case. The corresponding results for case 6 with a chevron initial perturbation are shown in Figs. 9(a)–9(d). For both first and second shocks, case 6 had bubble initial conditions similar to case 5. As a result, the developments of bubble and spike amplitudes and velocities show similar behavior to the sinusoidally perturbed case 5.

We summarize bubble and spike velocities from all the simulations listed in Table I in Figs. 10(a) and 10(b). For bubbles, the velocities are scaled with the potential flow model [Eq. (2)] evaluated with k_eff. Figure 10(a) shows that the scaling collapses bubble velocities from cases 1–6, and in agreement with Eq. (2). Similarly, we scale spike velocities with Eq. (8) using k_eff, leading to a collapse of all the data, in agreement with the proposed empirical model. As noted earlier, the short-term transients due to secondary instabilities do not follow this scaling, but these effects subside as the dominant mode asserts itself at late times. The dependence of the asymptotic spike velocities from the second shock on the scaled bubble amplitude k_effh₀ is shown in Fig. 11. The velocities from the FLASH simulations and models are scaled with the impulsive model prediction from Eq. (1). We find that over the range of double-shock conditions investigated here (Mikh0−−≫1), our empirical model is in very good agreement with simulation data. We also plot in Fig. 11 the terminal spike velocity at the end of the first shock (k_effh₀ = 0.12), which is predicted accurately by our empirical model [Eq. (8)]. In contrast, the empirical model of Ref. 16 over-predicts the first-shock data owing to the choice of the parameter ϕ_sp in Eq. (7).

In this section, we provide a possible approach to analyze experimental data from doubly shocked ejecta. A significant challenge in high-strain-rate experiments is that the interface properties immediately preceding the second shock are not available from measurements. We suggest a method to evaluate our empirical model for experimental conditions when the interface shape, particularly the bubble surface length Lbu−−, is not directly available. We then compare the model predictions with the recent double-shocked experimental campaign^24,25 conducted at LANL.

We briefly review the experimental configuration used in Refs. 24 and 25, and refer the reader to those articles for additional details. The experiments described there were based on a two-shockwave tool developed at LANL, in which a machined Sn target is driven by explosive loading from a HE PBX 9501 package coupled to TNT. The detonation constitutes the first shock incident on the target, while a reflected wave bounces off a back anvil to form the second shock. Diagnostics fielded in the experiments included lithium niobate (LN) piezoelectric probes that convert measured pressure signals from ejecta impact on the probe surface to mass measurements using a time-of-flight technique.²⁴ Laser Doppler velocimetry (LDV) probes were also used, and provide time-resolved line-of-sight velocity measurements associated with the ejecta, free surface, and receding bubble. The physics package included a circular Sn coupon, with machined grooves corresponding to a dominant wavelength λ ∼ 80 μm and half-amplitude h₀ ∼ 1.5 μm, with kh0−∼0.12 for the first shock. The experiments reported in Refs. 24 and 25 were at three different initial shock pressures, and in the following we discuss the two shots with the highest pressures, 24.5 GPa and 26.4 GPa, which were labeled shots 1755 and 1756, respectively. The third shot was conducted at 18.5 GPa, below the melt line for Sn, so the ejecta would have been in a solid state.²⁴

As discussed earlier, our approach for evaluating the effective wavelength of the mature bubble surface involves fitting the shape to a flycut profile using the independent parameters (hbu−−, Lbu−−, λ). Owing to the high strain rates and limited optical access in experiments, the pre-second-shock bubble amplitude and lateral length scale are typically not available. This is further complicated by the speculation²⁴ that the recompaction from the second shock causes some of the first-shock ejecta to detach from the bulk, leaving behind a truncated interface of unknown amplitude. Given these uncertainties in the interface profile, we adopt the following two-step approach to determine hbu−− and Lbu−− for evaluating the amplitude and the effective wavelength of the equivalent flycut.

The saturated bubble amplitude may be determined from the cumulative ejecta mass reported by the LN probes at the end of the experiments. The cumulative mass from the second shock is sourced by the bubble before second-shock arrival. Note that, in contrast to our numerical simulations, this assumes the bubbles in the experiment arising from the first shock eventually saturate to a finite amplitude and form the initial conditions for the second shock. This is a reasonable assumption, since the experiments are run to very late times kV₀t, and the bubble motion is likely hindered by viscous and other effects. The ejecta mass (from the second shock) may then be related to the saturated bubble (from the first shock) initial conditions using the model of Ref. 16:σejecta++=ρA++|hbu−−|,(18)where σejecta++ is the measured ejecta areal density from second shock in the experiments, and ρA++ is the post-second-shock density of Sn and can be determined from Table IV in Ref. 16 for a given incident shock pressure. Using Eq. (18) and reported values of σejecta++ measured in the experiments, we obtain the pre-second-shock bubble amplitudes for shots 1755 and 1756 as hbu−−∼ 5.6 and 8.4 μm, respectively.

To determine the bubble half-width Lbu−− , we assume that the bubbles have evolved self-similarly, and compute the limiting streamlines from the self-similar potential flow model of Ref. 11 to determine the intercept with the free surface. The equation for the streamlines in a bubble-attached coordinate system (x′, y′) is¹¹dy′dx′=−arctan(sin⁡y′cos⁡y′+e−x′)ab+12ln(1+2ex′⁡cos⁡y′+e2x′),(19)where a and b are dimensionless coefficients that must satisfy the flow boundary conditions. Equation (19) was derived in Ref. 11 in self-similar coordinates in which the bubble tip is fixed at −xa′, while the spike tip is at +∞. It can be used to compute the limiting streamline that skirts the bubble surface [originating from (−∞, π)], following which the solution is transformed back to the laboratory coordinates (x, y) using the mapping¹¹x→λ2π[x′−a⁡ln(tτ+)],y→λ2πy′.(20)Thus, we determine the bubble width as the intercept of the limiting streamline with the free surface (located at x = 0). The missing area is then computed using Eq. (15) for a flycut, and we obtain A_shape = 328 μm² and 587 μm² for shots 1755 and 1756, respectively. The corresponding effective wavelengths of the saturated bubble are 58 μm (keffhbu−−∼0.6) and 70 μm (keffhbu−−∼0.75) for shots 1755 and 1756, respectively. Finally, we also determine the bubble effective wavelength by directly integrating for the area under the limiting streamlines and obtain A_shape = 329 μm² and 591 μm² for the above cases. This confirms that the flycut shape approximation recommended here accurately captures the properties of the bubble, in agreement with the self-similar, potential flow theory.

In Fig. 11, we compare the predictions from different models (based on the above estimation of keffhbu−−) with ejecta velocity data from LDV pins in the LANL experiments.^24,25 Since there is some variation in the data between the pins for a given experiment, the experimental data points in Fig. 11 are plotted as rectangular symbols with the vertical extent spanning the range of reported values. As expected, the best agreement with the hydrodynamic models is observed for shot 1756, which had the highest shock pressures, ensuring the ejecta were in a molten state. For this case, all three models provide good agreement with data. The growth rates from shot 1755 were likely compromised by material strength effects, since the ejecta were likely in a mixed state. For these conditions, all the hydrodynamic models over predict the experimental data point. We also show in Fig. 11 the asymptotic ejecta velocities from the first shock from the LDV readout (kh0−=0.12). This data point straddles the prediction curves from our empirical model and the potential flow model of Ref. 36. Owing to the lower values of kh0− for the first shock, the empirical model of Ref. 16 over-predicts the experimental data.

We have proposed an approach for predicting ejecta velocities from metallic surfaces when they are driven by multiple shocks in sequence. The procedure is validated for the special case of two shocks, assuming each shock is strong enough to melt the material (on shock or on release²³). Our central hypothesis is the bubble cavity immediately prior to the second shock forms the initial conditions for the subsequent ejecta formation and is responsible for sourcing the ejecta mass. If the bubble has already reached a nonlinear state, then the initial conditions correspond to a nonsinusoidal surface. We find that a flycut approximation best describes the bubble shape at such late times. Our prescription for describing the subsequent dynamics is to use an effective wavelength λ_eff of such a flycut profile to describe the bubble shape in the models for bubble and spike velocities as well as ejecta mass. The approach has been validated for ejecta mass previously,³¹ and we have shown the applicability of this method to velocities in this paper.

Thus, our suggestion is to treat the pre-second-shock bubble initial condition as a flycut surface and compute the corresponding λ_eff for use in the velocity models. For bubbles, the potential flow theory of Ref. 10 has been evaluated using this strategy. For spikes, we have tested three contending models and have found that our empirical model provides good agreement with data, since it accounts for both nonlinearities as well as compressibility through the shock Mach number. A key aspect of our model is the recognition that these two effects act in concert through the M_ikh₀ term, so that, even at moderate shock Mach numbers, the growth rates can be affected if the initial amplitudes are nonlinear. This is certainly the case in experiments, which are often initialized with finite amplitudes, and especially true in double-shock experiments where the bubble has been allowed to develop to a fairly mature state before being reshocked. The incompressible, potential flow model of Ref. 36 is exact and has been included in our comparison to highlight the significance of the combined M_ikh₀ term even at low values of M_i. The proposed model contains no tunable coefficients and computes ejecta velocities without requiring information on the initial growth rate Vsp0, which can be ambiguous. We have tested these ideas by comparison with continuum numerical simulations using the FLASH code, for different initial amplitudes and shapes at the second shock. For all the cases studied here, our empirical model given by Eq. (8) provides the best agreement with data when combined with the approach described above.

We have also provided an approach for analyzing experimental data where there is often significant uncertainty in determining the interface properties before the second shock. To compute the effective wavelength, we suggest using a self-similar model such as that of Ref. 11 and determining the bubble width Lbu−− as the intercept of the calculated streamlines with the free surface. To constrain the bubble tip location, the bubble amplitude is estimated from the measured cumulative ejecta mass using a source model¹⁶ that assumes a saturated bubble at late times. The approach appears promising, and, for the experiments where the shock pressure exceeded the melt condition, the measured velocities are in agreement with models evaluated using this approach. However, the comparison with experiments provided here is preliminary, and a broader examination of experimental data is required. We close by noting that the models and approaches presented here will be of value in interpreting such experimental campaigns involving shocked targets.