Asteroid impacts are destructive and low-probability threats to the Earth. The numerical simulation is considered an applicable analysis tool in asteroid deflection programs. As a novel shock-capturing strategy, the space–time conservation element and solution element (CESE) method can reliably predict shock waves and mechanical behaviors under high pressure and large strain conditions. In this paper, based on an elastoplastic flow model and an updated CESE scheme, the laboratory-scale iron asteroid impacts are modeled numerically, and the multi-material boundary treatment and the interface tracing strategy are introduced. Under hypervelocity impacts of the projectile to the iron asteroid target, the construction and realization of morphologies of impact craters and the implantation of projectile material into the target are numerically calculated. Numerical results show that the crater diameter and depth increase with increasing impact velocity and with increasing temperature, which softens the target. Computational results are compared with experimental observations available in the open literature, and good agreement is found. Therefore, the CESE method is successfully extended for capturing the key features of laboratory-scale hypervelocity asteroid impacts.

## I. INTRODUCTION

Asteroid impacts with the Earth are evidenced by several impact craters around the world. The catastrophic climate change, the topological formation, and the dinosaur extinction are attributed to amounts of impact scenarios.^{1} The shocked quartz crystal is the benchmark geochemistry material generated by shock pressures, showing the planar deformation features (e.g., shock lamellae).^{2–4} Formation of impact craters is characterized with the uplifted ring and basin-shaped bottom.^{5,6} Experiments of shock compression, such as quartz impacts, validate the equation of state and the Hugoniot^{7} and accurately elucidate the impact dynamics in laboratory observations.^{8–11}

Numerical simulations can be adopted to analyze the asteroid impact dynamics for assessing potential damages from triggered earthquakes, tsunamis, and airbursts. The dynamic response and the stress wave propagation under hypervelocity impacts are nonlinear.^{12} Therefore, many kinds of numerical models are developed.^{13–15} The elastic and plastic constitutive, the equation of state, and the particle flow rule are proposed for parameterization under hypervelocity impacts.^{16} Currently, the adaptive Lagrangian,^{17} the ENO (non-oscillatory schemes) shock capturing,^{18} the Godunov scheme,^{19} the finite element method,^{20} the space–time conservation element and solution element method,^{21–23} and the two-step Eulerian solution scheme^{24} are valuable tools for modeling the damage and deformation in impact simulations. The Earth impact effects program^{25} is established with mathematical schemes for predicting the fate of impactors during atmospheric entry, thermal radiation, and the intensity of ground shaking. The crater formation and the segmentation process, which are induced by shock waves, are reconstructed numerically and agree well with the observed large-scale crater formation.^{6,26} Physical properties of the shocked rock trapped within craters are generally consistent with that observed from the *in situ* hypervelocity impact experiments.^{11} Numerical simulations of deflecting asteroids by manmade-projectile impacts are performed to estimate the potential asteroid impact effects.^{14,27}

Asteroid impacts at a hypervelocity involve complex physics of flows, such as large deformation and phase changes.^{28} Under high pressures and temperatures, both the asteroid and the target will lose their own shear strength resisting the shear damage and exhibit fluid behaviors. Supplemented with pressure and temperature below the critical points, the projectile and the target can retrieve their own solid properties and show elastoplastic behaviors.^{29} Transition from the solid phase to the fluid phase can be described by an elastoplastic flow model,^{30,31} which can be formulated with compressible fluid Euler equations^{17,21,23} with the solid constitutive involved. Moreover, both local and global physics conservations are important for simulating shock waves triggered by hypervelocity impacts.^{18,32} Meanwhile, the contact boundaries configured between the projectile and the target determine the deformation and motion of craters.^{17,19,21} Thus, the multi-material interface configuration^{23,24,33,34} is vital for describing the asteroid impact dynamics.

