STRESS AND FAILURE ANALYSIS OF RAPIDLY ROTATING ASTEROID (29075) 1950 DA

Near-Earth asteroid (29075) 1950 DA is currently considered to be among the most hazardous asteroids due to its close encounter to the Earth in 2880 (Giorgini et al. 2002). Busch et al. (2007) conducted comprehensive radar observations of this asteroid and derived two possible shape models. The first is a prograde model, which is a spherical shape, and the second is a retrograde model, which is an oblate shape. Farnocchia & Chesley (2014) reported that the orbital semi-major axis of 1950 DA has been changing due to the Yarkovsky effect and confirmed that the retrograde model was consistent with their analysis. Because of its spin period, 2.1216 hr, this object may be close to its structural failure point if this body is a rubble pile. If 1950 DA has no cohesion, i.e., zero shear resistance at a zero-normal stress, the bulk density should be higher than 3.5 g cm⁻³ to prevent the body from failing structurally (Busch et al. 2007). Using an advanced thermophysical model and archival WISE thermal-infrared data to derive a bulk density of 1.7 ± 0.7 g cm⁻³, Rozitis et al. (2014) indicated that 1950 DA should have cohesive strength to keep its current shape. They used Holsapple's limit analysis technique to derive the lowest cohesion for keeping the current shape of 1950 DA. The detailed development was described by Holsapple (2004, 2007). However, since the applied technique was based on taking the average of stress components over the whole volume, it discards crucial information about the disruption processes.

Disruption events of asteroids have been recently reported. For example, active asteroid P/2013 R3 was observed breaking into multiple components (Jewitt et al. 2014). From this fact, Hirabayashi et al. (2014) found that this asteroid has a level of cohesion ranging between 40 Pa and 210 Pa. Other disrupting bodies have also been observed. Active astroid P/2013 P5 is currently shedding discontinuous dust tails (Jewitt et al. 2013) and Sheppard & Trujillo (2014) discovered a new active asteroid, (62412) 2000 SY178, exhibiting a tail from its nucleus. These observations reveal that disruption events in our solar system are more common than previously known. Thus, an understanding of such disruption events provides better constraints on the formation of small bodies.

Holsapple (2008) developed a finite element model that takes into account plastic deformation characterized by the von Mises yield criterion, a pure-shear-dependent criterion, and applied it to the failure condition of a rotating, non-gravitating ellipsoid. We extend his model to a model that can take into account self-gravity and plastic deformation characterized by a shear-pressure yield criterion. For our computations we use a commercial finite element software, ANSYS, version 15.03. We investigate the internal condition of the 1950 DA retrograde model at the current spin period. Our analysis provides a more precise lower bound on the necessary cohesion to hold 1950 DA together in a stable state. More importantly, our stress analysis also shows that the preferred mode of failure is for the central region of the body to collapse along the rotation axis, which would cause the equator to expand outward. The analysis indicates that the body's failure state is not close to either surface material being shed or landsliding. Given that plastic flow will generally increase volume in an associated flow rule (Chen & Han 1988), this failure mode predicts that the central core of the asteroid may have a reduced density, which is a previously unpredicted state for such oblate, rapid rotators.

Information from radar observations can constrain an asteroid's near-surface bulk density, as a function of the material density and porosity (Magri et al. 2001). Using this technique, Busch et al. (2007) obtained the minimum surface material density as 2.4 g cm⁻³ and derived two compositional possibilities of the surface, enstatite chondrite and nickel-iron. On the other hand, the thermophysical model used by Rozitis et al. (2014) gave a bulk density of 1.7 ± 0.7 g cm⁻³. Here, we use the bulk density, derived by Rozitis et al. (2014), which also includes the estimated value by Busch et al. (2007).

The retrograde model of 1950 DA derived by Busch et al. (2007) is oblate and has an equatorial ridge (Figure 5 in their article). In the following discussion, the minimum, intermediate, and maximum moment of inertia axes are defined as the x₁, x₂, and x₃ axes, respectively. The mesh is constructed using the Busch et al. (2007) retrograde model. Each element is a 10 node tetrahedron. The numbers of elements and nodes are 5569 and 8595, respectively. The initial size of the shape is fixed to be constant and so the initial volume is 1.145 km³. The spin period is 2.1216 hr. It is assumed that 1950 DA is a principal axis rotator and is affected only by self-gravity and rotation. Here, the rotation pole is fixed along the x₃ axis. We define Young's modulus as 1.0 × 10⁷ Pa, a typical value for geological sand and gravel, and Poisson's ratio as 0.25, giving compressibility of the volume. A friction angle, ϕ, is fixed at 35°, the mean value of a friction angle of a geological material ranging from 30° and 40° (Lambe & Whitman 1969). Cohesion is kept as a free parameter and we will determine the lowest cohesion that prevents this object from failing structurally. These properties are given in Table 1. Here, we calculate finite element solutions for three different bulk density cases: 1.0 g cm⁻³, 1.7 g cm⁻³, and 2.4 g cm⁻³.

