Diurnal temperature variation as the source of the preferential direction of fractures on asteroids: theoretical model for the case of Bennu
DDiurnal temperature variation as the source of thepreferential direction of fractures on asteroids:theoretical model for the case of Bennu
D. Uribe a,b, ∗ , M. Delbo b , P.-O. Bouchard a , D. Pino Mu˜noz a a Centre de Mise en Forme des Mat´eriaux (CEMEF)MINES ParisTech, PSL Research UniversityCNRS UMR 7635, CS 10207 rue Claude Daunesse, 06904 Sophia Antipolis Cedex, France b Universit´e Cˆote d’Azur, Observatoire de la Cˆote d’AzurCNRS-Lagrange, CS 34229, 06304 - NICE Cedex 4, France
Abstract
It has been shown that temperature cycles on airless bodies of our Solar System can causedamaging of surface materials. Nevertheless, propagation mechanisms in the case of spaceobjects are still poorly understood. Present work combines a thermoelasticity model togetherwith linear elastic fracture mechanics theory to predict fracture propagation in the presenceof thermal gradients generated by diurnal temperature cycling and under conditions similarto those existing on the asteroid Bennu. The crack direction is computed using the maximalstrain energy release rate criterion, which is implemented using finite elements and the so-called Gθ method (Uribe-Su´arez et al. 2020. Eng. Fracture Mech. 227:106918). Using theimplemented methodology, crack propagation direction for an initial crack tip in differentpositions and for different orientations is computed. It is found that cracks preferentiallypropagate in the North to South (N-S), in the North-East to South-West (NE-SW) and inthe North-West to South-East (NW-SE) directions. Finally, thermal fatigue analysis wasperformed in order to estimate the crack growth rate. Computed value is in good agreementwith available experimental evidence. Keywords:
Asteroid Bennu, Thermoelastic model, Crack propagation Direction, EnergyRelease Rate, Thermal fatigue crack growth
The formation of fractures and their growth are key processes in man-made structures andmaterials. The large majority of the cases studied in the literature are on the aforementionedstructures, given their importance in the everyday life of humans. These processes also occurin natural objects, such as rocks, boulders, and cliffs [1, 2, 3, 4, 5, 6] and are documented ∗ Corresponding author.E-mail address: diego.uribe [email protected] (D. Uribe Su´arez). a r X i v : . [ a s t r o - ph . E P ] F e b n several objects of our solar system, including Earth [7, 8], Mars [6, 9], our Moon [10, 11],the nuclei of comets [12, 13, 14], asteroids [15, 16, 17, 18, 19], and meteorites [20].The source of the driving forces that provokes crack nucleation and propagation canbe very diverse, ranging from unloading of the pressure stresses under which certain rocksformed in the deep crust of Earth, tectonic stresses, rapid mechanical stresses from impactsand thermal stresses. In the latter case, the presence of water can enhance the crackingphenomena via the known freeze-thaw effect [see e.g. 21, and references therein]. The ef-fectiveness of thermal stresses in cracking rocks and other geological units in the absenceof water has been long debated: a famous laboratory experiment by Griggs [22] arguedagainst earlier claims that rocks could be fractured by temperature variations only, the pro-cess that in general and hereafter is called thermal cracking. However, thermal crackinggained momentum recently and came to great attention to planetary scientists thanks tonew measurements, modelling, and observations.More in details, it has been shown that temperature variations resulting from the cyclesbetween day and night can damage materials on airless bodies of our Solar System [for ref-erences on this topic, see the introductions of 23, 24]. This damaging process consists in thenucleation and growth of micro-fractures inside the material due to the mechanical stressesinduced by the diurnal temperature cycles. In general the mechanical stresses resulting fromtemperature gradients due to the day and night cycles are smaller than the strength of thematerial [see e.g. 20, 25, 26, 27, 23]. In this case the crack can still open and grow in aregime that is said to be sub-critical [28]; the material is increasingly damaged at each cycleand it is usually spoken of thermal fatigue. Eventually, the application of a large numberof cycles can produce important crack growth; the crack tip can reach a boundary of thematerial, such as a discontinuity or the edge of a rock; this will cause a rapid transition fromsub-critical to critical failure and lead to material failure [see e.g. 29, 20, for the case of as-teroids and meteorites]. For example, Liang et al. [30] studied volumetric stress distributionin an L6 ordinary chondrite’s microstructure subjected to thermal and mechanical loadingsthrough the combination of experiments and micromechanical models. It was found thatunder thermal cycling, the stress concentrates more uniformly along with particle interfaces.The authors interpret that thermal fatigue crack propagation could result in the debondingof particles from the surrounding matrix.Hence, on solar system bodies without an atmosphere thermal fatigue of surface rocks,in addition to the impact of micrometeorites, can eventually lead to rocks’ breakup andproduce regolith [29, 20], the latter being the layer of unconsolidated material that coversplanetary surfaces [31, 32, 33]. In the case of the near-Earth asteroid (101955) Bennu [16],Molaro et al. [34] propose that thermal cracking is also able to eject sub-cm-sized particleaway from the asteroid surface, thereby offering an explanation for the observed activityof this asteroid [35]. Furthermore, it is also proposed by several studies that macroscopicfractures, mass-wasting, and material breakdown on asteroids and cometary nuclei could beexplained as a consequence of thermal effects [29, 14, 36, 12, 19].For all the reasons above, thermal fatigue cracking is now considered a space weatheringmechanism. On the other hand, direct evidence of thermal cracking on asteroids (and comets)is still relatively scarce (but strongly growing) and the details of the process in terms ofspatial and temporal scales still poorly understood. One of the first studies that invokes thisphenomenon to explain certain in situ asteroid observations is the work of [29].Using images2btained by NASA’s Shoemaker mission of the surface of the asteroid (433) Eros, theseauthors noted boulders that appear to break and erode in place, producing fragments thatfill the inside of craters, creating characteristic “ponds” of regolith. Another observationalevidence, obtained from images of NASA’s OSIRIS-REx mission [37], is constituted by thedetection of exfoliation sheets on some of the boulders on the asteroid (101955) Bennu[19]. The thickness of the exfoliation sheets is consistent with the depth inside boulders atwhich thermoelastic simulations show stress concentration as a result of diurnal temperaturevariations [19]. However, it has been shown on Earth [38] and Mars [6] that one of themost diagnostic observations of thermal cracking induced by diurnal temperature variationsis a preferential meridional direction (north to south) of the fractures on surface rocks. Thereason is very simple: during the day the Sun moves in the sky from the east to the west.As a consequence, the temperature gradients are directed essentially in the same direction(west to east). Fractures mainly propagate in a direction perpendicular to that of maximumprincipal stress. Therefore it is expected that this direction of propagation is essentially fromthe north to the south, when they are driven by these diurnal temperature cycles. There arehints of predominance of fracture directed in the north to the south and in the north-westto the south-east on the boulders of the asteroid Bennu [39]. However, a modelling of thecrack propagation direction in conditions similar to those existing on Bennu is still lacking.The aim of this work is to provide theoretical foundation for analysis and interpretation offracture directions on small asteroids with properties similar to those of Bennu. For oursimulations, a thermo-mechanical model will be used to predict stresses on Bennu’s bouldersand an appropriate crack propagation method based on principles of fracture mechanics.The paper is organised as follows. In Section 2, the developed thermomechanical modelis presented. This section also reviews crack growth direction theory using an energeticapproach for linear thermoelasticity fracture mechanics problems. One simple example ispresented in Section 3 to show the accuracy of the proposed method when calculating thecrack propagation direction due to the presence of thermal gradients. The results are com-pared against observed crack propagation directions on asteroid Bennu. In section 4 thecomputed crack propagation directions are discussed. Finally some concluding remarks arepresented in Section 5. Primary goal of this work is to compute the direction of crack propagation due to thermalstrain in a geometry corresponding to a typical boulder on the surface of the asteroid (101955)Bennu, the target of NASA’s sample return mission OSIRIS-REx. Here it is made thehypothesis that most of the fractures observed on the surface of the boulders by [16], [18],and [39] are due to the growth of surface cracks. [34] also discuss stresses inside bouldersthat could cause cracking deep inside the rock mass.The geometry of our problem is schematised in Fig. 1, where it is shown a cubic-likeboulder extruding from the equator of the asteroid, of which only about half of the equatorialbelt is simulated.The mesh is divided in triangular facets. A thermophysical model [40] is used to calculatethe temperatures of all facets as a function of time, as described in the following sections.3 astWest North
Rotation direction
EastWestWest view East view a”a’ E a”a’W
Figure 1: Schematic representation of the cubic-like boulder extruding from the equatorof the asteroid, of which only about half of the equatorial belt is simulated to obtain thetemperature distribution on faces E and W .The temperatures of the boulder facets are then used in a thermoelastic model in order tocompute the strain and stresses as a function of the position in the boulder and time. Next,fracture mechanics theory is used to estimate the propagation direction of a tiny notch thatis placed on the horizontal (a’, a”) face of the cubic boulder. This is done by computing thestrain energy release rate , which is a well-known fracture mechanics parameter [41]. Theimplemented thermoelastic model as well as the used methodology for the computation ofthe energy release rate ( G ) are explained in the following sections. The first step of the presented method consists in using a well-established thermophysicalmodel [42, 40] to solve the one-dimensional heat diffusion problem. Temperature is calculatedas a function of time for all the surface elements of the mesh depicted in Fig. 1. Boundaryconditions are given by the variable day/night illumination including the shadows cast bythe local terrain of the mesh on itself, radiation of the heat in space, conduction in thesub surface, and mutual heating [43]. The physical parameters of the material used in thiswork are given in Table 1. These properties were taken from [20]. They correspond to theCarbonaceous Chondrite sub-type CM2 Murchison meteorite. This is considered to be agood analog of asteroids belonging to the C-complex broad spectroscopic class [44]. Theasteroid Bennu also belong to the C-complex [45].
Delbo et al. [39] observed and mapped cracks on boulders on the surface of the asteroid Bennuusing OSIRIS-REx images that can be approximated as parallel to the local surface of the The strain energy release rate represents the change of elastic strain energy per unit area of crackextension. P s 15,469.2 [1]Bulk modulus, K MPa 29,000 [2]Shear modulus, µ MPa 18,000 [2]Young modulus, E MPa 44,742.857 [2]Poisson’s ratio, ν ρ kg m − λ W mm − K − × − [2]Thermal expansion coefficient, α K − × − [2]Heat capacity, c J kg − K −
500 [2]Reference temperature, T ref K 250 [3]Paris pre-factor, C m [MPa √ m] − n × − [2]Paris exponent, n (cid:15) ), which can be decomposed into an elastic, plastic and thermal part, as it isshown in equation (1): (cid:15) ij = (cid:15) elasticij + (cid:15) plasticij + (cid:15) thermalij ; (1)where (cid:15) ij is the total strain tensor, (cid:15) elasticij , (cid:15) plasticij and (cid:15) thermalij are respectively, the elastic,plastic and thermal strain tensors. In this work it is assumed that there is no plastic strain( (cid:15) plasticij = 0). Therefore, equation (1) simplifies into: (cid:15) ij = (cid:15) elasticij + (cid:15) thermalij ; (2)where thermal strain is defined as: (cid:15) thermalij = α ∆ T δ ij , where ∆ T = T − T ref (3)In equation (3), α is the thermal expansion coefficient [ K − ], ∆ T is the difference betweenthe computed temperature and the reference temperature, latter one being the temperaturewhere there is no strain. Finally δ ij is the Kronecker delta ( δ ij = 0 for i (cid:54) = j, δ ij = 1 for i =j). It is assumed that the reference temperature is the average of the temperatures given inSection 3. It is also assumed an isotropic thermal expansion coefficient.In the presented method a weak thermomechanical coupling is assumed. Which meansthat the temperature is initially obtained from the heat problem and then introduced intothe mechanics computation. This is possible since the characteristic time scale of the thermalproblem is several orders of magnitude greater than the characteristic time scale of the crackpropagation problem. For the solution of a thermoelastic problem, first, the temperaturedistribution inside the body is computed by solving the heat transfer problem. Then, theresulting temperature distribution is input to the mechanical problem as an initial strain. Inelasticity, the Hooke’s law for homogeneous and isotropic materials is defined by Equation(4): 6 ij = λδ ij (cid:15) kk + µ ( (cid:15) ij + (cid:15) ji ) (4)where σ ij is the stress tensor and λ and µ are respectively, the Lam´e’s first and secondparameters. When thermal strain is included, Equation (4) becomes: σ ij = λδ ij ( (cid:15) kk − α ∆ T ) + 2 µ ( (cid:15) ij − α ∆ T δ ij ) (5) Crack propagation requires energy. The amount of energy released during the fractureprocess is known as energy release rate ( G ). Commonty, G is widely used in the literaturein order to find the crack propagation direction for a given configuration. To compute thecrack propagation direction, a criterion based on the azimuthal distribution of the energyrelease rate around the crack tip is used. Multiple crack propagation directions are testedand then select the one that maximises G . Namely, in the plane of the a-face, the azimuthof the maximum energy release rate is identified, as it is known that crack propagation willtake place in that direction. This method was tested and validated in Uribe-Su´arez et al.[47].The G -value is evaluated using all virtual and kinematically admissible crack lengthdisplacements. The direction of crack propagation ( θ ) can be determined by: (cid:0) dGdθ (cid:1) θ = θ = 0 , (cid:16) d Gdθ (cid:17) θ = θ ≤ . (6)In equation (6), according to the work done by Erdogan and Sih [48], there is a limit anglecorresponding to pure shear: θ = ± . ◦ . In this work, for the computation of the energyrelease rate ( G ), the numerical technique known as Gθ method [49] is used. The Gθ method’simplementation is quite simple and multiple extensions are available.The strain energy release rate is the decrease in the total potential energy ( w p ) duringa growth of crack area ( dA ). To determine the variation of the total potential energy incracked solid Ω, F ε is defined as an infinitesimal geometrical perturbation ε in the vicinityof the crack tip: F ε : R → R ∀ M ∈ Ω , F ε ( M ) = M ε = M + ε V ( M ) (7)where the virtual field V gives the location of each point of the perturbated solid using itsinitial position ( M ) before the perturbation. When the perturbation ε is sufficiently small,Destuynder et al. [49] showed that the stress field ( σ ) and the displacement field ( u ) on theperturbed configuration may be expressed as: σ ε = σ + εσ u ε = u + εu (8)7here σ and u are the first order variations of the stress and displacement fields duringthe infinitesimal perturbation ε on Ω. The total potential energy variation during a crackextension is then obtained when ε aims towards 0: dW p da = lim ε → W εp − W p ε (9)The virtual displacement field V representing the virtual kinematics of the crack has thefollowing properties: • V is parallel to the crack plan. • V is normal to the crack front. • The support of V is only needed in the vicinity of the crack. • (cid:107) V (cid:107) is constant in a defined region around the crack tip.When implementing the Gθ method, first, it should be defined two contours C and C around the crack tip. Those contours divide this region into three domains C int , C ring and C ext as it is shown in Fig. 3.Figure 3: Contours and domains used to computed G using the Gθ method.The virtual displacement field, V ( v , v ), is defined by equation (10). V = v = (cid:0) − ABAC (cid:1) cos( θ ) v = (cid:0) − ABAC (cid:1) sin( θ ) (10)Where O is the crack tip, B is an integration point belonging to the ring, A and C are theintersections between OB and C and C respectively (inside and outside contours of thering); and θ is the virtual direction of propagation measured with respect to the crack axis.Propagation direction is defined as a function of an angle θ . Hence the values of V in thethree domains are: • The norm of V in C int is constant and equal to 1.8 The norm of V in C ring varies continuously from 1 to 0. • The norm of V in C ext is equal to 0.When there is neither thermal strain nor load applied directly to the crack faces, thetotal potential energy may be expressed as: W p = 12 (cid:90) Ω σ ij U j,i d Ω − (cid:90) Ω f i U i d Ω (11)where σ ij is the stress field, U i is the displacement field, and f i the external loads. Whenan infinitesimal perturbation ε takes place, derivatives and integrals on the perturbed partcan be expressed using a first order “limited development” of operations related to the non-perturbed part. Due to this, the energy release rate may be expressed as: G = (cid:90) ring (cid:18) σ ij U j,k V k,i − σ ij U j,i V k,k (cid:19) dA ring (12)where V i is the virtual displacement field, U j,k is the gradient of the displacement field, V k,i is the gradient of the virtual displacement field, V k,k is the divergence of the virtualdisplacement field and A ring is the integration region. Additionally, as V i varies only inside C ring , the integration will only be performed over C ring ( A ring ). In Fig. 4- a , elementsbelonging to C ring are marked with a green dot. A discrete set of θ values will be usedfor a virtual crack propagation. The θ values are selected inside the range [ − ◦ , ◦ ] [48]with respect to the crack axis. For each one of the θ values, an associated G value will becomputed. In Fig. 4- a the virtual ring in which all the integration points are consideredis depicted. Fig. 4- b shows that G can be computed for each value of θ , which makesthe identification of the θ value that maximises G ( θ ) straightforward. Fig. 4- b shows thecurve G ( θ ) for a crack where its crack axis is equal to 90 ◦ (E), in this case the direction ofpropagation is equal to 120 ◦ (NW - SE).Figure 4: a) Ring of elements around the crack tip. b) G ( θ ) curve for the maximun energyrelease rate criterion for a crack where its crack axis equal to 90 ◦ .Furthermore, due to the presence of thermal strain in this work, the inclusion of additionalterms in the Gθ method is necessary. In elasticity according to Hooke’s law, the stress andstrain tensors can be defined as follows: 9 ij = C ijkl (cid:15) kl (13)where C ijkl is the fourth-order elastic constitutive tensor. Combining equations (13), (2) and(3) leads to the definition of the stress tensor taking into account thermal strain: σ ij = C ijkl ( (cid:15) kl − α ∆ T δ kl ) (14)When thermal strain is accounted for, equation (8) is still valid. Therefore, according toSuo and Combescure [50], Brochard and Suo [51], the equation for the Gθ method in thepresent of thermal strain is given: G = (cid:90) ring (cid:18) σ ij U j,k V k,i − σ ij ( U j,i − α ∆ T δ ij ) V k,k + ασ ii T ,j V j (cid:19) dA ring (15)where α is the thermal expansion coefficient, ∆ T (∆ T = T − T ref ) is the difference betweenthe computed temperature and the reference temperature, δ ij is the Kronecker delta and T ,j is the gradient of the temperature. A schematic representation of the sequential solution ofthe thermomechanical problem presented in this work is provided in Fig. 5.Figure 5: Schematic representation for the thermoelastic fracture analysis performed in thiswork. For linear elastic materials, K and G are uniquely related as follows [52]:10 = K I E (cid:48) + K II E (cid:48) + K III µ with E (cid:48) = E for plane stress E (cid:48) = E − ν for plane strain µ = E ν ) (16)where G is the energy release rate, K I , K II and K III are respectively the stress intensityfactors for mode I , mode II and mode III , E is the Young modulus, ν is the Poisson’s ratioand µ is the shear modulus. Due to the fact that present work deals with thermal cyclicstress (tension-compression), only mode I loading is assumed. Once the energy release rate( G ) is known, it is possible to compute the stress intensity factor ( K I ). Assuming planestrain conditions, a relation between G and K I can be stated from equation (16): K I = (cid:114) G E − ν (17)where K I is the stress intensity factor, G is the energy release rate, E is the Young modulusand ν is the Poisson’s ratio. Paris’s Law [53] allows to relate the rate at which the cracklength ( a ) grows as a function of the number of cycles ( N ) to the maximum variation ofstress intensity factor (∆ K I ) over the whole cycle. Paris’s law is presented in equation (18): dadN = C [∆ K I ] n (18)where a is the crack length [ m ], N the number of cycles [ − ], ∆ K I the range of stress intensityfactor [ P a m . ], C [ m ( P a m . ) − n ] and n [ − ] are material properties fitted to experimentalfatigue data.The beam presented in Fig. 2, also referred as the a-face in this work, is embedded in amesh that is different from the one showed in Fig. 1. The aforementioned mesh is an isotropicunstructured triangular one, that is refined in the neighborhood of the crack tip (Fig. 4- a ). The thermoelasticity problem is solved on this mesh through the finite element method(FEM). The developed thermoelastic model has been implemented in CimLib, a C ++ in-house finite element library developed at CEMEF [54]. Regarding the thermal problem, theavailable finite element framework in CimLib is a classic one where the variable temperature( T ) is solved through an implicit formulation [55]. The mechanical problem is solved usingan implicit formulation (mixed) with first-order elements using linear interpolation functionsfor both velocity and pressure. In this formulation both variables are P
1. This elementcombination, linear in both −→ v and p ( P /P inf − sup condition [56,57]. The latter, being a condition related with the mathematical convergence characteristicsof the finite element formulation. The inf − sup condition ensures the solvability, stabilityand optimality of the finite element solution and is crucial in establishing its quality [58].Through the inf − sup condition, a problem can be assessed in order to see if it is well-posedor not. A problem is called well-posed if: (i) a solution exists, (ii) the solution is unique and(iii) the solution depends continuously on the given data (in some reasonable topology) [59].To overcome this issue, an extra degree of freedom (DOF) in the center of each element isadded in −→ v . This extra DOF is condensed and thus written as a function of the remainingDOFs of the element. The use of this extra DOF is a well-known stabilization technique11ften called “bubble” or P + /P
1. This technique satisfies the inf − sup condition, and hasbeen implemented in the CimLib [60, 61]. In this work, the simulations begin by calculating temperatures for the mesh of Fig. 1 usingthe asteroid thermophysical model described in section 2. The asteroid is placed at a distanceof 1.12 au from the sun. This corresponds to the semi major axis of the orbit of the asteroidBennu. The resulting temperatures of the E- and W-faces are shown in Fig. 6. Next, thesetemperatures are input into the thermo-mechanical model. This second model predicts thestress over the domain of the a-face as a function of space and time. Since the presence ofcracks highly affects the stress field, different simulations for different crack positions andorientations for the same temperature fields are run.Figure 6: Temperatures T W (W-face) and T E (E-face) used as thermal boundary conditions.Namely, the propagation direction for an initial crack in different positions inside thedomain of the a-face is computed, as shown in Fig. 2. In addition, for each position of eachinitial crack, its orientation is varied. Namely, when the crack is attached to the W-face thecrack propagation direction is computed for different initial orientations with azimuth anglesranging from 46 ◦ to 134 ◦ . Azimuth angle is increasing clockwise (i.e. from the north to theeast) and is equal to zero when it is pointing vertical up (i.e. to the north). When the crackis attached to the E-face, the crack propagation direction is simulated for different initialcrack orientations having azimuth angles ranging from 226 ◦ to 314 ◦ . Finally, when the crackis attached to the north (top) side of the a-face (also called beam), the crack propagationdirection is simulated for different initial crack orientations with azimuth angles ranging12rom 136 ◦ to 224 ◦ . In all the described cases, the length and the width of the a-face werefixed to 100 mm and 4 mm respectively. These values were chosen after studying: (i) thecharacteristic time of the thermal problem as well as (ii) the minimum required width ofthe a-face in order to avoid the boundary effects on the results and (iii) to avoid a highcomputational cost. Due to the fact that there is no temperature gradient in the verticaldirection in the aforementioned simulations, a larger height of the domain does not affectinitial crack propagation direction. Similarly, this work does not intend to evolve the crackwith time.The results are presented using the so-called windrose diagrams, which can be consideredas circular histograms that represent the distribution of the computed crack propagationdirections. It should be noticed that for a rotation of 180 ◦ of the circular histograms (win-droses) these are identical. This is due to the fact that crack propagation direction (azimuthalangle) has to be independent of the beginning and ending point of a crack, i.e., for each oneof the computed crack propagation directions θ , in the histograms there are also count thedirection θ + 180 ◦ . For the scenario where crack propagation direction is computed for aninitial crack tip placed in different positions inside the domain of the a-face and for differentcrack orientations, the windroses are presented separately.Figure 7a shows the windrose diagram for the case of a crack attached to the W-face.Computed crack propagation directions are preferentially oriented in the North-West toSouth-East (NW-SE) direction, regardless of the initial different orientations of the crackstaken into account. Furthermore, a minimum number of computed crack propagation direc-tions are oriented East (E) to West (W). In Fig. 7b it is shown the windrose diagram for thecase of a crack attached to the E-face. In this case the preferential orientation of the com-puted crack propagation directions is in the North-East to South-West (NE-SW) direction.There are also few cracks propagating in a direction aligned East (E) - West (W). The lastcase, where the crack is attached to the north (top) is shown in Fig. 7c. According to thewindrose diagram, the preferential orientation of the computed crack propagation directionsis in the North to South (N-S) direction. (a) W-face (b) E-face (c) north (top) Figure 7: Windrose diagrams for different crack orientations and cracks attached to a) theW-face, b) the E-face and to c) the north.Finally in Fig. 8 it is presented a windrose diagram that gathers all the cases described13bove. In summary, the distribution of all the computed crack propagation directions arepreferentially oriented in a higher concentration in the North to South (N-S), in the North-West to South-East (NW-SE) and in the North-East to South-West (NE-SW) directions asit was already described.Figure 8: Windrose diagram gathering the cases where crack is attached to the W-face, tothe E-face and to the north (top) for different crack orientations.It is noted that the cracks simulated up to now have a length considerably smaller than theones observed on the boulders on the surface of Bennu [19, 16, 18, 17, 39]. In the following,the crack propagation direction is studied as a function of the initial crack size. From thedifferent configurations shown in Fig. 2, it was selected the case for the crack attached tothe E-face. For this, crack propagation direction is computed by placing the crack tip at 3different positions inside the domain of the a-face; for each one of this positions 3 differentlengths of the a-face were used; for each one of these cases the orientation of the crack wasvaried as it was done in the previous simulations.we varied the orientation of the crack as it was done in the previous simulations.Figure 9 shows that for a small crack length ( | sin azimuth | mm ), increasing the length of thea-face does not play an important role in the distribution of the computed crack propagationdirections. Computed crack propagation directions are preferentially oriented in the North-East to South-West (NE-SW) direction, with a few amount of computed crack propagationdirections going to the East (E) and to the West (W).In Fig. 10, it is possible to see that when increasing the length of the a-face for crackswhose lengths are | sin azimuth | mm , the distribution of the computed crack propagation direc-tions oriented in the North to South (N-S) direction decreases, while the cracks oriented inthe North-East to South-West (NE-SW) direction increase. It should also be noted that the14 a) a-face length: 200 mm (b) a-face length: 500 mm (c) a-face length: 1000 mm Figure 9: Windrose diagrams for a crack attached to the E-face when crack length equals to | sin azimuth | mm and the length of the beam was varied for different crack orientations.amount of cracks directed to the East (E) and the West (W) decreases when increasing thelength of the a-face from 500 mm to 1000 mm . (a) a-face length: 200 mm (b) a-face length: 500 mm (c) a-face length: 1000 mm Figure 10: Windrose diagrams for a crack attached to the E-face when crack length equalsto | sin azimuth | mm and the length of the beam was varied for different crack orientations.Finally, in Fig. 11 it is presented the windrose diagrams for a configuration where cracklength is equal to | sin azimuth | mm . Taking into account that the length of the a-face took thevalues of 200 mm , 500 mm and 1000 mm , for this last case, the configuration of the a-facewith a length equal to 200 mm was not simulated. Fig. 11 shows that when increasingthe length of the a-face from 500 mm to 1000 mm , a redistribution of the computed crackpropagation directions takes place. The initially cracks preferentially oriented in the North-East to South-West (NE-SW) direction changed their orientation into the North-West toSouth-East (NW-SE) direction.Figure 12 shows the stress field through the maximum principal stress for three crackswith different orientations: attached to the W-face and directed to the East, attached tothe north face of the boulder and directed to the South, and attached to the E-face of the15 a) a-face length: 500 mm (b) a-face length: 1000 mm Figure 11: Windrose diagrams for a crack attached to the E-face when crack length equalsto | sin azimuth | mm and the length of the beam was varied for different crack orientations.boulder and directed to the West. It is worth noting that in all the cases, the computedcrack propagation direction (dashed black line), the one that maximises G , is perpendicularto the direction in which the maximum principal stress takes place (dashed white line). Forall cases, the maximum stress near the crack tip over a temperature cycle exceeds some MPa,which is a significant fraction of the typical strength of carbonaceous chondrites and evencomparable to the strength of boulders on Bennu [62]. This indicates that cracks propagateunder these stresses. Although it should be noted that, theoretically, the stress at the cracktip is supposed to be almost infinity ( ∞ ), it depends a lot on the mesh size. The use ofquarter-point elements [63] is recommended for better accuracy of the stress field.To study whether the crack is propagating it can also be used the typical approach ofthermal fatigue and determine the crack growth length per cycle using the Paris’ law. To doso, the configuration in which the crack is attached to the W-face was chosen. Additionally,it was selected an orientation of the crack (azimuthal angle) equal to 90 ◦ , i.e. pointing to theeast. In this case, as the maximum energy release rate ( G max ) is known, and assuming planestrain conditions, it is possible to compute the maximum stress intensity factor ( K Imax ) usingequation (17). In this case, G max in the thermal cycle is equal to 5 . × − [ M P a . mm ],therefore, using equation (17), K I is equal to 4 .
95 [
M P a . mm . ]. Due to the fact that inthis work it is dealt with stresses generated by thermal cycles, the crack tip is subjected toboth tensile and compressive stresses over a full cycle. Thus, ignoring crack closure [64], thelowest stress intensity factor experienced by the crack tip is simply zero. According to this, itfollows that ∆ K I is equal to 4 .
95 [
M P a . mm . ]. For the case where constant crack growthis assumed, it can be computed that for the presented thermal cycle, the crack growth rateis 2 . × − [ mm . cycle − ], which is of the order of 0 . mm . yr − ] (or ∼ . m ] inthousand years), in good agreement with previous studies [20]. This indicates that cracks16igure 12: Maximum principal stress ( M P a ) for a crack attached to a) the W-face, b) theE-face and to c) the north (top). In all cases, the dashed black line represents the computedcrack propagation direction, while the dashed white line represents the axis associated withthe maximum principal stress.would propagate on Bennu’s boulders solely due to diurnal thermal stresses. It is cautionedhere that crack propagation is a non linear process, where the rate of propagation is amongother parameters a function of the position of the crack tip with respect to the temperaturegradient. This means that a full simulation of the crack growth from the beginning untilits size is comparable to the size of the hosting boulder should be carried out in order toestimate the time required to fracture the boulder, which is beyond the scope of this paper. We begin to notice that in Fig. 7a, 7b and 7c, the range of computed crack propagationdirections is preferentially concentrated in some sectors of the windrose diagrams. Thiscould be explained by the fact that when a crack orientation is defined, the range of possible17rack propagation directions, is determined by the limit angle corresponding to pure sheardefined in the work of [48]. This value is approximately ± ◦ with respect to the crack axisorientation. It should also be noted that, the crack was only placed at three different placesinside the domain, i.e. attached to the East, the West faces and attached to the North rim ofthe domain. This means that not all the angles belonging to the range [0 ◦ , ◦ ] were takeninto account. For example, in Fig. 7b, i.e., for the case of a crack attached to the E-face,the computed crack propagation directions range from 202 . ◦ to 292 . ◦ . From this figure, itis worth noting that most of the computed direction are oriented to the South-West (SW),and due to the already explained symmetry of the wind rose, the preferential orientationof the computed crack propagation directions is in the North-East to South-West (NE-SW)direction.When plotting the windrose diagrams for the three evaluated cases (crack attached tothe E-, W-face and the north rim of the domain) in the same graph, as shown in Fig. 8,the sector grouping the majority of the computed directions is the one oriented in the Northto South (N-S) direction, followed by the groups oriented in the North-West to South-East(NW-SE) and in the North-East to South-West (NE-SW) directions, respectively. Lookingat the results of the three simulated cases, it is worth noting that cracks on asteroid Bennupropagate mainly in the North-South (N-S) direction.In Fig. 12 it should be noticed that the computed crack propagation directions in thethree cases presented there, are perpendicular to the axis of the maximum principal stresses,as stated by the literature. This makes us feel confident about the results coming out fromthe coupled thermoelastic model with the linear elastic fracture mechanics approach that wasimplemented in this work. Additionally, it is worth mentioning that our model to computecrack propagation direction was validated in previous work [47].From Fig. 9 it can be noticed that for short crack lengths, the increasing of the lengthof the a-face does not highly affect the tendency of the distribution of the computed crackpropagation directions. Cracks mainly continue propagating oriented in the North-East toSouth-West (NE-SW) direction. When increasing the length of the a-face from 500 mm to1000 mm , a few amount of cracks changed its orientation into the North-West to South-East(NW-SE) direction. It tells us that even if the temperature gradient along the a-face changes,cracks of short length are not affected.According to Fig. 10, when crack length increases, a redistribution of the computeddirections takes place. It is possible to see two main groups, one oriented almost in theNorth to South (N-S) direction and another one oriented almost in the West to East (W-E)direction. Meanwhile, for the case when the crack length is very long compared with thelength of the a-face, as depicted in Fig. 11, a marked tendency is clearly observed. Computedcrack propagation directions are preferentially oriented in the North-West to South-East(NW-SE) direction, which is consistent with the observations performed by Delbo et al. [39].Taking into account that the duration of the thermal cycle used in this work is 4.3 h ,and assuming a constant crack propagation rate, the crack growth rate computed here,whose value is 2 . × − [ mm . cycle − ], is equivalent to a crack growth rate of about0 . mm . year − ]. This value is not quite different from the one measured by Delbo et al.[20] through laboratory experiments on two meteorites: a carbonaceous chondrite and anordinary chondrite. However, as stated in Delbo et al. [20], the crack propagation rate istypically a nonlinear function between the crack size, the maximum variation of the stress18ntensity factor over the whole cycle and material properties. However, considering thatthis work is a first attempt to describe crack growth due to thermal fatigue, the results arepromising. An approach to compute 2D crack propagation direction under the presence of thermal gra-dients has been presented. It combines thermoelasticity and linear elastic fracture mechanicstheories in order to compute crack propagation directions driven by diurnal temperature cy-cling and under conditions similar to those existing on the asteroid Bennu. Through thedifferent scenarios simulated in this work, the robustness and the accuracy of the proposedapproach in terms of crack propagation direction is shown.It is found that the distribution of the computed crack propagation directions for differentconfigurations on a simple domain simulating Bennu’s boulders, show a preferential directionfrom North to South (N-S) and from North-East to South-West (NE-SW) for shorter cracks.For longer cracks, the preferential direction is from North to South (N-S) and from North-West to South-East (NW-SE). This conclusion is supported by observations performed by[39].According to the results obtained by means of our thermomechanical model combinedwith a well-known fatigue model (Paris’s law) and assuming constant crack growth length, itis found that on asteroid Bennu cracks grow at rate approximately equal to 0 . mm . year − ].This is a value very similar to the one found on experiments perform in the work of Delboet al. [20], giving us a good insight of the promising approach proposed.However, it is cautioned that crack propagation is a highly non-linear process. In orderto estimate the time to fracture Bennu’s meter-sized rocks a full-blown simulation shall beperformed to calculate the crack growth rate at different crack growth stages. This is reallyimportant because crack growth rate highly depends on the geometry and on the stressintensity factor at the crack tip, and those will change as the crack tip evolves trough thedomain. Acknowledgements
We acknowledge support from Academies of Excellence on Complex Systems and Space, En-vironment, Risk and Resilience of the Initiative d’EXcellence ”Joint, Excellent, and DynamicInitiative” (IDEX JEDI) of the Universit´e Cˆote d’Azur and the Center for Planetary Origin(C4PO) A Material approach ( ). Wealso acknowledge support from The DIGIMU ANR industrial chair ( https://chaire-digimu.cemef.mines-paristech.fr/ ) as well as support from the French space agency CNES.
References [1] R.-H. Cao, P. Cao, H. Lin, X. Fan, C. Zhang, and T. Liu. Crack initiation, propagation,and failure characteristics of jointed rock or rock-like specimens: A review.
Advances n Civil Engineering , 2019:1–31, 02 2019. doi: 10.1155/2019/6975751.[2] A. Al-Mukhtar and B. Merkel. Simulation of the crack propagation in rocks usingfracture mechanics approach. Journal of Failure Analysis and Prevention , 15:90–100,01 2015. doi: 10.1007/s11668-014-9907-2.[3] M. Wang, F. Wang, Z. Zhu, Y. Dong, M. Mousavi Nezhad, and L. Zhou. Modellingof crack propagation in rocks under shpb impacts using a damage method.
Fatigue &Fracture of Engineering Materials & Structures , 42(8):1699–1710, 2019. doi: 10.1111/ffe.13012. URL https://onlinelibrary.wiley.com/doi/abs/10.1111/ffe.13012 .[4] B. K. Atkinson. Subcritical crack propagation in rocks: theory, experimental re-sults and applications.
Journal of Structural Geology , 4(1):41 – 56, 1982. ISSN0191-8141. doi: https://doi.org/10.1016/0191-8141(82)90005-0. URL .[5] Vastola, G. Asymmetric crack propagation near waterfall cliff and its influence on thewaterfall lip shape.
EPL , 96(4):49002, 2011. doi: 10.1209/0295-5075/96/49002. URL https://doi.org/10.1209/0295-5075/96/49002 .[6] M.-C. Eppes, A. Willis, J. Molaro, S. Abernathy, and B. Zhou. Cracks in Martianboulders exhibit preferred orientations that point to solar-induced thermal stress.
NatureCommunications , 6:6712, March 2015.[7] B. D. Collins and G. M. Stock. Rockfall triggering by cyclic thermal stressing of exfo-liation fractures.
Nature Geoscience , 9(5):395–400, May 2016.[8] B. D. Collins, G. M. Stock, M.-C. Eppes, S. W. Lewis, S. C. Corbett, and J. B. Smith.Thermal influences on spontaneous rock dome exfoliation.
Nature Communications , 9(1):762, February 2018.[9] H. Viles, B. Ehlmann, C. F. Wilson, T. Cebula, M. Page, and M. Bourke. Simulatingweathering of basalt on Mars and Earth by thermal cycling.
Geophysical ResearchLetters , 37(1):18201, September 2010.[10] O. Ruesch, E. Sefton-Nash, J. L. Vago, M. K¨uppers, J. H. Pasckert, K. Krohn, andK. Otto. In situ fragmentation of lunar blocks and implications for impacts and solar-induced thermal stresses.
Icarus , 336:113431, January 2020.[11] Y. Li, A. T. Basilevsky, M. Xie, and W.-H. Ip. Shape of boulders ejected from smalllunar impact craters.
Planetary and Space Science , 145:71–77, October 2017.[12] N. Attree, O. Groussin, L. Jorda, S. Rodionov, A. T. Auger, N. Thomas, Y. Brouet,O. Poch, E. K¨uhrt, M. Knapmeyer, F. Preusker, F. Scholten, J. Knollenberg, S. Hviid,and P. Hartogh. Thermal fracturing on comets. Applications to 67P/Churyumov-Gerasimenko.
Astronomy and Astrophysics , 610:A76, March 2018.2013] C. Matonti, N. Attree, O. Groussin, L. Jorda, S. Viseur, S. F. Hviid, S. Bouley,D. N´ebouy, A. T. Auger, P. L. Lamy, H. Sierks, G. Naletto, R. Rodrigo, D. Koschny,B. Davidsson, M. A. Barucci, J. L. Bertaux, I. Bertini, D. Bodewits, G. Cremonese,V. Da Deppo, S. Debei, M. De Cecco, J. Deller, S. Fornasier, M. Fulle, P. J. Guti´errez,C. G¨uttler, W. H. Ip, H. U. Keller, L. M. Lara, F. La Forgia, M. Lazzarin, A. Luc-chetti, J. J. Lopez Moreno, F. Marzari, M. Massironi, S. Mottola, N. Oklay, M. Pajola,L. Penasa, F. Preusker, H. Rickman, F. Scholten, X. Shi, I. Toth, C. Tubiana, andJ. B. Vincent. Bilobate comet morphology and internal structure controlled by sheardeformation.
Nature Geoscience , 12(3):157–162, February 2019.[14] M. El Maarry, N. Thomas, A. Gracia Bern´a, R. Marschall, A. Auger, O. Groussin,S. Mottola, M. Pajola, M. Massironi, S. Marchi, S. H¨ofner, F. Preusker, F. Scholten,L. Jorda, E. K¨uhrt, H. Keller, H. Sierks, M. A’Hearn, C. Barbieri, M. Barucci,J. Bertaux, I. Bertini, G. Cremonese, V. Da Deppo, B. Davidsson, S. Debei,M. De Cecco, J. Deller, C. G¨uttler, S. Fornasier, M. Fulle, P. Gutierrez, M. Hofmann,S. Hviid, W. Ip, J. Knollenberg, D. Koschny, G. Kovacs, J. Kramm, M. K¨uppers,P. Lamy, L. Lara, M. Lazzarin, M. J. Lopez, F. Marzari, H. Michalik, G. Naletto,N. Oklay, A. Pommerol, H. Rickman, R. Rodrigo, C. Tubiana, and J. Vincent. Frac-tures on comet 67p/churyumov-gerasimenko observed by rosetta/osiris.
GeophysicalResearch Letters , 42(13):5170–5178, 2015.[15] A. J. Dombard, O. S. Barnouin, L. M. Prockter, and P. C. Thomas. Boulders andponds on the asteroid 433 eros.
Icarus , 210(2):713 – 721, 2010. ISSN 0019-1035. doi:https://doi.org/10.1016/j.icarus.2010.07.006. URL .[16] D. S. Lauretta, D. N. DellaGiustina, C. A. Bennett, K. J. Becker, O. S. Barnouin,W. F. Bottke, C. Y. Drouet d’Aubigny, J. P. Dworkin, J. P. Emery, V. E. Hamilton,C. W. Hergenrother, M. R. M. Izawa, B. Rizk, K. J. Walsh, C. W. V. Wolner, O.-R.Team, H. C. Connolly, B. E. Clark, H. Campins, T. L. Becker, W. V. Boynton, C. Y. D.d’Aubigny, H. L. Enos, D. R. Golish, E. S. Howell, S. S. Balram-Knutson, M. C. Nolan,H. L. Roper, P. H. Smith, D. J. Scheeres, and H. H. Kaplan. The unexpected surfaceof asteroid (101955) Bennu.
Nature , 568(7):55–60, March 2019.[17] K. J. Walsh, E. R. Jawin, R. L. Ballouz, O. S. Barnouin, E. B. Bierhaus, H. C. Connolly,J. L. Molaro, T. J. McCoy, M. Delbo, C. M. Hartzell, M. Pajola, S. R. Schwartz,D. Trang, E. Asphaug, K. J. Becker, C. B. Beddingfield, C. A. Bennett, W. F. Bottke,K. N. Burke, B. C. Clark, M. G. Daly, D. N. DellaGiustina, J. P. Dworkin, C. M. Elder,D. R. Golish, A. R. Hildebrand, R. Malhotra, J. Marshall, P. Michel, M. C. Nolan,M. E. Perry, B. Rizk, A. Ryan, S. A. Sandford, D. J. Scheeres, H. C. M. Susorney,F. Thuillet, D. S. Lauretta, and O.-R. Team. Craters, boulders and regolith of (101955)Bennu indicative of an old and dynamic surface.
Nature Geoscience , 12(4):242–246,March 2019.[18] D. N. DellaGiustina, J. P. Emery, D. R. Golish, B. Rozitis, C. A. Bennett, K. N.Burke, R. L. Ballouz, K. J. Becker, P. R. Christensen, C. Y. Drouet d’Aubigny, V. E.21amilton, D. C. Reuter, B. Rizk, A. A. Simon, E. Asphaug, J. L. Bandfield, O. S.Barnouin, M. A. Barucci, E. B. Bierhaus, R. P. Binzel, W. F. Bottke, N. E. Bowles,H. Campins, B. C. Clark, B. E. Clark, H. C. Connolly, M. G. Daly, J. D. Leon, M. Delbo,J. D. P. Deshapriya, C. M. Elder, S. Fornasier, C. W. Hergenrother, E. S. Howell, E. R.Jawin, H. H. Kaplan, T. R. Kareta, L. Le Corre, J. Y. Li, J. Licandro, L. F. Lim,P. Michel, J. Molaro, M. C. Nolan, M. Pajola, M. Popescu, J. L. R. Garcia, A. Ryan,S. R. Schwartz, N. Shultz, M. A. Siegler, P. H. Smith, E. Tatsumi, C. A. Thomas, K. J.Walsh, C. W. V. Wolner, X. D. Zou, D. S. Lauretta, and O.-R. Team. Properties ofrubble-pile asteroid (101955) Bennu from OSIRIS-REx imaging and thermal analysis.
Nature Astronomy , 3:341–351, March 2019.[19] J. L. Molaro, K. J. Walsh, E. R. Jawin, R. L. Ballouz, C. A. Bennett, D. N. DellaG-iustina, D. R. Golish, C. Drouet d’Aubigny, B. Rizk, S. R. Schwartz, R. D. Hanna, S. J.Martel, M. Pajola, H. Campins, A. J. Ryan, W. F. Bottke, and D. S. Lauretta. In situevidence of thermally induced rock breakdown widespread on Bennu’s surface.
NatureCommunications , 11(1):2913–11, June 2020.[20] M. Delbo, G. Libourel, J. Wilkerson, N. Murdoch, P. Michel, K. T. Ramesh, C. Ganino,C. Verati, and S. Marchi. Thermal fatigue as the origin of regolith on small asteroids.
Nature , 508(7):233–236, April 2014.[21] K. Hall. Evidence for freeze-thaw events and their implications for rock weathering innorthern Canada.
Earth Surface Processes and Landforms , 29:43–57, January 2004.[22] D. T. Griggs. The Factor of Fatigue in Rock Exfoliation.
The Journal of Geology , 44(7):783–796, November 1936.[23] J. L. Molaro, S. Byrne, and L. J. L. Thermally induced stresses in boulders on airlessbody surfaces, and implications for rock breakdown.
Icarus , 294:247–261, September2017.[24] J. L. Molaro, S. Byrne, and S. A. Langer. Grain-scale thermoelastic stresses and spa-tiotemporal temperature gradients on airless bodies, implications for rock breakdown.
Journal of Geophysical Research: Planets , 120(2):255–277, February 2015.[25] B. Ravaji, V. Ali-Lagoa, M. Delbo, and J. W. Wilkerson. Unraveling the Mechanics ofThermal Stress Weathering: Rate-Effects, Size-Effects, and Scaling Laws.
Journal ofGeophysical Research: Planets , 124(1):3304–3328, December 2019.[26] C. El Mir, K. T. Ramesh, and M. Delbo. The efficiency of thermal fatigue in regolithgeneration on small airless bodies.
Icarus , 333:356–370, November 2019.[27] K. Hazeli, C. El Mir, S. Papanikolaou, M. Delbo, and K. T. Ramesh. The origins ofAsteroidal rock disaggregation: Interplay of thermal fatigue and microstructure.
Icarus ,304:172–182, April 2018.[28] B. K. Atkinson. Subcritical crack growth in geological materials.
Journal of GeophysicalResearch: Solid Earth , 89(B):4077–4114, June 1984.2229] A. J. Dombard, O. S. Barnouin, L. M. Prockter, and P. C. Thomas. Boulders and pondson the Asteroid 433 Eros.
Icarus , 210(2):713–721, December 2010.[30] B. Liang, J. Cuadra, K. Hazeli, and S. Soghrati. Stress field analysis in a stonymeteorite under thermal fatigue and mechanical loadings.
Icarus , 335:113381, 2020.ISSN 0019-1035. doi: https://doi.org/10.1016/j.icarus.2019.07.015. URL .[31] H. Yano, T. Kubota, H. Miyamoto, T. Okada, D. Scheeres, Y. Takagi, K. Yoshida,M. Abe, S. Abe, O. Barnouin-Jha, A. Fujiwara, S. Hasegawa, T. Hashimoto, M. Ishiguro,M. Kato, J. Kawaguchi, T. Mukai, J. Saito, S. Sasaki, and M. Yoshikawa. Touchdownof the hayabusa spacecraft at the muses sea on itokawa.
Science , 312(5778):1350–1353,2006.[32] J. Veverka, P. C. Thomas, M. Robinson, S. Murchie, C. Chapman, M. Bell, A. Harch,W. J. Merline, J. F. Bell, B. Bussey, B. Carcich, A. Cheng, B. Clark, D. Domingue,D. Dunham, R. Farquhar, M. J. Gaffey, E. Hawkins, N. Izenberg, J. Joseph, R. Kirk,H. Li, P. Lucey, M. Malin, L. McFadden, J. K. Miller, W. M. Owen, C. Peterson,L. Prockter, J. Warren, D. Wellnitz, B. G. Williams, and D. K. Yeomans. Imaging ofsmall-scale features on 433 eros from near: Evidence for a complex regolith.
Science ,292(5516):484–488, 2001.[33] N. Murdoch, P. S´anchez, S. R. Schwartz, and H. Miyamoto. Asteroid Surface Geo-physics. in Asteroids IV (P. Michel, et al. eds.) University of Arizona Press, Tucson. ,pages 767–792, 2015.[34] J. L. Molaro, C. W. Hergenrother, S. R. Chesley, K. J. Walsh, R. D. Hanna, C. W.Haberle, S. R. Schwartz, R. L. Ballouz, W. F. Bottke, H. J. Campins, and D. S. Lauretta.Thermal Fatigue as a Driving Mechanism for Activity on Asteroid Bennu.
Journal ofGeophysical Research: Planets , 125(8):e06325, August 2020.[35] D. S. Lauretta, C. W. Hergenrother, S. R. Chesley, J. M. Leonard, J. Y. Pelgrift, C. D.Adam, M. Al Asad, P. G. Antreasian, R. L. Ballouz, K. J. Becker, C. A. Bennett,B. J. Bos, W. F. Bottke, M. Brozovic, H. Campins, H. C. Connolly, M. G. Daly, A. B.Davis, J. de Le´on, D. N. DellaGiustina, C. Y. Drouet d’Aubigny, J. P. Dworkin, J. P.Emery, D. Farnocchia, D. P. Glavin, D. R. Golish, C. M. Hartzell, R. A. Jacobson,E. R. Jawin, P. Jenniskens, J. N. Kidd, E. J. Lessac-Chenen, J. Y. Li, G. Libourel,J. Licandro, A. J. Liounis, C. K. Maleszewski, C. Manzoni, B. May, L. K. McCarthy,J. W. McMahon, P. Michel, J. L. Molaro, M. C. Moreau, D. S. Nelson, W. M. Owen,B. Rizk, H. L. Roper, B. Rozitis, E. M. Sahr, D. J. Scheeres, J. A. Seabrook, S. H.Selznick, Y. Takahashi, F. Thuillet, P. Tricarico, D. Vokrouhlick´y, and C. W. V. Wolner.Episodes of particle ejection from the surface of the active asteroid (101955) Bennu.
Science , 366(6), December 2019.[36] V. Al´ı-Lagoa, M. Delbo’, and G. Libourel. Rapid temperature changes and the earlyactivity on comet 67p/churyumov-gerasimenko.
The Astrophysical Journal Letters , 810(2):L22, 2015. 2337] D. S. Lauretta, A. E. Bartels, M. A. Barucci, E. B. Bierhaus, R. P. Binzel, W. F. Bottke,H. Campins, S. R. Chesley, B. C. Clark, B. E. Clark, E. A. Cloutis, H. C. Connolly, M. K.Crombie, M. Delbo, J. P. Dworkin, J. P. Emery, D. P. Glavin, V. E. Hamilton, C. W.Hergenrother, C. L. Johnson, L. P. Keller, P. Michel, M. C. Nolan, S. A. Sandford,D. J. Scheeres, A. A. Simon, B. M. Sutter, D. Vokrouhlick´y, and K. J. Walsh. TheOSIRIS-REx target asteroid (101955) Bennu: Constraints on its physical, geological,and dynamical nature from astronomical observations.
Meteoritics & Planetary Science ,page 113, November 2014.[38] L. D. McFadden, M. C. Eppes, A. R. Gillespie, and B. Hallet. Physical weathering inarid landscapes due to diurnal variation in the direction of solar heating.
GeologicalSociety of America Bulletin , 117:161, n/a 2005.[39] M. Delbo, K. Walsh, J. Molaro, M. Al Asad, D. DellaGiustina, M. Pajola, C. Bennett,E. Jawin, R. Ballouz, S. Schwartz, B. Rizk, and D. Lauretta. Distribution of fractures onboulders on (101955) Bennu: OSIRIS-REx searching for evidence of thermal cracking. In
EPSC-DPS Joint Meeting 2019 , volume 2019, pages EPSC–DPS2019–176, September2019.[40] M. Delbo, M. Mueller, J. P. Emery, B. Rozitis, and M. T. Capria. Asteroid Thermo-physical Modeling. in Asteroids IV (P. Michel, et al. eds.) University of Arizona Press,Tucson. , pages 107–128, n/a 2015.[41] A. A. Griffith. VI. The phenomena of rupture and flow in solids.
Phil. Trans. R. Soc.Lond. A , 221(582-593):163–198, 1921.[42] J. R. Spencer, L. A. Lebofsky, and M. V. Sykes. Systematic biases in radiometricdiameter determinations.
Icarus , 78(2):337–354, April 1989.[43] B. Rozitis and S. F. Green. Directional characteristics of thermal-infrared beaming fromatmosphereless planetary surfaces - a new thermophysical model.
Monthly Notices ofthe Royal Astronomical Society , 415(3):2042–2062, August 2011.[44] F. E. DeMeo, R. P. Binzel, S. M. Slivan, and S. J. Bus. An extension of the Bus asteroidtaxonomy into the near-infrared.
Icarus , 202(1):160–180, July 2009.[45] V. E. Hamilton, A. A. Simon, P. R. Christensen, D. C. Reuter, B. E. Clark, M. A.Barucci, N. E. Bowles, W. V. Boynton, J. R. Brucato, E. A. Cloutis, H. C. Connolly,K. L. Donaldson Hanna, J. P. Emery, H. L. Enos, S. Fornasier, C. W. Haberle, R. D.Hanna, E. S. Howell, H. H. Kaplan, L. P. Keller, C. Lantz, J. Y. Li, L. F. Lim, T. J.McCoy, F. Merlin, M. C. Nolan, A. Praet, B. Rozitis, S. A. Sandford, D. L. Schrader,C. A. Thomas, X. D. Zou, D. S. Lauretta, and O.-R. Team. Evidence for widespreadhydrated minerals on asteroid (101955) Bennu.
Nature Astronomy , 3(4):332–340, March2019.[46] O. S. Barnouin, M. G. Daly, E. E. Palmer, R. W. Gaskell, J. R. Weirich, C. L. Johnson,M. M. Al Asad, J. H. Roberts, M. E. Perry, H. C. M. Susorney, R. T. Daly, E. B.Bierhaus, J. A. Seabrook, R. C. Espiritu, A. H. Nair, L. Nguyen, G. A. Neumann,24. M. Ernst, W. V. Boynton, M. C. Nolan, C. D. Adam, M. C. Moreau, B. Rizk, C. Y.Drouet d’Aubigny, E. R. Jawin, K. J. Walsh, P. Michel, S. R. Schwartz, R. L. Ballouz,E. M. Mazarico, D. J. Scheeres, J. W. McMahon, W. F. Bottke, S. Sugita, N. Hirata,S. i. Watanabe, K. N. Burke, D. N. DellaGiustina, C. A. Bennett, D. S. Lauretta, andO.-R. Team. Shape of (101955) Bennu indicative of a rubble pile with internal stiffness.
Nature Geoscience , 12(4):247–252, March 2019.[47] D. Uribe-Su´arez, P.-O. Bouchard, M. Delbo, and D. Pino-Mu˜noz. Numerical modelingof crack propagation with dynamic insertion of cohesive elements.
Engineering FractureMechanics , 227:106918, 2020.[48] F. Erdogan and G. C. Sih. On the crack extension in plates under plane loading andtransverse shear.
Journal of Basic Engineering , 85(4):519–525, 1963.[49] P. Destuynder, M. Djaoua, and S. Lescure. Quelques remarques sur la m´ecanique de larupture ´elastique.
J. M´eca. Th´eo. Appl. , 2(1):113–135, 1983.[50] X.-Z. Suo and A. Combescure. On the application of g ( θ ) method and its comparisonwith de lorenzi’s approach. Nuclear Engineering and Design , 135(2):207 – 224, 1992.ISSN 0029-5493. doi: https://doi.org/10.1016/0029-5493(92)90223-I. URL .[51] J. Brochard and X. Suo. Le taux de restitution d’´energie g en m´ecanique de la rupturenon-lin´eaire, formulation de la m´ethode gθ et description de la programmation danscastem 2000. Technical report, C.E.A., 1994. Rapport DMT/94-640.[52] T. Anderson. Fracture Mechanics: Fundamentals and Applications, Third Edition . Tay-lor & Francis, 2005.[53] M. Janssen, J. Zuidema, and R. Wanhill.
Fracture Mechanics . Spon Press, 2 nd edition,2004.[54] H. Digonnet, L. Silva, and T. Coupez. Cimlib: A Fully Parallel Application For Nu-merical Simulations Based On Components Assembly. In J. M. A. Cesar de Sa andA. D. Santos, editors, Materials Processing and Design; Modeling, Simulation and Ap-plications; NUMIFORM ’07 , volume 908 of
American Institute of Physics ConferenceSeries , pages 269–274, 2007.[55] A. J. Ryan, D. Pino Mu˜noz, M. Bernacki, and M. Delbo. Full-field modeling ofheat transfer in asteroid regolith: Radiative thermal conductivity of polydisperseparticulates.
Journal of Geophysical Research: Planets , 125(2), 2020. doi: 10.1029/2019JE006100. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JE006100 .[56] D. N. Arnold, F. Brezzi, and M. Fortin. A stable finite element for the stokes equations.
CALCOLO , 21(4):337–344, 1984.[57] F. Brezzi and M. Fortin.
Mixed and Hybrid Finite Element Methods . Springer series incomputational mathematics. Springer, 1991.2558] K.-J. Bathe. The inf–sup condition and its evaluation for mixed finite elementmethods.
Computers & Structures , 79(2):243 – 252, 2001. ISSN 0045-7949. doi:https://doi.org/10.1016/S0045-7949(00)00123-1. URL .[59] J. Hadamard. Sur les problemes aux derivees partielles et leur signification physique.
Princeton university bulletin , pages 49–52, 1902. URL https://ci.nii.ac.jp/naid/10030321135/en/ .[60] R. El Khaoulani.
Pr´ediction fiable de l’endommagement ductile par la m´ethode des´el´ements finis mixtes : endommagement non local et adaptation de maillage . PhDthesis, Ecole Nationale Sup´erieure des Mines de Paris, 2010.[61] T. S. Cao.
Modeling ductile damage for complex loading paths . PhD thesis, EcoleNationale Sup´erieure des Mines de Paris, 2013.[62] R.-L. Ballouz, K. Walsh, O. Barnouin, D. DellaGiustina, M. Al-Asad, E. Jawin, M. Daly,W. Bottke, P. Michel, C. Avdellidou, M. Delbo, R. Daly, E. Asphaug, C. Bennett,E. Bierhaus, H. C. Jr., D. Golish, J. Molaro, M. Nolan, M. Pajola, B. Rizk, S. Schwartz,D. Trang, C. Wolner, and D. Lauretta. Asteroid bennu’s boulders record near-earthhistory of impacts by millimeter- to centimeter-scale objects.
Nature , forthcoming.[63] R. S. Barsoum. On the use of isoparametric finite elements in linear fracture mechanics.
International Journal for Numerical Methods in Engineering , 10(1):25–37, 1976.[64] S. Suresh.
Fatigue of Materials . Cambridge University Press, 2 ndnd