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.

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 Hugoniot7 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 scheme24 are valuable tools for modeling the damage and deformation in impact simulations. The Earth impact effects program25 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 equations17,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 configuration23,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 boundary21,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.

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

(1)

where U=[ρ,ρu,ρv,e,sxx,syy,szz,sxy,εP,T]T, E=[ρu,ρu2σxx,ρuvσxy,ueσxxvσxy,usxx,usyy,uszz,usxy,uεP,uT]T, F=[ρv,ρuvσxy,ρv2σyy,veσyyuσxy,vsxx,vsyy,vszz,vsxy,vεP,vT]T, and S=[ηρu,ηρv2σxx+σzz,ηρuvσxy,ηveσxxvσxy,Ssxx,Ssyy,Sszz,Ssxy,SεP,ST]T.

Here, U is the conservative vector; E and F are the corresponding flux vectors in the x and y directions, respectively; S is the source term vector; ρ 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; sxx, sxy, and syy 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, η=ix, i = 0 for the plane strain problem and i = 1 for the axisymmetric one. The deviatoric stress tensor sij 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

(2)

where G is the shear modulus, ṡij is the substantial derivative of sij, Ḋij [=12(uixj+ujxi)] is the strain rate tensor, and Ωij [=12(uixjujxi)] is the spin tensor. ui and xi are the ith velocity and coordinate components, respectively. The source terms of Ssxx, Ssyy, Sszz, and Ssxy in Eq. (1) are given by

(3)

We denote as Ωxy=Ωyx=12(uyvx) the spin tensor components and as Φ=13(ux+iux+vy) 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

(4)

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

(5)

where A, B, C, m, and n are the material constants. ε̇* [=ε̇Pε̇0] is the relative equivalent plastic strain, and ε̇0=1.0s1 for the SI unite system. T*=TTrTmTr, where Tm and Tr 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 model41,42 reflects the thermal softening effect of materials and mainly results from the plastic work satisfying the following equation:

(6)

Here, we denote that T is the temperature and C is the specific heat. β represents the Taylor–Quinney parameter23,41,42 reflecting the fraction of converting mechanical power into thermal power and equals 0.9. ẆP[ε̇P32sijsij] is the plastic work rate, and ε̇P[=23ḊijPḊijP] is the plastic work. The source terms of the equivalent plastic strain and the plastic temperature in Eq. (1) are given by

(7)

The Mie–Grüneisen equation of state describes the hydrostatic pressure of the solid phase under high pressures23 and is represented in the form

(8)

Here, μ[=ρρ01]; s1, s2, and s3 are the solid material parameters; and α is the first-order volume correction to the Gruneisen parameter γ0.

We assume that the normal stress components are zero18,19,31 as the physical conditions on a free surface. The mirror method21,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. nx and ny 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

(9)

Let us formulate a relationship of stress components between the local coordinates and global Cartesian coordinates σ̂=ATσA, where σ̂=σnnσntσtnσtt, σ=σxxσxyσyxσyy, and A=nxnynynx. We obtain the components of the stress tensor in global coordinates at point P,

(10)

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 ν̂=ATν, where ν̂=[vnvt] and ν=[uv]. On the sliding interface, the velocity at point P in the global coordinates can be expressed in the form

(11)

provided that νn = uP,2nx + vP,2ny and νt = uP,1nyvP,1nx.

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

(12)

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

(13)

With regard to the multi-material interface treatment, the hybrid particle level-set function33 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 functions40 are used for defining different material interfaces. The governing equations for each material are solved separately in the region φ(x,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).

FIG. 1.

Schematics of the multi-material interface tracing and the free and contact boundary treatment, which are implemented through the use of the hybrid particle level-set function.

FIG. 1.

Schematics of the multi-material interface tracing and the free and contact boundary treatment, which are implemented through the use of the hybrid particle level-set function.

Close modal

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 form36,38,40

(14)

where H=U,E,F, SV is the boundary of an arbitrary close space–time region V, and ds = dσ · n, with dσ and n, respectively, being the area and the outward unit normal of a surface element on SV. Integrating the second-order Taylor expansion36,40 on the boundary of each CE, we determine the relationship between the current conservative vector and that in the next half time step,

