Computational Analysis of the Metal Free-Surface Instability , Fragmentation , and Ejecta under Shock

We conducted numerical simulations of the related processes of interface instability, tensile fragmentation, and jetting resulting from four kinds of typical macro defect perturbations (chevron, sine wave, rectangle, and square) on a Cu free surface under a reflected shock wave when Cu impacts a solid wall at a speed of 2.5 km/s and found that, for the chevron and sine wave cases, the ejecta velocities of the head are 6.28 and 5.88 km/s, respectively. Some parts of the inner material are in a tensile state without any fragmentation, which is observed only in the main body of the material owing to the tension effect. Furthermore, for the other two initial perturbations (rectangle and square), the highest ejecta velocities may even reach 9.14 and 9.59 km/s, respectively. Fragmentation caused by multilayer spallation can be observed on a large scale in the Cu main body, and there are granules in the front area of the ejecta but the degree to which fragmentation occurs is much less in the Cu main body and there is a notable high-speed, low-density granule area in the ejecta head. Finally, we present a detailed analysis of the spatial distribution of the granules, ejecta mass, pressure, temperature, and grid convergence.


Introduction
Metal interface instability occurs as a reaction to the effect of a shock wave or acceleration or to a shear load on the perturbed interface, and may afterward lead to fragmentation and mixing.This phenomenon, common in explosions and shock processes under extreme conditions (high temperature and high pressure), is jointly controlled by three scale factors, i.e., 1) the thermodynamics and the geometric boundary (macroscopic scale), 2) the initial perturbation of the metal interface (mesoscopic scale), and 3) material properties (microscopicmesoscopic scale).Metal interface instability is a typical multiscale, multiphysical, strongly coupled, nonequilibrium, and complex flow phenomena.
Related studies of the Rayleigh-Taylor (RT) instability of metal began in the 1970s.In particular, in the USA and Russia, theoretical studies, numerical simulations, and experimental investigations have been conducted.The prior research by Barnes, Blewett, McQueen, Meyer, and Venable [1] mainly focused on the metal interface instability growth of a flat plate driven by expansion of detonation products, as observed in the Los Alamos high-energy X-ray facility.By modeling with the two-dimensional elastic-plastic numerical hydrodynamics code MAGEE, they confirmed that amplitude growth was apparently governed mainly by the dynamic yield strength of the material.This pioneering experiment has been regarded as the standard experimental model in solid RT instability research.To confirm the theoretical analysis of Drucker [2], the additional experiments conducted by Barnes, Janney, and London et al. [3] showed that no more than semiquantitative agreement can be expected between the experiments and the results.However, the main conclusion reached by Drucker [2], that amplitude and wavelength determine the onset of the RT instability in a solid, has been proved in the experiments.Nevertheless, since Barnes, Blewett, McQueen, Meyer, and Venable's research [1], which was limited by the testing and diagnostic techniques, related experiments on the RT instability in a solid have not been improved by much.In 2006, Lindquist, Cavallo, and Lorenz et al. [4] studied the metal RT instability under high pressure and high strain rate, utilizing the same experimental configuration as did Barnes, and combining with numerical simulation, they eventually confirmed the applicability of the standard constitutive models: Steinberg Guinan (SG) and Preston-Tonks-Wallace (PTW).In recent years, owing to the rapid development of laser equipment and the construction of large laser-loaded devices, related laser-loaded experiments have been conducted.To study the effect of material strength on the metal RT instability under Mbar pressure, Weber, Kalantar, and Colvin et al. [5] have conducted a hohlraum experiment of an X-ray-driven metal layer using NOVA.Because of the X-ray preheating and metal-melting problem, these experiments failed to provide any definitive results.To eliminate the effect of preheating and metal melting, Edwards, Lorenz, and Remington et al. [6] developed a new laserdriven method for dynamic shockless loading of solid materials.Based on this platform, Lorenz, Edwards, and Jankowski et al. [7] [8] used the Omega laser to research metal RT instability under high pressure and high strain rate.In 2009, Park, Remington, and Beckeret et al. [9] also utilized the Omega laser to study RT instability under Mbar pressure.In their experiment, the target samples were driven quasi-isentropically.To solve the problems resulting from preheating and metal melting, plasma from a laser-driven polymer layers acts on the CH-based epoxy samples across the vacuum gap.In 2014, Smalyuk, Casey, and Clark et al. [10] firstly observed the hydrodynamic instability growth in indirectly driven implosions at the National Ignition Facility; their results are in good agreement with theoretical and simulation studies.Their work is the first step in experimentally studying hydrodynamic instability in a converging configuration.American and Russian researchers [11] [12] [13] conducted explosive and magnetic driving experiments to research material constitutive behavior by using perturbation growth and then developed a new material constitutive model to predict the results of experiments under high pressure and high strain rate.In the 1990s, using a high-speed camera, Chen, Zeng, and Wang et al. [14] performed experimental research on the interface instability of a alloy/water system in a cylindrical implosion configuration that was numerically simulated based on the hydrodynamics by Lin, Xia, and Zhang [15].Recently, Wang, Bai, and Cao et al. [16] have conducted experiments and numerical simulations of an aluminum flyer driven by an explosion and observed the sine perturbation growth at early time.
The theoretical study of the metal RT instability is limited by the material models used and also the lack of comparison with related experiments.Particularly, at the turbulent mixing stage, it is difficult to predict the development of the interface instability.Despite some progress made in experimental studies, understanding of metal instability in detail is still hindered by diagnostic techniques.Therefore, because the problems in algorithm precision, the constitutive model, and the equation of state (EOS) are still too tough to solve, simulation work stands out as especially worthwhile.Some effective software and codes, such as ABQUS, Lagrange codes FCI2, MAGEE, TOODY II [17], and Euler codes ALE3D and CHAP [18] and the two-dimensional (2D) Euler MHD program, have been designed and improved just for this purpose.By combining simulation and experiment, the evolution pattern of the metal interface instability is presented with time.The instability growth obtained indicates that the material strength and load history and the perturbation wavelength and amplitude have a vital influence on the RT instability.Specifically, the material strength restrains the growth and development of the interface instability, and only when the initial amplitude reaches or exceeds a specific value does the instability develop.
Under shock loading, macro defects perturbations on metal free surface may lead to ejecta and fragmentation.We numerically studied the growth of the instability on the metal interface with four kinds of typical macro defect (chevron, sine wave, rectangle, and square) perturbations and the dynamic characteristics of ejecta.

