Rotational Failure of Rubble-pile Bodies: Influences of Shear and Cohesive Strengths

In a granular assembly, each particle can interact with its surroundings through short-range interactions (e.g., mechanical contact forces, van der Waals forces) and long-range interactions (e.g., gravity, electrostatic forces).⁷ Within the high-efficiency parallel tree code framework, pkdgrav (Richardson et al. 2000; Stadel 2001), a soft-sphere model including four components in the normal, tangential, rolling, and twisting directions is used for computing particle contact forces (Schwartz et al. 2012; Zhang et al. 2017). Figure 1 presents the directions of the forces and torques acting on particle i generated by the contact with particle j. The linear spring-dashpot normal contact force ${{\boldsymbol{F}}}_{N}$ and the tangential stick-slip force ${{\boldsymbol{F}}}_{S}$ are given by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{N}=-{k}_{N}{\delta }_{N}\hat{{\boldsymbol{n}}}+{C}_{N}{{\boldsymbol{u}}}_{N}\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{S}=\min ({k}_{S}{{\boldsymbol{\delta }}}_{S}+{C}_{S}{{\boldsymbol{u}}}_{S},\,{\mu }_{S}| {{\boldsymbol{F}}}_{N}| {{\boldsymbol{\delta }}}_{S}/| {{\boldsymbol{\delta }}}_{S}| ),\end{eqnarray} \tag{ 1b }$$

which depend on the stiffness constants, k_N and k_S, the viscous damping coefficients, C_N and C_S, and the sliding friction coefficient μ_S.⁸ The variable δ_N is the mutual compression between two particles, and ${{\boldsymbol{\delta }}}_{S}$ is the sliding displacement from the equilibrium contact point. The unit vector $\hat{{\boldsymbol{n}}}$ gives the direction from the center of particle i to that of particle j. The dashpot force is linearly proportional to the normal and tangential relative velocities, ${{\boldsymbol{u}}}_{N}$ and ${{\boldsymbol{u}}}_{S}$. Recently, in order to capture behaviors of granular materials under quasi-static loadings, a spring-dashpot-slider rotational resistance model was implemented in pkdgrav (Zhang et al. 2017). The resulting rolling and twisting torques, ${{\boldsymbol{M}}}_{R}$ and ${{\boldsymbol{M}}}_{T}$, are given by