(15)

with the mth vector components of

(16)

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

(17)

where Nm can be Um, Em, or Fm. 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

(18)

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 code40 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 function40 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.

ALGORITHM 1.

Process flow diagram of the computational modeling that incorporates the definition of parameters of the constitutive model and the calculation of the hybrid particle level-set function along with the CESE algorithm of the second-order accuracy.

(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. 

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.

FIG. 2.

Mesh convergence study in Taylor bar impact calculation with a velocity of 227 m/s. Final shapes of 80 × 300 (black line), 120 × 500 (blue line), and 160 × 700 (orange line) meshes.

FIG. 2.

Mesh convergence study in Taylor bar impact calculation with a velocity of 227 m/s. Final shapes of 80 × 300 (black line), 120 × 500 (blue line), and 160 × 700 (orange line) meshes.

Close modal

Here, the computational region is (8 × 35 mm2) 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.

TABLE I.

Parameters for the Mie–Grüneisen equation of state.

Materialc0 (m/s)γ0s1s2s3A
Quartz47  5900 0.7 1.85 
Gibeon48  4050 1.69 1.92 
Materialc0 (m/s)γ0s1s2s3A
Quartz47  5900 0.7 1.85 
Gibeon48  4050 1.69 1.92 
FIG. 3.

Equivalent plastic strain distributions using the 160 × 700 mesh at the times of 20, 40, 60, and 80 µs.

FIG. 3.

Equivalent plastic strain distributions using the 160 × 700 mesh at the times of 20, 40, 60, and 80 µs.

Close modal
FIG. 4.

von Mises flow stress (Pa) distributions using the 160 × 700 mesh at the times of 20, 40, 60, and 80 µs.

FIG. 4.

von Mises flow stress (Pa) distributions using the 160 × 700 mesh at the times of 20, 40, 60, and 80 µs.

Close modal

Camacho and Ortiz17 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 Udaykumar50 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.

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.

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 mm3). The incident speeds are varied between 2940 and 5300 m/s. The initial bulk density is 2200 kg/m3 for the quartz and 7900 kg/m3 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 m3) 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 

TABLE II.

Parameters of the Johnson–Cook constitutive model. Note that symbol * means we set Tr the experimental temperature to simulate a low temperature target.

Materialρ (kg/m3)E (GPa)νA (MPa)B (MPa)nCmS (J/kg K)Tm (K)Tr (K)
Quartz48  2200 72 0.17 930 88 0.77 0.003 0.35 741 1670 298 
Gibeon46  7900 209 0.3 676 1015 0.53 0.027 0.52 473 1783 
Materialρ (kg/m3)E (GPa)νA (MPa)B (MPa)nCmS (J/kg K)Tm (K)Tr (K)
Quartz48  2200 72 0.17 930 88 0.77 0.003 0.35 741 1670 298 
Gibeon46  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.

FIG. 5.

Contours of the equivalent plastic strain at the times of 20, 40, 60, and 80 µs with the impact velocity of 5300 m/s for the projectile diameter of 6.35 mm, representing the AVGR experiment [shot ID180408].

FIG. 5.

Contours of the equivalent plastic strain at the times of 20, 40, 60, and 80 µs with the impact velocity of 5300 m/s for the projectile diameter of 6.35 mm, representing the AVGR experiment [shot ID180408].

Close modal
FIG. 6.

Impact crater morphology from the AVGR experiment [shot ID180408]. Left panel: top-view of the experimental crater topographic profile [reprinted with permission from Marchi et al., “Hypervelocity impact experiments in iron-nickel ingots and iron meteorites: Implications for the NASA Psyche mission,” J. Geophys. Res.: Planets 125, e2019JE005927 (2020). Copyright 2020 The American Geophysical Union]. Right panel: final crater topographic profile calculated. Black denotes the iron meteorite target, and gray denotes the quartz projectile.

FIG. 6.