Numerical Methods
Given a multimaterial, large-deformation, strong-shock physics problem, we improved our previous large eddy simulation code MVFT (multiviscous flow and turbulence) [19] [20] [21] by taking into consideration explosive detonation and the elastic-plastic behavior of the material.An Euler finite volume algorithm and code describing detonations and shocks dynamics with three-order accuracy has been developed using the two-step Euler algorithm with three-order accuracy.More specifically, the Lagrange and remap algorithms are applied to solve the mass, momentum, and energy conservation equations.In a single Lagrange step, the integration of the Lagrange control equation at one time step is obtained, the deformation of the material leads to deformation of the initial grid, and there is no material flow between the grids, whereas, in a single remapping step, the velocity and energy, deviation stress, and other state parameters in the Lagrange deforming grid are remapped back to the initial grid.The governing equations used are given as follows: where V is the control volume, u j is the velocity, S is the surface area of control volume, n  is the external normal line, s ij is the deviatoric stress tensor, P is the static pressure, and E is the total energy per unit mass.First, in our numerical simulation, a dimension splitting method is used to split Equation (1) into three one-dimensional problems, which are then solved using the piecewise parabolic method (PPM) method to perform the interpolation and reconstruction of the physical quantity in each grid.Owing to the lack of automatic monotonicity of PPM, there would be numerical oscillation at the discontinuity point, which leads to a decrease in the accuracy of the discontinuous solution.To restrict the numerical oscillation, we introduce a flow restrictor.In adopting this method, a monotonic limiter is utilized so that the calculation of δQ j is limited monotonically, as shown in the following: )( ) In a Lagrange step, the material strength model, artificial viscosity, and explosive detonation model are needed for calculation and, therefore, the Jones-Wilkins -Lee (JWL) equation of state is used for the explosive detonation and the Mie-Gruneisen equation of state is applied for the strength material.As for the strength model, we utilize the Steinberg-Guinan constitutive model.
In the Steinberg-Guinan constitutive model [22], terms for pressure, temperature, and strain rate are added to the elastic-plastic constitutive equation while the coupling effect of pressure and strain rate on flow stress has a characteristics of separating variables.Additionally, because the flow stress in the Steinberg-Guinan constitutive model relies on pressure, the material constitutive equation and the state equation are coupled, indicating stress hardening of the metal under a high pressure.The shear modulus and flow stress in the Steinberg-Guinan constitutive model can be expressed as ( ) ( ) ( ) ( ) where Y 0 signifies the yield strength in the initial state, G 0 is the shear modulus in the initial state, P and T are the pressure and temperature, ( ) respectively, the partial derivatives of shear modulus to the pressure and temperature in the initial state, A and α correspond to and n are material strain hardening parameters, ε is the strain, and is the material compression ratio.The form taken by the Steinberg-Guinan constitutive model is irrelevant to the strain rate, yet the requirement upon the strain rate is that it should be >10 5 s −1 .The reason for this restriction is that the effect of metal softening counteracts that of metal hardening, so this model can be used to describe multimaterial flow stress under high pressure, which is also the mostly used constitutive model under high pressure.

