Structure-dependent fragmentation risk of rubble-pile asteroids on low-impacts

Abstract

Deflection of stray asteroids by low-momentum impacts is an Earth’s key defense strategy. Using a proof-of-principle numerical model, we show that even such impacts on rubble pile asteroids (RPAs) may lead to fragmentation because of inherent internal stress chains. Such chains occur in three- and two-dimensional (2D) aggregates and we study this phenomenon in 2D. We first establish the existence of stress chains in our aggregates and then show that the dynamic shear stress and displacements are highest along the chains, increasing significantly the risk of fracturing there. Our simulations suggest that this phenomenon is independent of the aggregate’s particle size distribution. We conclude that future studies are needed to quantitatively assess fragmentation risk as a function of stress chain statistics in RPAs.

Defense against asteroid impacts is essential to Mankind’s survival. A key strategy is deflection of asteroids into safer trajectories by low-momentum projectile impacts [1]. This is often preferable to blasting them away, which could rain large fragments on Earth and cause more damage than the original object. A key consideration for the deployment of a projectile is informed choice of impact parameters: angle and momentum. In particular, it is important to determine impacts’ effects on cratering and local shattering [2, 3], on the range of damage propagation into the asteroid [4, 5], as well as on potential generation of long fractures. The extent of these effects depends not only on the impact parameters but also on the asteroid’s internal structure. To complicate matters, most asteroids are rubble-pile asteroids (RPAs), consisting of loosely aggregated particles and boulders [walsh2018rubble, jourdan2023rubble]. These fragment more easily than solid objects and the effects of impacts on them require more careful consideration. Crucially, loosely packed particle aggregates are known to transmit stresses non-uniformly as force chain networks [8, 9], whose characteristics depend on the aggregate’s internal structure. The more porous the aggregate the more non-uniform the stress field [10] and the more fragile the asteroid.

Discrete element method (DEM) studies of impact-induced seismic waves in RPAs have been useful to gain insight into such phenomena [5, 11]. However, previous studies have ignored the effects of non-uniform stress chains. Since the material density is higher along such chains than around them the response to impacts should propagate faster along such chains. This is expected to result in elevated shear stress and deformations, which should increase the risk of fractures and fragmentation. We test this expectation and show that stress chains indeed introduce such vulnerability. This suggests that structure-based fragmentation must be taken into consideration in RPA-deflection strategies.

Stress in granular aggregates is transmitted by stress chains in two and three dimensions. It is then instrumental to study this issue, as a proof of concept, in two dimensions (2D). Our 2D toy model illustrates clearly that stress chains modify significantly the response to even low-momentum impacts and increase the fracturing risk along them.

Our 2D models comprise rubble-pile clusters (RPCs) of particles aggregated under gravity. We used different size distributions for the aggregation, as detailed below. We then quantified their internal structures and stress fields, verifying that they exhibit stress chains. Next, we simulated impacts at, and away from, stress chains and studied the effects of the initial stress configurations and the dynamic stress and displacement responses.

For our numerical simulations, we used the in-house discrete-element code DEMBody [12, 13, 14]. It features nonlinear Hertz–Mindlin–Deresiewicz interparticle contact forces [15], determined by particle overlaps (typically less than $1\%$ particle sizes). The interactions include normal and tangential components, rolling and twisting torques, and cohesion [12]. The system’s planarity obviates interparticle twisting torques. We also neglected cohesion. Particle parameters, given in Table 1, were chosen to mimic the mechanical properties of real 3D RPA regolith particles and boulders. The shape parameter, $\beta$, is a statistical measure that models particle non-sphericity and controls the rotational resistance to relative motions in the rolling direction [16].

Table 1: Particle characteristic and contact parameters in our simulations.

Parameter	Symbol	Value
Young modulus	$Y$	$5\times 10^{9}$ Pa
Poisson’s ratio	$\nu$	$0.3$
Static friction coefficient	$\mu_{\mathrm{s}}$	$1.0$
Dynamic friction coefficient	$\mu_{\mathrm{d}}$	$1.0$
Restitution coefficient	$\varepsilon$	$0.4$
Shape parameter	$\beta$	$0.8$
Material density	$\rho$	$3.4\times 10^{3}$ kg/m3