Impact crater morphology from the AVGR experiment [shot ID180408]. Left panel: top-view of the experimental crater topographic profile [reprinted with permission from Marchi et al., “Hypervelocity impact experiments in iron-nickel ingots and iron meteorites: Implications for the NASA Psyche mission,” J. Geophys. Res.: Planets 125, e2019JE005927 (2020). Copyright 2020 The American Geophysical Union]. Right panel: final crater topographic profile calculated. Black denotes the iron meteorite target, and gray denotes the quartz projectile.

Close modal

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 dD = 0.35 ± 0.05, which is in accordance with the experimental one dD = 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

FIG. 7.

Equivalent plastic strain contours at the final time with impact velocity of (a) 3880 m/s, (b) 2940 m/s, (c) 5300 m/s, and (d) 5230 m/s for the projectile diameter of 3.175 mm. The maximum value locates at the axisymmetric axis on the rigid wall, and a large gradient of the equivalent plastic strain appears near the material interface.

FIG. 7.

Equivalent plastic strain contours at the final time with impact velocity of (a) 3880 m/s, (b) 2940 m/s, (c) 5300 m/s, and (d) 5230 m/s for the projectile diameter of 3.175 mm. The maximum value locates at the axisymmetric axis on the rigid wall, and a large gradient of the equivalent plastic strain appears near the material interface.

Close modal
FIG. 8.

Deformation manifested in the inlet at the final time with impact velocity of (a) 3880 m/s, (b) 2940 m/s, (c) 5300 m/s, and (d) 5230 m/s for the projectile diameter of 3.175 mm. Black denotes the iron meteorite target, and gray denotes the quartz projectile. The simulation is conducted up to 80 µs when the velocities of the projectile and target nearly equal zero.

FIG. 8.

Deformation manifested in the inlet at the final time with impact velocity of (a) 3880 m/s, (b) 2940 m/s, (c) 5300 m/s, and (d) 5230 m/s for the projectile diameter of 3.175 mm. Black denotes the iron meteorite target, and gray denotes the quartz projectile. The simulation is conducted up to 80 µs when the velocities of the projectile and target nearly equal zero.

Close modal
TABLE III.

Details on AVGR impact experiments47 and numerical results of CESE simulations. The crater diameter and depth are measured from the pre-impact surface.

AVGR shot IDCrater43 diam. (mm)Crater43 depth (mm)Simulated crater diam. (mm)Simulated crater depth (mm)Maximum plastic strainImpact43 velocity (km/s)Target43 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 IDCrater43 diam. (mm)Crater43 depth (mm)Simulated crater diam. (mm)Simulated crater depth (mm)Maximum plastic strainImpact43 velocity (km/s)Target43 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 
FIG. 9.

Depth to diameter ratio of head-on impacts.

FIG. 9.

Depth to diameter ratio of head-on impacts.

Close modal

Figure 10 shows the x-ray scanning section43 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 surface32 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 nucleation11,21,49 problems.

FIG. 10.

Left panel: x-ray scan of the impacted target of the AVGR experiment [reprinted with permission from Marchi et al., “Hypervelocity impact experiments in iron-nickel ingots and iron meteorites: Implications for the NASA Psyche mission,” J. Geophys. Res.: Planets 125, e2019JE005927 (2020). Copyright 2020 The American Geophysical Union]. It shows that the target is heavily fractured due to the shocks. Right panel: final pressure distribution with the impact velocity of 5300 m/s for the projectile diameter of 6.35 mm, representing the CESE calculated result.

FIG. 10.

Left panel: x-ray scan of the impacted target of the AVGR experiment [reprinted with permission from Marchi et al., “Hypervelocity impact experiments in iron-nickel ingots and iron meteorites: Implications for the NASA Psyche mission,” J. Geophys. Res.: Planets 125, e2019JE005927 (2020). Copyright 2020 The American Geophysical Union]. It shows that the target is heavily fractured due to the shocks. Right panel: final pressure distribution with the impact velocity of 5300 m/s for the projectile diameter of 6.35 mm, representing the CESE calculated result.

Close modal

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.

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).

The author has no conflicts to disclose.

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

