© C. Sunday et al. 2022

1 Introduction

Two recent small body missions have involved direct interactions with asteroid surfaces: the JAXA Hayabusa2 mission (Watanabe et al. 2017) and the NASA OSIRIS-REx mission (Lauretta et al. 2017). The Hayabusa2 spacecraft arrived at asteroid (162173) Ryugu in June 2018 and deployed several modules (a lander and two mini-hoppers) to the asteroid’s surface (Van Wal et al. 2018; Ho et al. 2021). In addition, the spacecraft successfully collected a surface sample (Morota et al. 2020) and performed an impact experiment to observe the creation of an artificial crater (Arakawa et al. 2020). From December 2018 to June 2021, the OSIRIS-REx spacecraft surveyed and mapped the surface of asteroid (101955) Bennu. In October 2020, it performed a successful touch-and-go maneuver to retrieve a sample from the asteroid’s surface (Bierhaus et al. 2018). Ryugu and Bennu, similar to most asteroids, are covered with boulders and a layer of loose granular material, referred to as regolith (Murdoch et al. 2015; Hestroffer et al. 2019). Granular materials display notoriously complex behaviors on Earth, and the fundamental contact laws used to describe particle-to-particle behavior remain the same for different gravity levels. However, the electrostatic and cohesive forces that are often ignored on Earth can become nonnegligible in low-gravity (Scheeres et al. 2010; Hestroffer et al. 2019), and the bulk granular system can experience regime changes under different external triggers (Brucks et al. 2007; Murdoch et al. 2017). The nonintuitive nature of granular flow in low-gravity, paired with frequent unknowns regarding target body surface materials, makes it difficult to plan and analyze spacecraft operations similar to those listed above. It also makes it very challenging to deduce the origins of the various surface features observed on asteroids.

The ultimate objective of most small body missions is to understand the formation and evolution of these celestial objects. The Hayabusa2 and OSIRIS-REx spacecrafts are not the first to have visited small bodies, nor will they be the last to conduct in situ surface experiments. JAXA’s upcoming Martian Moons eXploration mission (MMX) will deploy a rover and will land on the surface of Phobos (Michel et al. 2022; Usui et al. 2020). The spacecraft will be equipped with a system to collect and return a sample to Earth in 2029. The NASA DART mission, planned for launch in late November 2021, will perform an impact experiment on the small moon called Dimorphos, which is part of the binary asteroid (65803) Didymos (Cheng et al. 2018). Following DART, ESA will launch the Hera mission to observe the outcomes of the impact and to track the dynamical and physical changes of the asteroid system (Michel et al. 2018). The Hera mission includes two CubeSats that may attempt to land on the surface of Dimorphos and return information on the surface response using accelerometer measurements (Goldberg et al. 2019; Ferrari et al. 2021). A basic understanding of regolith dynamics is required to design and operate the next generation of impactors, landers, and sampling systems. More critically, scientists need reasonably accurate models to help analyze observations from tool-regolith interactions. With the right data, the act of impacting or landing on a granular surface can provide valuable insight regarding a surface’s friction, cohesion, and density properties, not to mention its response to external forces. This information is crucial for developing a robust timeline of the objects’ geological history, which is linked to that of the Solar System and its surface chronology.

In this study, we introduce a framework that can be used to analyze the outcome of slow impacts and interactions with granular materials, independent of gravity level. This model can therefore be applied to regolith-covered celestial bodies of any size. Several existing studies have investigated the link between low-velocity impact dynamics and landing on small body surfaces. One popular area of research focuses on rebound behavior and bouncing dynamics in reduced-gravity environments (Colwell & Taylor 1999; Colwell 2003; Çelik et al. 2021), while another investigates the mass transfer and ejecta patterns associated with slow collisions (Schwartz et al. 2014; Brisset et al. 2018, 2020). Several works have attempted to simulate and analyze specific lander and spacecraft interactions, such as those associated with the Hayabusa2 and the OSIRIS-REx missions (Maurel et al. 2018; Thuillet et al. 2018, 2020; Ballouz et al. 2021). Here, we investigate the general behavior of a spherical projectile impacting a bed of spherical grains at low speeds. In accordance with our previous experiments, we assume that the projectile comes to a rest without rebounding (Murdoch et al. 2017, 2021; Sunday et al. 2021).

The objective of this study is to develop a theoretical collision model that links a projectile’s final penetration depth to its impact velocity, the gravity level of the system, and the properties of the surface material. In Sect. 2.1, we present existing models and discuss why these solutions cannot be used to predict impact behavior for different gravitational accelerations. In Sects. 2.2 and 2.3, we introduce a revised collision model, and we reformulate the model in terms of the dimensionless Froude number, which takes into account both impact velocity and gravity. Then, we extend the work of Sunday et al. (2021) by performing impact simulations using the soft-sphere discrete element method (SSDEM) and the open-source code Chrono (Tasora et al. 2016; Sunday et al. 2020). In Sect. 3, we use the simulation results to show how collision behavior scales with impact velocity and gravity, and we assess the accuracy of the revised model. Finally, in Sect. 4, we demonstrate how collision behavior relates to the physical and mechanical properties of a surface material, and we discuss the implications for upcoming small body missions.

2 Method

2.1 Existing collision models

With the help of diverse experimental and numerical techniques, low-velocity collisions into granular materials have been studied extensively for different surface materials and projectile shapes (Omidvar et al. 2014; Katsuragi 2016). While some of this research focuses on the refinement of empirical scaling relationships (Uehara et al. 2003; de Bruyn & Walsh 2004), a significant portion of the existing work attempts to characterize projectile behavior in terms of Newton’s second law of motion (see Omidvar et al. 2014 for more information): $$m\ddot{z}=mg-f(z)-b(z)\dot{z}-h(z){\dot{z}}^2.$$(1)

In Eq. (1), m is the mass of the projectile, z is the depth of the projectile, g is gravity, and f(z), b(z), and h(z) are the coefficients associated with the drag forces that act on the projectile as it moves through a granular medium. f(z) is a quasi-static resistance force, b(z) ż a viscous resistance force, and h(z) ż² is an inertial resistance force. f(z) can either be constant or depth-dependent depending on the model, but b(z) and h(z) are often considered to be constant and depth-independent (see Table 1).

By neglecting or modifying the form of the drag terms in Eq. (1), the general force law can be simplified to better represent impact behavior within specific ranges of collision velocities. For example, the viscous and inertial drag terms are generally neglected from the model when $v_{\text{c}}<\sqrt{2gd}/10$, where v_c is the projectile’s impact velocity and d is the average grain diameter (Albert et al. 1999). Above this threshold, however, the behavior of the projectile becomes velocity-dependent. When the inertial term is neglected, Eq. (1) resembles the Bingham force model, where $m\ddot{z}=mg-f(z)-b(z)\dot{z}$. If the viscous term is neglected instead of the inertial term, Eq. (1) takes on the form of the Poncelet model, where $m\ddot{z}=mg-f(z)-h(z){\dot{z}}^2$ (Omidvar et al. 2014; Katsuragi 2016). Table 1 provides a summary of the most frequently cited models in reference to low-velocity granular impacts (i.e., impact where v_c < 5−6 m s⁻¹).

The models proposed by de Bruyn & Walsh (2004) and Clark & Behringer (2013) seem appropriate for impacts where v_c ≈ 1–6 m s⁻¹, but they exhibit nonphysical behavior for dense materials and for lower collision velocities. More specifically, these models fail to capture the nonzero initial sinkage of a projectile when v_c = 0. Replacing the constant friction term f_o with a depth-dependent f(z) term improves a model’s ability to predict collision behavior as v_c → 0. Katsuragi & Durian (2013), Lohse et al. (2004), and Kang et al. (2018) introduce models where f(z) varies linearly with depth, while Tsimring & Volfson (2005) assume that f(z) varies with the area of the projectile in contact with the grains. Alternatively, Brzinski III et al. (2013) assert that f(z) varies by the volume of the material displaced by the projectile. Pacheco-Vázquez et al. (2011) introduce yet another form of f(z) that takes into accountthe effect of container’s size on the pressure in the system (i.e., the Janssen effect, described in Katsuragi 2016). The models listed in Table 1 succeed to various degrees when used to predict the penetration depth or the drag force for collisions where v_c < 6 m s⁻¹. With the exception of Brzinski III et al. (2013), however, these references do not address how their findings apply to different gravity levels, despite the fact that g appears as a variable in the expressions.