The initial rubble-pile clusters (RPCs) were constructed by placing a nuclear four-particle aggregate at the origin and surrounding it with a cloud of $N\approx 3000$ particles inside an annulus around it (Fig. 2 in the supplemental material (SM) [17]). The particles aggregated under gravity as they ‘fell’ toward the nucleus. Particle motion was damped by a standard viscous term to dissipate the kinetic energy on the approach to a static state. The RPC was considered static when the total kinetic energy fell below $10^{-6}$ of the total initial gravitational potential energy. The density and porosity of the final RPC were controlled by the initial density of the particle cloud in the annulus and the particle size distribution (PSD). Mechanical stability was ascertained by checking that the bulk mean coordination numbers, excluding rattlers, were slightly larger than 3 [18]. Most of our investigations were carried out using an equal-numbers bi-disperse PSD: $d_{1}=5$ m and $d_{2}=10$ m, with and without few surface boulders. We ran further simulations with a power-law PSD, $P(d)\sim d^{-3}$, often observed in RPAs [19, 20, 21], which corroborated our conclusions.

The stress field in the pre-impact RPC is slightly off-symmetric because the particles retain small residual translational and rotational velocities, as well as residual contact torques, also known as rolling torques, caused by local elastic deformations. Particle stress fields were determined using [22, 23],

$$\sigma_{\alpha\beta,p}=\frac{1}{A_{p}}\sum_{p^{\prime}\in p}f_{p^{\prime}p,\alpha}r_{p^{\prime}p,\beta}+\frac{J_{p}}{A_{p}}\left(\dot{\omega}_{p}\varepsilon_{\alpha\beta 3}+\omega_{p}^{2}\delta_{\alpha\beta}\right)\ ,$$ (1)

where $p^{\prime}$ are particles contacting $p$, $f_{p^{\prime}p,\alpha}$ the $\alpha$-component of the force that $p^{\prime}$ exerts on $p$, $r_{p^{\prime}p,\beta}$ the $\beta$-component of the position vector to the contact between $p^{\prime}$ and $p$ (see Fig. 1 in the SM [17]), $J_{p}$ the moment of inertia of $p$, $\omega_{p}$ its angular velocity, $J_{p}\dot{\omega}_{p}$ the total torque moment on $p$, $\varepsilon_{\alpha\beta 3}$ the component of the Levi-Civita tensor, $\delta_{\alpha\beta}$ Kronecker’s delta, and $A_{p}$ the area associated with particle $p$. We determine $A_{p}$ using the quadron method [18, 24], reviewed in the SM [17]. The negative stress trace reflects the compressive interparticle forces. Using these local measures, we visualized the traces of particle stress tensor and identified stress chains to inform the following impact simulations.

The impact simulations on the static RPCs were designed to elucidate the effects of stress chains on the dynamic response. An impact is a velocity pulse to a particle or a surface boulder. The pulse magnitudes ranged from $1\times 10^{-3}$ m/s to $10$ m/s, corresponding to momenta from $1.78\times 10^{3}$ kg$\cdot$m/s to $1.78\times 10^{7}$ kg$\cdot$m/s. The low limit was chosen to enable observations of preferred propagation directions of stress and displacements, which get washed out by extended damage on more violent impacts. Impact locations were chosen judiciously to isolate the effects of stress chains and impact momenta.

For generality, we use non-dimensional lengths, time, velocities, and stress components, scaled by, respectively:

$$l_{0}\equiv\bar{d}\quad;\quad t_{0}\equiv\sqrt{\frac{l_{0}}{\bar{g}}}\quad;\quad v_{0}\equiv\frac{l_{0}}{t_{0}}\quad;\quad\sigma_{0}\equiv\frac{\bar{m}}{t_{0}^{2}}\ .$$ (2)

Thus, $t\to\tau=t/t_{0}$, $\bm{v}\to\bm{u}=\bm{v}/v_{0}$, $\sigma\to\tilde{\sigma}=\sigma/\sigma_{0}$, and lengths are in units of $l_{0}$. Here, $\bar{d}$ is the mean particle diameter, $\bar{g}=(\mathrm{G}/N)\sum_{i}\sum_{j\neq i}m_{j}/\left|\bm{r}_{j}-\bm{r}_{i}\right|^{2}$ is the pre-impact gravitational acceleration magnitude, averaged over all the RPC particles, and $\bar{m}$ is the mean particle mass. The size of the $N$-particle RPC is defined as its radius of gyration,