1.
W.
Benz
and
E.
Asphaug
, “
Catastrophic disruptions revisited
,”
Icarus
142
(
1
),
5
20
(
1999
).
2.
J. J.
Petrovic
, “
Mechanical properties of meteorites and their constituents
,”
J. Mater. Sci.
36
,
1579
1583
(
2001
).
3.
D. A.
Kring
,
A. R.
Hildebrand
, and
W. V.
Boynton
, “
Provenance of mineral phases in the Cretaceous-Tertiary boundary sediments exposed on the southern peninsula of Haiti
,”
Earth Planet. Sci. Lett.
128
,
629
641
(
1994
).
4.
A. E.
Gleason
,
C. A.
Bolme
,
H. J.
Lee
et al., “
Ultrafast visualization of crystallization and grain growth in shock-compressed SiO2
,”
Nat. Commun.
6
,
8191
(
2015
).
5.
M. K.
Shepard
et al., “
Radar observations and shape model of asteroid 16 Psyche
,”
Icarus
281
,
388
403
(
2017
).
6.
B.
den Hond
, “
New simulation supports Chicxulub impact scenario
,”
Eos
99
(27 April
2018
).
7.
S.
Root
,
J. P.
Townsend
, and
M. D.
Knudson
, “
Shock compression of fused silica: An impedance matching standard
,”
J. Appl. Phys.
126
,
165901
(
2019
).
8.
D.
Stoeffler
,
D. E.
Gault
,
J.
Wedekind
, and
G.
Polkowski
, “
Experimental hypervelocity impact into quartz sand: Distribution and shock metamorphism of ejecta
,”
J. Geophys. Res.
80
(
29
),
4062
4077
, (
1975
).
9.
A.
Nakamura
and
A.
Fujiwara
, “
Velocity distribution of fragments formed in a simulated collisional disruption
,”
Icarus
92
(
1
),
132
146
(
1991
).
10.
R. T.
Daly
and
P. H.
Schultz
, “
Projectile preservation during oblique hypervelocity impacts
,”
Meteorit. Planet. Sci.
53
(
7
),
1364
1390
(
2018
).
11.
N. M.
Masahiro
,
K.
Hayashi
,
J.
Nakagawa
, and
Y.
Ito
, “
Influence of temperature on crater and ejecta size following hypervelocity impact of aluminum alloy spheres on thick aluminum alloy targets
,”
Int. J. Impact Eng.
42
,
37
47
(
2012
).
12.
G. S.
Collins
,
H. J.
Melosh
, and
B. A.
Ivanov
, “
Modeling damage and deformation in impact simulations
,”
Meteorit. Planet. Sci.
39
(
2
),
217
231
(
2004
).
13.
D. D.
Durda
,
R.
Greenberg
, and
R.
Jedicke
, “
Collisional models and scaling laws: A new interpretation of the shape of the main-belt asteroid size distribution
,”
Icarus
135
(
2
),
431
440
(
1998
).
14.
T. P.
Remington
,
J. M.
Owen
,
A. M.
Nakamura
et al., “
Numerical simulations of laboratory-scale, hypervelocity-impact experiments for asteroid-deflection code validation
,”
Earth Space Sci.
7
,
e2018EA000474
(
2020
).
15.
S. D.
Raducan
,
T. M.
Davison
, and
G. S.
Collins
, “
Morphological diversity of impact craters on asteroid 2(16) Psyche: Insight from numerical models
,”
J. Geophys. Res.: Planets
125
,
e2020JE006466
, (
2020
).
16.
J.
Michael Owen
, “
A compatibly differenced total energy conserving form of SPH
,”
Int. J. Numer. Methods Fluids
75
(
11
),
749
774
(
2014
).
17.
G. T.
Camacho
and
M.
Ortiz
, “
Adaptive Lagrangian modeling of ballistic penetration of metallic targets
,”
Comput. Methods Appl. Mech. Eng.
142
,
269
301
(
1997
).
18.
H. S.
Udaykumar
,
L.
Tran
,
D. M.
Belk
, and
K. J.
Vanden
, “
An Eulerian method for computation of multimaterial impact with ENO shock capturing and sharp interfaces
,”
J. Comput. Phys.
186
(
1
),
136
177
(
2003
).
19.
I. S.
Menshov
,
A. V.
Mischenko
, and
A. A.
Serejkin
, “
Numerical modeling of elastoplastic flows by the Godunov method on moving Eulerian grids
,”
Math. Models Comput. Simul.
6
,
127
141
(
2014
).
20.
S.
Kobayashi
,
S. I.
Oh
, and
T.
Altan
,
Metal Forming and the Finite Element Method
(
Oxford University Press
,
1989
).
21.
Q.
Chen
,
J.
Wang
, and
K.
Liu
, “
Improved CE/SE scheme with particle level set method for numerical simulation of spall fracture due to high-velocity impact
,”
J. Comput. Phys.
229
(
19
),
7503
7519
(
2010
).
22.
K. X.
Liu
,
S. Y.
Wu
, and
Q. Y.
Chen
, “
Numerical analysis of hypervelocity impact problems with large deformations, high strain rates and spall fractures
,”
Key Eng. Mater.
535–536
,
441
444
(
2013
).
23.
J.
Wang
,
K.
Liu
, and
D.
Zhang
, “
An improved CE/SE scheme for multi-material elastic-plastic flows and its applications
,”
Comput. Fluids
38
(
3
),
544
551
(
2009
).
24.
C. E.
Anderson
, Jr.
,
C.
Sidney
,
P.
Rory
, and
Bigger
, “
Time-resolved penetration into glass: Experiments and computations
,”
Int. J. Impact Eng.
38
,
723
731
(
2011
).
25.
G. S.
Collins
,
H. J.
Melosh
, and
R. A.
Marcus
, “
Earth impact effects program: A web-based computer program for calculating the regional environmental consequences of a meteoroid impact on Earth
,”
Meteorit. Planet. Sci.
40
(
6
),
817
840
(
2005
).
26.
L. E.
Senft
and
S. T.
Stewart
, “
Modeling impact cratering in layered surfaces
,”
J. Geophys. Res.
112
,
E11002
, (
2007
).
27.
A. M.
Stickle
,
E. S. G.
Rainey
,
M. B.
Syal
et al., “
Modeling impact outcomes for the double asteroid redirection test (DART) mission
,”
Procedia Eng.
204
,
116
123
(
2017
).
28.
W. F.
Bottke
, Jr.
,
D.
Vokrouhlický
,
D. P.
Rubincam
, and
D.
Nesvorný
, “
The Yarkovsky and Yorp effects: Implications for asteroid dynamics
,”
Annu. Rev. Earth Planet. Sci.
34
,
157
191
(
2006
).
29.
M. L.
Wilkins
, “
Calculation of elastic-plastic flow
,” in
Methods in Computational Physics
, edited by
B.
Adler
,
S.
Fembach
, and
M.
Rotenberg
(
Academic Press
,
New York
,
1964
), Vol. 3, pp.
211
263
.
30.
M. B.
Tyndall
, “
Numerical modeling of shocks in solids with elastic-plastic conditions
,”
Shock Waves
3
,
55
66
(
1993
).
31.
A.
Kapahi
,
S.
Sambasivan
, and
H. S.
Udaykumar
, “
A three-dimensional sharp interface Cartesian grid method for solving high speed multi-material impact, penetration and fragmentation problems
,”
J. Comput. Phys.
241
,
308
332
(
2013
).
32.
N. N.
Smirnov
,
A. B.
Kiselev
, and
P. P.
Zakharov
, “
Numerical simulation of the hypervelocity impact of the ball and the spherical containment in three-material statement
,”
Acta Astronaut.
171
,
215
224
(
2020
).
33.
D.
Enright
,
R.
Fedkiw
,
J.
Ferziger
, and
I.
Mitchell
, “
A hybrid particle level set method for improved interface capturing
,”
J. Comput. Phys.
183
,
83
(
2002
).
34.
R. P.
Fedkiw
,
T.
Aslam
,
B.
Merriman
et al., “
A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)
,”
J. Comput. Phys.
152
(
2
),
457
492
(
1999
).
35.
S.-C.
Chang
, “
The method of space-time conservation element and solution element—A new approach for solving the Navier-Stokes and Euler equations
,”
J. Comput. Phys.
119
,
295
324
(
1995
).
36.
G.
Wang
,
H.
Zhu
,
Q.
Sun
et al., “
An improved CE/SE scheme and its application to dilute gas-particle flows
,”
Comput. Phys. Commun.
182
(
8
),
1589
1601
(
2011
).
37.
S.-C.
Chang
,
X.-Y.
Wang
, and
C.-Y.
Chow
, “
The space-time conservation element and solution element method: A new high-resolution and genuinely multidimensional paradigm for solving conservation laws
,”
J. Comput. Phys.
156
(
1
),
89
136
(
1999
).
38.
H.
Shen
,
C.-Y.
Wen
,
M.
Parsani
et al., “
Maximum-principle-satisfying space-time conservation element and solution element scheme applied to compressible multifluids
,”
J. Comput. Phys.
330
,
668
692
(
2017
).
39.
Y.
Jiang
,
C.-Y.
Wen
, and
D.
Zhang
, “
Space-time conservation element and solution element method and its applications
,”
AIAA J.
58
(
12
),
5408
5430
(
2020
).
40.
D. X.
Yang
, “
Numerical analysis of drop coalescence in cavity flow using CESE method
,”
AIP Adv.
10
(
5
),
055116
(
2020
).
41.
G. R.
Johnson
and
W. H.
Cook
, “
Fracture characteristics of three metals subjected to various strains strain rates, temperatures and pressures
,”
Eng. Fract. Mech.
21
(
1
),
31
48
(
1985
).
42.
G. R.
Johnson
and
W. H.
Cook
, “
A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures
,” in
Presented at the Seventh International Symposium on Ballistics
,
The Hague, The Netherlands
,
April 1983
.
43.
S.
Marchi
,
D. D.
Durda
,
C. A.
Polanskey
et al., “
Hypervelocity impact experiments in iron-nickel ingots and iron meteorites: Implications for the NASA Psyche mission
,”
J. Geophys. Res.: Planets
125
,
e2019JE005927
, (
2020
).
44.
E. F.
Dudarev
,
A. B.
Markov
,
G. P.
Bakach
et al., “
Shock-wave and spalling phenomena in ultrafine-grained and coase-grained (α + β) alloy Ti-Al-V treated by a nanosecond relativistic high-current electron beam
,”
AIP Conf. Proc.
1783
,
020047
(
2016
).
45.
D.
Veysset
,
J.-H.
Lee
,
M.
Hassani
,
S. E.
Kooi
,
E. L.
Thomas
, and
K. A.
Nelson
, “
High-velocity micro-projectile impact testing
,”
Appl. Phys. Rev.
8
,
011319
(
2021
).
46.
R.
Ogawa
,
A. M.
Nakamura
,
A. I.
Suzuki
, and
S.
Hasegawa
, “
Crater shape as a possible record of the impact environment of metallic bodies: Effects of temperature, impact velocity and impactor density
,”
Icarus
362
,
114410
(
2021
).
47.
Q. L.
Zeng
,
J. W.
McCauley
, and
K. T.
Ramesh
, “
A mechanism-based model for the impact response of quartz
,”
J. Geophys. Res.: Solid Earth
126
(
3
),
e2020JB020209
, (
2020
).
48.
A.
Chadegani
,
K. A.
Iyer
,
D. S.
Mehoke
, and
R. C.
Batra
, “
Hypervelocity impact of a steel microsphere on fused silica sheets
,”
Int. J. Impact Eng.
80
,
116
132
(
2015
).
49.
A. I.
Suzuki
,
C.
Okamoto
,
K.
Kurosawa
,
T.
Kadono
,
S.
Hasegawa
, and
T.
Hirai
, “
Increase in cratering efficiency with target curvature in strength-controlled craters
,”
Icarus
301
,
1
(
2018
).
50.
L. B.
Tran
and
H. S.
Udaykumar
, “
A particle-level set-based sharp interface Cartesian grid method for impact, penetration, and void collapse
,”
J. Comput. Phys.
193
,
469
510
(
2004
).