As a new shock-capturing strategy, the space–time conservation element and solution element (CESE) method, as a novel conservation algorithm, can accurately predict shock waves and material behaviors under wide ranges of pressure and temperature.^{23,35–38} The hybrid particle level-set function enables high-resolution interface tracing and accurately calculates interface deformation.^{23,33,39,40} A renewed treatment of multi-material contact boundary^{21,23} can maintain the reliable physics flux across the material interface. Until now, the novel CESE method has never been utilized to solve the laboratory-scale asteroid impacts including complex elastoplastic flows with solid features. In this paper, we demonstrate advantages of the CESE scheme, the hybrid particle level-set function, and the multi-material treatment for modeling hypervelocity laboratory-scale asteroid impacts. We consider the numerical modeling of hypervelocity impact experiments in iron meteorites. Primary attention is paid to the construction and realization of morphologies of impact craters (for instance, crater diameter, depth, and floor-rim shape) and the implantation of the projectile material into the target. With the use of the novel CESE method and an elastoplastic flow model, computational results are compared with the experimental data. The potential is found of what can be accomplished for extending the CESE scheme and the elastoplastic flow model to the laboratory-scale asteroid-impact computational physics.

## II. METHODOLOGY

### A. Elastoplastic flow model

An elastoplastic flow model can be expressed as the Euler equation neglecting external forces, heat conduction, and thermal diffusion,^{21–23}

where $U=[\rho ,\rho u,\rho v,e,sxx,syy,szz,sxy,\epsilon P,T]T$, $E=[\rho u,\rho u2\u2212\sigma xx,\rho uv$$\u2212\sigma xy,ue\u2212\sigma xx\u2212v\sigma xy,usxx,usyy,uszz,usxy,u\epsilon P,uT]T$, $F=[\rho v,\rho uv$$\u2212\sigma xy,\rho v2\u2212\sigma yy,ve\u2212\sigma yy\u2212u\sigma xy,vsxx,vsyy,vszz,vsxy,v\epsilon P,vT]T$, and $S=[\eta \rho u,\eta \rho v2\u2212\sigma xx+\sigma zz,\eta \rho uv\u2212\sigma xy,\eta ve\u2212\sigma xx\u2212v\sigma xy,Ssxx,$$Ssyy,Sszz,Ssxy,S\epsilon P,ST]T$.

Here, ** U** is the conservative vector;

**and**

*E***are the corresponding flux vectors in the**

*F**x*and

*y*directions, respectively;

**is the source term vector;**

*S**ρ*is the density;

*u*and

*v*are the

*x*and

*y*components of the velocity;

*e*is the total volume energy;

*p*is the hydrostatic pressure;

*σ*

_{xx},

*σ*

_{xy}, and

*σ*

_{yy}are the components of the stress tensor;

*s*

_{xx},

*s*

_{xy}, and

*s*

_{yy}are the components of the deviatoric stress tensor;

*ɛ*

^{P}is the equivalent plastic strain; and

*T*is the temperature generated through the plastic work. Here, $\eta =\u2212ix$,

*i*= 0 for the plane strain problem and

*i*= 1 for the axisymmetric one. The deviatoric stress tensor

*s*

_{ij}denotes the shear deformation and obeys Hook’s law and the plastic flow model. The hydrostatic pressure

*p*reflects the volume deformation and is determined by the thermodynamic equation of state (EOS). Hook’s law in the Lagrangian representation is formulated in the incremental form

where *G* is the shear modulus, $s\u0307ij$ is the substantial derivative of *s*_{ij}, $D\u0307ij$ [$=12(\u2202ui\u2202xj+\u2202uj\u2202xi)$] is the strain rate tensor, and Ω_{ij} [$=12(\u2202ui\u2202xj\u2212\u2202uj\u2202xi)$] is the spin tensor. *u*_{i} and *x*_{i} are the *i*th velocity and coordinate components, respectively. The source terms of $Ssxx$, $Ssyy$, $Sszz$, and $Ssxy$ in Eq. (1) are given by

We denote as $\Omega xy=\u2212\Omega yx=12(\u2202u\u2202y\u2212\u2202v\u2202x)$ the spin tensor components and as $\Phi =13(\u2202u\u2202x+iux+\u2202v\u2202y)$ the volume strain.

The relationship between the plastic stress and the strain satisfies the von Mises yield criterion and the Prandtl–Reuss flow rule. The von Mises yield criterion condition has the form