$$R_{\mathrm{g}}=\sqrt{\frac{\sum_{i=1}^{N}m_{i}r_{i}^{2}}{\sum_{i=1}^{N}m_{i}}}\ .$$ (3)

The time-dependent mean front of the propagating displacement response is

$$D=\frac{\sum_{i=1}^{N}H\left(\left|\bm{u}_{i}\right|-\left|\bm{u}_{\mathrm{th}}\right|\right)\left|\bm{r}_{i}-\bm{r}_{\mathrm{impact}}\right|}{R_{\mathrm{g}}\sum_{i=1}^{N}H\left(\left|\bm{u}_{i}\right|-\left|\bm{u}_{\mathrm{th}}\right|\right)}\ .$$ (4)

$H(x)$ is the Heaviside step function, $\left|\bm{u}_{\mathrm{th}}\right|=\alpha\left|\bm{u}_{\mathrm{impact}}\right|$ a velocity threshold, $\bm{u}_{i}$ the velocity of particle $i$, and $\bm{r}_{\mathrm{impact}}$ and $\bm{r}_{i}$ the position vectors of the impacted particle and particle $i$ respectively, with the origin at the RPC’s center of mass. $D$ depends somewhat on the choice of $\alpha$, but the results and conclusions are independent of it.

The PSDs used for the impact simulations were: bi-disperse (Fig. 1), bi-disperse and larger surface boulders (Fig. 3), and power-law (Fig. 4). In all the simulations, we first verified that gravity alone generates networks of stress chains and identified the strong chains before simulating impacts.

Real impacts take place at RPAs’ surfaces and we simulated such impacts. However, since our main aim was to study the effect of stress chains on damage pattern and propagation, we also mimicked internal impacts on stress chains. For thoroughness, those impacts were applied in different directions relative to the chain orientations. For comparison, we applied the same impact magnitudes to particles in weak-stress regions. The post-impact dynamics data were used to determine the correlations between the initial stress chains and the response propagation patterns and ranges.

Bi-disperse PSD: A typical example of the gravity-induced non-uniform stress distribution in a bi-disperse system is shown in Fig. 1(a). The occurrence probability of $\left|\mathrm{Tr}(\tilde{\sigma}_{0})\right|$ decays exponentially, as shown in Fig. 3 in the SM [17], and the probability that $\left|\mathrm{Tr}(\tilde{\sigma}_{0})\right|>100$ is $0.0073$. Impacts were applied to large ($M=1.78\bar{m}$) internal and surface particles in judiciously chosen locations. All impact magnitudes were $8.90\times 10^{3}$ kg$\cdot$m/s $=1.78\bar{m}v_{0}$, with $\left|\bm{u}_{\mathrm{impact}}\!\right|=1$. In Fig. 1, we show typical results of impact on stress chain-resident particle I. The impacts were applied in three angles to the chain direction: $0^{\circ}$, $60^{\circ}$, and $90^{\circ}$. We applied the same impact to particle II, residing in a weak-stress region at the surface. These impact simulations are labeled $a-d$, respectively.

Figure 1: An example of the bi-disperse packing and the impact tests. (a) A color map of the pre-impact Tr($\tilde{\sigma}_{0}$). Several impacts, labeled $a-d$ (yellow arrows), are applied at I and II (yellow circles). (b) A zoom on Tr($\tilde{\sigma}_{0}$) around location I, with the region around I divided into 8 sectors. (c) The particle displacement rate (blue-red) at $\tau=6.63\times 10^{-4}$, superposed on Tr($\tilde{\sigma}_{0}$) in grey. Particles with $\left|\bm{u}\right|\lesssim 0.01$ are invisible. (d) The bi-variate probability scatter map in the $\epsilon-\zeta$ plane at $\tau=6.63\times 10^{-4}$.

The effects of the pre-impact stress structure on the stress and displacement rate responses are illustrated in Figs. 1(b)-(d) and in Fig. 4 in the SM [17]. We present results for case $b$, in which the impact momentum is at $60^{\circ}$ to the chain orientation, but the results are typical across all our simulations. In Fig. 1(b), we show the particles carrying strong pre-impact stress chains ($\left|{\rm Tr}\left(\tilde{\sigma}_{0}\right)\right|\gtrsim 20$) and, to study the local response, we divide the region around the impacted particle I into eight sectors. In Fig. 1(c), we show the local particle displacement rate at $\tau=6.63\times 10^{-4}$ after the impact, superimposed on the pre-impact stress chains. Focusing on high stress chains, Fig. 1(d) shows the bi-variate probability of the normalized values of $\epsilon\equiv\left|{\rm Tr}\left(\tilde{\sigma}_{0}\right)\right|-\left|{\rm Tr}\left(\tilde{\sigma}_{0}\right)\right|_{\rm med}$ and $\zeta\equiv\left|\tilde{\sigma}_{12}\right|-\left|\tilde{\sigma}_{12}\right|_{\rm med}$, with $(.)_{\rm med}$ denoting the median value.

