Strength and energy consumption of inherently anisotropic rocks at failure

Using a discrete-element approach and a bonding interaction law, we model and test crushable inherently anisotropic structures reminiscent of the layering found in sedimentary and metamorphic rocks. By systematically modifying the level of inherent anisotropy, we characterize the evolution of the failure strength of circular rock samples discretized using a modified Voronoi tesselation under diametral point loading at different orientations relative to the sample's layers. We characterize the failure strength, which can dramatically increase as the loading becomes orthogonal to the rock layers. We also describe the evolution of the macroscopic failure modes as a function of the loading orientation and the energy consumption at fissuring. Our simulation strategy let us conclude that the length of bonds between Voronoi cells controls the energy being consumed in fissuring the rock sample, although failure modes and strength are considerably changing. We end up this work showing that the microstructure is largely afected by the level of inherent anisotropy and loading orientation.


Introduction
The failure strength of rocks is linked to the genesis and mineral composition and organization in space. In sedimentary or metamorphic rocks, the layering of minerals strongly affects joint family characteristics, bonds, and fissuring (i.e., an inherent geometrical anisotropy is imprinted in these materials by formation, tear and wear). Ultimately, the mechanical behavior of anisotropic rocks may strongly vary due to the layering properties and the loading orientation on rock samples [1][2][3][4][5]. Typically, the failure strength of rock cores presenting an anisotropic inner structure is tested under diametral point load, also known as the brazilian test in which the loading orientation is varied from totally aligned to the layering up to orthogonal. In general, the failure strength turns out to be maximal for the orthogonal loading and minimal for diagonal-like load orientation relative to the layering.
More recently, numerical methods, mostly based on the discrete-element approach, have been employed to model and analyze the behavior of anisotropic rock given its capabilities of reproducing failure and crushing (see, for instance, Refs. [5,6]). These methods have added valuable information concerning the modeling and failure mechanism of anisotropic brittle structures. Nonetheless, the microstructure characteristics, failure modes evolution, and energy consumption partitions by rupture modes, still lack an extended analysis. This is not a trivial modeling and analysis given that anisotropic rock structures need a numerical approach that allows samples to break in multiple irregular pieces of diverse shape and size, while mea- * e-mail: david.cantor@polymtl.ca suring energy consumption as fissures appears within the bulk. As we will show in this paper, this is possible by using a discrete-element approach and a modified Voronoi decomposition of samples, allowing us to create cells reminiscent of the mineral organization, and letting us finely control the level of inherent anisotropy. The outcomes of this research axis are large and may cover large scales for tunnel and other underground earthworks, up to the failure of granulates and conglomerates for powder technology. This paper is organized as follows. In Sec. 2, we describe the numerical method employed to build and test anisotropic rock samples. In Sec. 3, we present the characterization of the failure strength as a function of the level of inherent anisotropy and loading orientation. We then describe the failure modes evolution as a function of the loading orientation in Sec. 4. We also quantify the energy consumption at fissuring from the beginning of the loading up to the breakage in Sec. 5. Finally, in Sec. 6, we show how the loading affects the microstructure by comparing the intact connectivity of cells and that on the onset of failure. We finish this work with some conclusions and perspectives.