where *σ*_{Y} is the current von Mises flow stress. We suppose that the plastic yield stress is given by the Johnson and Cook model, which determines the effects of the isotropic strain hardening, the strain rate, and the thermal softening.^{20,23} Hence, we obtain

where *A*, *B*, *C*, *m*, and *n* are the material constants. $\epsilon \u0307*$ [$=\epsilon \u0307P\epsilon \u03070$] is the relative equivalent plastic strain, and $\epsilon \u03070=1.0s\u22121$ for the SI unite system. $T*=T\u2212TrTm\u2212Tr$, where *T*_{m} and *T*_{r} are, respectively, the melt temperature of the solid material and the room temperature. *E* is the elastic modulus, and *ν* is the Poisson ratio. The temperature in the Johnson–Cook constitutive model^{41,42} reflects the thermal softening effect of materials and mainly results from the plastic work satisfying the following equation:

Here, we denote that *T* is the temperature and C is the specific heat. *β* represents the Taylor–Quinney parameter^{23,41,42} reflecting the fraction of converting mechanical power into thermal power and equals 0.9. $W\u0307P$$[\u2248\epsilon \u0307P32sijsij]$ is the plastic work rate, and $\epsilon \u0307P$$[=23D\u0307ijPD\u0307ijP]$ is the plastic work. The source terms of the equivalent plastic strain and the plastic temperature in Eq. (1) are given by

The Mie–Grüneisen equation of state describes the hydrostatic pressure of the solid phase under high pressures^{23} and is represented in the form

Here, $\mu [=\rho \rho 0\u22121]$; *s*_{1}, *s*_{2}, and *s*_{3} are the solid material parameters; and *α* is the first-order volume correction to the Gruneisen parameter *γ*_{0}.

### B. Free and contact boundary conditions

We assume that the normal stress components are zero^{18,19,31} as the physical conditions on a free surface. The mirror method^{21,23} is adopted for configuring the free boundary conditions, and we use the settings *σ*_{nn,P} = −*σ*_{nn,I}, *σ*_{nt,P} = −*σ*_{nt,I}, and *σ*_{tt,P} = *σ*_{tt,I}. We suppose that point *I* is the mirror point of *P* of the free surface. The subscripts *I* and *P* denote the values at point *I* and *P*, respectively. The normal stress components *σ*_{nn} and *σ*_{nt} are zero, and the tangential stress component *σ*_{tt} can be nonzero on the free surface. *n*_{x} and *n*_{y} are the components of the unit normal vector in the *x* and *y* directions. Components of the stress tensor in the local coordinates at the mirror point *I* have the form

Let us formulate a relationship of stress components between the local coordinates and global Cartesian coordinates $\sigma \u0302=AT\sigma A$, where $\sigma \u0302=\sigma nn\sigma nt\sigma tn\sigma tt$, $\sigma =\sigma xx\sigma xy\sigma yx\sigma yy$, and $A=nx\u2212nynynx$. We obtain the components of the stress tensor in global coordinates at point *P*,

In terms of multi-material contact boundary conditions, the normal velocity and the normal stress are continuous across the material interface, whereas the tangent velocity and the tangent stress can discontinue. It is assumed that, in the local coordinates, the velocity vector is $\nu \u0302=AT\nu $, where $\nu \u0302=[vnvt]$ and $\nu =[uv]$. On the sliding interface, the velocity at point *P* in the global coordinates can be expressed in the form

provided that *ν*_{n} = *u*_{P,2}*n*_{x} + *v*_{P,2}*n*_{y} and *ν*_{t} = *u*_{P,1}*n*_{y} − *v*_{P,1}*n*_{x}.

Subscript 1 denotes the linearly extrapolated values from material 1, and subscript 2 denotes the actual velocity at point P in material 2. We therefore have components of the stress tensor in the local coordinates

and hence, the components of the stress tensor in the global coordinates