$$\begin{eqnarray}&&{{\boldsymbol{M}}}_{R}=\left\{\begin{array}{l}{k}_{R}{{\boldsymbol{\delta }}}_{R}+{C}_{R}{{\boldsymbol{\omega }}}_{R},\quad | {k}_{R}{{\boldsymbol{\delta }}}_{R}| \lt {M}_{R,\max }\\ {M}_{R,\max }{{\boldsymbol{\delta }}}_{R}/| {{\boldsymbol{\delta }}}_{R}| ,\quad | {k}_{R}{{\boldsymbol{\delta }}}_{R}| \geqslant {M}_{R,\max }\end{array}\right.\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&{{\boldsymbol{M}}}_{T}=\left\{\begin{array}{l}{k}_{T}{{\boldsymbol{\delta }}}_{T}+{C}_{T}{{\boldsymbol{\omega }}}_{T},\quad | {k}_{T}{{\boldsymbol{\delta }}}_{T}| \lt {M}_{T,\max }\\ {M}_{T,\max }{{\boldsymbol{\delta }}}_{T}/| {{\boldsymbol{\delta }}}_{T}| ,\quad | {k}_{T}{{\boldsymbol{\delta }}}_{T}| \geqslant {M}_{T,\max },\end{array}\right.\end{eqnarray} \tag{ 2b }$$

where k_R and k_T are the material stiffness along rolling and twisting directions, and C_R and C_T are the damping coefficients for rolling and twisting motions. ${{\boldsymbol{\omega }}}_{R}$ and ${{\boldsymbol{\omega }}}_{T}$ describe the relative rolling and twisting motions between two particles in contact. Similar to the tangential force model, ${{\boldsymbol{\delta }}}_{R}$ and ${{\boldsymbol{\delta }}}_{T}$ are the rolling and twisting angular displacements from the equilibrium contact point, respectively. When the magnitude of the rolling/twisting torque reaches a critical value, M_R,max or M_T,max, the surfaces of the two particles start to slide relative to each other along the rolling/twisting direction, whereupon the system enters into a dynamic friction condition. The motion of a particle can be obtained by integrating the total force and torque acting on it.

Zoom In Zoom Out Reset image size

The interactions between particles and wall boundaries with a wide range of geometries as well as the effect of gravity are also implemented in pkdgrav, allowing precise reproduction of laboratory experiments (see Schwartz et al. 2012 for details). The numerical approach has been validated through comparison with laboratory experiments, e.g., low-speed impact experiments (Schwartz et al. 2014), and has been successfully used to study various behaviors of granular systems, e.g., the Brazil-nut effect (Matsumura et al. 2014; Perera et al. 2016; Maurel et al. 2017) and avalanches (Yu et al. 2014). In addition, the parallel hierarchical tree code in pkdgrav provides a highly efficient way to carry out long-range interaction calculations for an N-body system, giving an ${ \mathcal O }(N\mathrm{log}N)$ scaling dependence, where N is the number of particles in a simulation (Richardson et al. 2000). With the current processing power, pkdgrav has the ability to simulate the evolution of tens of thousands to millions of particles; this has been well tested in a wide range of environments (e.g., Ballouz et al. 2015, 2017; Zhang et al. 2017; Schwartz et al. 2018).

Several cohesion models have been implemented previously in pkdgrav, but for a different cohesive mechanism, i.e., sintering forces, which can build very strong bonds giving rise to the formation of particle aggregates (see Richardson et al. 2009; Schwartz et al. 2013). Although the thermal processes on asteroids may facilitate the formation of sintering bonds between grains at the asteroids' surface, the van der Waals force is a more universal mechanism for bodies in a low-gravity, ultra-high-vacuum environment (Scheeres et al. 2010), which has been verified for lunar regolith (Mitchell et al. 1972).

The cohesive force introduced in this study is a combination of the granular bridge idea proposed by Sánchez & Scheeres (2014, 2016) and a hypothetical contact area characterized by a shape parameter, β, first introduced in the contact model of Jiang et al. (2013, 2015). With the assumption that there exist substantial micro-sized grains to cover larger boulders, this model can mimic the cumulative effect of the cohesive regolith between two large boulders without simulating each individual fine grain. The rotational resistance arising from the interstitial regolith is also taken into account via the hypothetical contact area.

As shown in Figure 1, taking particles i and j as boulder representatives and assuming that there are many fine particles distributed around the gaps between these two boulders to form a granular bridge, we can obtain the cumulative forces and torques transmitted through the granular bridge by integration over an effective contact area. The orange polygon shown in Figure 1 can serve as the effective contact area that contains all the points that are at the same distance from the two boulders' surfaces and are closer than these surfaces to any other boulders' surfaces. The area of this polygon can be accurately calculated using the Voronoi tessellation per pair of neighboring boulders in principle. To avoid this burdensome calculation, we introduce the shape parameter, β, to represent a statistical measure of the area where the interstitial grains are important and approximate this effective contact area using a rectangle:

$$\begin{eqnarray}&&{A}_{\mathrm{eff}}=4{(\beta R)}^{2},\end{eqnarray} \tag{ 3 }$$

where the effective radius R = r_ir_j/(r_i + r_j), and r_i and r_j are the radii of the corresponding boulders. The resulting cohesive force between the two boulders can be modeled as

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{C}={{cA}}_{\mathrm{eff}}\hat{{\boldsymbol{n}}},\end{eqnarray} \tag{ 4 }$$

where the interparticle cohesive tensile strength c reflects the physical properties of the interstitial regolith, e.g., porosity, size distribution, surface energy, and Hamaker constants (Seizinger et al. 2013; Sánchez & Scheeres 2014; Gundlach & Blum 2015).

The interstitial cohesive regolith will also strengthen the friction between boulders to resist relative rotational motion, commensurate with the shape effect (Li et al. 2011; Jiang et al. 2013). With a similar contact area characterized by the shape parameter β, the rolling/twisting stiffness and damping coefficients can be expressed as

$$\begin{eqnarray}{k}_{R} & = & {k}_{N}{(\beta R)}^{2},\quad {C}_{R}={C}_{N}{(\beta R)}^{2},\\ {k}_{T} & = & 2{k}_{S}{(\beta R)}^{2},\quad {C}_{T}=2{C}_{S}{(\beta R)}^{2}.\end{eqnarray} \tag{ 5 }$$

When the relative rolling/twisting rotational motion exceeds a critical threshold, the cohesive contact on the edge of the contact area begins to break, and the peak rolling/twisting torque is given as

$$\begin{eqnarray}{M}_{R,\max } & = & {\mu }_{R}\beta R| {{\boldsymbol{F}}}_{N}| ,\\ {M}_{T,\max } & = & {\mu }_{T}\beta R{\mu }_{S}| {{\boldsymbol{F}}}_{N}| \end{eqnarray} \tag{ 6 }$$

where μ_R and μ_T represent the static friction coefficients for rolling and twisting, respectively.⁹

In our implementation, the cohesive force described in Equation (4) is added to the calculation of the interaction between two particles in contact, and is ignored when the two particles are separated. By these means, the breakage and recovery of the cohesive effect of interstitial regolith can be captured. Although the behavior of individual interstitial grains is not physically simulated, our model can capture the resulting tensile strength and rotational resistances well, which is sufficient for our purpose of demonstrating the effect of cohesive strength on rubble-pile asteroids. In addition, since the effective contact area A_eff given by Equation (3) could be much larger than the actual contact area between two spherical particles (i.e., the dashed yellow circle shown in Figure 1), the resulting rotational resistance and cohesion can reduce excess particle rotations and restrict the particles' freedom to rotate. Therefore, the effects of surface roughness and the interlocking of non-spherical shapes found for real-world granular materials can be reflected to some extent in our model. In Section 3, a series of drained triaxial compression tests were carried out to validate the ability of our SSDEM model to simulate cohesive granular materials.

Before proceeding to lay out numerical simulations, it is necessary to clarify how we set the parameters used in our SSDEM model, including two stiffness constants, k_N and k_S, two damping coefficients, C_N and C_S, three friction coefficients, (μ_S, μ_R, μ_T), a shape parameter β, an interparticle cohesive tensile strength c, and an integration timestep Δt. The normal stiffness k_N and the timestep Δt are set to ensure that particle overlaps for each entire simulation do not exceed 1% of the minimum particle radius (in cases where it is important to match the sound speed of real materials, k_N can also be chosen to control the speed of energy propagation through a medium represented by soft-sphere particles; see Section 2.1 in Schwartz et al. 2012 for details). To keep normal and tangential oscillation frequencies equal to each other, the tangential stiffness k_S is often set to (2/7)k_N (see Section 2.6 in Schwartz et al. 2012, for details). C_N and C_S can be derived from the normal and tangential coefficients of restitution, ε_N and ε_S, which are set to 0.55 so that the granular system is subject to sufficient damping (these values are also reasonable for terrestrial rocks; Chau et al. 2002). μ_R and μ_T represent the capability of the interstitial regolith to resist breakage, which are set to 1.3 and 1.05 for a medium level, respectively (the derivation is similar to the cohesionless case; we refer the reader to Jiang et al. 2015, for details).

By these means, the dimension of the SSDEM parameter space is narrowed down to three, i.e., μ_S, β, and c, which control the sliding resistance strength, the rotational resistance strength, and the cohesive strength between particles in contact, respectively. In SSDEM simulations, these three parameters can serve as the adjustable material parameters, allowing the modeling of granular materials with various material properties, as revealed in the following section.

Figure 2 illustrates the simulation setup. The triaxial compression procedure is automatically performed using pkdgrav. In the beginning, a uniform confining pressure is applied to the specimen and then the vertical axial strain ε₁ is increased linearly in compression while the motions of the lateral walls are controlled by a servo mechanism that can maintain a constant confining pressure continuously acting on the specimen, i.e., σ₂ = σ₃ = constant.¹⁰ To meet quasi-static loading conditions, the top and bottom walls are moved toward each other at a constant strain rate of 3% per minute.

Zoom In Zoom Out Reset image size

The granular specimen used in our triaxial tests consists of ∼14,000 particles, with a −3-index power-law distribution in radius ranging from 6 mm to 18 mm and a material density of 3.0 g/cc. A two-stage procedure is used to prepare the specimen for representing the macroscopic homogeneity of a real sand sample (the same procedure is also applied to generate rubble-pile asteroid representatives in Section 4).¹¹ First, a granular assembly is created by placing particles with a predefined size distribution randomly in a spherical cloud and allowing the cloud to gravitationally collapse with highly inelastic collisions. Frictionless and cohesionless SSDEM parameters are used to minimize the void space between particles. Second, a portion with a predefined shape is extracted from the interior of the granular assembly produced in the first step.

The prepared specimen is isotropically compressed to a uniform confining pressure condition in a 340 mm × 340 mm × 720 mm box with six frictionless boundary walls (see Figure 2). Given that the particle arrangement is sensitive to the external gravity field, leading to a biased direction of the contact forces, compression processes (including the isotropic and triaxial) are conducted in a gravity-free condition. Figure 3 presents the distribution of contact orientations prior to vertical loading, showing that this specimen is virtually homogeneous in all directions. The packing efficiency of our tested sample is ∼0.62, which is within the typical range of terrestrial sands observed in laboratory tests (Cabalar et al. 2013).

Zoom In Zoom Out Reset image size

To reproduce the desired elastic mechanical property and match the sound speed of real sand, a normal stiffness constant k_N of 5.1 × 10⁶ kg s⁻² is chosen (a stiffer k_N normally leads to a higher Young's modulus), and the integration timestep Δt is set to 3.6 × 10⁻⁶ s so that particle overlaps for each entire simulation do not exceed 1% of the minimum particle radius. The effect of μ_S, β, and c is explored in Section 3.3.

The material strengths of a granular material can be derived from its failure stress state. As one of the most common yield criteria for geophysical granular materials, the Drucker–Prager failure criterion has been successfully applied to analyzing the failure of rubble-pile structures (e.g., Holsapple 2007; Sharma et al. 2009; Hirabayashi 2015), which is given as (Jaeger et al. 2009)

$$\begin{eqnarray}&&\sqrt{{J}_{2}}\leqslant C+{{sI}}_{1},\end{eqnarray} \tag{ 7 }$$

where I₁ is the first invariant of the Cauchy stress tensor, and J₂ is the second invariant of the deviatoric stress tensor, which can be written in terms of the principal stresses, σ₁, σ₂, and σ₃ (from largest to smallest), as I₁ = σ₁+σ₂+σ₃, ${J}_{2}\,=[{({\sigma }_{1}-{\sigma }_{2})}^{2}+{({\sigma }_{2}-{\sigma }_{3})}^{2}+{({\sigma }_{3}-{\sigma }_{1})}^{2}]/6$.

The failure envelope (Equation (7) with the equality) plotted in the I₁–$\sqrt{{J}_{2}}$ plane is a straight line with a slope of s and an intercept of C, in which s and C represent the shear and cohesive strengths of the granular material, respectively, as shown in Figure 4. Note that this cohesive strength C represents the bulk property of the granular material, which is different from the interparticle cohesive tensile strength c used in the cohesive force model (Equation (4)). We distinguish this large C as the macroscopic cohesive strength. Since the Drucker–Prager failure envelope is a smooth version of the Mohr-Coulomb failure envelope, the slope s can be expressed as a function of the commonly used friction angle ϕ, i.e., $s=2\sin \phi /[\sqrt{3}(3-\sin \phi )]$.

Zoom In Zoom Out Reset image size

The constants ϕ and C can be determined from results of triaxial compression tests. Generally, three or more tests should be performed to make an appropriate estimation. As shown in Figure 4, the stress states at failure (i.e., the state at which the deviator stress achieves its peak in triaxial compression tests) under different confining pressure are plotted in the I₁–$\sqrt{{J}_{2}}$ plane, and a best-fit straight line (i.e., the failure envelope) can be derived from these data points. Then the values of ϕ and C can be extracted from the failure envelope.

The macroscopic mechanical properties of a specimen are characterized by four parameters, i.e., the friction angle ϕ, the macroscopic cohesive strength C, the Young's modulus E, and the Poisson's ratio ν. ϕ and C can be estimated using the method introduced in Section 3.2, and E and ν can be obtained from the initial slopes in the linear-elastic regime of the stress-strain curve and the volumetric strain curve, respectively (Belheine et al. 2009),

$$\begin{eqnarray}&&E=\displaystyle \frac{{\rm{\Delta }}{\sigma }_{1}}{{\rm{\Delta }}{\varepsilon }_{1}\,},\quad \nu =\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{{\rm{\Delta }}{\varepsilon }_{v}}{{\rm{\Delta }}{\varepsilon }_{1}}\right).\end{eqnarray} \tag{ 8 }$$

We investigated the effects of the SSDEM parameters, μ_S, β, and c, on these macroscopic properties. To draw the failure envelope, triaxial tests for each parameter set were conducted under three confining pressures, i.e., σ₂ = σ₃ = 50 kPa, 100 kPa, and 200 kPa, respectively. As shown in Table 1, the value ranges of the Young's modulus E and the Poisson's ratio ν and their dependences on the confining pressure are consistent with laboratory experiments on real sands (e.g., for the Labenne medium sand, E ∼ 100 MPa, ν ∼ 0.27, and ϕ ∼ 36°; see Kolymbas & Wu 1990; Mestat & Berthelon 2001), indicating the validity of our SSDEM model.

3.3.1. Effect of Interparticle Friction Parameters: μ_S and β

Figure 5 shows the evolution of the deviator stress, σ₁–σ₃, the volumetric strain, ε_v, and the mechanical coordination number (CN, i.e., the mean number of contacts for particles with two or more contacts; particles with one or zero contact cannot transmit forces) of tested specimens during triaxial compression processes for different sets of parameters, μ_S and β, corresponding to simulations 1 to 8 in Table 1. At the beginning, the vertical axial loading compresses the specimen, leading to slight increases in the volumetric strain ε_v. Meanwhile, a biased orientation of the force chains along the vertical direction is formed to resist the axial loading, causing breakage in the contacts along other directions. As a result, the CN then decreases rapidly. As the vertical compression continues, the lateral walls move outward to maintain their equilibrium. Therefore, the initial volumetric contraction is followed by dilation. Particles also move outward to fill the space, weakening the force networks along the vertical direction and the ability to resist the axial loading. So, the growth rate of the deviator stress slows down and exhibits strain-softening behavior after the peak state is achieved, and a constant CN is maintained. The results show that the typical mechanical response of dense sands, e.g., the nonlinear stress-strain behavior and volumetric dilation behavior, can be well reproduced using our SSDEM model.

Zoom In Zoom Out Reset image size

In addition, some conspicuous features of the effect of μ_S and β can be observed under each confining pressure:

• 1.  
  The deviator stress as well as its peak value and the corresponding axial strain at the peak state all increase with μ_S and β.
• 2.  
  The maximum compaction in the volumetric strain and the corresponding axial strain increase with μ_S and β.
• 3.  
  The value of CN at the dilation stage decreases as μ_S or β increases.

The above findings indicate that higher friction resistance can improve the robustness of the contact force network and require more effort to disrupt its structure.

The shape parameter β seems to play a smaller role than the sliding friction coefficient μ_S does, especially for the cases of μ_S = 0.3. Figure 6(a) compares the peak stress variables and the corresponding failure envelopes in the I₁–$\sqrt{{J}_{2}}$ plane for different μ_S with various β. Generally, a higher friction angle ϕ can be achieved with a larger μ_S or β, which is physically expected. However, as shown in Figure 6(a) and Table 1, the effect of β on the friction angle ϕ declines when β becomes larger, and the maximum achievable friction angle is limited by the value of μ_S, e.g., ∼32° for μ_S = 0.3. This implies that the shear strength of a granular assembly is strongly controlled by the sliding frictional resistance between components. Once particles start sliding relative to each other, the rotational resistance, regardless of how large it is, cannot stop this process, and the shear strength will be weakened, resulting in the strain-softening behavior. As shown in Figure 5, this strain-softening behavior is more pronounced with larger μ_S, in agreement with previous simulations (e.g., Göncü & Luding 2013).

Zoom In Zoom Out Reset image size

3.3.2. Effect of the Interparticle Cohesive Tensile Strength c

The granular representations of rubble-pile asteroids with a given shape are constructed using the same procedure as described in Section 3.1. An oblate model with axis ratios α₁ of ∼0.9, α₂ of ∼1.0, a maximum diameter of ∼1 km, and a bulk density ρ_B of ∼2.4 g/cc is taken as the rubble-pile sample in this analysis, where α₁ = a₃/a₁, α₂ = a₂/a₁. The axis lengths a₁, a₂, and a₃ are given as the maximum extensions along the principal axes of inertia (from largest to smallest) of the rubble pile. The bulk density of a rubble pile is calculated as ρ_B = η_Bρ_p, where η_B is the internal packing efficiency (defined in Equation (9)) and ρ_p is the particle material density. This rubble-pile model consists of ∼18,000 particles, with a −3-index power-law distribution in particle radius ranging from 10 m to 30 m, and ρ_p of ∼3.3 g/cc.¹² The normal stiffness k_N and the timestep Δt are set to meet the maximum particle overlap requirement.

During a spin-up simulation, with given SSDEM parameters, the simulated rubble pile evolves under its self-gravity, following a spin-up path, which sets its spin period T to a prescribed value as a function of time t (see the top frame in Figure 7 for an example; only maximum principal-axis rotation is considered in this study, where z is the spin axis). The spin-up procedure is similar to our previous study (see Section 2.3.4 in Zhang et al. 2017, for details). In brief, in the beginning of a simulation, the rubble-pile body is allowed to settle down at a slow spin state T_start (T_start = 4 hr is used in this study), and then, the body is linearly spun up from T_start until global disruption occurs, after which we stop maintaining the spin period of the simulated body and let the body evolve freely under its own gravity.¹³

Zoom In Zoom Out Reset image size

Considering that the mechanical behaviors of granular materials have different interpretations when analyses are carried out at different scales (Andrade & Tu 2009), variables defined at three different scales are monitored throughout the spin-up process. At the global scale, the shape, characterized by α₁ and α₂, the internal packing efficiency, η_B, and the mechanical coordination number, CN, are used to capture the global properties of the simulated body. The internal packing efficiency η_B is the fraction of the volume in the internal structure of a rubble pile that is occupied by the constituent particles, which can be measured based on the Voronoi tessellation of this rubble pile (see Section 3.2 in Zhang et al. 2017, for details),

$$\begin{eqnarray}&&{\eta }_{B}=\displaystyle \frac{1}{{N}_{\mathrm{internal}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{internal}}}\displaystyle \frac{4\pi {r}_{i}^{3}}{3{V}_{i}},\end{eqnarray} \tag{ 9 }$$

where r_i is the radius and V_i is the volume of the Voronoi cell for the ith particle, and N_internal is the number of particles in the interior of the rubble pile.

The global dynamical behaviors of a granular material depend on the interaction of individual particles at the microscopic scale, which can be captured by the force chain distribution. The contact networks in the simulated rubble pile can be visualized by tracing the normal contact strength, S_N, i.e., the total normal contact force divided by the effective contact area, so that the contacts that reach the cohesion limit (i.e., the given value of the interparticle cohesive tensile strength c used in the simulation) are revealed.

Furthermore, to link the properties of the simulated rubble pile with the conventional continuum medium and present the local characteristics at the macroscopic scale, it is also useful to analyze the stress distribution in the rubble pile. The stress analyses are carried out on the scale of representative volume elements (RVEs; see Section 3.3.1 in Zhang et al. 2017, and references therein for details). The size of an RVE is often taken to be about 12 times the radius of the typical particle for stress analyses (Masson & Martinez 2000). We use a tree code to divide the rubble-pile body into several similar-sized sub-assemblies, i.e., RVEs. Each RVE contains 200 to 300 particles. The average stress tensor of a local region, e.g., the jth RVE, in a rubble pile can be assessed by homogenization and averaging methods,

$$\begin{eqnarray}&&{\bar{{\boldsymbol{\sigma }}}}_{j}^{\mathrm{RVE}}=\displaystyle \frac{1}{{V}_{\mathrm{RVE}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{RVE}}}{V}_{i}{{\boldsymbol{\sigma }}}_{i}=\displaystyle \frac{1}{{V}_{\mathrm{RVE}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{RVE}}}\displaystyle \sum _{k=1}^{{N}_{i,C}}{{\boldsymbol{x}}}^{i,k}\otimes {{\boldsymbol{f}}}^{i,k},\end{eqnarray} \tag{ 10 }$$

where N_RVE is the number of particles in this RVE and V_RVE is the total volume of this RVE. The Cauchy stress tensor ${{\boldsymbol{\sigma }}}_{i}$ for a single particle i is given as the summation of the dyadic product of the branch vector ${{\boldsymbol{x}}}^{i,k}$ that connects the particle center with the contact point and the corresponding contact force ${{\boldsymbol{f}}}^{i,k}$ for the kth contact. N_i,C is the number of contacts for the ith particle. The first invariant of the average stress tensor ${I}_{j,1}^{\mathrm{RVE}}$ and the second invariant of the deviatoric stress tensor ${J}_{j,2}^{\mathrm{RVE}}$ for the jth RVE can be determined from the averaged stress tensor ${\bar{{\boldsymbol{\sigma }}}}_{j}^{\mathrm{RVE}}$. Analyzing the distribution of these stress-state variables can help us interpret how the simulated rubble pile responds to the spin-up loads.

To understand the effect of the SSDEM parameters on the dynamical behaviors and material strengths of rotating rubble piles, multiple spin-up simulations were carried out in this study, as summarized in Table 2.

4.3.1. Effect of Interparticle Friction Parameters: μ_S and β

Figure 7 shows the evolution of the state variables of the simulated body during the spin-up process for different sets of parameters, μ_S and β, with c = 0 Pa. The top left panel of the Figure 12 animation presents the spin-up and disruption process in the case of μ_S = 0.5, β = 0.5, and c = 0 Pa. As expected, an increase in the sliding friction or rotational friction can postpone the initiation of failure. When the body has been spun up past its failure limit, its shape becomes more oblate, and the packing efficiency and coordination number gradually decrease. The evolution of these state variables exhibits similar failure behaviors, regardless of the value of μ_S or β. The finding that the failure mode does not show strong dependencies on the friction parameters is in agreement with the analytical theory by Hirabayashi (2015). The simulations of Sánchez & Scheeres (2016), using a different SSDEM code, show a contrary result, namely that different friction angles lead to different disruption patterns. A plausible explanation for this is that they achieved different friction angles by turning off the rolling friction and the sliding friction separately, which changes the intrinsic properties of the contact network, and so the dynamical characteristics are different. In reality, neither the rolling friction nor the sliding friction is negligible for natural granular materials. Therefore, the failure mode should not be strongly affected by the frictional strength if both kinds of friction are turned on.

We take the case of μ_S = 0.5 and β = 0.5 as a representation to detail the characteristics of a cohesionless rubble pile during the spin-up process. Figure 8 illustrates the distribution of the stress-state variables and the contact strength networks over a cross-section of the simulated rubble pile at the slow spin state (T = 4.00 hr), the moderate spin state (T = 2.90 hr), and the critical spin state (T = 2.56 hr, which is the critical spin period for this model; see Section 5 for details about how to determine the critical spin limit of a rubble-pile body using SSDEM). Figure 9 compares the distributions of contact orientations in the rubble pile at different spin states. As would be expected for this random configuration, the distributions of contact force orientations are uniform at slow spin (see the areas encompassed by the blue curves in Figure 9). Since the centrifugal forces become larger as the spin rate increases, the force chains parallel to the x–y plane gradually break or rearrange themselves in other directions, leading to an asymmetric distribution of contact orientations (Figures 8(g)–(i); also see the areas encompassed by the red and orange curves in Figure 9). As shown in the first row of Figure 8, increasing the centrifugal forces also diminishes the central pressure for this body. Meanwhile, given that the forces along the vertical direction caused by gravity remain constant, the body experiences larger shear stresses as the spin rate increases (second row of Figure 8), and the strong force chains gradually align parallel to the rotation axis (bottom row of Figure 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the evolution of the state variables follows the same trend and the failure mode is independent of μ_S or β, the failure spin of the rubble pile shows clear dependencies on the interparticle friction parameters. That is, the rubble pile can be stable at a higher spin rate as μ_S or β increase (see Table 2 for a summary of the critical spin periods). However, the effect of β is less significant than that of μ_S. When cohesion is not included, the shear failure of a rubble-pile structure is caused by microscale sliding, rolling, and twisting between particles in contact. The friction coefficient μ_S dominates the sliding resistance, while the shape parameter β dominates the rolling and twisting resistance in SSDEM simulations. Figure 10 shows the distributions of the orientations for the sliding/rolling/twisting contacts at the critical spin state for the case of μ_S = 0.5 and β = 0.5 (other cases show similar results). The sliding/rolling/twisting contact is defined as the contact mode for which the tangential force/rolling torque/twisting torque exceeds the critical value, μ_SF_N/M_R,max/M_T,max. As seen in Figure 10, the number of sliding contacts is much larger than that of the rolling or twisting contacts at the critical spin state. This implies that at the microscale of a spinning rubble pile, particle rearrangements are mainly caused by sliding movements between particles in contact, i.e., sliding rearrangement is the preferred mode of particle-level deformation. Therefore, the friction coefficient μ_S has a much greater impact on the structural stability of rubble piles. Our study reveals the complexity and importance of the interparticle contact mode and friction parameters in controlling the strength and deformation mode of a granular assembly. Furthermore, it is interesting to note that the distributions of the three contact modes have different preferred directions in the y–z and z–x plane. For example, since the mid-latitude regions are subject to the highest slope for such a shape, as revealed by Harris et al. (2009), sliding failure is most likely to develop in these regions (see the areas compassed by the blue curves in Figure 10).

Zoom In Zoom Out Reset image size

4.3.2. Effect of Interparticle Cohesive Tensile Strength c

Since the cohesive strength of the interstitial regolith arising from van der Waals forces could range from hundreds of pascals to thousands of pascals in asteroids (Sánchez & Scheeres 2014; Gundlach & Blum 2015), the interparticle cohesive tensile strength c was selected to be 100 to 1600 Pa for the cohesive SSDEM simulations presented in this section, in which μ_S and β were both set to 0.5. As shown in Figure 11, the evolution of the rubble pile progresses in different manners with changes in c. The existence of cohesion can improve the strength of the contact force network by increasing the coordination number and the packing efficiency. As c increases, the rubble-pile structure can maintain its shape and internal packing efficiency at a higher spin rate, and the amount of deformation that appears before disruption decreases, followed by a more abrupt disruption process.

Zoom In Zoom Out Reset image size

In addition, the style of fragmentation during the disruption process also shows a strong dependence on c, as shown in Figure 12 (also see the animation of Figure 12 for the spin-up and disruption processes with different c). Recall that the spin period of the simulated body is strictly constrained to the spin-up path until α₁ is reduced to half of its initial value, so the body is continuously subject to the spin-up loading even after its structure fails. Normally, for kilometer-size asteroids, the timescale of the YORP loading process exceeds 10⁴ years (Rubincam 2000), which is much longer than the dynamical time needed for a rubble pile to settle down in the absence of external forces. Therefore, in reality, the low-cohesion rubble pile is more likely to shed some surface materials or slightly deform to offset the extra spin, rather than completely break up. However, for the high-cohesion cases, since the deformation of the rubble pile before disruption is inhibited, the disaggregation behaviors of the high-cohesion rubble piles observed in the spin-up tests could represent the failure behaviors of real cohesive rubble-pile asteroids. Taking the cases of c = 0 Pa and c = 800 Pa as examples, we carried out two other runs without continuous rotational acceleration. In both cases, the body is spun up to a spin period below its critical value and is forced to stay at this spin period until α₁ is reduced by 0.1, after which the body is set to be free of external forcing. Figure 13 presents the outcomes of these two simulations, and the animation of Figure 13 records the evolution processes. As can be seen in the animation of Figure 13, instead of continuous deformation as observed previously (see the top left panel in the animation of Figure 12), landslides and mass shedding are the main failure features for the cohesionless body, while the failure behavior for the cohesive case is almost the same as that shown in the continuous spin-up test (see the bottom left panel in the animation of Figure 12). The distinguishing failure behaviors of the cohesive body could be used to reveal the level of cohesion in a disintegrating asteroid from observations.

Download figure:

Video Standard image High-resolution image

Download figure:

Video Standard image High-resolution image

As shown in Figure 12 (also see the bottom panels in the animation of Figure 12), when the cohesion is high enough, i.e., the cases of c = 800 and c = 1600 Pa, the rubble pile can directly split into multiple distinct fragments consisting of many particles. For this reason, the value of CN after global disruption in the high-cohesion cases can stay at a moderate constant value (see the bottom frame in Figure 11). Figure 14 shows the cumulative mass distributions of the fragments produced in the rotational fission events for different c. The size distribution of the resulting fragments depends on the cohesive forces that bond particles together. When cohesion is low, the rubble pile will disrupt into individual particles; however, when cohesion is high enough, the rubble pile will disaggregate into similar-size fragments. The disaggregating behaviors of cohesive rubble piles could be a plausible mechanism to form asteroid pairs (Pravec et al. 2010). The resulting fragments escape from each other at speeds on the order of ∼0.5 m s⁻¹, close to the mean pair-wise sky-plane velocity dispersion between fragments of the disintegrating asteroid P/2013 R3 (i.e., 0.33 ± 0.03 m s⁻¹; Jewitt et al. 2017). The properties of the fragments (e.g., size distribution, velocity dispersion) in the high-cohesion cases show many similarities to the observations of P/2013 R3, indicating that this asteroid may possess relatively high cohesion in its structure and experience rotational fission in a manner similar to the high-cohesion cases.

Zoom In Zoom Out Reset image size

The effect of cohesion, whereby deformation is diminished prior to disruption and the component size of the fissioned material is increased, was also observed in the numerical study of Sánchez & Scheeres (2016) using a different SSDEM code. However, unlike our high-cohesion cases, where the body is observed to axisymmetrically disaggregate into similar-size fragments, their simulations generally show that the rubble pile is disrupted by non-axisymmetric deformation and partial fission, even in the spherical case (see Table B.3 in Sánchez & Scheeres 2016). This behavior is expected when the simulated rubble pile possesses a non-uniform particle arrangement in the rotation axis direction, which leads to a non-uniform distribution of structural strength, but not when the body is axisymmetric (Holsapple 2008). The observed difference, therefore, may be due to the axial non-uniformity of the particle arrangement in the rubble-pile models used by Sánchez & Scheeres (2016). In addition, the discrepancy in the spin-up path and the numerical model between our study and theirs may also account for this difference.

Figure 15 shows the distributions of the stress-state variables and the contact strength networks over a cross-section of the simulated body with c = 800 Pa, at the slow spin state (T = 4.00 hr), the critical spin state for the cohesionless case (T = 2.56 hr), and the critical spin state for itself (T = 2.20 hr). The contact force networks consist of a significant number of cohesive force chains (i.e., the chains with S_N < 0), even in the slow spin state, which can slow down the breakage of force chains. As the spin rate is increased, some of the force chains, especially those parallel to the equatorial plane, reach the cohesion extremum, i.e., 800 Pa in magnitude for this case, implying that the corresponding particles are barely in contact with each other, and the connection could break as the spin-up process continues. As shown in the bottom frame of Figure 11, the amount of contacts that are lost before disruption dramatically decreases as c increases. Thus, failure in this high-cohesion body is likely to be stimulated when a sufficient number of cohesive force chains reach the cohesion extremum and break all together within a short time span.

Zoom In Zoom Out Reset image size

Comparing the case for T = 2.56 hr in Figure 15 with that in Figure 8, it is evident that higher cohesive strength can impose higher pressure and ease the shear stresses in the interior. Remarkably, the stress states of some surface regions transform from a compressive mode into a tensile mode as the spin rate increases (Figure 15(c)). The failure envelope analyses in Section 5.2 indicate that the failure point for this body occurs in the second quadrant in the ${I}_{1}-\sqrt{{J}_{2}}$ plane (see Figure 16). Therefore, the damage of the surface structure is mainly caused by tensile failure rather than shear failure in the high-cohesion case. As observed in the bottom panels in the animation of Figure 12, the fracture of the surface region can propagate into the interior region and stimulate the disaggregation of the whole body.

Zoom In Zoom Out Reset image size

According to the Drucker–Prager failure criterion, the failure condition of a granular assembly depends on both the shear stress and the pressures that act on it, so that the failure envelope can be drawn at several failure points under different confining pressures, as shown in the triaxial compression tests (e.g., see Figure 6). Normally, when compared to the material strength of a rubble-pile asteroid, the surface pressure, e.g., caused by the solar radiation pressure, is negligible, and the boundary of this body can be treated to be stress-free, which means the confining pressures are all zero. Therefore, instead of varying the surface pressure as done in the triaxial compression tests, we achieve various pressure conditions by changing the bulk density of the simulated rubble pile. For each SSDEM parameter set, (μ_S, β, and c), in addition to the nominal runs with the bulk density of 2.4 g/cc as used in Section 4, two other simulations with bulk densities of 1.8 g/cc and 3.0 g/cc were carried out to support the failure envelope analyses.

Since the stress state will vary throughout a rubble-pile structure (see Figure 8 for examples), the next step is to determine the failure stress state for each bulk density case. Previous studies show that the starting point of failure in a rubble pile always occurs locally in a region that is subject to the highest shear-pressure ratio (e.g., Hirabayashi 2015; Zhang et al. 2017). The RVEs can serve as the smallest element for the stress analyses. To measure how close an RVE is to its structural failure point, we introduce a new variable, the stress ratio λ, which for the jth RVE is given by

$$\begin{eqnarray}&&{\lambda }_{j}^{\mathrm{RVE}}=\sqrt{{J}_{j,2}^{\mathrm{RVE}}}/(C+{{sI}}_{j,1}^{\mathrm{RVE}}),\end{eqnarray} \tag{ 11 }$$

where C and s are the strength parameters used in the Drucker–Prager failure criterion (see Section 3.2). If this stress ratio exceeds 1, the jth RVE should experience plastic deformation. When the material is cohesionless, λ^RVE_j is proportional to the ratio of $\sqrt{{J}_{j,2}^{\mathrm{RVE}}}$ to ${I}_{j,1}^{\mathrm{RVE}}$. This implies that the part in a cohesionless rubble pile most sensitive to failure should possess the highest $\sqrt{{J}_{j,2}^{\mathrm{RVE}}}/{I}_{j,1}^{\mathrm{RVE}}$.

Figure 16(a) shows the stress-state variable distribution of a rubble pile's RVEs at the critical spin state in the I₁–$\sqrt{{J}_{2}}$ plane for the case of μ_S = 0.5, β = 0.5, and c = 0 Pa. It is evident that the shear stress and the pressure at the failure state of this rubble pile increase in proportion to the bulk density. The failure point is identified as the one that has the highest $\sqrt{{J}_{j,2}^{\mathrm{RVE}}}/{I}_{j,1}^{\mathrm{RVE}}$ for each bulk density case, i.e., the red symbols in Figure 16(a). Then, the values of ϕ and C can be determined from the best-fit failure envelope, i.e., the dashed line in Figure 16(a). We use this method to derive the strength parameters for all cohesionless cases. As shown in Table 2 (simulations 1 to 6), a higher friction angle ϕ can be achieved with a larger μ_S or β. Since there is no apparent cohesion in the material, the values of C are very close to zero. Furthermore, it is interesting to note that there are significant differences between the friction angle obtained from the spin-up tests (see Table 2) and that in the triaxial tests (see Table 1), indicating the high sensitivity of the properties of a granular material to the loading scenario and the particle arrangement.

When cohesion is present, the search for the failure point becomes complex due to the uncertainties of both C and s in Equation (11). Our triaxial simulation results showed that the friction angle is essentially independent of the interparticle cohesive tensile strength c. For simulations with the same μ_S and β, it is therefore proper to assume that their failure envelopes have the same slope. By translating the failure envelope obtained from the corresponding cohesionless case up and down, the failure point for each simulation can be identified as the one that has the largest intercept in the I₁–$\sqrt{{J}_{2}}$ plane.

Figure 16(b) shows the stress-state variable distribution of a rubble pile at the critical spin state for the simulation case of μ_S = 0.5, β = 0.5, and c = 800 Pa. The failure points determined using the approach introduced above for each bulk density case are highlighted by red symbols. The failure envelope for this cohesive case is then established through a linear least-squares procedure, which provides the friction angle ϕ and the cohesion C. Using this method, the values of ϕ and C were derived for all cohesive rubble piles in our simulations, as presented in Table 2. The small differences in estimates of the friction angle shown in the results once again confirm the independence between c and ϕ. Thus, as a practical matter, to save the effort of conducting two more spin-up tests, the value of C can be directly estimated from one failure point by assuming the friction angle is the same as the cohesionless case.

As shown in Table 2 (simulations 7 to 11), the macroscopic cohesive strength C is on the order of one percent of the interparticle cohesive tensile strength c and displays a nonlinear correlation with c. This is because a higher c can increase the number of cohesive connections and enhance the efficiency of cohesion, as discussed in Section 4.3.2. Our high-cohesion cases show that C of ∼30 Pa is enough to stimulate the disaggregating behavior of rubble-pile asteroids. Therefore, the observed active asteroid P/2013 P3 may possess cohesion close to or higher than this value. This argument is consistent with the velocity dispersion between fragments from observations, which indicates an effective cohesive strength ∼50–100 Pa (Jewitt et al. 2017).

Figure 17 plots the distributions of the stress ratios over the internal structure for the cases of c = 0 Pa and c = 800 Pa, with μ_S = 0.5, β = 0.5. The corresponding values of ϕ and C listed in Table 2 are used to calculate the stress ratios for each RVE. It is evident that increases in cohesion can postpone the onset of failure. As shown in Figure 17, the rubble piles are well below the failure condition at the slow spin state, and the stress ratios generally grow with the spin rate. The initial failure occurs at the surface near the equatorial region in both cases, but different failure behaviors are observed. A surface shedding-failure mode caused by shear failure of particles resting on the surface is confirmed for the cohesionless case (see the left panel in the animation of Figure 13), while a surface tensile-failure mode is identified for c = 800 Pa (see the right panel in the animation of Figure 13 and the last paragraph in Section 4.3.2). Our study indicates that in addition to the distribution of stress ratios, information on pressure states is also needed to properly predict the failure mode.

Zoom In Zoom Out Reset image size

The complete elastic solution of a gravitational spinning ellipsoid in equilibrium can be obtained by solving the elastic equations (where Hooke's law is used to describe the constitutive relation). Dobrovolskis (1982) presented the derivation in detail. Hirabayashi (2015) substituted the equilibrium stress solution into the Drucker–Prager failure criterion to analyze the failure condition and mode of a cohesive spherical body. Using a similar approach, we carry out stress analyses for a 1-km-diameter oblate body with an axis ratio α₁ = 0.9 at the failure spin for the cases with μ_S = 0.5, β = 0.5, and c = 0 Pa and μ_S = 0.5, β = 0.5, and c = 800 Pa (corresponding to C = 12.7 Pa), as shown in Figure 18. The quantity $\hat{{\rm{\Omega }}}$ is the spin rate normalized by $\sqrt{4\pi {\rho }_{B}G/3}$, where G is the gravitational constant. In both cases, the central region is subject to the highest pressure and the largest shear stress, and as the spin rate increases, the pressure decreases while the shear stress increases. The distribution characteristics of the stress field and the response to rotational acceleration of our simulation results are all consistent with the analytical solutions (e.g., see Figures 8 and 15). The magnitude of the pressure and shear stress field in the cohesionless case even shows quantitative consistency. For example, as shown in Figures 8(c) and (f), the maximum values of ${I}_{1}^{\mathrm{RVE}}$ and $\sqrt{{J}_{2}^{\mathrm{RVE}}}$ appear near the center, which are larger than 200 Pa and 40 Pa at the failure state for the simulated body, respectively. Similar phenomena are also shown in the results of the elastic theory (i.e., Figure 18).

Zoom In Zoom Out Reset image size

As characterized by the stress ratio λ (bottom row in Figure 18), the equatorial surface is the most sensitive part when the rubble-pile structure yields for C = 0 Pa, in agreement with our numerical results for c = 0 Pa (see Figure 17), while the central region is subject to the highest stress ratio in the high-cohesion case. This phenomenon was discussed in the analytical study of Hirabayashi (2015), where he found that for a homogeneous, cohesive body, the region subject to the highest relative shear stress migrates from the equatorial surface to the central region as the dimensionless spin rate increases to near 1. The inconsistency between our numerical results and Hirabayashi's theory in the high-cohesion case may come from the differences between the cohesion mechanism in the discrete element modeling and that in the continuum method. The cohesive strength can be completely controlled by a specific variable in continuum analyses (e.g., C in the Drucker–Prager failure criterion), while the distribution of the macroscopic cohesive strength cannot be perfectly uniform in a discrete element model even for a constant c due to heterogeneity of radial particle arrangements. As revealed by Zhang et al. (2017), the local packing efficiency and coordination number for the surface particles is much lower than that for the interior particles. Therefore, the local macroscopic cohesive strength is lower in the surface region for a rubble-pile model, resulting in a weaker surface shell. As a result, surface tensile failure is observed for high-cohesion bodies in the SSDEM simulations.

Overall, our simulation results show good consistency with the results of elastic theory in the cohesionless case, while the elastic theory fails to explain our findings for the cohesive model. In fact, in the course of solving the elastic equations, no information on the cohesion is taken into account. The value of C only affects the calculation of the stress ratio, and the distributions of the pressure and the shear stress are irrelevant to the cohesive strength. This treatment in the analytical method is therefore incorrect. As highlighted in our study (Section 4.3.2; also see Figure 15), the existence of cohesion changes the topology of the force network, where higher cohesive strength can impose higher pressure and ease the shear stress in the interior of a rubble pile. Furthermore, the usage of elastic theory implicitly requires an initial globally stress-free state, which is an unrealistic assumption for natural granular mediums in space (Holsapple & Michel 2006). Meanwhile, the stress-strain relation for a granular assembly is nonlinear generally (see Figure 5 for examples). Therefore, the actual stress distribution in a rubble pile (especially when cohesion is present) could be quite different from the results of elastic theory. Special care should be taken when using the elastic theory to predict the behavior of a cohesive rubble-pile body.

Since the state of a body when global failure will occur is independent of the loading history, the global-volume-averaging method can give an upper-bound condition on the critical spin rate (Holsapple 2008). Figure 19 shows the scaled critical spin rate for our 11 simulation cases with μ_S = 0.5 and β = 0.5 listed in Table 2 compared to the theoretical upper spin limits for various C derived by Holsapple (2007). The spin limits of the cohesionless cases are well described by the continuum theory, while the cases with nonzero c plot above the theoretical upper spin limit. For example, the model with c = 200 Pa, whose macroscopic cohesive strength C is determined to be 2.1 ± 0.1 Pa, can maintain its original shape at a spin rate higher than the theoretical spin limit derived for C = 5 Pa. This implies that our method for determining the material strengths of simulated rubble piles may underestimate the macroscopic cohesion. In fact, as discussed in Section 6.1, the surface region in a simulated rubble pile is naturally weaker than its interior. Given that the surface of the cohesive body fails first and the macroscopic cohesion is estimated according to the stress in this region, the obtained value of C is expected to be lower than the averaged cohesive strength of the whole body as used in the continuum theory.

Zoom In Zoom Out Reset image size