Validation of the Simulation Code
Through simulation of Barnes's three experiments [1] [3] on an 1100-0 Al flat plate accelerated by a detonation product, we confirmed the simulating capability of the high-fidelity detonation and shock dynamics calculation program.In these experiments, the initial perturbations of the Al flat plate were 1) initial wavelength 5.08 cm and initial amplitude 0.02 cm, 2) initial wavelength 2.54 cm and initial amplitude 0.01 cm, and 3) initial wavelength 4.8 cm and initial amplitude 0.01 cm.They utilized a P80 detonation plane-wave lens to trigger a shock initiation on a 3.81-cm-thick PBX-9404, after which the detonation products form an isentropic load for the Al plate through a vacuum gap 2.54 cm long.The parameters of the PBX-9404 JWL EOS used in the simulation are given in Table 1.The 1100-0 Al Mie-Gruneisen EOS and Steinberg-Guinan constitutive model parameters are listed in Table 2 and Table 3.
In the first experiment (wavelength 5.08 cm and amplitude 0.02 cm), it can be seen that, by a seminumerical linear analysis, when Y 0 = 0.055 GPa, the numerical simulation agreed well with the experiment data.Although the 2D calculations of our program agree well with the experiment (Y 0 = 0.075 GPa), there are differences between the above results and those from the 1100-0 Al Steinberg-Guinan constitutive model under the condition of Y 0 = 0.04 GPa.This is true especially when Y 0 = (0.325GPa) and (0.2 GPa), in which case the results significantly differ from those in ref. [1] [23], shown in Figure 1.World Journal of Mechanics  It is easily observed from Figure 1 that the amplitude is almost a constant owing to the large strength.When Y 0 declines to 0.04 GPa, the growth of the amplitude calculated by our program is greater than that in the experiment, which indicates that the influence of the strength on the amplitude growth is weak.For the case when Y 0 = 0.075 GPa, the simulation result agrees well with that from the experiment.However, the result for Y 0 = 0.055 GPa coincides with Barnes's further analytical results, which took account of the deviation of X-ray photographic magnification [3].Given the uncertainty deviation of the experiments, it is impossible to obtain a precise result for the initial yield strength by the comparison of numerical simulation and experiment.Nevertheless, the proper initial yield strength obtained from the simulation is supposed to illu-strate the applicability of the constitutive model.For the other two experiments, we conducted the simulation based on the initial yield strength Y 0 = 0.075 GPa of 1100-0 Al calibrated by using our program.
In the simulation of the experiment (with wavelength 2.54 cm and amplitude 0.01 cm), three kinds of grid resolution were adopted: ∆x = ∆y = 4, 8, and 15 μm.In Figure 2, simulation results are seen to agree well with experimental results and show grid convergence.The simulation results in Figure 3 (with wavelength 4.8 cm and amplitude 0.01 cm) indicate the obvious convergence trend., our calculation is in accordance with the experiment, and the grid convergence is also in the same trend.