Several other groups have performed experiments using either fluidized granular beds or drop-tower setups in order to better understand the role of gravity within the context of select collision models. Goldman & Umbanhowar (2008) use a simple drop-tower to investigate how collision behavior scales with projectile size, surface density, and gravity while comparing their results against several of the models in Table 1. Costantino et al. (2011) use a fluidized granular bed to evaluate gravity scaling for the f(z) term in the Tsimring & Volfson (2005) model. Altshuler et al. (2014), in comparison, reference the f(z) term in the Pacheco-Vázquez et al. (2011) model while studying settling behavior using a drop-tower setup. Murdoch et al. (2021) also use a drop-tower setup (described in Sunday et al. 2016), but they consider impact, not settling behavior, and they compare their experimental measurements to the model referenced in Clark & Behringer (2013).

The above studies highlight the limitations of applying the existing collision models to impacts in low-gravity environments. Even though the scaling relationships for gravity are still unclear, the shared and most notable issue with the current models is that the definition for a “low-velocity” collision does not account for gravity. Interestingly, granular materials have been shown to fluidize more readily in low-gravity (Murdoch et al. 2017), so the model that best describes a 1 m s⁻¹ collision on Earth is not likely to be valid for a 1 m s⁻¹ collision on a small body.

In the remainder of this section, we describe a collision model that is applicable for all gravity levels. We begin by introducing a revised model that describes impacts up to 10 m s⁻¹ under terrestrial gravity. Then, we reformulate the model to highlight the role that gravity plays in the collision process, and we define two distinct impact regimes that are characterized by both collision velocity and gravity. At the end of the section, we provide an overview of the simulation setup that was used to evaluate the revised model.

2.2 Revised collision model