In the current finite element model, if a stress state is below the yield condition, the behavior of a material follows linear elasticity. It is assumed that the behavior of plastic deformation does not include material hardening and softening (perfect plasticity). Body forces are defined at each node so that the total force acting on each element is equally split and apparently concentrates on its nodes. In nature, an asteroid rotates in free space, implying that there are no boundary conditions. However, in our numerical environment, the setting of body forces makes it difficult for the simulations to converge correctly because it could cause rigid motion, i.e., rotation and translation. To avoid this issue we artificially constrain 6 dof to eliminate such motion.

1950 DA is currently considered to be a rubble pile (Rozitis et al. 2014). For such a body, the mechanical behavior of its material depends on shear and pressure (Lambe & Whitman 1969; Jaeger & Cook 1972). In this study, we use the Drucker–Prager yield criterion, a smooth function in the principal stress space, which is given as

$$\begin{eqnarray} &&f (I_1, J_2) = \alpha I_1 + \sqrt{J_2} - k = 0, \end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray} I_1 &=& \sigma _1 + \sigma _2 + \sigma _3, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} J_2 &=& \frac{1}{6} [(\sigma _1 - \sigma _2)^2 + (\sigma _2 - \sigma _3)^2 + (\sigma _3 - \sigma _1)^2]. \end{eqnarray} \tag{ 3 }$$

σ₁, σ₂, and σ₃ are the principal components of the stress state, and α and k are free parameters that need to be fixed. We choose these parameters by the following consideration. Since the retrograde shape model of 1950 DA is oblate, only the stress components along the x₁ and x₂ axes are affected by the centrifugal force. In such a case, the relation of the stress components is described as σ₁ ∼ σ₂ > σ₃ (here, a negative value means compression, while a positive value indicates tension). Thus, since the actual stress state is near the compression meridian of the Mohr–Coulomb yield criterion, σ₁ = σ₂ ⩾ σ₃, we choose α and k such that the Drucker–Prager yield envelope touches the Mohr–Coulomb yield envelope at this meridian. These parameters are then given as (Chen & Han 1988)

$$\begin{eqnarray} \alpha &=& \frac{2 \sin \phi }{\sqrt{3} (3 - \sin \phi)}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} k &=& \frac{6 c \cos \phi }{\sqrt{(}3) (3 - \sin \phi)}, \end{eqnarray} \tag{ 5 }$$

where c is cohesion and ϕ is a friction angle.

Plastic deformation is assumed to be small. The constitutive law of plastic deformation is constructed based on an associated flow rule. This rule guarantees that the direction of plastic deformation is always perpendicular to the yield envelope. The plastic strain rate $\dot{\epsilon }^p_{ij}$, where i and j are indices, is described by the rate of an arbitrary coefficient λ times the partial derivative of the function f with respect to the stress tensor σ_ij, which is given as

$$\begin{eqnarray} &&\dot{\epsilon }^p_{ij} = \dot{\lambda }\frac{\partial f(\sigma _{ij})}{\partial \sigma _{ij}}. \end{eqnarray} \tag{ 6 }$$

Note that it is known that the behavior of a typical geological material would rather follow a non-associated flow rule than be characterized by an associated flow rule. However, since (1) choosing an associated flow improves the convergence of a finite element solution in the current model and (2) consideration of a non-associated flow rule requires additional free parameters (usually unknown), we use the associated flow rule described in Equation (6). We leave the application of a non-associated flow rule as future work.

A plastic solution depends on loading path. Although such a path represents evolutional history, in general it is usually unknown. Here, we assume that the gravitational and centrifugal forces acting on small elements ramp up linearly. This load step implies that the evolution of 1950 DA results from its accretion process due to a catastrophic disruption. The timescale of the accretion process due to a catastrophic disruption is considered to be negligibly short compared to the life time of an asteroid, probably much less than a year (Michel et al. 2001, 2002, 2004).