Results and Discussion
Figure 4 shows the four kinds of typical macro defect perturbations (chevron, sine wave, rectangle, and square) on a Cu metal free surface under the reflected shock wave produced from copper impacting the left solid wall at a high speed of 2.5 km/s.When the reflected shock waves reach the free surface, the reflected rarefaction wave in the material leads to the acceleration of the interface particles, as well as to tensile and unloading effects.At the same time, various tensile failures in the material are also under complex and interacting influences of the solid wall's reflected shock wave and rarefaction wave.Owing to the different initial perturbations on the interface, 0 p ρ ∇ ×∇ ≠ results in the interface insta- bility, fragmentation, and jetting.Numerical simulations are used to analyze quantitatively the difference in interface instability for four typical macro defects and to study the various distributing features of fragmentation and ejecta inside the material and around the interface.In our simulations, the Cu EOS and constitutive model parameters are given in Table 4 and Table 5.
Configuration of the calculation modes are shown in Figure 4; the computational domain is [0 cm, 8 cm] × [−1 cm, 1 cm], the length of Cu is 5 cm, and the boundary conditions on the left, top, and bottom edges are all the solid walls.The positions of the four defects (chevron, sine wave, rectangle, and square) are located at the center of the free surface, with the same area of 0.25 cm 2 , and the Cu mass of defect unit thickness is 2.1 g.For the chevron, the defect shape is an isosceles right-angled triangle with a 1-cm-long bottom edge.The sine wave amplitude and wavelength are 0.3927 and 2 cm, respectively.Additionally, the length and width of the rectangle are 0.7071 and 0.35355 cm, and the length of the square is 0.5 cm.
To verify grid convergence, we take the chevron model as an example and select two different sizes (25 and 12.5 µm). Figure 5(a) and Figure 5(b) demonstrate the grid convergence of the velocity history of the interface vertex and the velocity distribution on the symmetrical axis at t = 13.5 µs.For the coarse and fine grids, the peak velocities are 6.32 and 6.28 km/s, the starting times of the vertex are 5.845 and 5.854 µs, and the relative errors satisfying the convergence  condition are 0.64% and 0.15%, respectively.In Figure 5(b), it is obvious that the velocity distributions in the calculation area mostly satisfy the convergence condition, whereas the calculations of the quantity in the fragmentation area on the 4 to 5 cm horizontal axis display significant differences.The left panel of Figure 6 shows a partially enlarged view of Figure 5(b), from which the size of the fragmentation area can be obtained.The width of the fragmentation area for coarse grid is 2.925 mm, and that in the fine grid is clearly larger (5.812 mm), which means that the size of the fragmentation area is related to the computational grid.Moreover, the size of the fragment is also easy to determine from that of the velocity discontinuity's discrete space.The size of the World Journal of Mechanics large fragments is almost the same in the two different grids, and that of the smallest fragment even reaches one grid size, so the size of the smallest fragment is linked to the computational grid.The right panel of Figure 6, in contrast, shows the allelic distributions of the density in the coarse and fine grids at t = 13.5 µs, so the ejecta mass can be calculated (with the length of the horizontal axis exceeding 5.4 cm).Specifically, the ejecta masses are 0.799 and 0.737 g, respectively, and the difference is about 8.41%, suggesting that the ejecta mass is also closely linked to the computational grid.
The above findings reveal that the numerical calculations of the fragmentation are related to the size of the calculation grid.To eliminate the influence of the grid's mesh, we chose the same grid size, 12.5 µm in length, in the four models.Thus, the Euler mesh encompasses 12,800,000 points in the whole computational domain [0 cm, 10 cm] × [−1 cm, 1 cm].
To compare the various characteristics of the interface instability, fragmentation, and ejecta under the influence of the four typical defects, Figure 7 shows the material distribution at two moments (t = 10.5 and 13.5 µs), where the chevron ejecta shape is a nonuniform pole with a sharp-pointed nose.For the sine wave, the ejecta is also a nonuniform pole but with a mushroom head.The common characteristic in the former conditions is that hardly any distinct particle jet can be found.For the rectangle and square defects, the shared characteristic is that high-speed particles concentrate on the symmetrical axis of the ejecta head while they differ in that the rectangle has a conical hole shape in the front and the square is shaped like a large-size arc.
Table 6 lists the ejecta mass of the four defects when t = 10.5 and 13.5 µs.At 10.5 and 13.5 µs, the horizontal axis of the calculation domain is larger than 4.5 and 5.4 cm, respectively.From Table 6 we find that the ejecta mass of the chevron is 0.533 g, and those of the sine wave, rectangle, and square are, respectively, 0.581, 0.662, and 0.586 g.At t = 13.5 µs, the ejecta mass of the chevron is 0.737 g, and for cases the sine wave, rectangle, and square the ejecta masses are World Journal of Mechanics  0.803, 0.848, and 0.712 g, respectively.In summary, although all the ejecta masses increase with time, it is obvious that the law of ejecta mass and growth varies according to the kind of typical macro defect perturbation.The left panel of Figure 7 illustrates that the distributions of the fragmentation area near the interface located on both sides of the ejecta are basically the same, both being caused by the unloading of the free surface.Nevertheless, the interactions of the stretch formed from the meeting of the unloading wave in the material and the reflected wave from the solid wall leads to a totally different fragmentation area at t = 13.5 µs.For the chevron and sine wave, the fractured zones behind the metal main body are much larger than those from the other defects, which are less obvious.
For further analyze of the velocity and spatial distributions of the jet particles, Figure 8 and Figure 9 show the velocity distribution and density distribution around the ejecta head on the symmetrical axis at t = 13.5 µs. Figure 8 show that the ejecta head velocities of the chevron and sine wave cases are 6.28 and 5.88 km/s, respectively, and hardly any discrete velocity distribution can be detected in the front of the ejecta in either defect condition.Figure 9 shows that the World Journal of Mechanics low-density area can be found only in part of the ejecta head, indicating that there is tensile stress in the material, without any fragmentation.This explains why no discrete broken particles are found in the front of the ejecta and why only fragmentation can be found in the main body owing to the stretch.However, the ejecta velocities for the rectangle and square are even higher than in the former two conditions, and the highest velocities reach 9.14 and 9.59 km/s, respectively.Additionally, the discrete density distributions are distinctly presented in Figure 9, as are the density distributions.The density near the ejecta head is zero, indicating that most of the material fragments are the result of multilayer spallation.Therefore, in the rectangle and square conditions, discrete particles are mainly located in the front area of the ejecta, and the degree of fragmentation is minor in the main body.and 8.16 -8.27 cm.For the square, the particle distribution domains are 7.61 -8.30, 8.92 -9.19, 9.33 -9.44, and 9.59 -9.82 cm, so the conclusion may be drawn that the square undergoes more serious fragmentation than the rectangle.On a macro level, both ejecta heads of the square and rectangle are in a high-speed, low-density area in the particle state.At the same time, the fragmentations and the distribution are closely connected with the initial defect shape.To further understand the particles' thermodynamic state, Figure 10 shows the pressure and temperature distributions of the particles around the ejecta head.From Figure 10 we see that the pressure of the particles under the rectangle and square conditions is almost zero.The particle temperature has reached 1350 -1700 K, and at some points it even climbs to Cu's melting temperature (~1400 K).

