Numerical Simulations of Laboratory‐Scale, Hypervelocity‐Impact Experiments for Asteroid‐Deflection Code Validation

Asteroids and comets have the potential to impact Earth and cause damage at the local to global scale. Deflection or disruption of a potentially hazardous object could prevent future Earth impacts, but due to our limited ability to perform experiments directly on asteroids, our understanding of the process relies upon large‐scale hydrodynamic simulations. Related simulations must be vetted through code validation by benchmarking against relevant laboratory‐scale, hypervelocity‐impact experiments. To this end, we compare simulation results from Spheral, an adaptive smoothed particle hydrodynamics code, to the fragment‐mass and velocity data from the 1991 two‐stage light gas‐gun impact experiment on a basalt sphere target, conducted at Kyoto University by Nakamura and Fujiwara. We find that the simulations are sensitive to the selected strain models, strength models, and material parameters. We find that, by using appropriate choices for these models in conjunction with well‐constrained material parameters, the simulations closely resemble with the experimental results. Numerical codes implementing these model and parameter selections may provide new insight into the formation of asteroid families (Michel et al., 2015, https://doi.org/10.2458/azu_uapress_9780816532131‐ch018) and predictions of deflection for the Double Asteroid Redirection mission (Stickle et al., 2017, https://doi.org/10.1016/j.proeng.2017.09.763).


Introduction
We rely on computer modeling to assess how to deflect or disrupt a hazardous asteroid in order to prevent the possibility of it impacting Earth. Code validation is key to ensuring confidence in such simulation results and is important in advancing our understanding of impact processes on asteroids in general. Benchmarking hypervelocity laboratory experiments using well-characterized materials allows us to perform useful validation tests of such models.
In a collaboration between Lawrence Livermore National Laboratory and Kobe University, we compare simulations quantitatively with laboratory experiments, examining the fragmentation of rocky bodies due to hypervelocity impacts. The results are characterized by the fragment distribution in mass and velocity. We focus on the Nakamura & Fujiwara, 1991 two-stage light gas-gun experiments, wherein a spherical basalt target is impacted at an angle with a spherical nylon impactor (Nakamura & Fujiwara, 1991). We model these experiments using Spheral (Owen, 2010(Owen, , 2014Owen et al., 1998), an open-source adaptive smoothed particle hydrodynamics (ASPH) code available on GitHub, well suited to track stresses and strains during deformation of solids. Spheral works by solving three components simultaneously: the conservation ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2018EA000474
Key Points: • We investigated the accuracy of our code by comparing our simulation results to data from a 1991 hypervelocity experiment • The simulation results indicate that our code can produce results that closely resemble the experimental findings • This work provides insight into model and material parameter selections in our code, potentially applicable in other numerical models Several years after the 1991 gas-gun experiment, Benz and Asphaug (1994) used the experimental results from Nakamura and Fujiwara (1991) to validate their own smoothed particle hydrodynamics (SPH) code. At that time, the material parameters of the target had not yet been measured by Nakamura, leaving them to rely on simulation scans to find best fits to the experimental results. A decade after the work by Benz and Asphaug, the material parameters of the Yakuno basalt used in the original 1991 gas-gun experiments were measured. With the experimentally measured material parameters now available for code validation, we utilize these values in Spheral and investigate the sensitivity of the simulation results to material parameter, strength, and strain-model selection in the code.

Benchmarking the Experiment
The Nakamura and Fujiwara experiment (Nakamura & Fujiwara, 1991) shot a spherical nylon impactor into a spherical basalt target, with the impactor 0.7 ± 0.01 cm in diameter striking the target 6.0 ± 0.05 cm in diameter at a velocity of 3.2 km·s −1 . The impactor hit the target at an angle incident to the surface of 30° ( Nakamura & Fujiwara, 1991) (Figure 1). This single experiment is the subject for all comparative numerical work in this paper. All simulations are carried out in a full three-dimensional geometry.
Our code-validation investigation focuses on our constitutive model, which includes separate strength and strain models. We keep the equation of state and damage model fixed for all simulations, implementing the Tillotson equation of state for both the basalt target (Benz & Asphaug, 1999) and nylon impactor (Benz & Asphaug, 1994), as well as a tensor generalization of the Benz-Asphaug damage model (Benz & Asphaug, 1994). Additionally, we employ a friends-of-friends algorithm to identify the ASPH points that make up the individual fragments still connected by competent (undamaged) material in order to characterize fragments in the simulations.
We begin by comparing the simulation results of the damage morphology for the two strain-model selections available in Spheral to the high-speed photographs taken during the experiment, as well as the recovered fragment data from the experiments. We then investigate two different strength models in Spheral to inform users of the merits and limitations of each model. Finally, we present the simulation results for a combination of different material-parameter and strength-model selections and compare to the experimental data.

Sensitivities in Spheral
Sufficient computational resolution for the basalt spheres is first determined by assessing at what resolutions total damage within the target begins to plateau. We ran seven different scans of resolution, increasing the particles per centimeter for each successive run by 30-100%. The results are presented in Figure 2. At low resolution, damage decreases with increasing resolution, until converging at resolutions above approximately 50 particles across the diameter of the target. This trend is consistent with previous work using Spheral, which demonstrated damage can be overestimated with insufficient resolution (Benz & Asphaug, 1995;Bruck Syal, Rovny, et al., 2016). For this work, we subsequently adopt a resolution of 150 particles across the diameter of the basalt sphere, corresponding to a total number of~1.75 million particles within the target volume, requiring 512 processors per simulation.

Strain Models
As a part of this study, we compare two models of strain available in Spheral. A strain model is used by a brittle-damage model to determine when the stress in a piece of material (represented by an ASPH particle) exceeds a local critical strength threshold and begins to accumulate damage. Essentially, the strain model bridges the strength and damage models within the SPH framework. The Benz-Asphaug strain model (Benz & Asphaug, 1995) is used in many SPH codes, and the tensor generalization of this algorithm is available in Spheral, where E is the Young's modulus, S αβ is the deviatoric stress, and P is the pressure. The other strain model we consider follows the progression of the deviatoric stress (S αβ ) in the absence of plastic yielding: where μ i is the shear modulus. This is known as the "pseudo-plastic strain model," as it is intended to mimic the progression of plastic strain in the material. We find that the Benz-Asphaug strain model produces a characteristic spall wall and undamaged inner "core," as shown in Figure 3b (the right-hand panel illustrates the cross section of the target to its left, where the spall wall is green and intact core is blue), which was observed in the experiments ( Figure 3a). Note that the photograph shown on the right in Figure 3a is not the core fragment from the photograph shown to the left but rather an example core fragment from the replication of the experiment. In contrast, the pseudo-plastic strain model produces no spall wall, and, in addition, the damage propagates throughout the entire sphere (Figure 3c), unlike the experiment. We therefore choose to employ the Benz-Asphaug strain model in Spheral for further modeling of the experiment.

Strength Models and Parameter Selection
We investigate two specific strength models available for use in Spheral. The first is a pressure-dependent strength model (Collins et al., 2004), which is widely used for geologic materials simulated in shock physics and hydrodynamic codes. It describes how the shear strength of a rock increases from Y 0 (zero pressure) to Y m (the von Mises plastic limit at large confining pressures): where Y i is the yield strength of the rock at pressure P and μ i is the coefficient of internal friction (Collins et al., 2004;Lundborg, 1968). This strength model requires experimentally determined material values for μ i ,Y 0 , and Y m . Since these values are not available for the specific Yakuno basalt used in Nakamura and Fujiwara (1991), we use estimates of μ i ,Y 0 , and Y m for similar materials from previous studies. A summary of the basalt parameters used in our simulations is shown in Table 1. While the von Mises strength of basalt may vary between samples, a value of 3.5 GPa has been widely used in small-body fragmentation studies (Benz & Asphaug, 1999;Pierazzo et al., 2005;Senft & Stewart, 2007); hence, we utilize Y m = 3.5 GPa in our simulations. There is more variability present in the range of measured Y 0 values (Table 1), so we examine the sensitivity of our simulation results to the selected Y 0 for values of 66 MPa (Schultz, 1993), 130 MPa (Grady & Hollenbach, 1979), 600 MPa (basalt experiments by Stickle A. at the Applied Physics Laboratory, Johns Hopkins University), and 1 GPa by examining the largest-remaining target-fragment mass (M L ), also referred to as the "core." Y 0 equal to 1 GPa is not an experimentally measured parameter but rather is chosen for scanning purposes to approach the von Mises strength of basalt. These values of shear strength at zero pressure are well above the tensile strength of basalt, ranging between 11-30 MPa (Housen, 2009;Schultz, 1995) and measured by Nakamura et al. (2007) to be Y tensile = 19 MPa for the specific Yakuno basalt used in their 1991 experiment. We find that using Y 0 = 66 MPa and Y 0 = 130 MPa (Figures 4b and 4c) does not produce an M L that resembles the "core" recovered from the experiment (Figure 4a) but rather results in complete shattering of the basalt sphere into small fragments and dust. Expanding the strength investigation, we also compare the fragment masses and velocities from the simulation results to the experiment, using the pressure-dependent and constant-strength models. We find that using Y 0 = 600 MPa in the pressure-dependent strength model produces fragments with large spans in velocities 5.5-80 m/s (Figure 5a. Utilizing a value of Y 0 = 1 GPa slightly decreases the fragment velocity dispersion ( Figure 5b). When using Y m = 3.5 GPa and constant strength, the fragment velocity dispersion decreases further (Figure 5c), narrowing to a range of 18-50 m/s (not including the largest fragment). The lower values of Y 0 compared to Y m in the pressure-dependent strength model allow the material to have a wide range of strengths at varying pressures, resulting in a significant variance of failure strengths within the target. This

10.1029/2018EA000474
Earth and Space Science results in the wide span of fragment velocities seen in Figure 5a. By decreasing the range between Y 0 and Y m , the target is described by a more homogeneous strength; thus, the fragment velocities have a narrower distribution.
This analysis of the fragment masses and velocities shows that utilizing a yield strength larger than 600 MPa, where values are close or equal to the von Mises strength, eliminates much of the strength variance within the target so that the fragment velocities and masses closely resemble the experimental trend. Furthermore, the shear strength at zero pressure (Y 0 ) in the strength model should be greater than 600 MPa for the basalt material since this experimental value was determined by a dynamic tensile strength measurement. Because shear strength usually exceeds the tensile strength for a given material, values for the shear strength of the basalt should be larger than the tensile strength measurement of 600 MPa and the experimentally determined tensile strength of this specific Yakuno basalt of 19 MPa (Nakamura et al., 2007). In addition, our results are consistent with findings from Jutzi (2015), where the cohesion value (Y 0 = 100 MPa) for pumice utilized in their SPH code was also orders of magnitude larger than the measured tensile strength (σ t = 2.5 MPa) of the material, similar to our findings. We find that selecting a shear strength greater than 1 GPa produces an undamaged inner "core" resembling the M L recovered from the experiment shown in Figure 4, as well as fragment velocities that best resemble the experimental trend ( Figure 5).

Weibull Parameters
The deformation dynamics of a brittle solid depend upon its material properties. The two-parameter Weibull distribution is a statistical model used to describe flaws inherent in geologic materials, where N is the number density of flaws within a solid with failure strains less than ε, m is a shape parameter, and k is a scale parameter (Jaeger & Cook, 1979;Weibull, 1939). The m and k parameters are both experimentally measured for a given volume of material at a specific strain rate.
In 1991, when Nakamura and Fujiwara completed the hypervelocity gas-gun experiments on the Yakuno basalt targets, the m and k parameters were initially not measured. Three years later in 1994, Benz and Asphaug used the same set of experiments to help validate their SPH code. However, without the m and k values available, they performed Weibull parameter scans to find best fits to the experimental data. Their scans illuminated a relationship between sample volume, V, and the shape and scale parameters, where the simulation results did not change for values of m between 7 and 16 if ϕ = ln (kV)/m = 8.33 (±0.2). Using Weibull parameters m = 8.5 and k = 5.0 × 10 28 cm −3 for the basalt target (where ϕ = 8.33), Benz and Asphaug demonstrated their SPH code could produce an undamaged "core" matching the morphology and mass of M L from the experiments (Benz & Asphaug, 1994). In 2007, Nakamura teamed up with Michel and Setoh to measure the Weibull parameters of the original Yakuno basalt used in the 1991 experiments at different loading rates, determining a range of m values between 15-17 and k values between 10 51 -10 59 cm −3 (Nakamura et al., 2007). These experimental Weibull parameter measurements satisfy the dimensionless ratio proposed by Benz and Asphaug for catastrophic impacts, where ϕ ranges from 8.24 to 8.30 using the experimental m and k values, well within the 8.33 ± 0.2 bounds.
Employing the experimentally determined Weibull parameters for the Yakuno basalt in Spheral, we compare the cumulative mass distribution of fragments from the simulations to the experiment, using m = 17.2 and k = 3.43 × 10 59 cm −3 (Nakamura et al., 2007). For both the pressure-dependent and constant-strength models, we find that the simulation results underestimate M L when m = 17.2 and k = 3.43 × 10 59 cm −3 (Figure 6a). Conversely, implementing the Benz and Asphaug Weibull parameters, m = 8.5 and k = 5.0 × 10 28 cm −3 , overestimates the mass of the "core" (Figure 6b). It is likely that by varying the random seeding of flaws, M L will change, but the overall fragment distribution is expected to stay the same. Therefore, the comparison between the simulations and experimental results for this specific investigation is focused on the cumulative number distribution of the intermediate and smaller fragments.
As illustrated in Figure 6, using the experimentally measured Weibull parameters ( Figure 6a) produces a better fit to the measured fragment distribution than using the values determined through scans by Benz and Asphaug (Figure 6b) when comparing the cumulative number of intermediate and smaller fragments. Our results from Figure 6 indicate that using a shear strength at zero pressure greater than 1 GPa and Weibull parameters m = 17.2 and k = 3.43 × 10 59 for the basalt target yields the closest simulation results to the experimental findings. The pressure-dependent behavior is a more realistic model for rocks, although utilizing the von Mises strength model in Spheral performs similarly well in our simulations. However, it should be noted that it is likely that the tensile strength, which is controlled by the Weibull parameters, is playing a larger role than the shear strength in determining the final fragment distribution, which could explain why our shear strength-driven model slightly underproduces cumulative fragment number for intermediate fragment masses and overproduces the cumulative fragment number at the smallest fragment masses.

Conclusion
Potentially hazardous asteroids represent a low-probability but high-consequence risk. We depend on large-scale hydrodynamic codes to simulate mitigation options, requiring confidence that the codes are correct. Simulations are vetted through code validation by benchmarking against relevant laboratory-scale impact experiments. Comparisons between our numerical simulations and data from the Nakamura and Fujiwara (1991) basalt sphere impact experiment have provided new guidance for how to select appropriate strength and damage parameters for brittle, rocky materials for use in the Spheral hydrodynamics code.
We find that the Benz-Asphaug strain model is adequate to match the observed damage patterns in the experiment (preserved inner core and spallation at the antipodal from impact), whereas the pseudo-plastic strain model overestimates damage within the target. We also find that either a constant-strength model (Y m = 3.5 GPa) or a pressure-dependent strength with a relatively high value for Y 0 provides a reasonable fit to the fragment distribution. As discussed in section 3.3, our choice of such a large Y 0 in the pressure-dependent strength model may reflect a need to be closer in value to the von Mises strength, Y m = 3.5 GPa, in order to better approximate this particular basalt sample.
Finally, we determine that laboratory-based measurements of the Weibull parameters for the Yakuno basalt used in the experiment (Nakamura et al., 2007) provide a better fit to the fragment data than prior numerically determined Weibull parameters (Benz & Asphaug, 1994). It should be noted that these Weibull parameters were measured at strain rates well below the strain-rate regime of the 1991 experiment.
Nevertheless, our validation work shows that in this case, utilizing the experimentally measured Weibull parameters, albeit at a lower strain-rate regime, produces simulation results that compare well to the experimental fragment data, highlighting the value of performing careful characterization measurements on the target materials.
Ongoing efforts to prepare for the National Aeronautics and Space Administration directed Double Asteroid Redirection Test (DART) mission, a rare opportunity to field an impact deflection experiment (Cheng et al., 2018), have also emphasized the importance of code validation and verification (Stickle et al., 2016). This space mission will be the first demonstration of the kinetic-impactor technique. In 2022, a spacecraft will impact the secondary within the Didymos binary asteroid system at hypervelocity. Ground-based occultation measurements of Didymos will provide estimates of the resultant momentum transfer to the secondary (Cheng et al., 2018). We have preliminary indications, utilizing the findings from this work, that estimate a smaller momentum transfer than previously calculated for the DART impact. In a recent 2019 DART study, we compared our Spheral results using the pressure-dependent model and shear strength values used in this work to other codes. Although the initial conditions are unknown at this time for the Didymos system, which includes not knowing the material of the asteroid intended for deflection, basalt was the material of choice for the study. Results showed that the selection of strength model, and its parameters, had a substantial effect on the predicted crater size and momentum enhancement, more significant than the variation between codes (Stickle et al., 2019). This finding highlights the importance of model vetting in our code to have the needed confidence in our simulation results to correctly design a modeling plan for the upcoming DART mission.
In summary, we find that selecting the Benz-Asphaug strain model, a pressure-dependent strength model with Y 0 > 1 GPa, and Weibull m = 17.2 for basalt produced simulation results that closely resembled the experimental data. We recommend this approach for Spheral users and for consideration for use in other numerical models.