To derive the plastic deformation mode at the lowest cohesion that can keep the current shape of 1950 DA, we conduct an iteration method described as follows. First, we choose a value of cohesion and calculate a plastic solution that corresponds to it. If this solution only includes elastic states everywhere, we choose a lower value for cohesion and recalculate the stress solution. We iterate this process until we obtain a solution in which a majority of the internal structure reaches the yield condition. To go beyond such a condition, one usually encounters computational issues and so we terminate our iteration. Thus, the solution described in this study is as close as possible to the numerical computation limit of plastic analysis.

We implement a three-dimensional finite element model of 1950 DA (the retrograde model) on ANSYS, version 15.03. The purpose of this paper includes visualization of plastic deformation in the interior of an asteroid. To do so, we use the ratio of the current stress state to the yield stress, so-called stress ratio (Kohnke 2009). Although this parameter does not show the plastic state during a loading path, (in order words, it does not show the history of the plastic state), it is useful to visualize the final state. Using the derivation given by Kohnke (2009), we write the current stress state, the so-called equivalent stress, as

$$\begin{eqnarray} &&\sigma _e = \alpha I_1 + \sqrt{J_2}, \end{eqnarray} \tag{ 7 }$$

where I₁ and J₂ are obtained using the current stress state. The yield stress is given as

$$\begin{eqnarray} &&\sigma _y = k, \end{eqnarray} \tag{ 8 }$$

where k is given in Equation (5). Then, the stress ratio is obtained as

$$\begin{eqnarray} &&N = \frac{\sigma _e}{\sigma _y}. \end{eqnarray} \tag{ 9 }$$

If N = 1, an element in the body is in a plastic state.

Figure 1 shows a plastic solution for the case of a bulk density of 1.0 g cm⁻³. Figure 1(a) gives the total deformation vectors of the shape in meters. Figure 1(b) shows the stress ratio over the cross section along the x₁ and x₃ axes. The cohesion for this case is 75 Pa. The stress ratio around the center region is over 0.99, while the region with a stress ratio of 0.9 is much wider than that of 0.99. This indeed indicates that plastic deformation occurs in the center of the interior first. This plastic deformation comes from the fact that the strong centrifugal force acts in the horizontal plane, while only the gravitational force acts along the x₃ axis. Thus, the internal core experiences strong shear, equivalent to half of a difference between the principal components of the stress state, leading to a large scale of plastic flow. Specifically, the stress components in the x₁ and x₂ axes may be strongly tensile due to the fast spin period, while the component in the x₃ axis is always compressive. For this reason, plastic deformation occurs in the internal core and so the shape becomes more oblate than the current shape. On the other hand, Figure 1(c) indicates the stress ratio on the surface. The stress ratio of almost all the surface area is below 0.8, meaning that the near-surface region is below the yield condition. This implies that for the case of uniform material distribution, the central core of 1950 DA is the most sensitive part to structural failure, while the surface is not. Based on an associated flow rule, plastic flow must always be accompanied by an increase in volume (Chen & Han 1988). Hence, since the mass is conserved over the granular flow, the porosity of the internal core should be lower than that of other regions.

Zoom In Zoom Out Reset image size

We calculate finite element solutions for three different bulk density cases: 1.0 g cm⁻³, 1.7 g cm⁻³, and 2.4 g cm⁻³. Figure 2 gives the variation of the lowest cohesion with respect to different bulk densities by the solid line and the dark markers. The resolution of the lowest cohesion is 2 Pa. The shadow area shows the cohesion area in which the original shape can remain, while the white region indicates the prohibited cohesion for the existence of 1950 DA. The lowest cohesion ranges from 75 Pa to 85 Pa. It is also found that the lowest cohesion increases as the bulk density becomes higher. This comes from the fact that as the bulk density increases, higher body forces induce higher shear and normal stresses, requiring higher cohesion for keeping the current shape.

Zoom In Zoom Out Reset image size

The sensitivity of the internal core of an oblate body to failure is a new insight about the internal structure of an asteroid. Earlier studies using a volume averaging method (e.g., Holsapple 2007; Hirabayashi & Scheeres 2014) could only derive an upper bound condition for structural failure of an asteroid, while the present technique can consider a detailed deformation mode as well as the precise condition for it. For the case of 1950 DA, since the spin state is considered to be near its failure condition, we find that at the lowest cohesion for keeping its current shape, the internal core fails structurally, while the surface region does not. Because of this, horizontal outward deformation occurs in that plastic region and pushes the outer region to build the equatorial ridge. This result implies that if cohesion of the material of 1950 DA were constant over the whole volume, the formation of the equatorial ridges would result from plastic deformation of the internal core and not from a landslide. The present failure model for 1950 DA predicts that the core of the asteroid should be "underdense." This is the opposite of the prediction if the surface regolith flows to the equator to form the bulge—then the bulge is underdense. This is a distinction that can be specifically tested for future asteroid exploration missions such as the OSIRIS-REx mission, which should be able to distinguish whether the core is underdense or overdense relative to the rest of asteroid (101955) Bennu (Nolan et al. 2013; Takahashi & Scheeres 2014).