With regard to the multi-material interface treatment, the hybrid particle level-set function^{33} can determine the calculated grid points in or off the material domain.^{21,23,40} The interface tracing and the free and contact boundary treatments are illustrated in Fig. 1. Here, the level-set functions^{40} are used for defining different material interfaces. The governing equations for each material are solved separately in the region $\phi (x\u20d7,t)<0$ of each level-set function. After that, the level-set functions are reevaluated in terms of the material velocity fields. The contact boundary points control the interaction between two material blocks. Meanwhile, the free boundary points describe the motion of the free boundary. A mesh point is considered as a free boundary point, given that it presents a boundary point of one material and does not locate within the other material. A mesh point is defined as a contact boundary point, provided that it remains a boundary point of a specific material and locates in the other material. With the help of the mirror method,^{21} the causal transformation of the local and global coordinates is performed, as described in Eqs. (9)–(13).

### C. Space–time conservation element and solution element method

The space–time conservation element and solution element (CESE) method has been developed to give high accurate solutions for hyperbolic conservation laws.^{35} The basic idea of the CESE scheme is the global and local flux conservation in a space–time domain. Both independent physics variables and their derivatives are treated as entities and solved simultaneously.^{36–38} The CESE method consists of the following procedures:^{39,40} (i) element discretization, (ii) expansion into the time dimension in the Euclidian space, (iii) definition of the solution element (SE) along which the flow variables are approximated by the Taylor series and solution elements advance through time, and (iv) definition of the conservation element (CE), in which the integral form of flow equations is integrated on the space–time domain (CE).

By virtue of the Gauss divergence theorem, neglecting source terms, we can rewrite expression (1) in the form^{36,38,40}

where $H=U,E,F$, $SV$ is the boundary of an arbitrary close space–time region *V*, and *d*** s** =

*d*

**·**

*σ***, with**

*n**d*

**and**

*σ***, respectively, being the area and the outward unit normal of a surface element on $SV$. Integrating the second-order Taylor expansion**

*n*^{36,40}on the boundary of each CE, we determine the relationship between the current conservative vector and that in the next half time step,

with the *m*th vector components of

Herein, subscripts *i* and *j* and superscript *n* denote the sequence number of space–time mesh in the *x*, *y*, and *t* directions, respectively. Symbols ±^{a} and ∓^{a} assume that one of them makes a positive operator, while the other keeps negative.^{38} Operator ±^{a} is independent of ±^{b}. Caret ^ represents a new defined function through integration and satisfies the equation

where *N*_{m} can be *U*_{m}, *E*_{m}, or *F*_{m}. To obtain the value of $Umi,jn+1/2$, it is required to identify the second derivative $Umxxi,jn+1/2$ and $Umyyi,jn+1/2$ in the spatial direction of the current space–time point (*i*, *j*, and *n*). In the case of $Umxxi,jn+1/2$ and $Umxyi,jn+1/2$, the differential formulas are given by

Hence, the value of $Ui,jn+1/2$ can be found with the use of Eq. (15). Given that the continuity of variables in SE at the junction point,^{35,38,40} we obtain the *x*, *y* directional derivatives of conservative vector in the next time step.^{39} The source vectors in Eq. (1) can be treated by the linear approximation.^{21,23,40}

In this work, we updated the previous computer code^{40} based on the CESE method of the first-order accuracy. A novel computer code has been developed with the second-order accuracy [as described in Eq. (15)], and the diagram of the computational programming is shown in Algorithm 1. Meanwhile, the previous computer code of the hybrid particle level-set function^{40} is employed to perform the multi-material interface tracing and to calculate the evolution of the free and contact boundaries. With regard to the multi-material interface treatment, the hybrid particle level-set function can determine the calculated grid points in or off the material domain.^{21,23,40} It is worthy of noticing that the updated computer code introduces the second-order accuracy CESE scheme to the calculation of the elastoplastic flow model, and currently, we do not consider any damage mechanics models.

(1) Initialize. |

{Define global mesh grids, material parameters, level-set functions, and physics variables on grids} |

(2) Define the time step in terms of the CFL number. |

(3) Set up the boundary of each material with the level-set functions. |

(4) Solve the governing equation of each material using the CESE scheme. |

(5) Reevaluate the temporal evolution of the level-set functions through the flow fields. |

(6) Solve the evolution of the Lagrangian particle equation. |

(7) Update the level-set functions according to the non-mass particles. |

(8) Reinitialize the level-set functions. |

(9) Rectify the level-set functions by virtue of the non-mass particles. |

(10) Output the calculated results. |