This figure, as well as the evidence presented in Fig. 4 in the SM [17], illustrates several points. Firstly, even though the impact is at $60^{\circ}$ to the chain, the displacement rate is more pronounced along it than in the impact direction. Secondly, the stress propagates faster along the chain than in the impact direction. Thirdly, the post-impact dynamic shear stress is also highest along the initial stress chains. Since under compression only mode II fracture is possible in 2D [25] then it is the local shear stress magnitude that could open a fracture. This causal relation between shear stress and fracture risk makes the latter observation significant. The supplemental video “bidisp_abcd.mp4” also supports these conclusions.

In Fig. 2, we show the dynamic response front, $D$, in the different sectors and its dependence on the angle between the impact momentum and the chain direction. Particle I is on a stress chain that extends into sector 5, then curves into sector 4, and continues into sector 3. Figure 2(a) confirms that $D$ indeed propagates fastest into sector 5. At $\tau\approx 5\times 10^{-4}$, $D$ increases dramatically in sector 4, corresponding to the moment when the stress response, carried by the chain, entered sector 4. In contrast, $D$ remains low in sectors 6, 7, and 8 until this time, owing to initial low stress there, and in sector 1, which is opposite to the impact direction. Later, $D$ increases in sector 6 and then in 7, reflecting the eventual arrival of the stress response to these sectors along the chain path rather than directly through the medium. Figure 2(b) also shows that the stress response depends on the angle between the impact and initial chain directions – it is strongest when they are aligned (case $a$) and weakest when they are perpendicular (case $c$).

Figure 2: The dynamic response. (a) The evolution of $D$ in the eight sectors around I in impact test $b$, with $\left|\bm{u}_{\mathrm{th}}\right|=0.01$. (b) The evolution of the total value of $D$ in impact tests $a-d$.

The jumps in Fig. 2(a) make it possible to deduce the response propagation speed along the initial chains and potentially relate this speed to the initial stress, but this calculation is outside the scope of this study.

Consistent with the above, we see in Fig. 2(b) that the initial response ($\tau\lesssim 2\times 10^{-3}$) to an impact in the weak-stress region and also on a surface particle (case $d$) propagates much more slowly than in case $a$ and somewhat slower than in case $b$. Its speed is comparable to that in case $c$, in which the impact direction is also into a weak-stress region. The values of $D$ satisfy $a>b>c\approx d$ at all times. Past the maxima, $D$ in case $d$ is slightly lower than in $c$. In Fig. 5 in the SM [17], we show the displacement rate, $\left|\bm{u}\right|$, for these four cases at $\tau=6.63\times 10^{-3}$, superposed on the initial stress chains. A visual comparison of the dynamics following impacts $a-d$ is available in the supplemental video “bidisp_abcd.mp4”.

Bi-disperse PSD with surface boulders: The bi-disperse RPCs are useful for understanding the effects of stress chains, but real impacts take place at the surface. While the near-surface of a bi-disperse aggregate can only support weak stress due to gravity, surfaces of real RPAs are scattered with massive boulders. These can give rise to stress chains whose strengths depend on the boulder’s and RPA’s masses [14]. To assess the risk posed by these chains, we added to the static bi-disperse RPCs three surface boulders of diameters $d=20,40,50$ m. The boulders were released with zero velocity from a distance comparable to the RPC’s $R_{\mathrm{g}}$. In the example we show, strong stress chains formed only under the two more massive boulders (Fig. 6(a) in the SM [17]). We ran several simulations of impacts applied to these boulders and show the (typical) results of an impact applied to the $40$ m boulder ($M=85.57\bar{m}$). The impact’s momentum was $2.26\times 10^{4}$ kg$\cdot$m/s $\approx 4\bar{m}v_{0}$, with $\left|\bm{u}_{\mathrm{impact}}\right|=0.047$, and it was directed toward the RPC’s center of mass. In Fig. 3, we show the particle displacement rate (blue-red) at $\tau=1.71\times 10^{-3}$ after impact, superposed on the initial stress chains (grey-scale). We verified that, in all these simulations, the displacement rate is highest along the strong stress chains, bypassing the weak-stress regions (e.g., Fig. 3(a)), and so is the shear stress response. The latter is clearly observed in Fig. 3(b) and this conclusion is further supported by Fig. 6(b) in the SM [17] and by the supplemental video “boulder_Tr0S12switch_1010.mp4”.