From the current finite element model, it is unlikely that surface failure occurs. As seen in Figure 1(c), the surface region is always below the yield condition. This implies that although the centrifugal force exceeds the gravitational force, the surface region does not fail structurally and so the surface particles will not fly off because of its cohesion. Here, we examine whether or not a small boulder is stable in the case when the cohesion obtained above represents the material strength. This condition can be obtained by considering the force balance on a boulder at the equatorial region. Here, to simplify this analysis, we assume that 1950 DA and the boulder are spherical bodies. The masses of 1950 DA and the boulder are denoted as M and m, respectively. The radii of 1950 DA and the boulder are defined as R and r, respectively. The bulk density of 1950 DA and the material density of the boulder are given as ρ_B and ρ_G, respectively. Also, ω is the current spin rate of 1950 DA. In order for the boulder to separate, the centrifugal force acting on the boulder should exceed the gravitational force, which is given as

$$\begin{eqnarray} &&m (R + r) \omega ^2 \ge \frac{G M m}{(R+r)^2} + f_{{\rm cohesion}}. \end{eqnarray} \tag{ 10 }$$

The second term on the right-hand side indicates a cohesion force. Here, using the results in Sánchez & Scheeres (2014), we simply describe this term as ∼2cπr², where c represents cohesion and 2πr² is the apparent area of the boulder touching the surface. We assume that the radius of the boulder is very small compared to that of 1950 DA, i.e., r ≪ R. By choosing ρ_B  ∼  2.4 g cm⁻³, an upper bound for the bulk density of 1950 DA, ρ_G ∼ 6 g cm⁻³, from the fact that 1950 DA may be composed of a nickel-iron or enstatite chondritic composition (Busch et al. 2007), and R = 649 m, the mean radius of 1950 DA, we obtain the lowest cohesion to retain the boulder on the surface as

$$\begin{eqnarray} c_{{\rm min}} &\sim & \frac{2}{3} \rho _G r R \omega ^2 - \frac{8 \pi }{9} \rho _G \rho _B G r R, \\ &=& 0.015 \: {\rm [Pa/m]} \: r, \nonumber \end{eqnarray} \tag{ 11 }$$

where r is the radius of the boulder. Thus, even if the size of the boulder is on the order of a few hundred meters, the lowest cohesion to retain the boulder is on the order of a few pascals. This implies that since the necessary cohesion that prevents 1950 DA from failing structurally is ∼80 Pa, if this is representative of the cohesion within the regolith, any embedded boulders will be stable (Sánchez & Scheeres 2014). This simplified analysis also supports the primary result in this study that the internal core is more sensitive to structural failure than the surface region.

Finally, we emphasize our contributions. Rozitis et al. (2014) used Holsapple's limit analysis technique (Holsapple 2007) to derive the minimum cohesion of 1950 DA as 44 − 76 Pa. However, since this technique was based on the upper bound theorem (Chen & Han 1988), the cohesion derived by this technique indicates the lowest limit for the necessary condition only. Our finite element model is capable of giving a more precise failure condition and possible failure modes of an asteroid. We note that for this asteroid only a lower bound on strength can be given, as it is not currently seen to be failing. This is different from the study by Hirabayashi et al. (2014) who obtained both lower and upper limits for the cohesion of P/2013 R3 from its disruption event. In addition, we clearly show that the failure mode for this body is for the interior to fail. This leads to a new prediction for the central core of the body to be underdense. Future analysis will explore this in more detail.

M.H. deeply appreciates Eleanor Matheson for her dedicated review of grammar in the current paper. M.H. wishes to thank Keith A. Holsapple at University of Washington and Carlos A. Felippa at University of Colorado for their useful advice about development of a finite element model for plastic deformation. M.H. also acknowledges Ben Rozitis at University of Tennessee for constructive discussions about their technique used in Rozitis et al. (2014). Finally, M.H. appreciates Paul Sánchez for useful discussion about the interpretation of our study. This research was supported by grant NNA14AB03A from NASA's SSERVI program.