(11) If the final calculation time is approached, GOTO (12). Otherwise GOTO (2). |

(12) Stop. |

(1) Initialize. |

{Define global mesh grids, material parameters, level-set functions, and physics variables on grids} |

(2) Define the time step in terms of the CFL number. |

(3) Set up the boundary of each material with the level-set functions. |

(4) Solve the governing equation of each material using the CESE scheme. |

(5) Reevaluate the temporal evolution of the level-set functions through the flow fields. |

(6) Solve the evolution of the Lagrangian particle equation. |

(7) Update the level-set functions according to the non-mass particles. |

(8) Reinitialize the level-set functions. |

(9) Rectify the level-set functions by virtue of the non-mass particles. |

(10) Output the calculated results. |

(11) If the final calculation time is approached, GOTO (12). Otherwise GOTO (2). |

(12) Stop. |

### D. Method evaluation and validation

As a typical example of a two dimensional simulation, we consider the problem of a normal collision of a cylindrical rod with a rigid wall, called the Taylor impact problem, which in a similar statement is performed in Refs. 17, 19, and 23. During the collision, the kinetic energy is irreversibly converted into an internal one.^{19} The equivalent plastic strain of the rod is obtained in the form of a thickening in the contact area. It is assumed that the Taylor impact problem is an axisymmetric one. In view of cylindrical coordinates, the computational domain is rectangular and represents a half section of the rod along the radius.^{19} The rigid wall is treated as the fixed boundary in the *x* direction, and the rest boundaries are considered as free surfaces.^{23}

A grid independence analysis is investigated using 80 × 300, 120 × 500, and 160 × 700 meshes. Figure 2 depicts the final configurations at the time of 80 *µ*s for all the meshes considered. It is seen that the solution converges with grid refinement, and grid-independent results are approached for the finest mesh, 160 × 700.

Here, the computational region is (8 × 35 mm^{2}) with 160 × 700 meshes corresponding to Δ*x* = Δ*y* = 0.05 mm. We set the elastic modulus of 117 GPa, the Poisson ratio of 0.35, the tangent hardening modulus of 100 MPa, and the yield strength of 400 MPa. The linear isotropic hardening constitutive model is considered. The parameters of the Mie–Grüneisen equation of state are summarized in Table I. Figures 3 and 4 show the deformation, the equivalent plastic strain, and the von Mises flow stress of the copper rod at the times of 20, 40, 60, and 80 *µ*s. At time instant *t* = 80 *µ*s, an almost complete conversion of the kinetic energy into the internal one is found, which indicates a stop of the copper rod. The final length, the maximum radius, and the maximum equivalent plastic strain are 21.5, 7.24, and 2.73 mm, respectively. It is found that the maximum value appears at the axisymmetric axis on the rigid wall, and a large gradient of the equivalent plastic strain exists near the rigid wall. These results are in good agreement with the findings of Refs. 17, 19, and 50.

Camacho and Ortiz^{17} demonstrated the deformed configurations with adaption and full coarsening, which are calculated by the adaptive Lagrangian modeling at the times of 20, 40, and 80 *µ*s (as shown in Fig. 8 of Ref. 17). Tran and Udaykumar^{50} illustrated the equivalent plastic strain contours at different times using the 80 × 400 mesh, which are performed with an Eulerian, sharp interface, and Cartesian grid method along with the ENO scheme (as shown in Fig. 9 of Ref. 50). It is clear that these results agree well with the findings of Refs. 17 and 50.

## III. MODELING LABORATORY-SCALE ASTEROID IMPACTS

As motivated by the NASA Psyche mission,^{15} we consider the numerical modeling of hypervelocity impact experiments in iron meteorites.^{43} Primary attention is paid to the construction and realization of morphologies of impact craters (for instance, crater diameter, depth, and floor/rim shape) and the implantation of projectile material into the target. With the use of the novel CESE method and the elastoplastic flow model, computational results are compared with the experimental data.

### A. Morphology of the impact craters