There is strong evidence to support the hypothesis that the quasi-static friction force in Eq. (1) varies with depth when a projectile first impacts a surface, but that the friction force eventually saturates and becomes constant if the initial impact velocity of the projectile is high enough (Brzinski III et al. 2013; Tsimring & Volfson 2005). Here, we assume that f(z) becomes constant at a depth of z₁ and that $$m\ddot{z}=\{\begin{array}{ll}mg-f_o\frac{z}{z_1}-h{\dot{z}}^2\hfill & \text{when}z\le z_1\hfill \\ mg-f_o-h{\dot{z}}^2\hfill & \text{when}z>z_1,\hfill \end{array}$$(2)

where f_o and h are constants. We then define two collision regimes based on z₁ and the final penetration depth of the projectile, z_stop. In collision regime I, z_stop ≤ z₁ and f(z) = f_oz∕z₁ for the entire duration of the collision. In collision regime II, z_stop > z₁, and f(z) is depth-dependent when z ≤ z₁ but constant when z > z₁. Figure 1 illustrates the two phases of the collision during regime II.

Expressions for the penetration depth in collision regime I, z_stop,I, and the penetration depth in collision regime II, z_stop,II, can be obtained by recasting Eq. (1) in terms of kinetic energy K and solving the linear ordinary differential equation with the correct boundary conditions for each regime. Following the method detailed in Appendix A, $$z_{\text{stop},\text{I}}=\frac{m}{2h}\text{ln}\left[\frac{h{v}_{\text{c}}^2-mg-\frac{m{f}_o}{2h{z}_1}}{f_o\frac{z_{\text{stop},\text{I}}}{z_1}-mg-\frac{m{f}_o}{2h{z}_1}}\right],$$(3)

and $$z_{\text{stop},\text{II}}=\frac{m}{2h}\text{ln}\left[\frac{h{v}_{\text{c}}^2-mg+\frac{m{f}_o}{2h{z}_1}\left(\text{exp}\left(\frac{2h{z}_1}{m}\right)-1\right)}{f_o-mg}\right].$$(4)

Equation (3) is valid when z_stop,I ≤ z₁ (collision regime I) and Eq. (4) is valid when z_stop,II > z₁ (collision regime II). If v_c = 0, then z_stop,I = z_o, or the initial sinkage of the projectile. If z_o is a known value, then Eq. (3) can be rearranged to calculate z₁ as follows: $$z_1=\frac{f_o{z}_o}{mg}\left(\frac{\text{exp}\left(\frac{2h{z}_o}{m}\right)}{\text{exp}\left(\frac{2h{z}_o}{m}\right)-1}\right)-\frac{f_o}{2hg}.$$(5)

As will be demonstrated in Sect. 3.1, the above model succeeds in reproducing two important observation from previous works: (1) the projectile has a nonzero penetration depth when v_c = 0, and (2) the quasi-static friction force is constant for higher collision velocities. The key difference between Eq. (4) and existing solutions (i.e., those found in Clark & Behringer 2013), is that here, we consider that the kinetic energy of the projectile is $\frac{1}{2}m{v}_1^2$ when f(z) becomes constant, not $\frac{1}{2}m{v}_{\text{c}}^2$. The velocity term v₁ corresponds to the velocity of the projectile at depth z₁.

Fig. 1 Force models applied to collision regime II. The quasi-static drag force acting on the projectile f(z) varies linearly when 0 < z < z1 but is constant when z1 < z < zstop,II. m is the mass of the projectile, g is gravity, z is the depth of the projectile, z1 is the depth where f(z) becomes constant, zstop,II is the final penetration depth in collision regime II, and fo, and h are constants.

2.3 Collision model based on Froude scaling

When studying low-velocity impacts, it is common to adopt the hypothesis that the f(z) term in Eq. (1) originates from hydrostatic pressure and Coulomb friction such that f(z) ∝ ρ_ggAz, where A is the cross sectional area of the projectile and ρ_g is the bulk density of the surface material (Tsimring & Volfson 2005; Katsuragi 2016). On the other hand, the inertial drag force, h(z)v², can be expressed as the momentum transfer between the projectile and the surface material, where h(z) ∝ ρ_gA (Katsuragi 2016). At first glance, f(z) has a clear dependence on gravity, while h(z)v² does not.

With the above relationships, Katsuragi (2016) shows that the general force model can be broken down into five key dimensionless Π-groups using a, v, z, g, D, ρ_p, ρ_g, and μ. The ρ_p term is the density of the projectile, and μ is a coefficient related to the internal friction angle of the granular material. One of the dimensionless groups, v² ∕gD, corresponds to the Froude number Fr. The Froude number is the ratio of the inertial to gravitational forces in a system and is a powerful dimensionless parameter for demonstrating how certain processes scale with gravity. For instance, Brucks et al. (2007) use Fr to study granular flow in a rotating drum for varied gravity levels, Housen & Holsapple (2003) use a term resembling Fr to develop crater size scaling relationships, Hilton & Tordesillas (2013) and Faug (2015) use Fr to study drag through granular materials, and Wright et al. (2020) reference Fr in the context of ricochets on asteroid surfaces.

The exact expression of Fr varies by application and is different in all of the above works. As in Wright et al. (2020), we define the impact Froude number as $$\text{Fr}=\frac{v_{\text{c}}}{\sqrt{gD}},$$(6)

where D is the projectile diameter. The expressions in Sect. 2.2 can be reformulated in terms of Fr and by substituting Eq. (6) into Eqs. (3) and (4). Since f(z) is a function of gravity, we assume that $f_o=f_o^{\prime }g$, and we rewrite the penetration depth equations for regimes I and II as $$z_{\text{stop},\text{I}}^{\prime }=\frac{m}{2h^{\prime }}\text{ln}\left[\frac{h^{\prime }{\text{Fr}}^2-m-\frac{m{f}_o^{\prime }}{2h^{\prime }{z}_1^{\prime }}}{f_o^{\prime }\frac{z_{\text{stop},\text{I}}^{\prime }}{z_1^{\prime }}-m-\frac{m{f}_o^{\prime }}{2h^{\prime }{z}_1^{\prime }}}\right]$$(7)

and $$z_{\text{stop},\text{II}}^{\prime }=\frac{m}{2h^{\prime }}\text{ln}\left[\frac{h^{\prime }{\text{Fr}}^2-m+\frac{m{f}_o^{\prime }}{2h^{\prime }{z}_1^{\prime }}\left(\text{exp}\left(\frac{2h^{\prime }{z}_1^{\prime }}{m}\right)-1\right)}{f_o^{\prime }-m}\right].$$(8)

Here, $z_{\text{stop},\text{I}}^{\prime }=z_{\text{stop},\text{I}}/D$, $z_{\text{stop},\text{II}}^{\prime }=z_{\text{stop},\text{II}}/D$, $z_1^{\prime }=z_1/D$, and h′ = hD. $z_{\text{stop},\text{I}}^{\prime }$, $z_{\text{stop},\text{II}}^{\prime }$, and $z_1^{\prime }$ are dimensionless parameters while $f_o^{\prime }$ and h′ have the units of mass. In the reformulated expressions, the explicit gravity term disappears.

The transition from collision regime I to collision regime II, denoted as Fr_t, occurs when $z_{\text{stop},\text{I}}^{\prime }=z_1^{\prime }$. When Fr > Fr_t, the penetration depth exceeds z₁ and the collision takes place in regime II. Making the appropriate substitutions and rearranging Eq. (7) gives the following expression for Fr_t: $${\text{Fr}}_t=\sqrt{\frac{1}{h^{\prime }}\left[\left(f_o^{\prime }-m-\frac{m{f}_o^{\prime }}{2h^{\prime }{z}_1^{\prime }}\right)\text{exp}\left(\frac{2h^{\prime }{z}_1^{\prime }}{m}\right)+m+\frac{m{f}_o^{\prime }}{2h^{\prime }{z}_1^{\prime }}\right]}.$$(9)

2.4 Numerical simulations and data analysis

To assess the validity of the proposed collision model, numerical simulations of a spherical projectile impacting glass beads at different collision velocities and gravity levels were conducted using the soft-sphere discrete element method (SSDEM) and the Multicore module in the open-source code Chrono (Tasora et al. 2016; Sunday et al. 2020). For all tests, a 1 kg, 10 cm diameter projectile was dropped into a cylindrical container filled with 1 ± 0.05 cm diameterbeads. The container had a nominal diameter of 45 cm and was filled to a height of 22 cm. The beads were initially mixed and settled under terrestrial gravity conditions. Then, for the low-gravity tests, the beads were permitted to relax and re-settle under the desired gravity level. The same surface sample was used for all of the simulations, but the initial position of the projectile was randomized within a radius of two particle diameters at the beginning of each test to add variability to the impact configuration. The bead properties and the simulation parameters were selected by calibrating the simulation results against terrestrial-gravity impact experiments. A detailed description of the code, the simulation setup, the boundary conditions, and the simulation parameters can be found in Sect. 2.2 of Sunday et al. (2021).

In order to quantify the projectile’s collision behavior for different gravity levels, four collision parameters were extracted from each simulated impact: (1) the projectile’s collision velocity v_c, (2) the projectile’s peak acceleration a_peak, (3) the projectile’s final penetration depth z_stop, and (4) the total duration of the collision t_stop. Figure 2 shows the raw data curves that were obtained for a typical impact simulation. In this example, the projectile hits the granular bed with a collision velocity of 1 m s⁻¹. The projectile’s acceleration a, velocity v, and penetration depth z as a function of time t are shown in Figs. 2a–c, respectively. The results for a test case where g = 9.81 m s⁻² are shown in blue, and the results for a case where g = 0.1 m s⁻² are shown in red. The definition of the peak acceleration, the collision duration, and the penetration depth parameters are indicated bythe dashed lines and the text annotations on the plot. The exact method and criteria that was used to identify a_peak, z_stop, and t_stop can be found in Sect. 2.3 of Sunday et al. (2021). Additional information regarding the filtering and treatment of the acceleration data for this particular study can be found in Appendix B. Next, in Sect. 3, we compare the collision parameters for over 150 impact simulations with collision velocities ranging from 0 to 10 m s⁻¹.

Fig. 2 Acceleration, velocity, and position profiles for example impact simulations where g = 9.81 m s−2 (blue) and g = 0.1 m s−2 (red). For both test cases, the collision velocity of the projectile is 1 m s−1. The projectile’s approximate collision velocity vc, peak acceleration apeak, penetration depth zstop, and collision duration tstop are indicated by the circles, dashed lines, and text annotations on the plot.

3 Results

3.1 Collision behavior by gravity level

In this section, we compare impact behavior for three different gravity levels: g = 9.81 and 10 m s⁻² (terrestrial gravity), g = 0.1 m s⁻², and g = 0.001 m s⁻². The 0.1 m s⁻² gravity case is lower than but comparable to the gravity found on asteroid (4) Vesta, and the 0.001 m s⁻² gravity case is on the same order of magnitude as the gravities found on bodies such as asteroid (433) Eros and the moons of Mars, Phobos and Deimos (Murdoch et al. 2015; Murchie et al. 2015). Figure 3 shows the peak acceleration, the penetration depth, and the collision duration for the three gravity cases as a function of collision velocity. The measured peak acceleration in plot (a) has a quadratic dependence, but the data is shown on a log-log scale in order to better compare the full range of values.

As illustrated in Fig. 3, the penetration depth of the projectile and the duration of the collision clearly increase when gravity level decreases. This observation is consistent with the experimental results reported in Murdoch et al. (2021). Interestingly, it appears as though the peak acceleration has the same trend for all three gravity cases. It is difficult to draw conclusions from Fig. 3a, however, since there are so few overlapping test cases. The decision to test different collision velocities for different gravity levels was intentional, because the test cases were selected to cover a Froude number range of 0–10, not a collision velocity range of 0–10 m s⁻¹. Following Eq. (6), Fr = 10 for a 10 cm diameter projectile when g = 10 m s⁻² and v_c = 10 m s⁻¹. For the same Froude number, the projectile will only have an impact velocity of 0.1 m s⁻¹ when g = 0.001 m s⁻². The reason for comparing the collision behavior in terms of Froude number becomes apparent when inspecting Fig. 4.

In Fig. 4, the peak acceleration, the penetration depth, and the collision duration are expressed as dimensionless groups and are plotted with respect to the impact Froude number. Figure 4a is shown on a log-log scale for direct comparison with Fig. 3a. When z_stop is normalized by the projectile diameter D and plotted against Fr, the penetration depth collapses onto a single curve for all gravity cases. When t_stop is normalized by $\sqrt{D/g}$ and plotted against Fr, the collision duration also collapses for all gravity levels. The normalized collision duration appears to be slightly lower for the g = 0.001 m s⁻² case as compared to the other gravity cases. However, this offset is likely because the simulations for the g = 0.001 m s⁻² case were terminated earlier than desired (see Appendix B). The only collision parameter that does not collapse by Fr is the peak acceleration measurement. Several possible explanations for this discrepancy are presented in Sect. 3.4.

The curves on Fig. 4 represent the model fit and model predictions produced by the collision model in Sect. 2.3. The constants $f_o^{\prime }$, h′, and $z_1^{\prime }$ were determined by performing a nonlinear regression and least-squares fit using Eq. (8) and the penetration depth data. The fit, illustrated by the solid curve in Fig. 4b, was obtained by considering data only within a Froude range of 3–10, because the collisions within this range are assumed to take place completely within regime II. For the given data set, we found that $z_1^{\prime }$ = 0.4, $f_o^{\prime }$ = 4.07 kg, h′ = 1.66 kg. These results indicate that f(z) becomes constant at a depth slightly less than the radius of the sphere. The significance of the $f_o^{\prime }$ and h′ terms will be discussed in Sect. 4.1. The transition from collision regime I to collision regime II is calculated using Eq. (9), and based on the given fit parameters, occurs when Fr_t ≈ 1.5.

The dotted line in Fig. 4b represents the model prediction for the penetration depth in regimes I and the beginning of regime II. In Fig. 4b, the prediction was obtained using the fit parameters and Eqs. (7) and (8). The dotted line in Fig. 4c shows the model prediction for the collision duration in both regimes, as determined using the fit parameters and numerical integration of Eq. (1). Though the model slightly under-predicts the collision duration, it captures the correct form of penetration behavior at the lowest Froude numbers and is a significant improvement of the model used for a comparable study in Murdoch et al. (2021). Unlike the model in Murdoch et al. (2021), the revised model results in a nonzero penetration depth when Fr = 0 (i.e., when v_c = 0) and shows that the collision duration increases as expected when v_c → 0 (Goldman & Umbanhowar 2008).

Fig. 3 Collision behavior for a 1 kg, 10 cm diameter spherical projectile as a function of collision velocity vc for different gravity levels, where apeak is the peak acceleration, zstop is the penetration depth, and tstop is the collision duration. The measured peak acceleration in plot (a) has a quadratic dependence, but the data is shown on a log-log scale in order to better compare the full range of values.

Fig. 4 Collision behavior for a 1 kg, 10 cm diameter spherical projectile as a function of the Froude number for different gravity levels, where apeak is the peak acceleration, zstop is the penetration depth, tstop is the collision duration, g is gravity, and D is the projectile diameter. The solid line in plot (b) represents the least-squares fit to Eq. (8) when 3 < Fr < 10. The dotted line in plot (b) represents the model prediction based on Eqs. (7) and (8), and the dotted line in plot (c) represents the model prediction obtained through numerical integration of Eq. (1). Plot (a) is shown on a log-log scale for direct comparison with Fig. 3a.

Fig. 5 Evolution of the quasi-static friction force f(z) and the inertial drag force h(z)v2 for a 1kg, 10 cm diameter spherical projectile impacting glass beads under terrestrial gravity. Fd is the magnitude of the drag force. The vertical dashed lines represent zh and zstop, where zh is the depth where f(z) = h(z)v2 and zstop is the final penetration depth of the projectile.

3.2 Drag forces by collision regime

Murdoch et al. (2017) hypothesize that the inertial drag term dominates collision behavior at decreasing collision velocities as gravity level decreases. To test and better understand this concept, we compare the magnitude of the quasi-static friction force against the inertial drag force throughout each collision and calculate the percentage of the collision where the inertial drag term dominates the interaction. Since the friction drag f(z) depends on the depth of the projectile, the inertial drag component h(z)v² is always superior to the friction component at the start of a collision. As the collision progresses, however, the friction component increases until a depth z₁, where it remains constant until the projectile comes to a rest, and the inertial drag component decreases to zero as the projectile’s velocity decreases. The two drag components intersect at some depth z_h, where z_h is dependent on f_o, h, and v_c. The ratio of z_h to z_stop gives the portion of the projectile’s total travel distance where the inertial drag term exceeds the friction term. Figure 5 provides example drag profiles for the cases where g = 9.81 m s⁻² and v_c = 1.5 m s⁻¹ (Fr = 1.5), v_c = 2.2 m s⁻¹ (Fr = 2.2), and v_c = 3.5 m s⁻¹ (Fr = 3.5). The vertical dashed lines identify z_h and z_stop with respect to f(z) and h(z)v². In Fig. 5, the profiles for the quasi-static drag force, f(z), were calculated using the relationship f(z) = f_o when z ≤ z₁ and f(z) = f_oz∕z₁ when z > z₁. From Sect. 3.1, f_o = 39.9 kg m s⁻² and z₁ = 0.04 m. The profiles for the inertial drag force, h(z)v², were calculated using the relationship h(z)v² = hv². From Sect. 3.1, h = 16.6 kg m⁻¹. The velocity of the projectile v as a function of depth was found using Eq. (A.9) when z ≤ z₁ and Eq. (A.16) when z > z₁ (see Appendix A).

From z = 0 to z = z_h, the inertial drag term is superior to the friction term. From z = z_h to z = z_stop, the friction term is superior to the inertial drag term. As Fr increases, so does z_h and z_stop. The trend for z_h∕z_stop as a function of Fr is shown in Fig. 6. Here, the vertical dashed line indicates Fr_t, or the transition from collision regime I to collision regime II. The horizontal dashed line shows where the inertial drag term exceeds the friction term for more than 50% of the collision. In regime I, the friction drag dominates the collision, while in regime II, the inertial drag is more important. For the given material and fit parameters, z_h ∕z_stop ≈ 0.5 when Fr = Fr_t (see drag profile in Fig. 5a). Between Fr = 1.5 and Fr = 3.5, the impact occurs within regime II and the inertial drag term dominates the interaction, but the inertial term intersects the friction component when f(z) = f_oz∕z₁ (see example in Fig. 5b). The change in curvature around Fr = 3.5 in Fig. 6 corresponds to the moment where the inertial drag term begins intersecting the friction term at f(z) = f_o (see example in Fig. 5c).

Through this analysis, we see that granular interactions occur within an inertial regime sooner (i.e., at lower collision velocities) on small bodies than they do on Earth. Here, the transition from collision regime I to collision regime II, or the transition from a friction dominated system to an inertial system, occurs when v_c = 1.56 m s⁻¹ on Earth, but when v_c = 0.16 m s⁻¹ at g = 0.1 m s⁻², and when v_c = 0.015 m s⁻¹ at g = 0.001 m s⁻². It is importantto note that the transition Froude number, as well as the relative importance of the quasi-static and inertial drag terms, heavily depend on the properties of the granular surface material. For example, Katsuragi & Blum (2017) study impact dynamics into low-density dust agglomerates and find that friction drag dominates the majority of the collision, even for impact velocities of 3.2 m s⁻¹. The influence of different material properties on collision dynamics will be discussed further in Sect. 4.

Fig. 6 Ratio of a collision where the inertial drag force h(z)v2 is superior to the quasi-static friction force f(z) for a 1 kg, 10 cm diameter spherical projectile impacting glass beads at different Froude numbers. zh is the depth where the f(z) = h(z)v2 and zstop is the final penetration depth of the projectile. The vertical dashed line indicates Frt, or the transition from collision regime I to regime II, and the horizontal dashed line shows where the inertial drag term exceeds the friction term for more than 50% of the collision.

Fig. 7 Snapshots of two different impact simulations when g = 10 m s−2 and vc = 10 m s−1 (Fr = 10) and z = zstop. L is the container diameter and D is the projectile diameter. In image a, the 1 kg, 10 cm diameter projectile is dropped into a 31.5 cm diameter container. In image b, the same projectile is dropped into a 45 cm diameter container. The particles are colored by their initial radial distance from the center of the container.

3.3 Boundary effects by Froude regimes

Seguin et al. (2008) study the influence of container size on penetration depth for impacts up to 3 m s⁻¹ and find that the lateral container walls affect collision behavior when the distance between the projectile and the container walls is less than one projectile diameter. If the projectile is dropped too close to a wall, its penetration depth will be less than if dropped onto an unconfined surface. Goldman & Umbanhowar (2008) report similar findings for a collision velocity range of 1–5 m s⁻¹. Based on the Froude analysis presented in the previous section, it seems logical that boundary effects should be considered as a function of both collision velocity and gravity level, not just collision velocity.

To explore the relationship between boundary effects, collision velocity, and gravity, we performed impact simulations using three different container sizes: L = 3.15D, L = 4.5D, and L = 6.0D. L and D are the container and projectile diameters, respectively. Figure 7 provides snapshots of two different impact simulations when Fr = 10 and z = z_stop. The particles are colored by their initial radial distance from the center of the container. When the projectile is dropped into the smaller container, it causes particles to lift across the entire surface of the granular bed, even near the container walls. When the projectile is dropped into the larger container, the particles at the outermost edge of the container remain more or less unaffected by the projectile’s motion. In fact, the particles within two grain diameters of the wall have an average velocity of about 2.7 cm s⁻¹ at the end of the test with the 31.5 cm diameter container, but only 0.5 cm s⁻¹ at the end of the test with the 45 cm diameter container. Both visually and analytically, it seems as though the size constraints proposed by Seguin et al. (2008) are not applicable for impacts within the upper range of the tested Froude numbers.

Figure 8 shows how the projectile’s penetration depth trends with the Froude number for different container sizes. The filled symbols correspond to the results for the smallest container size (L = 3.15D) and the open symbols correspond to the results for the largest container size (L = 4.5D). Due to the large computation time, simulations were not conducted for the smallest container size and the g = 0.001 m s⁻² case when Fr > 2. The verticaldashed line in Fig. 8 indicates the transition from collision regime I to collision regime II, and the plot inset provides a zoomed-in view of the data where Fr ≤ Fr_t. Two different “large container” sizes were tested for the g = 10 m s⁻² gravity case,L = 4.5D and L = 6.0D. The penetration depth measurements for both of these cases are represented by the same open blue symbol, because the results are nearly identical for the two container sizes.

Expressing the collision behavior in terms of the Froude number helps demonstrate how boundary effects are sensitive to both collision velocity and gravity. Similar to Seguin et al. (2008) and Goldman & Umbanhowar (2008), we see that the projectile’s penetration depth is artificially low when the container diameter is too small. In collision regime I, a minimum separation of one projectile diameter between the projectile and container walls is sufficient to avoid influence from the walls. This regime includes impacts up to 1.5 m s⁻¹ in terrestrial gravity, but only 0.16 m s⁻¹ when g = 0.1 m s⁻² and 0.016 m s⁻¹ when g = 0.001 m s⁻². In comparison, a container diameter of at least 4.5D is required to avoid influence from the walls in collision regime II. We consider that the 4.5D size is sufficiently large because the penetration depth when L = 4.5D and L = 6.0D is the same for the terrestrial gravity case. As discussed in Sect. 3.2, collision regime I is characterized by friction drag while collision regime II is dominated by inertial drag. Since particles are mobilized more easily during impacts in regime II, it follows that a larger wall-projectile separation distance is required to avoid boundary effects.

Fig. 8 Penetration depth as a function of the Froude number for a 1 kg, 10 cm diameter projectile impacting glass beads under different gravity levels. zstop is the penetration depth, vc is the collision velocity, g is gravity, D is the projectile diameter, and L is the container diameter. The filled symbols correspond to the smallest container size (L = 3.15D) and the opensymbols correspond to the largest container sizes (L = 4.5D, or when g = 10 m s−2, L = 4.5D and L = 6.0D). The vertical dashed line indicates Frt, or the transition from collision regime I to collision regime II, and the plot inset provides a zoomed-in view of the data where Fr ≤ Frt.

3.4 Peak acceleration measurements

In Sect. 3.1, we see that the penetration depth and the collision duration collapse when expressed as a function of the Froude number. However, the peak acceleration does not. To highlight the differences between the trends, Fig. 9 shows the acceleration, the velocity, and the position profiles for collisions at three different gravity levels when Fr = 6.3. The profiles are normalized by gravity and projectile diameter as shown. The inset of Fig. 9a provides a closer view of the acceleration behavior directly following the impact.

The variations in the peak acceleration measurements might occur for several reasons. The fluctuations in the acceleration profiles are due in large part to the coarse size of the grains and roughness of the surface layer (Goldman & Umbanhowar 2008). To introduce variability between simulation trials, the initial position of the projectile was randomized within a radius of two grain diameters before the start of each test. Some of the differences in the peak measurements can therefore be attributed to changes in the impact point between the projectile and the surface material. The filtering method that is used to process the data also changesthe peak measurements (see Appendix B). Lastly, it is worth noting that the bulk density of surface material varies slightly between the three gravity cases (ranging from approximately 1.36 g cm⁻³ when g = 0.001 m s⁻² to 1.46 g cm⁻³ when g = 10 m s⁻²). It is likely that the peak acceleration is affected by the surface density, but more work is required to truly understand the relationship between the peak value, the gravity level, and the surface properties (see Sects. 4.1 and 5).

The sensitivity of the acceleration profile to grain size, material properties, and processing method makes it difficult to extract meaningful information from the peak acceleration parameter. An additional challenge of using a_peak to analyze collision behavior is that a clear relationship between a_peak and v_c has yet to be determined. Murdoch et al. (2021) derive an expression for the peak acceleration as a function of collision velocity using a model where f(z) = f_o and h(z) = h. However, the solution assumes that the peak acceleration occurs right at the moment of impact. Figure 10 illustrates how the arrival of the peak changes with the Froude number, and by extension, the collision velocity and the gravity level.

At the higher end of the tested Froude range (Fr ≥ 6.3), the peak acceleration occurs within a few milliseconds of the initial impact. By comparison, the arrival of the peak is delayed with respect to the initial impact for mid-range Froude values (Fr ≈ 3–4.5). At the lowest Froude numbers (Fr ≤ 1), the peak becomes almost indistinguishable from the other fluctuations in the profile. Goldman & Umbanhowar (2008) show how the general form of the acceleration profile also changes with surface material. All existing collision models, including this one, fail to reproduce the complex details of the acceleration profile. At the same time, however, they succeed in predicting macroscopic collision properties such as the penetration depth and the collision duration. For the above reasons, the penetration depth appears to be a more reliable metric for interpreting collision behavior than the peak acceleration.

Fig. 9 Example acceleration, velocity, and position profiles for a 1 kg, 10 cm diameter projectile impacting glass beads at three different gravity levels when Fr = 6.3. The profiles are normalized by gravity g and projectile diameter D as shown. The inset of plot a highlights the fluctuations in the acceleration profiles directly following the impact.

Fig. 10 Example acceleration, velocity, and position profiles for a 1 kg, 10 cm diameter projectile impacting glass beads at three different Froude values when g = 10 m s−2. The profiles are normalized by gravity g and projectile diameter D as shown. The inset of plot a highlights how the timing of the peak acceleration measurement changes with the impact Froude number.

4 Discussion

4.1 Interpreting penetration depth measurements

An undeniable challenge for studying small body surfaces is the lack of opportunity for extensive in situ exploration. The surfaces of asteroids (25143) Itokawa, (162173) Ryugu, and (101955) Bennu have been examined using hopping rovers, impact experiments, and/or touch-and-go sampling maneuvers (Kawaguchi et al. 2008; Lauretta et al. 2017; Watanabe et al. 2017). These operations have produced incredibly valuable, but limited data related to spacecraft-regolith interactions. As a result, many in the field would like to know how to constrain the material properties of a surface using a single accelerometer or force profile.

As discussed in Sect. 2.3, the quasi-static friction and inertial drag terms in Eq. (1) are proportional to hydrostatic pressure and momentum transfer, respectively. This means that f(z) and h(z)v² have physical links to the size and mass of the projectile and the bulk density and the friction angle of the surface material. If f(z) scales with Coulomb friction and hydrostatic pressure, then f_o = μρ_ggz₁A. If h(z)v² scales with momentum transfer and is proportional to the friction coefficient μ, then h = μρ_gA (Katsuragi 2016). Empirical analysis of experimental data has shown that f(z) and h(z)v² also scale by two dimensionless factors, α and η, and that $f_o\propto {({\rho }_{\text{g}}{\rho }_{\text{p}})}^{1/2}$, where ρ_p is the density of the projectile (Katsuragi & Durian 2013). The scaling for f_o and h then becomes $$f_o=\alpha \mu {\left({\rho }_{\text{p}}{\rho }_{\text{g}}\right)}^{1/2}gA{z}_1,$$(10)

and $$h=\eta \mu {\rho }_{\text{g}}A.$$(11)

Katsuragi & Durian (2013) find that α = 8.0 and η = 2.7 for a model where f(z) = kz and h(z) = h, where k and h are constants. In this study, Eqs. (10) and (11) can be used to calculate α and η. Using the fit parameters found in Sect. 3.1 and recalling that $f_o^{\prime }=f_o/g$, h′ = hD, and $z_1^{\prime }=z_{1}D$, we find that α ≈ 19.6 and η ≈ 3.7. The revised model from Sect. 2.2 can then be reformulated as $$m\ddot{z}=\{\begin{array}{ll}mg-19.6\mu {\left({\rho }_{\text{p}}{\rho }_{\text{g}}\right)}^{1/2}gAz-3.7\mu {\rho }_{\text{g}}A{\dot{z}}^2\hfill & \text{when}z\le z_1\hfill \\ mg-19.6\mu {\left({\rho }_{\text{p}}{\rho }_{\text{g}}\right)}^{1/2}gA{z}_1-3.7\mu {\rho }_{\text{g}}A{\dot{z}}^2\hfill & \text{when}z>z_1.\hfill \end{array}$$(12)

The scaling parameters α and η were determined using the mean bulk density of all simulated material beds (ρ_g ≈ 1.4 g cm⁻³) along with the expression μ = tanϕ, where ϕ is the approximate repose angle for the 10 cm diameter glass beads (ϕ ≈ 22°) (Sunday et al. 2021).

The scaled model given by Eq. (12) is a preliminary solution because the scaling expressions for f_o and h are still being debated in literature. More work is required to verify how f_o and h scale as a function of ρ_g, ρ_p, ϕ, and D, and Eq. (12) needs to be tested for different impactor sizes and shapes as well as different surface materials. We introduce the preliminary model in this section because it allows us to demonstrate how one can infer certain material properties from a projectile’s penetration depth. Figure 11 illustrates how the penetration depth varies with the bulk density and the internal friction angle of the surface material. In this analysis, we consider impacts for three Froude values, Fr = 0, Fr = 3, and Fr = 10, and we vary either the friction angle of the material (solid black line in Fig. 11) or the bulk density of the material (dashed blue line in Fig. 11).

The case where Fr = 0 corresponds to the initial sinkage of the projectile, z_o. The sinkage varies minimally within the given range of the ϕ and ρ_g values. By contrast, the penetration depth can more than double for the same range of the ϕ and ρ_g values if Fr = 10. The higher the collision Froude number becomes, the more the surface properties influence the projectile’s behavior. Curves similar to those shownin Fig. 11 can be constructed for a given landing event and used to constrain the bulk density or the friction angle of a surface material, as long as the impactor size, gravity, impact velocity, and penetration depth are known.

Fig. 11 Penetration depth for a 1 kg, 10 cm diameter projectile into surface materials with different friction angles ϕ and bulk densities ρg. The solid black lines represent the case where ρg = 1.4 g cm−3 and ϕ is varied. The dashed blue lines represent the case where ϕ = 22° and ρg is varied.

4.2 Application to small body missions

Over the past few years, the Hayabusa2 and OSIRIS-REx spacecrafts have been visiting and sampling the near-Earth asteroids Ryugu and Bennu, respectively. In the coming decade, the Martian Moons eXploration (MMX) mission will perform a Phobos sample return and a Deimos fly-by, and the upcoming DART and Hera missions will contribute to the first asteroid deflection test and rendezvous with the binary asteroid Didymos and its secondary, Dimorphos. The bulk densities of Ryugu and Bennu are measured to be about 1.19 g cm⁻³ (Lauretta et al. 2019; Watanabe et al. 2019), while the densities of Phobos and Didymos are estimated to be higher and around 1.85 g cm⁻³ and 2.17 g cm⁻³, respectively (Murchie et al. 2015; Naidu et al. 2020). The internal friction angles for the different surface materials are either assumed or are based on modeling and simulation results. The Environment Requirement Document (ERD) produced by the MMX rover team assumes that the regolith layer on Phobos could have a friction angle ranging anywhere 30–50° (S. Tardivel, priv. comm., 2021). Similarly, the friction angles on Ryugu, Bennu, and Didymos are thought to range from 30–40° (Barnouin et al. 2019; Zhang et al. 2021). Using the estimated properties for the different materials, we can predict how far a 1 kg, 10 cm diameter projectile might sink into the surface of each body. Figure 12 shows the estimated penetration depth for the sphere when dropped onto the surfaces of Bennu, Phobos, and Didymos at different Froude numbers. The case for a glass bead surface material is also provided for reference.

The shaded regions in Fig. 12 highlight the transition from collision regime I to collision regime II. In collision regime I (the gray shaded region on the plot), the penetration depth changes only slightly for the different material cases. The difference inthe penetration depth becomes much more pronounced for collisions in regime II (the blue shaded region on the plot). As expected, the highest penetration depth is associated with the lowest friction material (glass beads), and the lowest penetration depth occurs for impacts into the material with the highest bulk density and highest internal friction angle (Didymos and Phobos). It should be noted that the presented collision model has only been tested for collision regime I and when Fr ≤ 10 in collision regime II (the dark blue shaded region on the plot). Impact cases where Fr > 10 have not yet been tested (the light blue shaded region on the plot), but the prediction is thought to provide a reasonable or a lower-boundestimate for the penetration depth in this range. As the impact Froude number increases, the collision behavior might enter a third regime where the presented model is no longer valid (see Sect. 5 for additional discussion).

Table 2 provides an example of the impact Froude numbers that are associated with various in situ operations on past and upcoming small body missions. The Froude number is calculated using the estimated collision velocity v_c and the effective projectile diameter D′, where $D^{\prime }=2{(A/\pi )}^{1/2}$ (Kang et al. 2018). The MASCOT, MMX, and Hera landers are cubic in shape (Scholten et al. 2019; Sedlmayr et al. 2020; Goldberg et al. 2019), and for the purpose of this study, are assumed to land on the largest face of the cube. The OSIRIS-REx and MMX samplers are taken to be closed cylinders, even though they are hollow in the center (Ballouz 2017; Usui et al. 2020).

The surface operations that occurred during the Hayabusa2 and OSIRIS-REx missions likely took place within the untested region of collision regime II (where Fr > 10). The future Hera CubeSat landing will also occur in the untested region, but the MMX rover will have an impact Froude number of approximately 10, and the MMX sampling operations could take place completely within collision regime I. The impact Froude numbers for these events give an indication of how much resistance the different systems will encounter during their respective surface interactions. The outcome of the actual operations, however, are still highly dependent on the given surface terrain. Ballouz et al. (2021) simulate the touchdown of the OSIRIS-REx sampling mechanism on Bennu and revise the model by Katsuragi & Durian (2013) for Bennu-level gravity. The authors investigate the influence of friction, packing fraction and cohesion on the sampling operation and observe that certain types of surface materials result in a rebound of the lander. The OSIRIS-REx TAGSAM penetrated deep into the surface of Bennu (Lauretta & Team 2021). On the other hand, the MASCOT lander appears to have hit a boulder when it touched down on the surface of Ryugu and consequently rebounded (Scholten et al. 2019). The model described in this work is appropriate for granular surfaces that can be treated as a near continuum. Additional testing is required to understand the extent to which the model applies to cohesive surfaces or surfaces with a large dispersion of grain sizes.

Fig. 12 Predicted penetration depth for a 1 kg, 10 cm diameter spherical projectile impacting different granular materials, including a bed of glass beads and surface materials similar to what might be found on Phobos, asteroid (101955) Bennu, and asteroid (65803) Didymos. zstop is the penetration depth, D is the projectile diameter, vc is the collision velocity, g is the gravity level, ρg is the bulk density of the material, and ϕ is the internal friction angle of the material. The shaded regions on the plot show the transition from collision regime I (gray) to collision regime II (blue). The presented model has been tested against numerical and laboratory experiments in this study for collisions in regime II where Fr ≤ 10 (dark blue shaded region), but not for collisions where Fr > 10 (light blue shaded region).

5 Conclusions

In this study, we introduce a revised collision model that can be used to predict the final penetration depth of a projectile during slow granular impacts. The revised model resembles the Poncelet model because it includes both quasi-static and inertial drag components. The inertial drag force is quadratic with velocity, and we assume that the drag coefficient isconstant throughout the entire collision. Then, we assume that the quasi-static friction component is depth-dependent during the initial part of a collision, but that the friction force becomes constant after some threshold depth z₁. If the projectile comes to a rest before reaching z₁, then the impact occurs within collision regime I. If the projectile comes to a rest after passing through z₁, then the impact occurs within collision regime II. The penetration depth in regimes I and II can be calculated as a function of impact velocity using Eqs. (3) and (4).

In order to highlight the influence that gravity has on collision behavior, we reformulate the expressions for the penetration depth in regimes I and II in terms of the dimensionless Froude number Fr, where $\text{Fr}=v/{(gD)}^{1/2}$. The new solutions, given by Eqs. (7) and (8), include three parameters that are dependent on the properties of the projectile and the surface material: f_o, h, and z₁. We treat f_o, h, and z₁ as fit parameters, and we determine their values using simulation data for impacts under three different gravity levels and a Froude number ranging from 3 to 10. Then, we use the proposed model to predict the penetration depth for impacts where Fr ≤ 3 and the collision duration for impacts where 0 ≤Fr ≤ 10. By applying this method, we show that the revised model predicts collision behavior well for impacts where Fr < 10. In addition, we see that the penetration depth and collision duration measurements collapse onto a single curve when expressed as a function of the Froude number.

With the proper scaling, the presented framework can be used to analyze single landing or impact events from small body missions. It should be noted that it is difficult to use the peak acceleration measurement to assess collision behavior, because the relationship between the peak measurement, the surface material, and the collision velocity is still poorly understood. However, accelerometer or force sensor data can be used to obtain estimates for penetration depth and collision duration, which can in turn be used to determine certain surface properties, such as the internal friction angle or the bulk density of the material. More work is needed to verify the scaling relationships for the friction and inertial drag coefficients, and the expressions for f_o and h should be developed as a function of cohesion, grain size, grain size distribution, grain shape, internal friction angle, surface density, and projectile shape.

The presented model has not yet been tested for lower gravity levels or higher Froude numbers, because these test cases require larger containers, and consequently, higher computation times. However, collision behavior when g < 0.001 m s⁻² is expected to follow the revised collision model and the Froude scaling presented in this work. For higher Froude number impacts, the projectile will likely enter a third collision regime where additional forms of energy dissipation must be taken into account (see Omidvar et al. 2014 for an overview of the different types of energy loss that occur for quasi-static through supersonic impacts). An objective of future work will be to assess the validity of the presented model for impacts where Fr > 10 and to identify the transition point between the second and third collision regimes.

Acknowledgements

C.S. is jointly funded by the Centre National d’Études Spatiales (CNES) and the Institut Supérieur de l’Aéronautique et de l’Espace (ISAE) under a PhD research grant. All authors acknowledge support from CNES. P.M. acknowledges funding from Academies of Excellence: Complex systems and space, environment, risk, and resilience, part of the IDEX JEDI of the Université Côte d’Azur, and Y.Z. acknowledges funding from the Doeblin Federation and the “Individual grants for young researchers” program, also part of the IDEX JEDI of the Université Côte d’Azur. P.M. and Y.Z. acknowledge funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 870377 (project NEO-MAPP). Simulations with Chrono were performed using HPC resources from CALMIP under grant allocation 2020-P20020.

Appendix A Derivation of penetration depth expressions

We begin by considering the general force law $$m\ddot{z}=mg-f(z)-h(z){\dot{z}}^2,$$(A.1)

where m is the mass of the projectile, g is gravity, f(z) is the quasi-static friction force, h(z) is the inertial drag term, and z is the depth of the projectile. As demonstrated in Clark & Behringer (2013), Eq. (A.1) can be recast in terms of kinetic energy K such that $$\frac{dK}{dz}=mg-f(z)-\frac{2}{m}h(z)K,$$(A.2)

where $K=\frac{1}{2}m{v}^2$ and v is the velocity of the projectile. Eq. (A.2) is a linear ordinary differential equation with the following solution: $$K(z)=K_p(z)\left(K_o+\Phi (z)\right).$$(A.3)

In Eq. (A.3), K_o is the initialkinetic energy of the projectile, $$K_p(z)=exp\left(-{\displaystyle \int \frac{2}{m}}h(z^{\prime })d{z}^{\prime }\right),$$(A.4)

and $$\Phi (z)={\displaystyle \int \frac{mg-f(z^{\prime })}{K_p(z^{\prime })}}d{z}^{\prime }.$$(A.5)

The projectile impacts the surface with a collision velocity of v_c and comes to a rest at z_stop. We assume that the inertial drag term is constant and that h(z) = h throughout the entire collision. The quasi-static friction force, however, takes on one of two forms depending on the projectile’s depth. At first, f(z) increases linearly with depth such that f(z) = kz. Then, f(z) reaches a maximum value f_o and remainsconstant until the end of the collision. The transition between the depth-dependent and the constant friction states occurs at some depth z₁, and this transition separates what we refer to as collision regimes I and II.

If the collision velocity is small enough such that z_stop ≤ z₁, then the collision takes place completely within collision regime I. In this regime, f(z) = kz, where k = f_o∕z₁. f_o is a constant that depends the properties of both the surface material and the projectile. By evaluating Eqs. (A.4) and (A.5) from 0 to z with h(z) = h and f(z) = f_o∕z₁, K_p (z) and Φ(z) can be expressed as $$K_p(z)=exp\left(-\frac{2h}{m}z\right)$$(A.6)

and $$\Phi (z)=\frac{m}{2h}\left[\left(mg-f_o\frac{z}{z_1}+\frac{m{f}_o}{2h{z}_1}\right)exp\left(\frac{2h}{m}z\right)-mg-\frac{m{f}_o}{2h{z}_1}\right].$$(A.7)

Substituting Eqs. (A.6) and (A.7) into Eq. (A.3) yields $$\begin{array}{l}K(z)=\frac{m}{2h}[\left(\frac{2h}{m}{K}_o-mg-\frac{m{f}_o}{2h{z}_1}\right)exp\left(-\frac{2h}{m}z\right)+\hfill \\ mg-f_o\frac{z}{z_1}+\frac{m{f}_o}{2h{z}_1}].\hfill \end{array}$$(A.8)

When $K_o=\frac{1}{2}m{v}_c^2$, $$\begin{array}{l}v(z)=(\frac{1}{h}[\left(h{v}_c^2-mg-\frac{m{f}_o}{2h{z}_1}\right)exp\left(-\frac{2h}{m}z\right)+\hfill \\ {mg-f_o+\frac{m{f}_o}{2h{z}_1}])}^{1/2}.\hfill \end{array}$$(A.9)

The above expression for velocity as a function of depth is valid when z ≤ z₁. Eq. (A.8) can also be simplified using the correct boundary conditions to obtain an expression for z_stop,I, or the penetration depth in collision regime I. In regime I, $K_o=\frac{1}{2}m{v}_c^2$ and K(z_stop,I) = 0, so $$z_{\text{stop},\text{I}}=\frac{m}{2h}\text{ln}\left[\frac{h{v}_c^2-mg-\frac{m{f}_o}{2h{z}_1}}{f_o\frac{z_{\text{stop},\text{I}}}{z_1}-mg-\frac{m{f}_o}{2h{z}_1}}\right].$$(A.10)

If the collision velocity is large enough such that z_stop > z₁, then the collision occurs within regime II. In collision regime II, $$f(z)=\{\begin{array}{ll}kz\hfill & \text{when}0\le z\le z_1\hfill \\ f_o\hfill & \text{when}{z}_1<z\le z_{\text{stop},\text{II}},\hfill \end{array}$$(A.11)

where z_stop,II is the penetration depth in regime II and k = f_o∕z₁. Since f(z) at the beginning of regime II has the same form as f(z) in regime I, the velocity of the projectile at depth z₁, or v₁ can be determined using Eq. (A.8). During this phase of the collision, $K_o=\frac{1}{2}m{v}_c^2$ and $K(z_1)=\frac{1}{2}m{v}_1^2$, so $$v_1=\sqrt{\frac{1}{h}\left[\left(h{v}_c^2-mg-\frac{m{f}_o}{2h{z}_1}\right)exp\left(-\frac{2h{z}_1}{m}\right)+mg-f_o+\frac{m{f}_o}{2h{z}_1}\right]}.$$(A.12)

The expression for K(z) during the second phase of the collision must be obtained by evaluating Eqs. (A.4) and (A.5) from z₁ to z. When h(z) = h and f(z) = f_o, $$K_p(z)=exp\left(-\frac{2h}{m}\left(z-z_1\right)\right),$$(A.13) $$\Phi (z)=\frac{m}{2h}\left[\left(mg-f_o\right)exp\left(\frac{2h}{m}\left(z-z_1\right)\right)+f_o-mg\right],$$(A.14)

and $$K(z)=\frac{m}{2h}\left[\left(\frac{2h}{m}{K}_o+f_o-mg\right)exp\left(-\frac{2h}{m}\left(z-z_1\right)\right)+mg-f_o\right].$$(A.15)

Then, $K_o=\frac{1}{2}m{v}_1^2$, and $$\begin{array}{l}v(z)=(\frac{1}{h}[\left(h{v}_c^2-mg-\frac{m{f}_o}{2h{z}_1}\right)exp\left(-\frac{2h}{m}z\right)+\hfill \\ {mg-f_o+\frac{m{f}_o}{2h{z}_1}exp\left(-\frac{2h}{m}\left(z-z_1\right)\right)])}^{1/2}.\hfill \end{array}$$(A.16)

The above expression for velocity as a function of depth is valid when z > z₁. The penetration depth in regime II, z_stop,II, is determined by applying the correct boundary conditions to Eq. (A.15). During the second phase of the collision, $K_o=\frac{1}{2}m{v}_1^2$, where v₁ is given by Eq. (A.12), and K(z_stop,II) = 0, so $$z_{\text{stop},\text{II}}=\frac{m}{2h}\text{ln}\left[\frac{h{v}_c^2-mg+\frac{m{f}_o}{2h{z}_1}\left(exp\left(\frac{2h{z}_1}{m}\right)-1\right)}{f_o-mg}\right].$$(A.17)

Appendix B Post-processing of acceleration data

In this study, the projectile’s acceleration information is provided as a direct data output from the impact simulations. The data is reported at a high sampling frequency of 10 kHz. In Sunday et al. (2021), we calibrated the simulation parameters using terrestrial gravity impact experiments. In order to compare the peak acceleration measurements from the impacts experiments and simulations, we filtered the acceleration data using using a second order Butterworth filter with a cut-off frequency of 500 Hz (see Sunday et al. (2021) for more information). This filtering method resulted in a reasonable match between the different measurements because the impact simulation and experiment were conducted under the same gravity level. When comparing impact behavior under different gravity levels, however, the same processing technique cannot be used because the collision time scale varies with gravity. Rather than filtering the data in the time domain t as was donein Sunday et al. (2021), we filtered the data in this study using the normalized time domain t^*, where $t^{*}=t{(g/D)}^{1/2}$. Fig. B.1 illustrates why it is important to consider the collision time scale when comparing impacts under different gravity levels. The blue, red, and green lines on the plot represent the acceleration profiles for impacts where Fr = 6.3 and g = 10, 0.1 and 0.001 m s⁻² respectively.Fig. B.1 (a) shows the acceleration profiles when the data is filtered in the time domain t, and Fig. B.1 (b) shows the acceleration profiles when the data is filtered in the normalized time domain t^*.

As seen in Fig. B.1, the peak acceleration measurements change drastically with the processing method. Based on the frequency of the fluctuations in the acceleration data, it seems more appropriate to filter the data in the normalized time domain. Though not shown, it should be noted that the filtering method does not influence the velocity or depth profiles. However, the collision duration, just like the peak acceleration, should be determined in the normalized domain rather than the standard domain. In Sunday et al. (2021), the end of the collision was identified as the moment when the projectile’s velocity fell below 1 mm s⁻¹. In this study, the end of the collision is identified as the moment when the normalized velocity v^* falls below 1x10⁻², where $v^{*}=v{(gD)}^{1/2}$. If v^* = 1x10⁻², then v = 10, 1, and 0.1 mm s⁻¹ when g = 10, 0.1 and 0.001 m s⁻² respectively. The end criteria for the lowest gravity level is difficult to achieve within a reasonable computation time, so the impact simulations with the g = 0.001 m s⁻² gravity case often terminated prematurely, when v ≈ 0.05 mm s⁻¹.

Fig. B.1 Example acceleration profiles for a 1 kg, 10 cm diameter projectile impacting glass beads at three different gravity levels when Fr = 6.3. The acceleration data is normalized by gravity g, and the timeparameter t is normalized by g and the projectile diameter D. Plot (a) gives the acceleration profiles when the data is filtered in the time domain t, and plot (b) gives the profiles when the data is filtered in the normalized time domain t* where $t^{*}=t{(g/D)}^{1/2}$.

Appendix C List of symbols

Table C.1 lists of all of the symbols that are referenced in the main body of this text.

References

All Tables

All Figures

Fig. 1 Force models applied to collision regime II. The quasi-static drag force acting on the projectile f(z) varies linearly when 0 < z < z1 but is constant when z1 < z < zstop,II. m is the mass of the projectile, g is gravity, z is the depth of the projectile, z1 is the depth where f(z) becomes constant, zstop,II is the final penetration depth in collision regime II, and fo, and h are constants.

In the text

Fig. 2 Acceleration, velocity, and position profiles for example impact simulations where g = 9.81 m s−2 (blue) and g = 0.1 m s−2 (red). For both test cases, the collision velocity of the projectile is 1 m s−1. The projectile’s approximate collision velocity vc, peak acceleration apeak, penetration depth zstop, and collision duration tstop are indicated by the circles, dashed lines, and text annotations on the plot.

In the text

Fig. 3 Collision behavior for a 1 kg, 10 cm diameter spherical projectile as a function of collision velocity vc for different gravity levels, where apeak is the peak acceleration, zstop is the penetration depth, and tstop is the collision duration. The measured peak acceleration in plot (a) has a quadratic dependence, but the data is shown on a log-log scale in order to better compare the full range of values.

In the text

Fig. 4 Collision behavior for a 1 kg, 10 cm diameter spherical projectile as a function of the Froude number for different gravity levels, where apeak is the peak acceleration, zstop is the penetration depth, tstop is the collision duration, g is gravity, and D is the projectile diameter. The solid line in plot (b) represents the least-squares fit to Eq. (8) when 3 < Fr < 10. The dotted line in plot (b) represents the model prediction based on Eqs. (7) and (8), and the dotted line in plot (c) represents the model prediction obtained through numerical integration of Eq. (1). Plot (a) is shown on a log-log scale for direct comparison with Fig. 3a.

In the text

	Fig. 5 Evolution of the quasi-static friction force f(z) and the inertial drag force h(z)v2 for a 1kg, 10 cm diameter spherical projectile impacting glass beads under terrestrial gravity. Fd is the magnitude of the drag force. The vertical dashed lines represent zh and zstop, where zh is the depth where f(z) = h(z)v2 and zstop is the final penetration depth of the projectile.
In the text

Fig. 6 Ratio of a collision where the inertial drag force h(z)v2 is superior to the quasi-static friction force f(z) for a 1 kg, 10 cm diameter spherical projectile impacting glass beads at different Froude numbers. zh is the depth where the f(z) = h(z)v2 and zstop is the final penetration depth of the projectile. The vertical dashed line indicates Frt, or the transition from collision regime I to regime II, and the horizontal dashed line shows where the inertial drag term exceeds the friction term for more than 50% of the collision.

In the text

Fig. 7 Snapshots of two different impact simulations when g = 10 m s−2 and vc = 10 m s−1 (Fr = 10) and z = zstop. L is the container diameter and D is the projectile diameter. In image a, the 1 kg, 10 cm diameter projectile is dropped into a 31.5 cm diameter container. In image b, the same projectile is dropped into a 45 cm diameter container. The particles are colored by their initial radial distance from the center of the container.

In the text

Fig. 8 Penetration depth as a function of the Froude number for a 1 kg, 10 cm diameter projectile impacting glass beads under different gravity levels. zstop is the penetration depth, vc is the collision velocity, g is gravity, D is the projectile diameter, and L is the container diameter. The filled symbols correspond to the smallest container size (L = 3.15D) and the opensymbols correspond to the largest container sizes (L = 4.5D, or when g = 10 m s−2, L = 4.5D and L = 6.0D). The vertical dashed line indicates Frt, or the transition from collision regime I to collision regime II, and the plot inset provides a zoomed-in view of the data where Fr ≤ Frt.

In the text

	Fig. 9 Example acceleration, velocity, and position profiles for a 1 kg, 10 cm diameter projectile impacting glass beads at three different gravity levels when Fr = 6.3. The profiles are normalized by gravity g and projectile diameter D as shown. The inset of plot a highlights the fluctuations in the acceleration profiles directly following the impact.
In the text

	Fig. 10 Example acceleration, velocity, and position profiles for a 1 kg, 10 cm diameter projectile impacting glass beads at three different Froude values when g = 10 m s−2. The profiles are normalized by gravity g and projectile diameter D as shown. The inset of plot a highlights how the timing of the peak acceleration measurement changes with the impact Froude number.
In the text

	Fig. 11 Penetration depth for a 1 kg, 10 cm diameter projectile into surface materials with different friction angles ϕ and bulk densities ρg. The solid black lines represent the case where ρg = 1.4 g cm−3 and ϕ is varied. The dashed blue lines represent the case where ϕ = 22° and ρg is varied.
In the text

Fig. 12 Predicted penetration depth for a 1 kg, 10 cm diameter spherical projectile impacting different granular materials, including a bed of glass beads and surface materials similar to what might be found on Phobos, asteroid (101955) Bennu, and asteroid (65803) Didymos. zstop is the penetration depth, D is the projectile diameter, vc is the collision velocity, g is the gravity level, ρg is the bulk density of the material, and ϕ is the internal friction angle of the material. The shaded regions on the plot show the transition from collision regime I (gray) to collision regime II (blue). The presented model has been tested against numerical and laboratory experiments in this study for collisions in regime II where Fr ≤ 10 (dark blue shaded region), but not for collisions where Fr > 10 (light blue shaded region).

In the text

Fig. B.1 Example acceleration profiles for a 1 kg, 10 cm diameter projectile impacting glass beads at three different gravity levels when Fr = 6.3. The acceleration data is normalized by gravity g, and the timeparameter t is normalized by g and the projectile diameter D. Plot (a) gives the acceleration profiles when the data is filtered in the time domain t, and plot (b) gives the profiles when the data is filtered in the normalized time domain t* where $t^{*}=t{(g/D)}^{1/2}$.

In the text