[orcid=0000-0002-9724-8744]

[1]

1]organization=Department of Mechanical Engineering, University of Rochester,
street=235 Hopeman Building, P.O. Box 270132,, city=Rochester, state=NY,
postcode=14627, country=USA

2]organization=Colorado Center for Astrodynamics Research, University of
Colorado Boulder, street=3775 Discovery Dr.,, city=Boulder, state=CO,
postcode=80303, country=USA

3]organization=Department of Physics and Astronomy, University of Rochester,
street=206 Bausch and Lomb Hall, P.O. Box 270171,, city=Rochester, state=NY,
postcode=14627, country=USA

# Sub-surface granular dynamics in the context of oblique, low-velocity
impacts into angular granular media

Peter M. Miklavčič<EMAIL_ADDRESS>[  [  [  Paul Sánchez
<EMAIL_ADDRESS>Esteban Wright<EMAIL_ADDRESS>Alice
C. Quillen<EMAIL_ADDRESS>Hesam Askari<EMAIL_ADDRESS>
###### Abstract

Oblique, low-velocity impacts onto non-terrestrial terrain are regular
occurrences during space exploration missions. These are not only a necessary
component of landing and sampling maneuvers, but can also be used as impact
experiments to reveal characteristics of the interacting surfaces. We conduct
two-dimensional discrete simulations to model such impacts into a bed of
triangular grains. Finite element method provides the basis for simulation,
enabling the angular grain geometry. Our findings re-create the three classes
of impact behavior previously noted from experiments: full-stop, rollout, and
ricochet (Wright, Quillen, South, Nelson, Sánchez, Siu, Askari, Nakajima and
Schwartz, 2020). An application of Set Voronoi tessellation assesses packing
fraction at a high resolution, revealing how grains shift relative to each
other during an impact event. We also assess how packing fraction at the point
of impact influences different impact behavior types. Calculation of Von Mises
strain distributions then reveal how grains shift relative to the overall
system, leading to the notion of the ‘skin zone’. Intuition would suggest that
the region of perturbed grains would grow deeper with higher velocity impacts,
results instead show that increasing velocity may evoke a change in the
grains’ dissipative response that dispatches energy predominantly laterally
from the impact site instead of deeper into the bed. Finally, we consider how
sub-surface response could link with impactor dynamics to deepen our
understanding of oblique, low-velocity impact events, one day helping to
improve mission outcomes.

###### keywords:

collisional physics cratering impact processes regoliths

Finite Element Methods support modeling of impacts into discrete triangular
grains

Simulations reproduce experimental ricochet, roll-out, and full-stop behaviors
of projectile

Visualization of local grain packing reveals breadth of granular disturbance
from impact

Strain distributions reveal zones of perturbed media around the crater

Higher velocity impacts spur increasingly lateral granular dissipative
mechanisms

##  1 Introduction

In recent decades, navigation of non-terrestrial terrain has presented a new
front to the field of granular mechanics. Such environments differ from Earth
because of different levels of gravity, low atmospheric pressure, and
materials that differ from those commonly found on Earth. Asteroids are one
such environment of interest. The population of Near-Earth-Asteroids within
reach of spacecraft are theorized to be composed entirely of rubble (Burns,
1975; Pravec and Harris, 2000; Walsh, 2018). These asteroids form over time by
accretion of fragments of rock held together by self-gravity, a process that
has spurred simulation research of its own (Leinhardt et al., 2000; Bagatin et
al., 2001; Michel et al., 2002; Richardson et al., 2005; Sánchez and Scheeres,
2011, 2014; Sánchez, 2015; Ferrari et al., 2019; Ferrari and Tanga, 2020;
Michel et al., 2020; Sánchez et al., 2021). Missions to these asteroids have
returned images and even samples of the surface asteroid material, offering
insight into the harsh environments that spacecraft must inevitably endure
(Yano et al., 2006; Tsuchiyama et al., 2011; Nakamura et al., 2011; Watanabe
et al., 2017).

Low velocity impacts on non-terrestrial surfaces have been subject of limited
interest in the literature. Nonetheless, future exploration missions - to
moons, planets, or asteroids - will likely encounter landings upon a rubble or
regolith surface, making study of such impacts important for mission
preparedness. The MASCOT lander touching down on asteroid Ryugu underlines
this importance, as it has become the subject of simulations studying the
impact at the landing site (Thuillet, Florian et al., 2018, 2021). Other
literature studies the populations of craters observed on a planet’s or
asteroid’s surface (Gou et al., 2018) focusing on such metrics as the depth-
to-diameter ratio (Noguchi et al., 2021). Morphometry of craters have been
connected through various scaling research (Goldman and Umbanhowar, 2008;
Dowling and Dowling, 2013; Çelik et al., 2022). The space shuttle and
parabolic test flights have also supported experimental campaigns for studying
low-velocity normal impacts into regolith simulant (Brisset et al., 2018).
Oblique impacts have been less well studied experimentally, though recently
Wright et al. studied impacts of a marble into sand and gravel as a function
of impact angle and velocity, classifying resultant behavior into categories
dependent upon whether the projectile ricocheted off the media or not (Wright
et al., 2020, 2021). This body of work provides valuable information about the
surface and overall projectile behavior in low velocity impact conditions, but
the dynamics of sub-surface grains and how these correspond to the impact
dynamics are, as yet, rarely studied.

Discrete element simulations offer great potential for studying sub-surface
dynamics, though such methods are constrained by how accurately and
efficiently each grain and its local interactions are modeled. A standard
approach is to assume a spherical grain geometry which enables streamlined
algorithms for contact detection and response (Cundall and Strack, 1979).
Nonetheless, real grain geometries are never perfectly round. This could not
be more true for regolith, where grains are rougher than Earthen quartz sands
due to an absence of erosive forces such as water and atmosphere. Regolith
also behaves differently than quartz sands. Under cyclic loading conditions,
regolith has been shown to continually deform with subsequent loadings,
contrasting to quartz sands that exhausts all plastic deformation within the
first cycle (Sandeep et al., 2019). Given that quartz sand is often modeled
using simple discrete spheres, it stands to reason that a model well-suited
for simulating regolith or any irregular, coarse grain could be improved by
adopting a different grain geometry. A growing segment of literature has
exemplified this. The study of railroad ballast - the coarse gravel piled
beneath railroad ties - adopts polyhedral grains in discrete simulation (Lobo-
Guerrero and Vallejo, 2006; McDowell and Li, 2016; Suhr and Six, 2017; Li and
McDowell, 2018). Recent models of rubble-pile asteroids have adapted to
simulate polyhedral grains instead of standard spherical grains (Ferrari et
al., 2019; Ferrari and Tanga, 2020; Sánchez et al., 2021). Also, poly-
ellipsoidal grains are shown to better replicate the experimental results of a
rover wheel digging in regolith simulant compared to spherical grains (Knuth
et al., 2012).

Driven by the trend towards modeling angular grains, the sparsity of relevant
impact literature, and the potential wealth of knowledge in observing sub-
surface granular dynamics, we use Finite Element Method (FEM) - via commercial
software Abaqus - to discretely simulate a two-dimensional angular grain
system impacted by a disk-like projectile. We call this combination of finite
element and discrete element methodologies the Discrete Finite Element Method
(DFEM). Specifically, we consider media composed of mono-disperse triangular
grains, a simplistic angular geometry without complexities that may introduce
numerical instabilities or make models computationally inefficient. The space-
filling nature of the triangular grains also supports our later study of grain
packing fraction, since a much wider range of packing density is observed with
triangles instead of circles. The two-dimensional basis streamlines data
extraction and simplifies calculation procedures for obtaining sub-surface
packing and strain metrics. These simulations offer access to bulk data for
each element of the virtual system such as coordinate data or nodal fields.
The packing of grains within a granular system shows how they re-arrange in
response to a dynamic event (Arévalo et al., 2010; Clark et al., 2012; Jia et
al., 2012; Kondic et al., 2012; Xu et al., 2014; Zhao et al., 2017). Such
packings are observable experimentally with X-ray tomography (Fu et al., 2006;
Moreno-Atanasio et al., 2010; Weis and Schröter, 2017; Reimann et al., 2017),
though simulations allow the same perspective but with reduced hassle. We will
show this using an adapted form of Voronoi tessellation (Voronoi, 1908). Sub-
surface strain distribution caused by the surface impacts are also studied
using methods applied to other granular systems in the literature (Chen et
al., 2011; Guo and Zhao, 2016; Luo et al., 2017). In summation, our methods
enable a deeper level of system analysis that is unattainable with
experimental methods.

The application of DFEM - FEM for discrete simulations - is present in the
literature, albeit sparse. In the simulation of a two-dimensional shear cell,
Kabir et al. (2008) used FEM to model discrete triangular, circular, and
square grains with multiple elements per grain, visualizing force chains
developed within the cell. Some angular DEM simulations have also used FEM as
an intermediary between simulation steps to facilitate breakage of grains
during modeling of rockfill (Bagherzadeh Kh. et al., 2011; Raisianzadeh et
al., 2018, 2019; Seyyedan et al., 2021; Nitka and Tejchman, 2020). The field
of multi-particle finite element method (MPFEM) simulates powders discretely
with individual meshes to observe powder compaction under load (Güner et al.,
2015; Zhang et al., 2015).

In this paper, we present a computationally efficient method of modeling
angular grains in a DEM scheme based on FEM formulation called Discrete Finite
Element Method (DFEM) that allows the study of surface impacts and reveals the
role of sub-surface properties on the outcome of impact events. Particularly,
we present how local packing fraction can result in variations in impact
characteristics and the implications of impact events on crater shape and
depth as well as the affected area around the impact site. Section 2 will
review the simulation design, procedures for conducting impact tests, and the
ranges of impact conditions evaluated. Section 3 will review modeling results.
Finally, Section 4 will summarize findings and re-iterate their broader
implications. Appendix I showcases additional simulation work that verifies
the performance of the Abaqus platform against LAMMPS DEM software (Thompson
et al., 2022). Appendix II demonstrates the potential of nodal outputs from
discrete models for viewing stress distributions through granular media.

##  2 Simulation methods and setup

Monodisperse two-dimensional discrete triangles are the base grain geometry
modeled. This geometry represents a step towards the space-filling and rough
character of naturally-occurring regolith grains while prioritizing model
efficiency with a simple geometry. The finite-elements platform means that
only one element is required per grain when formulated in Abaqus, minimizing
the risk of improper deformation modes. The simplicity will also assist post-
processing algorithms later on.

The simulation is comprised of a circular impactor and the triangular grains.
The impactor has a diameter of 16.15mm - the same as the marble used by Wright
et al. (2020) \- and the edge length of the right-triangle grain is 3mm,
giving an impactor-to-grain size ratio of 5.38. Both entities are modeled as
elastic bodies with the material properties shown in Table 1. In this work,
the system is under Earth-like gravity to improve simulation efficiency
(grains settle faster, reducing the overall number of increments to solve) and
enable comparison to experimental literature for qualitative verification of
model performance. However, this simulation setup accommodates adjustments to
gravitational acceleration to study impacts in low gravity which will be
explored in future.

To set up the system for initial grain settling, a lattice grid is positioned
above the rigid bin that acts as a container for granular material. A
triangular grain is created and randomly positioned and oriented within each
lattice cell. An external MATLAB code generates the vertex coordinates for all
grains and writes those coordinates to a text file to be referenced by the
Abaqus input file which will control the bed-settling stage of the simulation
and results in a randomized initial configuration. The initial patterning of
grains and their randomized orientation is shown in the left pane of Figure 1.
Each grain was meshed as one single three-node linear plane strain element
(type CPE3 in Abaqus). The impactor was discretized into thirty-two, four-node
bi-linear plane strain quadrilateral elements with reduced integration (type
CPE4R in Abaqus).

Table 1: Material properties used in DEM simulation. | grain density | Young’s | Poisson
---|---|---|---
| ($kg/m^{3}$) | modulus ($GPA$) | ratio
disk | 2500 | 4.5 | 0.2
grains | 2600 | 7.0 | 0.2

An automatically chosen, flexible time increment was used to execute the
simulations. This selection is a feature of Abaqus where a stable time
increment is computed as a factor of the shortest time for a dilatational wave
to cross an element in the model. To form the granular bed, the grains are
allowed to fall freely and settle under Earth’s gravity into the rigid bin.
The tangential interaction between all of the grains as well as the impactor
are defined using a friction coefficient of $\mu=0.84$, appropriate for coarse
gravel (Wright et al., 2021). Normal interactions are dictated by the material
properties of the involved bodies. (Hibbitt, 1998) A ramped acceleration-based
shake is applied to the settled granular bed to further settle the system. The
details of the vertical shake procedure are included in Table 2. Once gravity
is returned to 9.81 $m/s^{2}$, the simulation is allowed to run from 1.1s to
1.3s such that the total energy - calculated by Abaqus as the sum of strain
energy, kinetic energy, viscous dissipation, frictional dissipation, and
internal heat energy - stabilizes to a constant value.

Table 2: Details of the gravity-based shake applied during the settling simulation to further settle the granular bed. Negative accelerations are directed downwards. sim time ($s$) | 0 | 0.3 | 0.55 | 0.8 | 1.1 | 1.4
---|---|---|---|---|---|---
gravity ($m/s^{2}$) | -9.81 | -9.81 | +5.0 | -12.0 | -9.81 | -9.81

Figure 1: The granular bed is formed from an initial sparse patterning of
randomly oriented and shifted grains (left) that are allowed to fall and
settle under gravity (right).

The dimensions of the granular bed are approximately 0.32 m long and 0.16 m
deep. The depth is intentionally about 10 times the diameter of the impactor.
The length is such that, for a 30° impact at 10 m/s, the peak reaction force
felt by the rigid container is approximately 5% of the peak reaction force
felt by the projectile. This indicates the vast majority of energy is
dissipated within the granular material and that effects of the rigid
container on projectile response are minimal.

An impact is simulated by accelerating the projectile, shown in starting
position about 2.2cm above the bed in Figure 2a, with a one millisecond
impulse to achieve desired impact velocity and approach angle. The impactor is
also subjected to Earth-like gravity. It is assumed that the impact angle and
velocity is minimally impacted by gravity acting on the projectile; for the
lowest velocity and shallowest impacts, which are most affected, gravity was
observed to change impact angle by only a couple of degrees and velocity by a
small fraction of a meter-per-second.

An example snap-shot of each impact behavior type is shown in Figure 2. A
ricochet event occurs when the impactor departs the initial crater, breaking
contact with the granular bed (see Figure 2b). The roll-out event occurs when
the impactor rolls out of the initial impact crater but does not break contact
with the granular surface (see Figure 2c). The full-stop impact behavior
occurs when the impactor comes to a complete stop in the initial impact crater
(see Figure 2d).

(a) Initial configuration of disc relative to grains

(b) Ricochet impact event

(c) Roll-out impact event

(d) Full-stop impact event

Figure 2: The initial configuration of impactor and granular bed system in the
DFEM model with impact angle convention (a) and final frame grabs of the three
impact types (b-d). Panels (b) through (d) do not display the full depth of
the granular bed for space-saving in print.

We carry out two batches of impact simulations. Section 3.1 reviews the macro
behavior of two-hundred eighty six impact simulations with velocities ranging
from 1 m/s to 7 m/s and angles between 20° and 70° (measured from horizontal).
These ranges are selected to be similar to those studied by Wright et al.
(2020, 2021) to ensure a full span of impactor behaviors. Section 3.2 studies
eighteen simulations; nine at an impact angle of 30° ranging from 2 m/s to 10
m/s and nine at 60° over the same velocity range. The majority of impact
simulations are conducted across twelve cores to maintain a roughly one-hour
wall time. A subset of simulations were run on four, six, and twelve cores and
each model returned identical results, confirming that the number of cores
does not introduce artefacts to the impact model. To present this work on the
same plane as Wright et al. (2020), we also report the dimensionless Froude
number of each impact, calculated in terms of projectile velocity $v$, radius
$r$, and gravitational acceleration $g$ by the expression $v/\sqrt{gr}$.

##  3 Results & Discussions

###  3.1 Impact behavior mapping

In this section, we carry out analysis of the simulations in order to
understand how impactor behavior evolves with impact angle and velocity.
Denoting each set of initial conditions and their corresponding post-impact
behavior as a distinct marker, Figure 3a shows the impact behavior map
obtained from the setup in Figure 2a. Red squares represent ricochet impact
events (ric), blue diamonds represent roll-out events (ro), and black circles
represent full-stop events (fs). Red diamonds represent a hybrid
ricochet/roll-out event (ro/ric), where the impactor may have become partially
airborne after contact as it moves away from the impact crater. Blue squares
represent a hybrid full-stop/roll-out event (fs/ro), where the impact clearly
does not settle in the basin of the initial crater but does not completely
escape it either.

Inspection of Figure 3a reveals several features of interest. In the lower
right region at higher velocities and shallower impact angles, the ricochet
impact response dominates. In the upper left region at lower velocities and
steeper impact angles, full-stop behavior dominates. This behavior also
appears in a band across the top of the figure, for the steepest impact angle
of 70°. In-between the full-stop zones and ricochet zones, the roll-out events
appear as a transitional behavior. These regions and their relative positions
are consistent with experimental findings by Wright et al. (2020), though it
is noted that quantitative agreement is neither expected nor observed due to
numerous differences in system specifications such as grain-to-projectile size
ratio and the triangular grain geometry. Experimental work is also three
dimensional whereas simulations are constrained to two dimensions, which may
greatly impact comparisons.

(a)

(b)

Figure 3: Impact event classification as a function of impact angle and impact
velocity for (a) granular bed AA from Figure 2a and (b) granular bed BB.

In Figure 3a, out-of-place roll-out behaviors are noticed at the steepest
impact angles and highest velocities. This observation spurs the question, how
stable is this behavior map? To answer, we re-generated the granular bed of
Figure 2a using identical parameters and initial patterning. This produces a
‘new’ granular bed, denoted bed BB (and renaming the original granular system
bed AA), that has the same granular properties with marginally different
packing characteristics. We repeated the same slew of impact tests on the
newly generated bed and classified the impact behaviors, plotting the results
in Figure 3b.

There are clear differences between the plots of Figure 3. The stochastic
roll-out boundary region observed in Figure 3a is no longer, replaced by a
much more regular roll-out region. The ricochet dominated region in Figure 3b
is reduced, with ricochets no longer presenting at steeper impact angles.
Nonetheless, there are similarities between the two figures that should not go
amiss. Full-stop behaviors still dominate in the upper left regions while
ricochets still dominate in the lower right. The ‘y-intercept’ of the boundary
between full-stop and ricochet is still close to 20°.

Wright et al. (2021) suggest the variation between the two plots is expected
due to the small impactor-to-grain size ratio (5.38). They suggest that
stochasticity in impact behaviors is anticipated for lower ratios whereas
stable and repeatable impact behaviors are produced for higher ratios.
Additional to their argument, we theorize that the angular nature of the
grains is a factor, exacerbating the effect that the new packing configuration
has on the projectile’s motion; looser packings allow easier deformation and
regular behaviors whereas tighter packings are more interlocked and firm,
producing stochastic responses. Tighter packings may also emphasize effects of
surface topography which is highly variable due to the angular grain geometry.

Regardless, revealing the systemic factors that cause these variations cannot
be fully assessed by visual observation of the impacts alone. In the next
section, an algorithm for evaluating packing fraction is presented and the
variations between the plots of Figure 3 are re-visited in a quantifiable
manner.

###  3.2 Sub-surface quantification

####  3.2.1 Packing fraction analysis

In this section, post-processing of simulation data from impacts angled at 30°
and 60° reveal information about the sub-surface granular dynamics that occur
during a surface impact. This is a meaningful study because such dynamics
contribute to the behavior of the impactor due to the dissipation of energy
through grains. Such pathways of dissipation are affected by how grains are
positioned relative to their neighbors and the incoming impulse from the
impactor. For instance, loose packing of grains could lead to compaction and
slippage between grains and a gradual dissipation of energy by friction,
creating an overall cushioning effect on an impactor. Tighter grain packings
may lead to direct transmission of forces along force chains, producing a
stouter impactor behavior reminiscent of solid-on-solid impact.

Impactor behaviors and sub-surface granular mechanics can be tied together by
considering changes in grain packing density and its distribution. Typically
packing analysis is done by evaluating a system’s overall change in packing
density, though discrete simulation method can provide a better understanding
of such change by considering its spatial distribution.

The void mapping procedure stems originally from Point Voronoi tessellation
(Voronoi, 1908), a method that can take a random collection of points and
divide the surrounding space into convex polygonal regions belonging to each
point. Point Voronoi tessellation may also be used for collections of uniform
circles (An et al., 2021) or lightly adapted with weighting techniques for
poly-disperse circles (Morley and Wilson, 2020).

However, Point Voronoi tessellation is insufficient for mapping void space
between more complex geometries because centroids of particles are no longer
equidistant from the edges. Therefore, tessellated regions may cut off parts
of the parent grain, causing inaccuracies in packing fraction calculation.
This limitation inspired the Set Voronoi tessellation method (Schaller et al.,
2013; Dong et al., 2016; Zhang et al., 2021), an evolution of the Point
Voronoi procedure for irregular geometries that uses additional discretization
and unification steps. In this work, MATLAB’s built-in Point Voronoi function
(simply named ‘voronoi’) has been merged with the additional steps of Set
Voronoi tessellation to permit accurate division of void space among the
triangular grains. The details of the process for two dimensions are described
below. With some modification, such a process should be achievable in three-
dimensions, but this is beyond the scope of the present paper.

  1. 1.

Raw positional data is read in and rearranged into a more accessible format.

  2. 2.

Each grain is shrunk to improve the accuracy of the Voronoi tessellation
procedure. Relative to the centroid of each grain, vertex coordinates are
adjusted to reduce the size of the grain by a certain factor.

  3. 3.

The edges of grains are discretized. Based on a certain length increment,
additional points are created between the vertices and added to the overall
point-set for each grain.

  4. 4.

A set of points bounding the full granular bed is generated to constrain the
Voronoi tessellation and maintain accuracy at the edges of the granular bed.

  5. 5.

Combining the data points belonging to the grains with the points belonging to
the boundary, the whole structure is passed to MATLAB’s built-in ‘voronoi’
function which completes the first step of the tessellation.

  6. 6.

Resulting polygonal cells created for each point are united together according
to their parent grains. The merger of these polygons produces the final
‘spatial cell’ of the granular bed that a particular grain owns.

  7. 7.

Finally, a local packing fraction ($\phi_{local}$) is estimated for each grain
by dividing the area of the grain itself ($A_{grain}$) by the area of its
spatial cell ($A_{cell}$) such that $\phi_{local}=A_{grain}/A_{cell}$.

(a) Original grain data

(b) Discretized edges

(c) Point Voronoi tessellation

(d) Merged final cells

(e) Merged final cells

Figure 4: The tessellation procedure for calculation of grain-level packing
fraction.

Figure 4a depicts a few grains positioned based on data exported from
simulations. Figure 4b depicts these grains after they are shrunken and
discretized with additional points along its boundary. Figure 4c shows the
result of MATLAB’s Voronoi tessellation utility which divides up the space in-
between points to spatial allotment closest to a particular point. Figure 4d
illustrates how the post-processing script merges the polygons belonging to
the larger grain into one larger polygon. Figure 4e shows the final result of
the tessellation procedure, with a polygonal cell that represents the region
of the granular bed belonging to a grain. The original grains are overlaid in
Figure 4e to show the two regions which are used for computing packing
fraction. Each color in Figure 4e corresponds to the core that calculated that
Voronoi cell, as parallel processing was used to expedite post-processing.

The procedure outlined allows for calculation of grain-level packing fraction
and such a result is shown in Figure 5a, depicting a homogeneous distribution
of overall packing density in the granular bed from Figure 2a before an
impact. A histogram of packing fraction measurements from Figure 5a is shown
in the black bars of Figure 5b. A neat bell curve is produced, indicating a
well-settled granular bed.

(a)

(b)

Figure 5: Evaluation of pre-impact granular bed quality using local packing
fraction estimates. A colored packing density plot shows the raw local packing
fraction measurements (a), and a histogram (b) illustrates the distribution of
these packing fraction measurements.

Our Voronoi-based packing analysis can re-visit the impact behavior maps of
Figure 3 and explain the differences between the two plots. Figure 5a depicts
packing measurements for the bed AA associated with Figure 3a. A similar
collection of packing measurements is also collected for bed BB, from Figure
3b, though not pictured. For each impact angle, the approximate point of
impact on the granular bed is determined and the impact site defined as the
area of grains within one projectile radius. Using the packing fraction
measurements for grains within each impact region, the mean packing fraction
of the impact site is calculated. Figure 6 presents this data for both
granular beds.

For the data of bed AA (red line in Figure 6), the mean packing fraction at
impact sites is shown to increase at steeper impact angles, reaching a maximum
tight packing of 0.83. This supports the stochastic behavior observed at the
upper right side in Figure 3a since it implies that the projectile is hitting
a much firmer surface. Therefore, the possibility of a full stop is lessened
and there is a higher tendency to display ricochet or roll-out events. In
contrast, bed BB data (blue line) displays looser packing at the steeper
impact sites, dropping as low as 0.66 locally. This indicates the projectile
is hitting a softer granular media that allows greater intrusion, explaining
the observation of more full-stop behavior in the upper right of the map in
Figure 3b. Finally, the greatest similarity between the plots of Figure 3 is
observed at the shallower impact angles. Figure 6 also upholds this, since the
mean site packing fraction is found to be fairly consistent ($\approx$0.74 -
0.76) between both versions of the granular bed below 45 °.

The mechanism causing firmer packings to produce stochastic behaviors is not
definitively known. For a hard and flat surface, instead of a granular
surface, one would still expect repeatable behavior from an angled impact. We
suggest that the reason stochasticity arises in our model is because of the
way large grains ‘rough-up’ the surface topography. As the projectile hits a
range of different points on the surface, it may encounter a wide range of
angled flat faces or even corners. While this condition is also true for
looser packings, the firmer packings’ diminished ability to dissipate energy
through grain sliding exacerbates the topographic effect, producing the
stochasticity we observe.

Figure 6: Pre-impact, mean packing fraction of impact site versus impact
angle.

We can also use the packing fraction analysis to assess how the granular bed
is re-configured globally and locally by the surface impact. Originally, to
demonstrate the reconfiguration, plain packing fraction measurements after
impact were assessed in the same manner shown in Figure 5a. In Figure 5b, the
after-impact histogram for a 10 m/s, 30° impact scenario is overlaid on the
pre-impact data with the red bars. The shift is slight and uninformative from
this perspective. The poor visualization also suggests that holistic packing
fraction for the entire bed, as would be typically measured with fluid
displacement or bulk density, would also fail to capture any grain
restructuring.

Instead, the per-grain change in packing fraction (relative packing density)
was found to be a better metric. Figure 7 presents packing density plots of
each velocity case for 60° impacts. Results for 30° simulations can be found
in Figure 8. For each grain, the packing fraction post-impact is subtracted
from its pre-impact value, so a positive change in packing fraction indicates
grain patterns tightening whereas a negative change implies loosening. The
upper and lower limits of the colorbar indicate a positive or negative change
in packing fraction of approximately 8.7%, which is taken from the standard
deviation of pre-impact packing fraction measurements in Figure 5b. The plots
of Figure 7 show how the breadth of disturbance changes with higher velocity
impacts.

Considering the 2 m/s density plot (Figure 7a), it is shown that only grains
immediately around the impact crater reconfigure asymmetrically, within about
a one to two impactor diameter radius. Disturbance is contained to this zone,
and surrounding surface and sub-surface grains remain relatively undisturbed.

Figure 7e shows that for an impact at 6 m/s the density plot has a region of
disturbed grains that has clearly grown with a tripling of the impactor
velocity. Interestingly, the right side of the crater shows a zone of grains
in which porosity has increased compared to other regions around the crater
where grains have packed closer together. Velocities equal to or greater than
6 m/s also show ejecta that has been strewn along the surface, evident from
the trail of overall shifted grains extending along the surface beyond the
crater.

Finally, at 10 m/s (Figure 7i), all aspects of the granular disturbance are
amplified. The region of disturbed grains around the crater is roughly about
four to five times the diameter of the impactor, the scattered ejecta layer
along the surface is much more prominent and several grains deep, and there
are even pockets of shifted grains deeper within the granular bed outside of
the primary disturbed zone. This latter observation is attributed to
structural instabilities within the granular bed that hold up for lower
velocities but begin to collapse for higher velocities and indicates that a
pressure pulse has propagated through the medium.

These findings intuitively demonstrate how the region of disturbed grains
increases steadily with impactor velocity in both the 30° and 60° impact
cases. The asymmetry of the sub-surface disturbance and ejecta spread is
another way in which 30° and 60° results may be compared. At the highest
velocities of 8 to 10 m/s, the immediate granular disruption underneath the
crater is more asymmetrical in the 30° case, though the more pronounced ejecta
spread of 60° overshadows this. Oddly, the same conclusion does not hold at
the lowest velocities of 2 to 4 m/s where sub-surface asymmetry actually
appears more pronounced for the 60° cases. For instance, at 2 m/s, the 60°
impact produces a much deeper and clearly more asymmetrical disturbance
compared to 30°. These observations suggest that the symmetry of the granular
disturbance is not solely influenced by the angle of the impact but also
affected by the nature of impact. In the full-stop events of 60° impacts, much
more energy is injected into the granular bed, amplifying the asymmetry at
high velocities. The 30° impacts are predominantly ricochets, so while higher
velocities yield increasingly glancing impacts, less energy is injected into
the bed and less asymmetry is observed.

(a) v = 2 m/s, Fr = 5.0

(b) v = 3 m/s, Fr = 7.5

(c) v = 4 m/s, Fr = 9.9

(d) v = 5 m/s, Fr = 12.4

(e) v = 6 m/s, Fr = 14.9

(f) v = 7 m/s, Fr = 17.4

(g) v = 8 m/s, Fr = 19.9

(h) v = 9 m/s, Fr = 22.4

(i) v = 10 m/s, Fr = 24.9

Figure 7: Colored relative packing density plots illustrating the local
packing fraction changes after impacts of various velocities at 60° from the
horizontal. Velocity and Froude number of each impact are reported.

(a) v = 2 m/s, Fr = 5.0

(b) v = 3 m/s, Fr = 7.5

(c) v = 4 m/s, Fr = 9.9

(d) v = 5 m/s, Fr = 12.4

(e) v = 6 m/s, Fr = 14.9

(f) v = 7 m/s, Fr = 17.4

(g) v = 8 m/s, Fr = 19.9

(h) v = 9 m/s, Fr = 22.4

(i) v = 10 m/s, Fr = 24.9

Figure 8: Colored relative packing density plots illustrating results for 30°
impacts.

####  3.2.2 Strain distributions

Grain reconfiguration is not necessarily accompanied by a change in packing
fraction, so other perspectives should be sought out. Therefore, we compute
the strain on each grain using displacement fields. A large-strain Eulerian
formulation (Lai et al., 2010), shown in Equation 1 in summation index
notation, shows how we calculate the strain tensor from the displacement
fields of the grains $u$.

$\epsilon_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)-\frac{1}{2}\frac{\partial
u_{m}}{\partial x_{i}}\frac{\partial u_{m}}{\partial x_{j}}$ (1)

Here, it is implied that we sum over index $m$.

To collapse the strain tensor for each grain to a single representative value,
the von Mises distortion
111urlhttps://dianafea.com/manuals/d944/Analys/node405.html is calculated as
follows

$\boldsymbol{\epsilon_{\text{von
Mises}}}=\sqrt{\frac{2}{3}\>\>\left(\boldsymbol{\epsilon}^{\prime}:\boldsymbol{\epsilon}^{\prime}\right)}.$
(2)

where $\boldsymbol{\epsilon}^{\prime}$ represents the deviatoric strain
tensor, with components given by
$\epsilon_{ij}^{\prime}=\epsilon_{ij}-\frac{1}{3}\>\delta_{ij}\>\epsilon_{ij}$
and $\delta_{ij}$ represents the Kronecker delta function. The ‘:’ symbol
represents a sum over both indices,
$\epsilon^{\prime}:\epsilon^{\prime}=\sum_{ij}\epsilon^{\prime}_{ij}\epsilon^{\prime}_{ij}$.

To visualize the strained region around the impact crater, a two-color contour
plot is created, the threshold of which is set at a strain level for the von
Mises distortion $\epsilon_{\text{von Mises}}=3\times 10^{-4}$. The colors,
then, imply that white region has permanently deformed while the black region
has deformed within elastic limit. This value is based on a quasi-static flow
strain measurement calculated from separate simulation. Using the same
granular bed, a 5 cm horizontal rigid plate is positioned just above the
surface. A displacement boundary condition pushes the plate downward at a slow
rate of 0.25 mm/s for 1 cm, breaking through the granular surface. Reading the
reaction force against the plate reveals the time at which the surface yields
and then post-processing the grain displacement data at that frame gives the
local strain corresponding to the yield, found to be $3\times 10^{-4}$. Figure
9 displays the strain contours for the 60° impact cases. Figure 10 shows
results for 30°. Within these plots, the white region bound by red lines is
the granular material strained beyond the quasi-static flow strain threshold.
The black region is granular media that is not strained beyond that point. The
upper red line marks the surface of the granular media, the lower red line
follows the contour boundary.

Figure 9a shows the distortion contour plot for the lowest velocity 2 m/s
impact. As was shown with the relative packing fraction results, the strained
region of grains is largely confined to the impact crater. A moderate velocity
level of 6 m/s in Figure 9e demonstrates a larger region of strained grains
with a noticeable strained layer to the right of the impact crater, indicating
ejecta along the surface. In Figure 9i at 10 m/s, all features of the 6 m/s
case are further amplified except, seemingly, for the depth of the perturbed
region. Otherwise, the layer of surface-scattered material is much deeper as
is the span of the directly strained region around the impact site. Figure 9
also illustrates the asymmetry in substrate response caused by the oblique
impact mentioned in the packing fraction analysis.

(a) v = 2 m/s, Fr = 5.0

(b) v = 3 m/s, Fr = 7.5

(c) v = 4 m/s, Fr = 9.9

(d) v = 5 m/s, Fr = 12.4

(e) v = 6 m/s, Fr = 14.9

(f) v = 7 m/s, Fr = 17.4

(g) v = 8 m/s, Fr = 19.9

(h) v = 9 m/s, Fr = 22.4

(i) v = 10 m/s, Fr = 24.9

Figure 9: Contour plots of Von Mises strain distribution in the granular bed
after impacts of various velocities at 60° from the horizontal. Velocity and
Froude number of each impact are reported. Red lines bound the white strained
region of material, denoted the ‘skin zone’. Black indicates granular media
that is not strained beyond the quasi-static flow strain threshold.

The observation that depth of the strained region did not continue to grow for
higher velocity impacts at 60° spurred further study of the strain contours.
The enclosure of the two red lines, marking the strained media, is denoted the
‘skin zone’. The skin zone, therefore, represents the region of grains that
have been permanently strained beneath the surface as well as the layer of
ejected grains that have settled atop the surface. The difference in height
between corresponding points of upper and lower boundaries produces a ‘skin
depth’ measurement.

(a) v = 2 m/s, Fr = 5.0

(b) v = 3 m/s, Fr = 7.5

(c) v = 4 m/s, Fr = 9.9

(d) v = 5 m/s, Fr = 12.4

(e) v = 6 m/s, Fr = 14.9

(f) v = 7 m/s, Fr = 17.4

(g) v = 8 m/s, Fr = 19.9

(h) v = 9 m/s, Fr = 22.4

(i) v = 10 m/s, Fr = 24.9

Figure 10: Contour plots of Von Mises strain distribution for 30° impacts.

Figure 11 shows skin depth measurements versus horizontal position for 30° and
60° impacts. For easier viewing, the results are smoothed with a moving
average and a vertical shift to align the leftmost points. The smoothing helps
to view the trend in skin zone geometry with increasing velocity impacts. In
Figure 11a, showing 30° impacts, a fairly straightforward trend is seen where
the skin zone deepens and widens in a steady manner as velocity increases.
Contrast to this the 60° impacts in Figure 11b. Starting from the blue line,
representing the 2 m/s case, a single sharp peak is shown with a general
return to zero moving further to the right. The next five impact cases, up to
about 6 m/s, show that same sharp peak continue to grow accompanied by a
levelling out of the zone moving to the right. However, for 7 m/s impacts and
above, the peak no longer grows and instead reduces in height and widens out
as lateral strain appears to dominate over strains deeper in the bed.

The observed trend of peaks in the skin zone depth measurements for 60°
suggests a potential change in the response of the granular media to the
surface impacts of increasing velocity. Up until a point, the strained region
does grow deeper with increasing velocity. However, around 7 m/s for this
specific scenario, the peak shortens with deeper permanently deformed regions
spanning further to the left and right of the impact crater. We have not been
able to identify this type of trend in literature though it suggests that the
energy dissipation mechanism that arrests projectile motion may change with
increasing velocity, leading to a more lateral plastic deformation.

(a) 30° impact data

(b) 60° impact data

Figure 11: Plots of skin depth - the difference in height between
corresponding points on the upper and lower boundaries shown in versus
horizontal position.

To further investigate this trend’s dependence on velocity, the raw skin depth
measurements are averaged to create a data point for mean skin depth at
corresponding impact velocity. The data are displayed in Figure 12, where mean
skin depth (MSD) measurements for 30° impacts and 60° impacts are marked by
red squares and blue circles, respectively. Additionally, for each simulation,
measurements of the transient crater depth (TCD) - the deepest point to exist
in the crater during the impact event - are also displayed for 30° and 60°
impacts as red and blue dashed lines.

From the four data sets of Figure 12, some interesting discussions arise. The
MSD data from both 30° and 60° simulations roughly coincide each other until
around 9 m/s, despite predominantly ricochet behavior in the 30° impact
simulations and full-stop behavior in the 60° cases. For instance, at the 5m/s
case, the mean skin depth of the 30° impact is near identical to 60°. The TCD
data, on the other hand, show that crater depths are generally greater for
across the velocity range for 60° impacts than for 30° impacts. Furthermore,
the TCD data for 60° also closely aligns with the MSD data for 60°. This last
observation suggests that for such an angle in these specific conditions,
measuring crater depth at the surface directly reveals the mean depth of the
disturbed grain region beneath. The analysis of skin depth reveals that for
steeper impacts, the trends of depth of the skin zone is expected to break at
smaller velocities since the granular bed can perform as a monolithic solid at
larger depths from the free surface.

Figure 12: Plots of Mean Skin Depth (MSD) and measured Transient Crater Depth
(TCD) for 30° (red) and 60° (blue) impacts.

##  4 Conclusion

We have presented simulations that study low-velocity oblique impacts into a
non-spherical granular media. Citing trends in granular dynamics literature
toward modeling angular grain geometries as a better representation of non-
terrestrial regolith, this work has developed discrete 2D simulations based on
finite elements, a technique called the Discrete Finite Element Method (DFEM),
that models the simplest possible polygonal grain geometry – a triangle.
Maintaining Earth-like gravity to enable comparison to experimental
literature, impacts at various angles and velocities of a circular projectile
into a bed of triangular grains have been presented. Our results demonstrate
the three classes of impact behavior and reproduce the velocity and impact
angle dependent division between ricochet and roll-out behavior, recently
observed in experiments (Wright et al., 2020).

We find that visualizing the per-grain local packing fraction variance within
8.7% reveals a sub-surface perturbed region that is lopsided with respect to
the impact point. Local porosity variations are focused in front of the impact
site by a distance that is similar in size to the crater radius. The region in
front of the impact can be compacted while substrate nearer the surface
increases in porosity. By strain analyses, we introduce the notion of ‘skin
zone’, or the region beneath the surface that is permanently deformed beyond
the quasi-static yield threshold. Measurements of this skin zone have
demonstrated trends that suggest granular dissipative mechanisms change when
faced with increasing impact velocity,especially at steeper impact angles.

There is more work to do. Mono-disperse triangular grains are not reflective
of any naturally occurring material and their use here was to ease development
of the DFEM model. Various gravity conditions, extending to micro-gravity,
will be studied. Polydispersity, another defining characteristic of regolith,
will also need to be considered to bring simulations closer in line with the
environments we wish to study. The finite elements basis also permits many
more perspectives beyond the grain coordinate and displacement based metrics
we have evaluated in this paper, such as nodal data. Outputting fields such as
nodal stress direct from the impact model will help us to visualize force
networks that propagate throughout granular media. An example of this output
is shown in Appendix II for a 3.5 m/s impact at 45°, as a showcase for future
efforts. While detailed analysis of stress contours is beyond the scope of
this paper, it will be a prime focus of the next steps in this modeling work.

Developing a comprehensive understanding of the phenomena entailed in low-
velocity impacts would support space missions targeting asteroid and planetary
exploration in a myriad of ways. Our work demonstrates the importance of
packing fraction in determining the behavior of an impactor, an especially key
realization for granular media such as regolith with very high packing
potential. Through our strain analyses, we have shown the existence of and
measured the depth to which granular media is disrupted, a finding that, when
reversed, informs on how deep the composition of a granular bed must be
understood to anticipate the behavior of a surface impact. Ultimately, these
conclusions could have ripple effects into other areas of research beyond just
impact events. Improved simulations of regolith media and better
virtualization of its behavior could aid development of surface vehicles,
methods for terrain cultivation and terraforming, drilling equipment, or in-
situ building materials.

## Acknowledgements

This material is based on work partially supported by NASA grant 80NSSC21K0143

## References

  * An et al. (2021) An, X., Dong, K., Yang, R., Zou, R., Yu, A., 2021\.  On the relationships between structural properties and packing density of uniform spheres.  Powder Technology 388, 139–148.  URL: https://www.sciencedirect.com/science/article/pii/S0032591021003661, doi:https://doi.org/10.1016/j.powtec.2021.04.079.
  * Arévalo et al. (2010) Arévalo, R., Zuriguel, I., Maza, D., 2010.  Topology of the force network in the jamming transition of an isotropically compressed granular packing.  Phys. Rev. E 81, 041302\.  URL: https://link.aps.org/doi/10.1103/PhysRevE.81.041302, doi:10.1103/PhysRevE.81.041302.
  * Bagatin et al. (2001) Bagatin, A., Petit, J.M., Farinella, P., 2001.  How many rubble piles are in the asteroid belt?  Icarus 149, 198–209.  URL: https://www.sciencedirect.com/science/article/pii/S001910350096531X, doi:https://doi.org/10.1006/icar.2000.6531.
  * Bagherzadeh Kh. et al. (2011) Bagherzadeh Kh., A., Mirghasemi, A., Mohammadi, S., 2011.  Numerical simulation of particle breakage of angular particles using combined dem and fem.  Powder Technology 205, 15–29.  URL: https://www.sciencedirect.com/science/article/pii/S0032591010003815, doi:https://doi.org/10.1016/j.powtec.2010.07.034.
  * Brisset et al. (2018) Brisset, J., Colwell, J., Dove, A., Abukhalil, S., Cox, C., Mohammed, N., 2018\.  Regolith behavior under asteroid-level gravity conditions: low-velocity impact experiments.  Progress in Earth and Planetary Science 5, 73.  URL: https://doi.org/10.1186/s40645-018-0222-5, doi:10.1186/s40645-018-0222-5.
  * Burns (1975) Burns, J.A., 1975.  The angular momenta of solar system bodies: Implications for asteroid strengths.  Icarus 25, 545–554.  URL: https://www.sciencedirect.com/science/article/pii/0019103575900354, doi:https://doi.org/10.1016/0019-1035(75)90035-4.
  * Chen et al. (2011) Chen, Q., Andrade, J.E., Samaniego, E., 2011.  Aes for multiscale localization modeling in granular media.  Computer Methods in Applied Mechanics and Engineering 200, 2473–2482.  URL: https://www.sciencedirect.com/science/article/pii/S0045782511001630, doi:https://doi.org/10.1016/j.cma.2011.04.022.
  * Clark et al. (2012) Clark, A.H., Kondic, L., Behringer, R.P., 2012.  Particle scale dynamics in granular impact.  Physical review letters 109, 238302–238302.
  * Cundall and Strack (1979) Cundall, P.A., Strack, O.D.L., 1979\.  A discrete numerical model for granular assemblies.  Géotechnique 29, 47–65.  URL: https://doi.org/10.1680/geot.1979.29.1.47, doi:10.1680/geot.1979.29.1.47, arXiv:https://doi.org/10.1680/geot.1979.29.1.47.
  * Dong et al. (2016) Dong, K., Wang, C., Yu, A., 2016.  Voronoi analysis of the packings of non-spherical particles.  Chemical Engineering Science 153, 330--343.  URL: https://www.sciencedirect.com/science/article/pii/S0009250916303785, doi:https://doi.org/10.1016/j.ces.2016.07.013.
  * Dowling and Dowling (2013) Dowling, D.R., Dowling, T.R., 2013\.  Scaling of impact craters in unconsolidated granular materials.  American journal of physics 81, 875--878.
  * Ferrari et al. (2019) Ferrari, F., Lavagna, M., Blazquez, E., 2019.  A parallel-GPU code for asteroid aggregation problems with angular particles.  Monthly Notices of the Royal Astronomical Society 492, 749--761.  URL: https://doi.org/10.1093/mnras/stz3458, doi:10.1093/mnras/stz3458, arXiv:https://academic.oup.com/mnras/article-pdf/492/1/749/31731024/stz3458.pdf.
  * Ferrari and Tanga (2020) Ferrari, F., Tanga, P., 2020\.  The role of fragment shapes in the simulations of asteroids as gravitational aggregates.  Icarus 350, 113871\.  URL: https://www.sciencedirect.com/science/article/pii/S0019103520302517, doi:https://doi.org/10.1016/j.icarus.2020.113871.
  * Fu et al. (2006) Fu, X., Dutt, M., Bentham, A.C., Hancock, B.C., Cameron, R.E., Elliott, J.A., 2006.  Investigation of particle packing in model pharmaceutical powders using x-ray microtomography and discrete element method.  Powder Technology 167, 134--140.  URL: https://www.sciencedirect.com/science/article/pii/S0032591006002348, doi:https://doi.org/10.1016/j.powtec.2006.06.011.
  * Goldman and Umbanhowar (2008) Goldman, D.I., Umbanhowar, P., 2008\.  Scaling and dynamics of sphere and disk impact into granular media.  Physical review. E, Statistical, nonlinear, and soft matter physics 77, 021308--021308.
  * Gou et al. (2018) Gou, S., Yue, Z., Di, K., Liu, Z., 2018.  A global catalogue of ceres impact craters greater than 1 km and preliminary analysis.  Icarus 302, 296--307.  URL: https://www.sciencedirect.com/science/article/pii/S0019103517303585, doi:https://doi.org/10.1016/j.icarus.2017.11.028.
  * Guo and Zhao (2016) Guo, N., Zhao, J., 2016.  3d multiscale modeling of strain localization in granular media.  Computers and Geotechnics 80, 360--372.  URL: https://www.sciencedirect.com/science/article/pii/S0266352X16300027, doi:https://doi.org/10.1016/j.compgeo.2016.01.020.
  * Güner et al. (2015) Güner, F., Ömer Necati Cora, Sofuoğlu, H., 2015.  Numerical modeling of cold powder compaction using multi particle and continuum media approaches.  Powder Technology 271, 238--247.  URL: https://www.sciencedirect.com/science/article/pii/S0032591014009061, doi:https://doi.org/10.1016/j.powtec.2014.11.008.
  * Hibbitt (1998) Hibbitt, K.a., 1998.  ABAQUS: Theory manual.  Hibbitt, Karlsson, and Sorensen, Pawtucket, R.I.
  * Jia et al. (2012) Jia, T., Zhang, Y., Chen, J., 2012.  Simulation of granular packing of particles with different size distributions.  Computational Materials Science 51, 172--180.  URL: https://www.sciencedirect.com/science/article/pii/S0927025611004368, doi:https://doi.org/10.1016/j.commatsci.2011.07.044.
  * Kabir et al. (2008) Kabir, M.A., Lovell, M.R., Higgs, C.F., 2008.  Utilizing the explicit finite element method for studying granular flows.  Tribology Letters 29, 85--94.  URL: https://doi.org/10.1007/s11249-007-9285-y, doi:10.1007/s11249-007-9285-y.
  * Knuth et al. (2012) Knuth, M., Johnson, J., Hopkins, M., Sullivan, R., Moore, J., 2012.  Discrete element modeling of a mars exploration rover wheel in granular material.  Journal of Terramechanics 49, 27--36.  URL: https://www.sciencedirect.com/science/article/pii/S0022489811000619, doi:https://doi.org/10.1016/j.jterra.2011.09.003.
  * Kondic et al. (2012) Kondic, L., Fang, X., Losert, W., O’Hern, C.S., Behringer, R.P., 2012.  Microstructure evolution during impact on granular matter.  Physical review. E, Statistical, nonlinear, and soft matter physics 85, 011305--011305.
  * Lai et al. (2010) Lai, W.M., Rubin, D., Krempl, E., 2010.  Introduction to Continuum Mechanics.  Butterworth-Heinemann/Elsevier.
  * Leinhardt et al. (2000) Leinhardt, Z.M., Richardson, D.C., Quinn, T., 2000.  Direct n-body simulations of rubble pile collisions.  Icarus 146, 133--151.  URL: https://www.sciencedirect.com/science/article/pii/S001910350096370X, doi:https://doi.org/10.1006/icar.2000.6370.
  * Li and McDowell (2018) Li, H., McDowell, G.R., 2018\.  Discrete element modelling of under sleeper pads using a box test.  Granular Matter 20, 26\.  URL: https://doi.org/10.1007/s10035-018-0795-0, doi:10.1007/s10035-018-0795-0.
  * Lobo-Guerrero and Vallejo (2006) Lobo-Guerrero, S., Vallejo, L.E., 2006\.  Discrete element method analysis of railtrack ballast degradation during cyclic loading.  Granular Matter 8, 195\.  URL: https://doi.org/10.1007/s10035-006-0006-2, doi:10.1007/s10035-006-0006-2.
  * Luo et al. (2017) Luo, T., Ooi, E., Chan, A., Fu, S., 2017.  The combined scaled boundary finite-discrete element method: Grain breakage modelling in cohesion-less granular media.  Computers and Geotechnics 88, 199--221.  URL: https://www.sciencedirect.com/science/article/pii/S0266352X17300824, doi:https://doi.org/10.1016/j.compgeo.2017.03.012.
  * McDowell and Li (2016) McDowell, G.R., Li, H., 2016\.  Discrete element modelling of scaled railway ballast under triaxial conditions.  Granular Matter 18, 66\.  URL: https://doi.org/10.1007/s10035-016-0663-8, doi:10.1007/s10035-016-0663-8.
  * Michel et al. (2020) Michel, P., Ballouz, R.L., Barnouin, O.S., Jutzi, M., Walsh, K.J., May, B.H., Manzoni, C., Richardson, D.C., Schwartz, S.R., Sugita, S., Watanabe, S., Miyamoto, H., Hirabayashi, M., Bottke, W.F., Connolly, H.C., Yoshikawa, M., Lauretta, D.S., 2020.  Collisional formation of top-shaped asteroids and implications for the origins of ryugu and bennu.  Nature Communications 11, 2655\.  URL: https://doi.org/10.1038/s41467-020-16433-z, doi:10.1038/s41467-020-16433-z.
  * Michel et al. (2002) Michel, P., Tanga, P., Benz, W., Richardson, D.C., 2002\.  Formation of asteroid families by catastrophic disruption: Simulations with fragmentation and gravitational reaccumulation.  Icarus 160, 10--23.  URL: https://www.sciencedirect.com/science/article/pii/S0019103502969484, doi:https://doi.org/10.1006/icar.2002.6948.
  * Moreno-Atanasio et al. (2010) Moreno-Atanasio, R., Williams, R.A., Jia, X., 2010.  Combining x-ray microtomography with computer simulation for analysis of granular and porous materials.  Particuology 8, 81--99.  URL: https://www.sciencedirect.com/science/article/pii/S1674200110000210, doi:https://doi.org/10.1016/j.partic.2010.01.001.
  * Morley and Wilson (2020) Morley, D.O., Wilson, M., 2020\.  Voronoi diagrams in quasi-2d hard sphere systems.  Journal of Statistical Mechanics: Theory and Experiment 2020, 093201.  URL: https://doi.org/10.1088/1742-5468/aba7af, doi:10.1088/1742-5468/aba7af.
  * Nakamura et al. (2011) Nakamura, T., Noguchi, T., Tanaka, M., Zolensky, M.E., Kimura, M., Tsuchiyama, A., Nakato, A., Ogami, T., Ishida, H., Uesugi, M., Yada, T., Shirai, K., Fujimura, A., Okazaki, R., Sandford, S.A., Ishibashi, Y., Abe, M., Okada, T., Ueno, M., Mukai, T., Yoshikawa, M., Kawaguchi, J., 2011\.  Itokawa dust particles: A direct link between s-type asteroids and ordinary chondrites.  Science 333, 1113--1116.  URL: https://www.science.org/doi/abs/10.1126/science.1207758, doi:10.1126/science.1207758, arXiv:https://www.science.org/doi/pdf/10.1126/science.1207758.
  * Nitka and Tejchman (2020) Nitka, M., Tejchman, J., 2020\.  Comparative dem calculations of fracture process in concrete considering real angular and artificial spherical aggregates.  Engineering Fracture Mechanics 239, 107309.  URL: https://www.sciencedirect.com/science/article/pii/S0013794420308912, doi:https://doi.org/10.1016/j.engfracmech.2020.107309.
  * Noguchi et al. (2021) Noguchi, R., Hirata, N., Hirata, N., Shimaki, Y., Nishikawa, N., Tanaka, S., Sugiyama, T., Morota, T., Sugita, S., Cho, Y., Honda, R., Kameda, S., Tatsumi, E., Yoshioka, K., Sawada, H., Yokota, Y., Sakatani, N., Hayakawa, M., Matsuoka, M., Yamada, M., Kouyama, T., Suzuki, H., Honda, C., Ogawa, K., Kanamaru, M., ichiro Watanabe, S., 2021\.  Crater depth-to-diameter ratios on asteroid 162173 ryugu.  Icarus 354, 114016\.  URL: https://www.sciencedirect.com/science/article/pii/S0019103520303791, doi:https://doi.org/10.1016/j.icarus.2020.114016.
  * Pravec and Harris (2000) Pravec, P., Harris, A.W., 2000\.  Fast and slow rotation of asteroids.  Icarus 148, 12--20.  URL: https://www.sciencedirect.com/science/article/pii/S0019103500964820, doi:https://doi.org/10.1006/icar.2000.6482.
  * Raisianzadeh et al. (2018) Raisianzadeh, J., Mirghasemi, A.A., Mohammadi, S., 2018.  2d simulation of breakage of angular particles using combined dem and xfem.  Powder Technology 336, 282--297.  URL: https://www.sciencedirect.com/science/article/pii/S0032591018304479, doi:https://doi.org/10.1016/j.powtec.2018.06.006.
  * Raisianzadeh et al. (2019) Raisianzadeh, J., Mohammadi, S., Mirghasemi, A.A., 2019.  Micromechanical study of particle breakage in 2d angular rockfill media using combined dem and xfem.  Granular Matter 21, 48\.  URL: https://doi.org/10.1007/s10035-019-0904-8, doi:10.1007/s10035-019-0904-8.
  * Reimann et al. (2017) Reimann, J., Vicente, J., Brun, E., Ferrero, C., Gan, Y., Rack, A., 2017.  X-ray tomography investigations of mono-sized sphere packing structures in cylindrical containers.  Powder Technology 318, 471--483.  URL: https://www.sciencedirect.com/science/article/pii/S0032591017304072, doi:https://doi.org/10.1016/j.powtec.2017.05.033.
  * Richardson et al. (2005) Richardson, D.C., Elankumaran, P., Sanderson, R.E., 2005.  Numerical experiments with rubble piles: equilibrium shapes and spins.  Icarus 173, 349--361.  URL: https://www.sciencedirect.com/science/article/pii/S0019103504003124, doi:https://doi.org/10.1016/j.icarus.2004.09.007.
  * Sánchez and Scheeres (2011) Sánchez, P., Scheeres, D.J., 2011\.  SIMULATING ASTEROID RUBBLE PILES WITH a SELF-GRAVITATING SOFT-SPHERE DISTINCT ELEMENT METHOD MODEL.  The Astrophysical Journal 727, 120\.  URL: https://doi.org/10.1088/0004-637x/727/2/120, doi:10.1088/0004-637x/727/2/120.
  * Sandeep et al. (2019) Sandeep, C.S., Marzulli, V., Cafaro, F., Senetakis, K., Pöschel, T., 2019.  Micromechanical behavior of dna-1a lunar regolith simulant in comparison to ottawa sand.  Journal of Geophysical Research: Solid Earth 124, 8077--8100.  URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JB017589, doi:https://doi.org/10.1029/2019JB017589, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2019JB017589.
  * Schaller et al. (2013) Schaller, F.M., Kapfer, S.C., Evans, M.E., Hoffmann, M.J., Aste, T., Saadatfar, M., Mecke, K., Delaney, G.W., Schröder-Turk, G.E., 2013.  Set voronoi diagrams of 3d assemblies of aspherical particles.  Philosophical Magazine 93, 3993--4017.  URL: https://doi.org/10.1080/14786435.2013.834389, doi:10.1080/14786435.2013.834389, arXiv:https://doi.org/10.1080/14786435.2013.834389.
  * Seyyedan et al. (2021) Seyyedan, S.M., Mirghasemi, A.A., Mohammadi, S., 2021.  Numerical simulation of direct shear test on granular materials composed of breakable angular particles: A dem-xfem approach.  Powder Technology 391, 450--466.  URL: https://www.sciencedirect.com/science/article/pii/S0032591021005623, doi:https://doi.org/10.1016/j.powtec.2021.06.038.
  * Suhr and Six (2017) Suhr, B., Six, K., 2017.  Parametrisation of a dem model for railway ballast under different load cases.  Granular Matter 19, 64\.  URL: https://doi.org/10.1007/s10035-017-0740-7, doi:10.1007/s10035-017-0740-7.
  * Sánchez (2015) Sánchez, P., 2015.  Asteroid evolution: Role of geotechnical properties.  Proceedings of the International Astronomical Union 10, 111–121.  doi:10.1017/S1743921315008583.
  * Sánchez et al. (2021) Sánchez, P., Renouf, M., Azéma, E., Mozul, R., Dubois, F., 2021.  A contact dynamics code implementation for the simulation of asteroid evolution and regolith in the asteroid environment.  Icarus 363, 114441\.  URL: https://www.sciencedirect.com/science/article/pii/S0019103521001238, doi:https://doi.org/10.1016/j.icarus.2021.114441.
  * Sánchez and Scheeres (2014) Sánchez, P., Scheeres, D.J., 2014\.  The strength of regolith and rubble pile asteroids.  Meteoritics & Planetary Science 49, 788--811.  URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/maps.12293, doi:https://doi.org/10.1111/maps.12293, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/maps.12293.
  * Thompson et al. (2022) Thompson, A.P., Aktulga, H.M., Berger, R., Bolintineanu, D.S., Brown, W.M., Crozier, P.S., in ’t Veld, P.J., Kohlmeyer, A., Moore, S.G., Nguyen, T.D., Shan, R., Stevens, M.J., Tranchida, J., Trott, C., Plimpton, S.J., 2022.  LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Comp. Phys. Comm. 271, 108171\.  doi:10.1016/j.cpc.2021.108171.
  * Thuillet, Florian et al. (2018) Thuillet, Florian, Michel, Patrick, Maurel, Clara, Ballouz, Ronald-Louis, Zhang, Yun, Richardson, Derek C., Biele, Jens, Tatsumi, Eri, Sugita, Seiji, 2018.  Numerical modeling of lander interaction with a low-gravity asteroid regolith surface - application to mascot on board hayabusa2.  A&A 615, A41.  URL: https://doi.org/10.1051/0004-6361/201832779, doi:10.1051/0004-6361/201832779.
  * Thuillet, Florian et al. (2021) Thuillet, Florian, Zhang, Yun, Michel, Patrick, Biele, Jens, Kameda, Shingo, Sugita, Seiji, Tatsumi, Eri, Schwartz, Stephen R., Ballouz, Ronald-Louis, 2021.  Numerical modeling of lander interaction with a low-gravity asteroid regolith surface - ii. interpreting the successful landing of hayabusa2 mascot.  A&A 648, A56.  URL: https://doi.org/10.1051/0004-6361/201936128, doi:10.1051/0004-6361/201936128.
  * Tsuchiyama et al. (2011) Tsuchiyama, A., Uesugi, M., Matsushima, T., Michikami, T., Kadono, T., Nakamura, T., Uesugi, K., Nakano, T., Sandford, S.A., Noguchi, R., Matsumoto, T., Matsuno, J., Nagano, T., Imai, Y., Takeuchi, A., Suzuki, Y., Ogami, T., Katagiri, J., Ebihara, M., Ireland, T.R., Kitajima, F., Nagao, K., Naraoka, H., Noguchi, T., Okazaki, R., Yurimoto, H., Zolensky, M.E., Mukai, T., Abe, M., Yada, T., Fujimura, A., Yoshikawa, M., Kawaguchi, J., 2011.  Three-dimensional structure of hayabusa samples: Origin and evolution of itokawa regolith.  Science 333, 1125--1128.  URL: https://www.science.org/doi/abs/10.1126/science.1207807, doi:10.1126/science.1207807, arXiv:https://www.science.org/doi/pdf/10.1126/science.1207807.
  * Voronoi (1908) Voronoi, G., 1908.  Nouvelles applications des paramètres continus à la théorie des formes quadratiques. deuxième mémoire. recherches sur les parallélloèdres primitifs.  Journal für die reine und angewandte Mathematik 1908, 198--287.
  * Walsh (2018) Walsh, K.J., 2018.  Rubble pile asteroids.  Annual Review of Astronomy and Astrophysics 56, 593--624.  URL: https://doi.org/10.1146/annurev-astro-081817-052013, doi:10.1146/annurev-astro-081817-052013, arXiv:https://doi.org/10.1146/annurev-astro-081817-052013.
  * Watanabe et al. (2017) Watanabe, S.i., Tsuda, Y., Yoshikawa, M., Tanaka, S., Saiki, T., Nakazawa, S., 2017\.  Hayabusa2 mission overview.  Space Science Reviews 208, 3--16.  URL: https://doi.org/10.1007/s11214-017-0377-1, doi:10.1007/s11214-017-0377-1.
  * Weis and Schröter (2017) Weis, S., Schröter, M., 2017\.  Analyzing x-ray tomographies of granular packings.  Review of Scientific Instruments 88, 051809.  URL: https://doi.org/10.1063/1.4983051, doi:10.1063/1.4983051, arXiv:https://doi.org/10.1063/1.4983051.
  * Wright et al. (2021) Wright, E., Quillen, A.C., Sanchez, P., Schwartz, S.R., Nakajima, M., Askari, H., Miklavcic, P., 2021.  Ricochets on asteroids ii: Sensitivity of laboratory experiments of low velocity grazing impacts on substrate grain size arXiv:2109.06249.
  * Wright et al. (2020) Wright, E., Quillen, A.C., South, J., Nelson, R.C., Sánchez, P., Siu, J., Askari, H., Nakajima, M., Schwartz, S.R., 2020.  Ricochets on asteroids: Experimental study of low velocity grazing impacts into granular media.  Icarus 351, 113963\.  URL: https://www.sciencedirect.com/science/article/pii/S001910352030333X, doi:https://doi.org/10.1016/j.icarus.2020.113963.
  * Xu et al. (2014) Xu, Y., Padding, J., Kuipers, J., 2014.  Numerical investigation of the vertical plunging force of a spherical intruder into a prefluidized granular bed.  Physical review. E, Statistical, nonlinear, and soft matter physics 90, 062203--1/10.
  * Yano et al. (2006) Yano, H., Kubota, T., Miyamoto, H., Okada, T., Scheeres, D., Takagi, Y., Yoshida, K., Abe, M., Abe, S., Barnouin-Jha, O., Fujiwara, A., Hasegawa, S., Hashimoto, T., Ishiguro, M., Kato, M., Kawaguchi, J., Mukai, T., Saito, J., Sasaki, S., Yoshikawa, M., 2006\.  Touchdown of the hayabusa spacecraft at the muses sea on itokawa.  Science 312, 1350--1353.  URL: https://www.science.org/doi/abs/10.1126/science.1126164, doi:10.1126/science.1126164, arXiv:https://www.science.org/doi/pdf/10.1126/science.1126164.
  * Zhang et al. (2021) Zhang, C., Zhao, S., Zhao, J., Zhou, X., 2021.  Three-dimensional voronoi analysis of realistic grain packing: An xct assisted set voronoi tessellation framework.  Powder Technology 379, 251--264.  URL: https://www.sciencedirect.com/science/article/pii/S0032591020310032, doi:https://doi.org/10.1016/j.powtec.2020.10.054.
  * Zhang et al. (2015) Zhang, Y.X., An, X.Z., Zhang, Y.L., 2015.  Multi-particle fem modeling on microscopic behavior of 2d particle compaction.  Applied Physics A 118, 1015--1021.  URL: https://doi.org/10.1007/s00339-014-8861-x, doi:10.1007/s00339-014-8861-x.
  * Zhao et al. (2017) Zhao, S., Zhang, N., Zhou, X., Zhang, L., 2017.  Particle shape effects on fabric of granular random packing.  Powder Technology 310, 175--186.  URL: https://www.sciencedirect.com/science/article/pii/S0032591016309743, doi:https://doi.org/10.1016/j.powtec.2016.12.094.
  * Çelik et al. (2022) Çelik, O., Ballouz, R.L., Scheeres, D.J., Kawakatsu, Y., 2022\.  A numerical simulation approach to the crater-scaling relationships in low-speed impacts under microgravity.  Icarus 377, 114882\.  URL: https://www.sciencedirect.com/science/article/pii/S0019103522000045, doi:https://doi.org/10.1016/j.icarus.2022.114882.

## Appendix I: Verifying the model performance

Table 3: Parameter setup for LAMMPS simulation. *Note: The projectile density is reduced to ensure the same projectile-to-grain mass ratio as for the Abaqus DEM simulation. projectile | projectile | grain | grain | Young’s | friction
---|---|---|---|---|---
diam. ($m$) | density* ($kg/m^{3}$) | diam. ($m$) | density ($kg/m^{3}$) | modulus ($GPa$) | ceoff
0.01615 | 541.7 | 0.002394 | 2600 | 7.0 | 0.84

Finite element method is not often applied to model discrete grains. A primary
difference between the triangular grain and a DEM spherical grain is the
number of nodes per discrete element that must be managed. Spherical DEM only
has one node and a fixed radius, allowing for high optimization of simulation
codes and efficient algorithms for contact and friction controls. In contrast,
the right triangle grains have three nodes and flat sides allowing for more
complex inter-grain contact to arise. Because of this, the contact algorithms
employed by Abaqus must be more versatile. This is examined using two parallel
simulations. The first simulation is of an impact into a bed of spheres,
modeled in Abaqus and thus dependent on Abaqus’ contact algorithm. The second
simulation tests the same scenario, but with the model built in LAMMPS per the
parameters in Table 3. LAMMPS, an open source molecular dynamics software
(Thompson et al., 2022), was selected as it is a widely used discrete granular
systems solver.

In both simulations, the impactor connects with the surface at an angle of 60°
from the horizontal and a velocity of 4 m/s.

Figure 13: Initial and final frames of LAMMPS impact simulation. The particle
representing the projectile is circled in red.

Figure 14: Initial and final frames of Abaqus impact simulation.

Figures 13 and 14 depict the initial and final frames of the impact
simulations carried out in LAMMPS and Abaqus, respectively. Qualitatively, it
is observed that both simulations result in a clear full-stop of the impactor.
The width and depth of both craters are similar and the ‘mounds’ of material
at the leading edges of the craters are of comparable size. Despite the minute
differences in how the LAMMPS and Abaqus models are formulated in terms of
material properties and contact controls, it is clear that both platforms
performed similarly.

In terms of computational efficiency, LAMMPS was found to perform faster than
both Abaqus DEM and DFEM models. Considering bed-settling and a single impact
modeling duration, LAMMPS simulations finished in about 15 minutes on a single
core while Abaqus DEM finished in about 5.5 hours on a single core. In
addition, the DFEM simulations with three times as many nodes finished in
about seven hours on 12 cores. However, Abaqus allows a wider range of
tailoring in model setup and simulation parameters compared to LAMMPS, which
justifies its use despite slower performance.

Figure 15: Plot of kinetic energy versus time for the entire simulation (top)
and zoomed in (bottom).

To strengthen the argument that Abaqus’ contact performs similar to LAMMPS’, a
comparison can also be made using kinetic energy. For each simulation, the
magnitude of the translational velocities of the impactor were exported and
read into a script to calculate kinetic energy throughout the simulation. This
approach allows dispersion of the impactor’s kinetic energy to be plotted
together, shown in Figure 15. The top plot shows the all kinetic energy
calculated through the simulation and the bottom plot shows the zoomed in
region between 0.01s and 0.08s where most kinetic energy dispersion into the
granular bed occurs. First, the slightly varied topography of each simulation
is evident by the deviation of kinetic energy trends. Since the granular beds
in each simulation were generated independently under gravity, the exact
‘pathways’ through which kinetic energy is drained are not identical during
the impacts. However, in both simulations, the kinetic energy settles to zero
at almost the exact same point in between 0.06s and 0.07s. This observation
implies that the average rate by which kinetic energy is drained into the
granular bed is very similar in both cases. Slight discrepancies noticed over
the course of the dispersion are attributed to the slight difference in
material property and contact definition in Abaqus and LAMMPS models. To
conclude, these findings support the performance of Abaqus’ contact algorithm
in the discrete modeling application.

## Appendix II: Plotting stress distribution throughout a bed of angular
grains under surface impact

Abaqus reports the stress experienced within each grain. We use a 3-noded
linear triangular element (CPE3 triangle element) and therefore all gradient
fields including stress and strain are constant within the element. This data
provides a way to look at force propagation throughout the granular bed
underneath an impact. Figure 16 shows this stress contour for several frames
during an impact at 3.5 m/s and 45°. The von Mises stress for each element is
plotted, normalized against the Young’s modulus prescribed to the grains (see
Table 1).

Figure 16a is the first frame of contact between projectile and bed. A clear
pathway of forces is seen from the site of impact to the lower right corner of
the bed at an angle similar to 45°, which corresponds to the angle of the
impact. Since this impact case produces a ricochet, Figure 16d shows the
stress field immediately before the projectile departs the impact site. The
spiderweb of forces has disappeared, though local stresses around the
projectile are still high, likely driving the projectile’s launch from the
surface. Figure 16b and Figure 16c show intermediate stress distributions seen
during impact.

(a) 0.006s, first contact

(b) 0.009s

(c) 0.012s

(d) 0.015s, beginning ricochet

Figure 16: A visualization of stress propagating throughout the granular bed
under a 3.5 m/s, 45° surface impact. Reported values are von Mises stress
normalized against the Young’s modulus of the grain material (7.0 GPa).

Stress distributions around the impact site will be a central component of our
future impacts work. Defining the characteristics of force chain networks and
their correspondence with surface impact dynamics will undoubtedly yield a
deeper understanding of low-velocity impact dynamics.