Numerical model
We use the discrete-element method known as contact dynamics (CD) [7,8] to simulate breakable brittle rock cores. The CD method is a discrete approach in which rigid bodies interact with solid contact, friction, and cohesion. In particular, the CD method is an implicit method in which the equations of motion are applied to a collection of bod- ies in contact for which velocities and interaction forces are computed on a time-stepping scheme.
In order to build anisotropic structures, we use a Voronoi tessellation on circular samples, so the resulting cells' geometry translates the inherent anisotropy level of a rock. Given that a random Voronoi tessellation produces cells of varied shapes and sizes, we modify the tessellation procedure to finely control the level of anisotropy. First, we employ a centroidal Voronoi tessellation in which iteratively we reuse the cells' centroid as seeding for next meshings. By doing so, we can obtain cells of similar geometries. Then, we elongate all the cells in a given direction so we can compute their average shape ratio η = h/L, with h and L the average small and larger dimensions. We produce a set of samples in which η varied from one to six in steps of 1 (see Fig. 1).
To provide mechanical strength to the assembly of cells, we used the bonded-cell method approach previously used in the simulation and analysis of crushable grains [9][10][11]. By adding a bonding law at the interfaces between cells (note that these interfaces are exclusively edge-edge contacts), the cells' interfaces can resist tensile and shearing stresses. This needs the definition of two values of cohesion C n and C t , that, for the sake of simplicity, we fix to the same value of 10 kPa. Additionally, we preset a debonding distance δ c necessary to break a cell-cell bond. In other words, a fissure between two cells can only be broken if once the stress threshold C n or C t is reached, they present a gap or relative sliding equals to δ c . The introduction of this debonding distance is not arbitrary. It allows us to represent the progressive debonding of two cells, and to characterize the type of material we are analyzing. Following fracture mechanics theory, we then choose a typical value of surface energy for silicate materials γ = 50J/m 2 [12], and then the debonding distance is simply defined as δ c = 2γ/C n .
Note that this construction strategy can generate configurations with varying number of cells depending on the anisotropy level. Several studies have used the cells' strategy to reproduce fragmentation and have spotlighted the role of the number of cells on the mechanical response. As a general observation, adding cells reduces the final strength of the samples, which is in agreement with the fact that more interactions between cells adds more failure potential paths. Nonetheless, to make the different samples comparable, we set the total length of cohesive interactions to be constant.
To test these samples, we used the diametral point load in which two rigid platens gradually increasing the applied force. We also varied the orientation of the load θ respect to the layering of the cells, as shown in the inset of Fig.  2 Finally, we perform our simulations by including the bonded-cell method in the simulation platform LMGC90, developed at the Université de Montpellier [13,14].

Failure strength
The failure strength of the samples is found as the loading reaches a critical force value F c , for which the structure can no longer bear load and collapses. By using the tensile failure criterion for circular geometries, we can characterize the failure strength σ c as with d the diameter of the sample. Figure 2 gathers the values of crushing strength as a function of the level of inherent anisotropy and loading orientation normalized by the cohesive strength C n . We can remark that for a non anisotropic structure η = 1, the strength remains near C n . Once some degree of anisotropy is introduced in the cells' configuration, the strength substantially varies. First, for loading orientations in the range [0 • , 70 • ] the strength is lower than C n and very similar between values of η. For orientations beyond 75 • , the failure strength rapidly increases and reaches higher values for larger η. This behavior is consistent with laboratory testing as well as the minimum failure strength for an orientation around 30 • after simple stress considerations [15]. In the following, we explore the mechanisms behind this strong difference in failure strength among samples and loading orientations.

Failure modes
The first qualitative characteristic we can highlight is the evolution of failure modes as the loading orientation varies. Figure 3 presents screenshots of the samples with η = 2 for increasing θ. For a loading orientation matching the cell's elongation (i.e., θ = 0 • ), the fissures are mostly vertical and percolating between the walls. For such a case, we can deduce that tensile stresses in the orthogonal direction to the loading were responsible of this failure. For loading orientations θ = 30 • and θ = 60 • , the failure paths are different and diagonal respect to the application point or force. The relative displacement between cells suggest that shearing is a major mechanism for debonding of cells. Finally, for the case θ = 90 • , we can observe that the fissuring paths are more complex, presenting larger amounts of damage near the contact with the loading plates, nonetheless a vertical series of cracks may suggest that tensile stresses were also important for this failure. This observations are in agreement with the drop of failure strength for anisotropic structures presenting diagonal fissuring trajectories. Indeed, mobilization of cells end up being traveling shorter distances for such cases. Nonetheless, the rapid increase of strength for more orthogonal-like loading conditions is impossible to grasp from the visualizations.