Conclusions
By considering four kinds of typical macro defect (chevron, sine wave, rectangle, and square) perturbations on a Cu free surface, our research in the present work World Journal of Mechanics and ejecta all originate from the initial interface defect and are also associated closely with the shape of the defect.According to the above findings, in the chevron and sine wave cases, the ejecta mass velocities of the head are 6.28 and 5.88 km/s, respectively.Some parts of the inner material are found in a tensile state without any fragmentation, which appears only in the main body of the metal owing to the tension effect.Additionally, for the other two initial perturbations (rectangle and square), the highest ejecta mass velocities are 9.14 and 9.59 km/s, respectively.Fragmentation appears in the large range of shaped pole owing to the multilayer spallation.There is a granule area with high-speed and low-density in the ejecta head.However, the degree of fragmentation is much lower in the main body of Cu.Overall, the jet masses under the four defect conditions vary greatly.Moreover, their time-dependent eject mass follow different law.
Although the Euler method we utilized to calculate the material failure under the stretch strength has met the convergence requirements for conserved quantities, the sizes of the fragmentation area and of the smallest particle are still related to the computational grid size.Our future work will be directed to largescale computing under a micrometer and to submicron grid size.With the combination of the dynamic experimental data, the best computational size is likely to be obtained.Therefore, our simulation ability will be developed from macro to micro scale.Moreover, because of the complex thermodynamic processes during loading and unloading, establishing a correct constitutive model and EOS will be crucial to improving the reliability of out numerical simulations, and considerably more work remains to be done to study the physical factors that influence the interface instability, fragmentation, and ejecta.

Figure 1 .
Figure 1.Experimental and numerical simulation amplitudes when the initial amplitude is 0.02 cm and the wavelength is 5.08 cm.

Figure 2 .
Figure 2. Experimental and numerical simulation amplitudes when the initial amplitude is 0.01 cm and the wavelength is 2.54 cm.

Figure 3 .
Figure 3. Experimental and numerical simulation amplitudes when the initial amplitude is 0.01 cm and the wavelength is 4.8 cm.

Figure 4 .Figure 5 .
Figure 4. Four typical initial macro defects (chevron, sine wave, rectangle, and square) on the Cu free surface.
Figure 9(c) and Figure 9(d) demonstrate the distribution of the ejecta particles in the rectangle and square conditions.The particles in the rectangle are distributed mainly in the domains of 7.47 -8.06 World Journal of Mechanics

Figure 9 .
Figure 9. Density distributions around the ejecta head on the asymmetrical axis in the four cases at t = 13.5 µs: (a) chevron; (b) sine wave; (c) rectangle; and (d) square.

Figure 10 .
Figure 10.Pressure (a) and temperature (b) distributions around the ejecta head on the asymmetry axis at t = 13.5 µs.

Table 6 .
Jet mass under the four defects at 10.5 and 13.5 µs.