We assume that the composition of iron meteorites is equivalent to the fine octahedrite-structured Gibeon iron, and the quartz represents an impactor selection. Under experimental conditions,^{43} a spherical quartz projectile with a diameter of 3.175/6.35 mm perpendicularly (head-on) penetrates a planar iron meteorite target at a hypervelocity. The iron meteorite target is cubic with (40 × 40 × 40 mm^{3}). The incident speeds are varied between 2940 and 5300 m/s. The initial bulk density is 2200 kg/m^{3} for the quartz and 7900 kg/m^{3} for the Gibeon meteorite. The Johnson–Cook constitutive model and the Mie–Grüneisen equation of state are adopted to model the quartz projectile and the Gibeon target. The model parameters are shown in Tables I and II. The impact process is modeled by the axisymmetric problem. The computational domain is (0.03 × 0.06 m^{3}) with 300 × 600 meshes, which corresponds to Δ*x* = Δ*y* = 0.1 mm. The left boundary is treated as the axisymmetric one, and the other three sides are the free surface boundaries. The calculation is complemented until an almost complete conversion of kinetic energy into the internal one is obtained,^{19,21,23} which means a stop of the projectile. The initial temperature of the target is similar to the experimental conditions.^{43}

Material . | ρ (kg/m^{3})
. | E (GPa)
. | ν
. | A (MPa)
. | B (MPa)
. | n
. | C
. | m
. | S (J/kg K)
. | T_{m} (K)
. | Tr (K)
. |
---|---|---|---|---|---|---|---|---|---|---|---|

Quartz^{48} | 2200 | 72 | 0.17 | 930 | 88 | 0.77 | 0.003 | 0.35 | 741 | 1670 | 298 |

Gibeon^{46} | 7900 | 209 | 0.3 | 676 | 1015 | 0.53 | 0.027 | 0.52 | 473 | 1783 | * |

Figure 5 shows the calculated equivalent plastic strain and the deformation at the times of 20, 40, 60, and 80 *µ*s for the AVGR experiments [shot ID180408]. Figure 6 shows the impact crater morphology of the AVGR experiment [shot ID180408] and the final shape calculated by the CESE method.

Figures 7 and 8 show the final deformed configuration, e.g., morphology of impact craters, with various impact velocities. The crater deformation, the excavation, and the implantation of projectile material into the target are observed. As shown, the spherical quartz projectile is squashed. The target mass is pushed upward and finally ejects out by the shock pressures to form the crater rims. Of particular interest are the crater diameter, depth, floor, and rim shape. The calculated crater morphology is characterized with the sharp, raised rim and the deep cavity, which is in agreement with the finding of the experiments.^{43} The crater diameter and depth, which are measured from the pre-impact surface, are summarized in Table III. With regard to numerical simulations of the five head-on impact experiments, we obtain the mean depth/diameter ratio *d*_{D} = 0.35 ± 0.05, which is in accordance with the experimental one *d*_{D} = 0.39 ± 0.04. As shown in Fig. 9, the depth to diameter ratio of head-on impacts between the experimental data and the numerically simulated results has been found to be highly correlated. The depth/diameter ratio significantly depends on the impact velocity and is higher at room temperature than at low temperatures.^{46} The crater diameter and depth increase with increasing impact velocity and with increasing temperature, which softens the target. We believe crater shape as a possible record of the impact environment of metallic bodies.^{46} The numerical simulations prove that the Johnson–Cook constitutive model can provide a reasonable interpretation for the laboratory-scale iron meteorite impact at the high strain rate and large strain.^{43,46,49}

AVGR shot ID . | Crater^{43} diam. (mm)
. | Crater^{43} depth (mm)
. | Simulated crater diam. (mm) . | Simulated crater depth (mm) . | Maximum plastic strain . | Impact^{43} velocity (km/s)
. | Target^{43} temperature (K)
. |
---|---|---|---|---|---|---|---|

180 215 | 6.9 | 2.5 | 6.5 | 2.0 | 2.39 | 3.88 | 296 |

180 218 | 5.3 | 1.8 | 5.4 | 1.6 | 1.86 | 2.94 | 109 |

180 222 | 6.9 | 3.0 | 7.1 | 2.7 | 2.86 | 5.30 | 116 |

180 401 | 8.0 | 3.2 | 7.3 | 3.1 | 3.64 | 5.23 | 292 |