Energy consumption
We can analyze the consumption of energy by partitions being consumed in tensile failures and shear failures. The total energy E c being consumed at fissuring can be split as with c * n and c * t , the number of bonds broken in traction or shearing, and l c * n and l c * t their corresponding lengths. Figure 4 presents the evolution of E n and E t during the loading of samples with η = 1, 3, and 6, up to the  moment of failure (t c ) and for loading orientation θ = 0. The data is normalized by the nominal energy consumption for a vertical straight failure path (i.e., 2γd). In effect, the ratio E c /2γd translates the tortuosity of the failure path in the cell configuration. As a matter of fact, fixing the total bonding length between cells to be constant produces failures consuming similar energy amounts independently of η. Note that this is the failure mechanism at the cells' scale and is not necessarily related to the 'apparent' macroscopic failure mode of the assembly. The trend of these curves also proves that no burst of energy occurs at failure. Conversely, the cohesive bonds are progressively broken (i.e., fissuring) up to a critical energy for which the particle collapses. This fact is in agreement with the failure mechanics theory and energy criticality at failure.
In Fig. 5 we observe the total energy consumed by mode (traction or shearing) as a function of the loading orientation and η. We see that, in general, more fissuring occurs in shearing mode than in traction. A slight increase of E n and E t also seems to occur for increasing loading orientation θ. Nonetheless, the variation is relatively small compared with the order of magnitude of the energy, and the fact that the ratio between energies remains roughly the same for all the values of η and θ.

Cells' connectivity
Using the coordination number, defined as Z = 2N c /N cl , with N c the number of bonds between cells, and N cl the number of cells, we compute the average number of neighboring interactions per cell. This parameter makes part of the signature of mineral organization within geomaterials.
Note that a subtle difference in coordination number can be found between the intact coordination number (i.e., at the beginning of the loading) and the coordination number instants prior failure due to fissuring. As we rather characterize the onset of failure as the state bearing σ c , let us consider the cohesive bonds on the onset of failure N * c as the effective number of interactions, so Z * = 2N * c /N cl . Figure 6 presents the averaged coordination number on the onset of failure as a function of η. We observe that the connectivity decreases as θ and η grow. This behavior is due to the fact that orthogonal loading respect to the cell's main orientation can reach larger forces through column-like structure. In the meantime, larger values of η imply fewer contacts within the assembly, so the lost of a given number of bonds imply a stronger drop of Z * . Note that a centroidal Voronoi tessellation tends to produce hexagonal-like configurations of cells, then a value of Z 6 is expected for intact samples. Finally, the important variation of Z * with η and θ suggests that the number of bonds participating in the load transmission can largely change at the onset of failure, implying strong variation in the microstructural characteristics producing σ c at the macroscopic scale.

Conclusions
In the frame of the contact dynamics method and the bonded-cell approach, we built and tested crushable samples recalling inherent anisotropic geometries often found in sedimentary or metamorphic rocks. This approach let us define a level of anisometry η of the internal structure of circular samples divided in cells via a modified Voronoi tessellation. By applying a diametral loading, we were able to characterize the failure strength as a function of η and the loading orientation θ. Our results, in agreement with experimental observations and fracture mechanics theory, showed that the strength decreases for orienta-tions θ < 70 • once an inherent anisotropic structure is included in the samples. But as θ gets values beyond 75 • , the failure strength rapidly increases with η.
Despite these variations of failure strength, our simulations showed that fixing the total length of cohesive bonds between cells determines the amount of energy being consumed at fissuring, which, ultimately, has not relation with the resistance. Then, another microstructural elements should be at the origin of the failure strength. We explore this microstructure using the coordination number Z. We observed that this parameter it strongly varies with η and θ. But, counterintuitively, it seems that stronger inherent anisotropic structures also are those who present lower connectivity of their components. Future research should focus on the force transmission mechanisms and should finely quantify the inherent and induced anisotropy characteristics to better asses the sources of the macroscopic observations from the cells and bonds scale.