Figure 3: (a) Particle displacement rate (blue-red) at $\tau=1.71\times 10^{-3}$, superposed on Tr($\tilde{\sigma}_{0}$). (b) The bi-variate probability scatter map in the $\epsilon-\zeta$ plane at $\tau=1.71\times 10^{-3}$.

Power-law PSD: To ascertain that our results are not limited to bi-disperse PSDs, we ran similar simulations with a power-law PSD, $P(d)\sim d^{-3}$. The power was chosen, following observations in real RPAs [19, 20, 21]. We chose $5\leq d\leq 15$ m to ensure the same RPC mass as in the bi-disperse case for a fair comparison. Fig. 4 depicts an example of the response to an impact on a particle ($M=1.78\bar{m}$) in a stress chain. The impact momentum was $8.93\times 10^{3}$ kg$\cdot$m/s $=1.78\bar{m}v_{0}$ ($\left|\bm{u}_{\rm impact}\right|=1$) and it was directed along the chain and perpendicular to a nearby one (Fig. 7(a) in the SM [17]). In all the simulations of these systems, an example of which is shown in Fig. 4(a), we observed the same enhanced propagation of the dynamic response in the chain direction. Superposing the initial stress chains, the displacement rates at $\tau=1.00\times 10^{-3}$ clearly illustrate this phenomenon. In places, the high displacement rate circumvents regions of weak stress and is even perpendicular to the impact direction as it follows the chain direction. Fig. 4(b) also shows that, by and large, high post-impact shear stress locations coincide with pre-impact strong stress chains. Further support is provided by the product Tr($\tilde{\sigma}_{0}$)$\cdot$ $\left|\tilde{\sigma}_{12}\right|$, shown in Fig. 7(b) in the SM [17], which is highest along the pre-impact chains.

Figure 4: Superposition of Tr($\tilde{\sigma}_{0}$) (grey-scale) and particle displacement rate (blue-red) at $\tau=1.00\times 10^{-3}$. (b) The bi-variate probability scatter map of particles in the parameter space formed by $\epsilon$ and $\zeta$ at $\tau=1.00\times 10^{-3}$.

We ran further simulations of RPCs with the same power-law PSD and the same particle size range as the bi-disperse range, which are unavoidably less massive. Those results are shown in Fig. 8 in the SM [17] and also support our conclusions that, along stress chains, the response propagates further and faster and the shear stress is generically highest.

To conclude, we set out to determine the effects of stress chains on the response of RPAs to low-momentum impacts, aimed only to change their trajectories. For proof of principle, we constructed a 2D model of random pile clusters (RPCs) and showed that pre-impact stress chains modify significantly the internal dynamic responses to impact. Specifically, because the dynamic shear stress and the displacement responses propagate faster and further along the denser stress chains, these fields are also locally highest there. This enhancement weakens as the relative angle between the impact direction and the orientation of the impacted stress chain increases. We have verified that these effects occur in both bi-disperse and power-law PSDs, suggesting that, at least within the particle size ranges we studied, they are not sensitive to the PSD.

Since fracturing is sensitive to shear stress and displacements, our results illustrate that stress chains increase significantly the risk of fragmentation along pre-existing chains even at low impact momenta. This risk must be taken into consideration when attempting to merely deflect RPAs. Ways to assess this risk should include investigations into the relations between the statistical characteristics of the stress chain network (strength and spatial distributions) and the magnitudes of the dynamic shear stress and displacements. Such studies are downstream of this work and undergoing investigations in this direction will be reported eventually.

R.B. and C.H. conceived the idea for the study; C.H., R.B., P.R.K., and Y.Y. developed the idea; C.H. conducted the simulations and the numerical analysis; C.H. and R.B. wrote the paper; P.R.K and Y.Y. reviewed and revised the paper.