180 408 | 14.5 | 6.0 | 15.4 | 5.6 | 4.15 | 5.30 | 122 |

AVGR shot ID . | Crater^{43} diam. (mm)
. | Crater^{43} depth (mm)
. | Simulated crater diam. (mm) . | Simulated crater depth (mm) . | Maximum plastic strain . | Impact^{43} velocity (km/s)
. | Target^{43} temperature (K)
. |
---|---|---|---|---|---|---|---|

180 215 | 6.9 | 2.5 | 6.5 | 2.0 | 2.39 | 3.88 | 296 |

180 218 | 5.3 | 1.8 | 5.4 | 1.6 | 1.86 | 2.94 | 109 |

180 222 | 6.9 | 3.0 | 7.1 | 2.7 | 2.86 | 5.30 | 116 |

180 401 | 8.0 | 3.2 | 7.3 | 3.1 | 3.64 | 5.23 | 292 |

180 408 | 14.5 | 6.0 | 15.4 | 5.6 | 4.15 | 5.30 | 122 |

### B. Discussion of cracking and spallation

Figure 10 shows the x-ray scanning section^{43} of the damaged target of the AVGR experiments [shot ID180408] and the final pressure distribution simulated by the CESE method. As expected, the shock wave and the rarefaction wave are generated in the course of impact, which may lead to radial distribution of cracks around the crater under the large strain.^{43} The rarefaction waves generated at the free surface affect the spall diameter and the total crater volume.^{49} In the AVGR experiments [shot ID180408], it is observed that the cracking and spallation initiate from the crater and propagate toward the edge of the surface during the excavation (as shown in Fig. 10, left panel). The spalling phenomenon is attributed to the interaction of the unloading wave with the rarefaction wave reflected from the free surface.^{19,21,49} As usual, an increase in the surface velocity results from the shock wave arrival to the near surface of the target, and the shock wave reflection from the surface leads to the appearance of a rarefaction wave.^{19} Subsequently, following the shock wave, an incident unloading wave arrives at the surface^{32} and results in a decrease in the surface velocity. Hence, the damage of the target near the rear surface is associated with the interaction of the opposing rarefaction and the unloading wave.^{44–46} We leave to future work to further assess the crack growth nucleation^{11,21,49} problems.

## IV. CONCLUSIONS

With the use of the elastoplastic flow model and the CESE scheme, we numerically simulate the laboratory-scale iron asteroid impacts and characterize the crater morphology generated by stress waves, which are dependent on the yield criterion, the strain model, and the equation of state. The crater deformation, the excavation, and the projectile implantation are numerically realized in the course of impact. The crater morphology has the sharp, raised rim and the deep cavity. The calculated depth/diameter ratio of impact craters is consistent with the experimental one. Numerical results show that the crater diameter and depth increase with increasing impact velocity and with increasing temperature, which thermally softens the target. We believe crater shape as a possible record of the impact environment of metallic bodies. Moreover, the target spallation is discussed, which results from the interaction of the unloading wave with the rarefaction wave reflected from the free surface. Numerical simulations confirm that the Johnson–Cook constitutive model can provide a feasible explanation for the laboratory-scale iron meteorite impact at the high strain rate and large strain. Hence, the shock wave and the crater deformation during laboratory-scale asteroid impacts are successfully simulated by the elastoplastic flow model and the high-accuracy CESE scheme. These results might provide new insights into laboratory-scale asteroid impact simulations. It notes that the updated computer code introduces the second-order accuracy CESE scheme to the calculation of the elastoplastic flow model, and, currently, we do not consider any damage mechanics models.

## ACKNOWLEDGMENTS

The author is grateful to senior Professor De-Liang Zhang (Institute of Mechanics, Chinese Academy of Sciences) for his support and encouragements, especially for his discussion on developing the novel computer code of the second-order accuracy CESE method.

This work was funded by the National Natural Science Foundation of China (Grant No. 41874113) and research grants from the National Institute of Natural Hazards, Ministry of Emergency Management of China (Grant Nos. ZDJ2018-20 and ZDJ2019-15).

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

## DATA AVAILABILITY

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

## REFERENCES

_{2}