David Cantor
email: david.cantor@tpolymtl.ca

Carlos Ovalle
email: carlos.ovalle@polymtl.ca

Emilien Azéma
email: emilien.azema@umontpellier.fr

Microstructural origins of crushing strength for inherently anisotropic brittle materials

Keywords: fabric, anisotropy, failure strength, Weibull statistics, fragmentation energy, discrete element method

We study the crushing strength of brittle materials whose internal structure (e.g., mineral particles or grains) presents a layered arrangement reminiscent of sedimentary and metamorphic rocks. Taking a discrete-element approach,

we probe the failure strength of circular-shaped samples intended to reproduce specific mineral configurations. To do so, assemblies of cells, products of a modified Voronoi tessellation, are joined in mechanically-stable layerings using a bonding law. The cells' shape distribution allows us to set a level of inherent anisotropy to the material. Using a diametral point loading, and systematically changing the loading orientation with respect to the cells' configuration, we characterize the failure strength of increasingly anisotropic structures. This approach lets us reproduce experimental observations regarding the shape of the failure strength curve, the Weibull modulus, failure

Introduction

The mechanical behavior of many solids can be tracked down to the level of molecules, defects, and dislocations that may cause stress concentrations and yielding of the material. However, the microstructural level can be equally or more important for characterizing their failure strength. Imagine, 5 for instance, arrangements produced by the genesis, layering, and stratification of a solid such as sedimentation, rock metamorphism, or even crystal growth. In these cases, the mechanical properties are not determined at the molecular level but instead at the microstructural level [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF].

We focus on the failure strength of brittle materials, which may include 10 rocks, soil grains, ceramics, and even ice. In the case of rock masses and rock aggregates, for example, it is well known that the mechanical behav-ior depends on mineralogy and grain characteristics (e.g., size and shape distribution of minerals), matrix level of cementation, joint characteristics, and fissuring [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF]. Many of these geometrical attributes are also called fabric or microstructure for geological materials. If any of those fabric properties present a preferred orientation or organization in space, then the material can be considered inherently anisotropic. Many studies have focused on the quantification of the level of inherent fabric anisotropy in rocks or on the impact of the loading orientation on the ultimate strength of such layered materials [START_REF] Hoek | Fracture of Anisotropic Rock[END_REF][START_REF] Oda | Fabric tensor for discontinuous geological materials[END_REF][START_REF] Amadei | Importance of anisotropy when estimating and measuring in situ stresses in rock[END_REF][START_REF] Chen | Determination of Deformability and Tensile Strength of Anisotropic Rock Using Brazilian Tests[END_REF][START_REF] Zhang | Estimating the intensity of rock discontinuities[END_REF][START_REF] Karakul | Empirical models and numerical analysis for assessing strength anisotropy based on block punch index and uniaxial compression tests[END_REF][START_REF] Khanlari | An Experimental Investigation of the Brazilian Tensile Strength and Failure Patterns of Laminated Sandstones[END_REF][START_REF] Roy | Effect of Heat Treatment and Layer Orientation on the Tensile Strength of a Crystalline Rock Under Brazilian Test Condition[END_REF][START_REF] Pouragha | Micromechanical correlation between elasticity and strength characteristics of anisotropic rocks[END_REF].

In experiments, cylindrical cores are often used to characterize the failure strength of inherently anisotropic rocks under diametrical point loading (commonly called Brazilian test). In those tests, the orientation θ of the applied force is gradually varied with respect to the orientation of the internal layering (see Fig. 1). For such a circular geometry, axial symmetry is found along the layering orientation; then, the failure strength can be fully characterized by varying θ in the range [0 • , 90 • ]. Note as well that the rock cores and test configurations are chosen, so the inherently anisotropic configuration is homogeneous as possible along the length of the cylinders (i.e., the mineral configuration can be considered transversely isotropic).

For rocks not presenting an inherent anisotropy, the failure strength is independent of the loading orientation θ. However, for highly layered rocks such as slate, schist or shale, the failure strength largely varies with θ in a 'U ' shape with minimal strength for an orientation around θ ≃ 25 • , and increasing strength as θ → 0 • or 90 • [START_REF] Amadei | Importance of anisotropy when estimating and measuring in situ stresses in rock[END_REF][START_REF] Karakul | Empirical models and numerical analysis for assessing strength anisotropy based on block punch index and uniaxial compression tests[END_REF][START_REF] Pouragha | Micromechanical correlation between elasticity and strength characteristics of anisotropic rocks[END_REF][START_REF] Saroglou | A modified Hoek-Brown failure criterion for anisotropic intact rock[END_REF][START_REF] Saeidi | A modified failure criterion for transversely isotropic rocks[END_REF][START_REF] Garagon | Predicting the strength anisotropy in uniaxial compression of some laminated sandstones using multivariate regression analysis[END_REF][START_REF] Xu | Discrete element modeling of transversely isotropic rocks with non-continuous planar fabrics under Brazilian test[END_REF]. As shown early by Hoek in 1964, the 'U ' shape is consistent with Griffith's theory of brittle fracture of materials, in which the maximum stress at the tip of the crack triggers the propagation of a fissure once a critical amount of energy is added to the system. The rock microstructure is then capable of signing the failure modes and patterns. This in turn significantly affects the stress-strain relations at a macroscopic level (i.e., stiffness, hardening/softening, and strength) [START_REF] Hurley | In situ grain fracture mechanics during uniaxial compaction of granular solids[END_REF][START_REF] Marinelli | Anisotropic breakage mechanics: From 38 stored energy to yielding in transversely isotropic granular rocks[END_REF].

It is important to clearly distinguish between inherent and induced anisotropies.

While the first is defined here as a property of the fabric, the latter refers to anisotropies arising from that primary structure, such as the joint distribution in space. Although the anisotropies mentioned up this point are only related to geometrical properties, they can also arise from the loading configuration (i.e., stress-induced anisotropies). Indeed, a detailed description of both inherent and induced anisotropies are key elements for understanding the behavior of brittle materials, as we will show in this paper.

Improving our understanding of anisotropic geological materials will help us better address problems involving rock and grain fragmentation, such as railway ballast design [START_REF] Indraratna | Advanced rail geotechnology-ballasted track[END_REF][START_REF] Lim | The application of Weibull statistics to the strength of railway ballast[END_REF], rockfill dam design [START_REF]Mechanical properties of rockfill dams[END_REF][START_REF] Ovalle | The effect of size on the strength of coarse rock aggregates and large rockfill samples through experimental data[END_REF][START_REF] Ovalle | Data Compilation from Large Drained Compression Triaxial Tests on Coarse Crushable Rockfill Materials[END_REF], rock tunneling processes [START_REF] Favier | Statistics of the deformation of the front of a tunnel-crack propagating in some inhomogeneous medium[END_REF][START_REF] Pindra | Geometrical disorder of the fronts of a tunnel-crack propagating in shear in some heterogeneous medium[END_REF], mining waste dumps construction [START_REF] Bard | Waste Rock Behavior at High Pressures: Dimensioning High Waste Rock Dumps[END_REF], surface subsidence [START_REF] Brzesowsky | Compaction creep of sands due to time-dependent grain failure: Effects of chemical environment, applied stress, and grain size[END_REF], slip stability analysis of fault gouges [START_REF] Main | Stress corrosion constitutive laws as a possible mechanism of intermediateterm and shortterm seismic quiescence[END_REF], filling rock mass discontinuities [START_REF] Sammis | The kinematics of gouge deformation[END_REF], the geological formation of glacial till [START_REF] Turcotte | Fractals and Fragmentation[END_REF], confined comminution [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF][START_REF] Ovalle | Population balance in confined comminution using a physically based probabilistic approach for polydisperse granular materials[END_REF], weathering and environmental degradation effects [START_REF] Zhang | Breakage mechanics for granular materials in surface-reactive environments[END_REF], the failure localization in progressively deteriorated disordered structures [START_REF] Berthier | Forecasting failure locations in 2-dimensional disordered lattices[END_REF], etc.

In this paper, we use bi-dimensional discrete-element modeling to study the failure strength of circular samples that have an inherent anisotropic configuration under varied loading orientation. In Sec. 2, we introduce our numerical strategy based on the contact dynamics and the bonded-cell methods, and the sample construction and testing procedures. In Sec. 3, we characterize the failure strength of two-inherent anisotropic structures showing a good agreement with experimental observations. We then analyze the failure strength variability in terms of Weibull's statistics and failure mode evolution. Section 4 focuses on a fine description of the microstructure in terms of fabric connectivity, force transmission, and inherent and induced anisotropies. In Sec. 5, we develop a theoretical analysis that allows us to discover the microstructural origins of the crushing strength in terms of the level of inherent anisotropy and loading orientation. This analytical approach based on the granular stress tensor and its harmonic decomposition linking microstructure and the macromechanical response. Finally, we conclude with a summary and perspectives.

Numerical modeling

Inherently anisotropic materials are challenging to characterize given the complex and multiscale properties that minerals, grains, bonds, and fissures can present in space. Numerical approaches have proven successful at an-alyzing these materials because they are able to reproduce complex failure mechanisms under controlled geometries. Some of these approaches use, for example, finite-elements [START_REF] Sulem | Finite element analysis of the indentation test on rocks with microstructure[END_REF][START_REF] Amir Reza Beyabanaki | Nodal-based three-dimensional discontinuous deformation analysis (3-D DDA)[END_REF], discrete-element methods with bonded bodies [START_REF] Potyondy | A bonded-particle model for rock[END_REF][START_REF] Cho | A clumped particle model for rock[END_REF][START_REF] Lan | Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading[END_REF][START_REF] Scholtès | A DEM model for soft and hard rocks: Role of grain interlocking on strength[END_REF][START_REF] Kazerani | Effect of micromechanical parameters of microstructure on compressive and tensile failure process of rock[END_REF][START_REF] Gao | The application of a modified Voronoi logic to brittle fracture modelling at the laboratory and field scale[END_REF], splitting or replacing mechanisms [START_REF] Cantor | Split-Cell Method for grain fragmentation[END_REF][START_REF] Ciantia | An approach to enhance efficiency of DEM modelling of soils with crushable grains[END_REF][START_REF] Gladkyy | DEM simulation of polyhedral particle cracking using a combined MohrCoulombWeibull failure criterion[END_REF][START_REF] Iliev | Evolution of fragment size distributions from the crushing of granular materials[END_REF], or coupled discretefinite element strategies [START_REF] Mahabadi | An example of realistic modelling of rock dynamics problems: FEM/DEM simulation of dynamic brazilian test on Barre Granite[END_REF][START_REF] Bagherzadeh Kh | Numerical simulation of particle breakage of angular particles using combined DEM and FEM[END_REF][START_REF] Guo | A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media[END_REF][START_REF] Ma | Combined FEM/DEM modeling of triaxial compression tests for rockfills with polyhedral particles[END_REF][START_REF] Tan | Brazilian Tests on Transversely Isotropic Rocks: Laboratory Testing and Numerical Simulations[END_REF].

Among these approaches, the discrete-element method (DEM) has become increasingly popular for dealing with fragmentation due to its versatility in reproducing grain fissuring, crushing, and many experimental observations [51? -53]. However, some modeling strategies employ circular particles to represent grains and blocks [START_REF] Cheng | Crushing and plastic deformation of soils simulated using DEM[END_REF][START_REF] Lin | A combined method to model grain crushing with dem[END_REF][START_REF] Tavares | An effective sphere-based model for breakage simulation in dem[END_REF], which does not capture the complex variability of fragments' shapes and sizes. Other studies use a 'replacement' method in which bigger grains are substituted by a set of smaller bodies once a criterion is reached [START_REF] Ciantia | An approach to enhance efficiency of DEM modelling of soils with crushable grains[END_REF][START_REF] De Bono | Discrete element modelling of onedimensional compression of cemented sand[END_REF][START_REF] Zhou | Discrete element modeling of particle breakage considering different fragment replacement modes[END_REF], at the expense of missing mass conservation or creating local over-stresses at the replacement instant. Finally, energy consumption is not traceable when using circular bodies, or ad-hoc parameters are necessary to estimate the fragmentation energy. While these approaches have enabled the exploration of certain mechanisms of rock and grain failure, no clear mapping between the variability of strength, failure modes, and the microstructure has been found for inherently anisotropic materials. A correct simulation of these materials requires a model in which bodies can break into irregular and size disperse fragments while simultaneously controlling the inherent anisotropy level. As we show in the next section, these conditions can be met in 2D simulations using irregular convex polygons.

Construction of inherently anisotropic samples

We build circular samples composed of smaller bodies called cells using a Voronoi decomposition of a unitary circle. This procedure generates an assembly of N cl adjacent cells that we 'glue' using a cohesive bonding law. This approach, known as the bonded-cell method (BCM), has been used in numerous studies of the mechanical behavior of crushable granular materials, both in 2D [START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF] and 3D [START_REF] Cantor | Three-dimensional bondedcell model for grain fragmentation[END_REF][START_REF] Orozco | Discrete-element model for dynamic fracture of a single particle[END_REF][START_REF] Huillca | Modeling size effect on rock aggregates strength using a DEM bonded-cell model[END_REF].

A random Voronoi tessellation normally creates a disordered distribution of cell shapes and sizes. In order to control the cell's geometry (and, in effect, the inherent anisotropy), we alter the initial tessellation in two steps. First, we iteratively rebuild the Voronoi tessellation using the centroids of previous tessellation seedings to produce similar cells. This approach is also called centroid tessellation [START_REF] Du | Convergence properties of the lloyd algorithm for computing the centrodial voronoi tessellations[END_REF]. Then, the cells are elongated along a given direction and an anisometry level is estimated using the average aspect ratio of the cells η = h/L, with h and L being the average short and long dimensions, respectively. This anisometry represents the inherent anisotropy configuration of the minerals in our model. We produced a set of samples with η = [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF][START_REF] Chen | Determination of Deformability and Tensile Strength of Anisotropic Rock Using Brazilian Tests[END_REF] in steps of 1 (see Fig. 3). Additionally, perturbations to the initial setting of the tessellation enabled us to have slightly different cell arrangements. For statistical representativeness, we built five different configurations for each value of η.

In order to give mechanical strength to the assembly of cells, we define a normal and tangential cohesion at the bonds (i.e., cell-cell interactions), C n and C t , respectively. C n prevents the interactions between cells from separating due to tensile stresses, while C t provides resistance against sliding.

We also preset a debonding distance δ c needed to effectively break a cohesive bond. By choosing a typical value of surface energy density for silicate minerals γ = 50 J/m 2 [START_REF] Jones | Engineering materials[END_REF], we can then determine the separation threshold as δ c = 2γ/C n , following fracture mechanics theory. Note that our model allows us to independently define the tensile and shear bonding strength, but for simplicity we set C n = C t . A detailed analysis of the combined effect of varying C n and C t can be found in Ref. [START_REF] Cantor | Three-dimensional bondedcell model for grain fragmentation[END_REF].

The critical rupture energy that a bond needs to break is thus E c = 2γl c , with l c the length of the interaction. Once E c is reached, the cohesive bond is removed, simulating a fissuring event. These fissures are considered dry frictional surfaces, with µ being the coefficient of friction that we set to 0.4 (see Fig. 2 for a schematic representation of the bonding law).

In addition, numerical studies have explored the effect of the number of cells on the failure strength of brittle materials, showing that an increased number of cells lowers the failure strength [START_REF] Nguyen | Bonded-cell model for particle fracture[END_REF][START_REF] Huillca | Modeling size effect on rock aggregates strength using a DEM bonded-cell model[END_REF]. However, it was recently shown that the scalability of failure strength is not simply linked to the number of cells, but more importantly to the length of bonding interactions [START_REF] Orozco | Discrete-element model for dynamic fracture of a single particle[END_REF][START_REF] Cantor | Strength and energy consumption of inherently anisotropic rocks at failure[END_REF]. Thus, to make the tests comparable, samples must present the same potential surface energy among the different values of η despite presenting a different number of cells. So, in our tests, the samples have the same total length of bonds.

Contact dynamics

The contact dynamics (CD) method is a discrete-element approach in which rigid bodies interact via non-smooth laws [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Dubois | The Contact Dynamics method: A nonsmooth story[END_REF], i.e., impacts are transmitted on an implicit time-stepping scheme. At the end of each time- step, particles' velocities and contact forces can therefore be simultaneously computed without requiring force-overlapping laws (i.e., no regularization of the contact law is needed). This allows the CD method to be unconditionally stable and capable of employing larger time-steps than in alternative smooth approaches. For details on implementation of the contact dynamics method, see Refs. [START_REF] Renouf | A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media[END_REF][START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF].

In two-dimensional simulations, three main interactions can occur between convex bodies: vertex-vertex, vertex-edge, edge-edge (see Fig. 4).

Vertex-vertex interactions are rare and unstable, so they are discarded from the computation and analysis. For edge-edge interactions, it is necessary to consider two contact points to correctly resolve the contact mechanics; however, only the resultant force is important, rendering the loci of the two contact points irrelevant. For the interaction detection and classification, we use the shadow-overlap method [START_REF] Saussine | Modelling ballast behaviour under dynamic loading. Part 1: A 2D polygonal discrete element method approach[END_REF], which creates a separating plane between two touching bodies via an iterative procedure. Updated body positions, velocities, and interaction forces are governed by the equations of motion and the cohesive bonding law we previously defined.

Finally, it is worth mentioning that our simulations were performed using the CD method on the free and open-source platform LMGC90 [START_REF] Dubois | 10e colloque national en calcul des structures[END_REF][START_REF] Dubois | LMGC90 wiki page[END_REF].

Test procedures

Samples are initially set matching the loading orientation with the cells' preferred direction, so θ = 0 • . Then, we apply a gradually increasing vertical force F using rigid platens up to the failure. To avoid dynamic perturbations during loading, we make sure that, over a time-step, a load increment is very small before C n d, with d being the diameter of the samples. We systemati- cally vary θ in the range [0 • , 90 • ] in steps of 5 • , i.e., 18 different orientations (see Fig. 1). Finally, we test five different configurations of η for each one of the angles for statistical representativeness. Videos of the tests can be found at the following link https://youtu.be/N66c_4crwlM. We present the averaged results of the total 540 simulations we performed.

Macroscopic observations

Failure strength

Our samples are able to reproduce a brittle material behavior by supporting a load that gradually increases up to a critical value that triggers the collapse of the assembly (see Fig. 5(a)). The critical force at failure F c allows us to characterize strength using the vertical stress at failure, defined as point load on cylindrical samples produces more subtle variations of failure strength [START_REF] Khanlari | An Experimental Investigation of the Brazilian Tensile Strength and Failure Patterns of Laminated Sandstones[END_REF] The increasing disagreement in θ c with respect to the theoretical estimation also suggests that the inherent anisotropy deeply modifies the stress configuration within the samples.

σ yy = F c d . (1)
For the case without inherent anisotropy, we can follow classical rock me-220 chanics testing and use the maximum tensile stress criterion to characterize the strength of our samples [START_REF] Hiramatsu | Determination of the tensile strength of rock by a compression test of an irregular test piece[END_REF][START_REF] Jaeger | Failure of rocks under tensile conditions[END_REF]. So the expression σ c = 2F c /πd allows one to deduce the maximum tensile stress at the center of the circular sample.

Using that equation, we find that, on average, σ c /C n ≃ 1. This result allows us to make two important observations: 1) our model correctly scales the indeed triggering the failure. However, experimental testing characterizing inherently anisotropic materials should be aware that the tensile stresses are not necessarily at the origin of failure.

Statistical variability of strength

To better understand the variability of failure strength, we analyze our results in terms of 'survival' using the Weibull probability distribution. This approach assumes that a sample's probability of not presenting failure P s (i.e., the survival) depends on the applied stress σ as

P s = exp { - ( σ σ 0 ) m } , (2)
with σ 0 being a reference stress for which P s = 1/e ≃ 37%. The exponent m is known as the Weibull modulus and is associated with the sharpness of the probability distribution. As m increases, so does the slope of the distribution, meaning that the failure strength is focused on a given value.

Conversely, as m decreases, the stress range within which the particle may break broadens. If we combine all values of failure strength -not distinguishing between η and θ -we find that m ≃ 1.5, which falls within a typical range of values for rocks or grains that are not necessarily anisotropic. This observation suggests that laboratory tests in which no special attention is given to the degree of inherent anisotropy or loading orientation can gather a wide variety of material characteristics. Such a simplified approach could thus be misleading and limit the predictability of the material failure strength.

Macroscopic failure modes and energy consumption

Microstructural analysis

As previously mentioned, rock microstructure is often related to mineral or grain size and shape distribution, joint spacing/density, fissuring, and etc. All of these geometrical characteristics clearly affect the failure strength. However, the microstructure cannot only be reduced to its geometrical aspects. Accounting for the connectivity between cells and the force transmission are key elements behind the macroscopic mechanical behavior.

In order to do this, we first need to define a framework of analysis. We have two possibilities when dealing with adjacent cells. First, the interaction frame in which the bond forces are defined as f = f n n + f t t, with n being the normal unit vector perpendicular to the contact line, and t being the tangential unit vector.

We can also define the inter-center vectors between cells, also called branch vectors, as ℓ = ℓ n n + ℓ t t, with ℓ n and ℓ t being the normal and tangential components. These branches let us define a second frame in which the unit vector n ′ is defined along ℓ, and t ′ is the tangential unit vector (see Fig. 8) [START_REF] Azéma | Stress-strain behavior and geometrical properties of packings of elongated particles[END_REF]. In this frame, the forces between the cells are written

as f = f ′ n n ′ + f ′ t t ′
, with f ′ n and f ′ t -the radial and ortho-radial forces, respectively -acting between the centers of the cells. Finally, the branch in this frame is simply written as ℓ = ℓn ′ , with ℓ being the length of the branch vector. For convenience, we used the branch frame for the following microstructural analysis.

As an illustration, Fig. 9 (top) presents the branch network and the force network (bottom) with lines whose thickness is proportional to the intensity of the force at the interactions. A visual inspection shows how the geometry of the cells dramatically modifies both networks. As η increases, the branch network becomes more irregular and the force chains more diffuse within the volume.

η = 1 η = 3 η = 6
Bonds Network Forces Network Traction forces are displayed in blue and compression forces in red. These screenshots are taken for cases in which the loading orientation is θ = 0.

Geometrical description 4.1.1. Connectivity

We can characterize the connectivity between cells by using the coordination number Z. This parameter shows the average number of neighboring interactions per cell as Z = 2N c /N cl , with N c being the number of bonds between cells and N cl the number of cells. There is, however, a subtle difference in coordination number between the intact state (i.e., at the beginning of the loading) and the coordination number we are able to compute instants before failure. Both values differ, as fissuring removes cohesive bonds previously capable of bearing force. We characterize the onset of failure as the state bearing σ yy , so let us consider the cohesive bonds at the onset of failure N * c as the effective number of interactions, so Z * = 2N * c /N cl . Figure 10 presents the averaged coordination number at the onset of failure as a function of η and θ. We observe that the connectivity decreases as θ and η increase. In other words, Z * varies conversely to σ yy showing that inherently anisotropic materials can bear larger stresses despite the fact that cells are less connected. The evolution of Z * provides a counterintuitive picture of the effect of θ and η on the microstructural properties at failure. However, as previously shown in Fig. 8, the branch and force networks also carry a strong anisotropic character that calls for a higher-order analysis accounting for their distribution in space.

Branch orientations

We can define the density probability distribution P c of branch vector orientations as

P c (Θ ′ ) = N * c (Θ ′ ) N * c (3)
with N * c (Θ ′ ) being the number of branches pointing at angle Θ ′ at the onset of failure. The inset of Fig. 11(b) presents the angular distributions P c (Θ ′ ) with symbols for three different values of inherent anisotropy and loading orientations θ = 0 • and θ = 90 • . We can see that when η = 1 the distribution remains almost circular, highlighting the fact that the bond network is nearly isotropic and independent of the assembly's rotation. Conversely, the distributions for η > 1 present preferential orientations matching the preferred orientation of the cells.

These angular distributions can also be described using truncated Fourier series, as

P c (Θ ′ ) = 1 2π {1 + a ′ c cos 2 (Θ ′ -Θ ′ c )} , (4)
with Θ ′ c being the preferential orientation of the distribution and a ′ c its anisotropy level, i.e., the branch vector orientation anisotropy. Note that a ′ c = 0 means a circular distribution P c , in which bonds are equally presented in all orientations Θ ′ . Conversely, as a ′ c increases, more bonds present a preferential orientation in space. Although we could fit Eq. ( 4) to our measures to find a ′ c and Θ ′ c , we prefer to use the fabric tensor defined as [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular material[END_REF]]

F ij = ∫ π 0 P c (Θ ′ )n ′ i (Θ ′ )n ′ j (Θ ′ )dΘ ′ , (5)
with n ′ = {cos Θ ′ , sin Θ ′ }. Equation ( 5) lets us define the anisotropy of branch orientations as We can observe that a ′ c is close to zero for η = 1, exposing the isotropic character of the branch network and its independence before θ. For η > 1, a ′ c can reach larger values as high as ≃ 1.1 for η = 6. There is also a slight drop in a ′ c occurring for values of θ > 70 • . This phenomenon, combined with the drop of Z * we observed before, suggests that contacts are mostly lost in the minor orientation of the fabric tensor as η increases.

a ′ c = 2(F 1 -F 2 ), with F 1 > F 2 being
In Fig. 11(b), we present the evolution of Θ ′ c as a function of the loading orientation θ. For η = 1, Θ ′ c is irrelevant given that a ′ c ≃ 0, so it is omitted from the plot. But as soon as η > 1, Θ ′ c decreases as π/2θ which matches the cells' orientation. Along with P c , we can also characterize the angular branch length distribution ⟨ℓ⟩(Θ ′ ). This distribution can be computed as a function of Θ ′ as

⟨ℓ⟩(Θ ′ ) = 1 N * c (Θ ′ ) ∑ c∈δΘ ′ ℓ c , (6)
with ℓ c being the length of the branches pointing at small intervals of angular orientation δΘ ′ . The inset in Fig. 12 As with branch orientations, the angular evolution of branch lengths can be described using the expression

⟨ℓ⟩(Θ ′ ) = ⟨ℓ⟩ {1 + a ′ ℓ cos 2 (Θ ′ -Θ ′ ℓ )} , (7)
with Θ ′ ℓ being the preferential orientation, and a ′ ℓ the level of anisotropy. In the insets of Fig. 12(b), we show that the branch length distributions become more anisotropic as η increases, and that the longest branches predominately point in the same direction as the cells are pointing.

To find the branch length anisotropy, it is convenient to build the branch tensor, which is defined in an integral form as [80]

H ℓ ij = ∫ π 0 ⟨ℓ⟩(Θ ′ )n ′ i (Θ ′ )n ′ j (Θ ′ )dΘ ′ . ( 8
)
Note that this integral is computed in the range [0, π] given the periodic evolution of the angular distributions in that interval. We can then compute the branch length anisotropy as

a ′ ℓ = 2(H ℓ 1 -H ℓ 2 )/(H ℓ 1 + H ℓ 2 )
, with H ℓ 1 and H ℓ 2 being the eigenvalues of the tensor, so

H ℓ 1 > H ℓ 2 .
The same construction allows us to compute Θ ′ ℓ as the major principal orientation of H ℓ using the same approach as with the fabric tensor. In the inset of Fig. 12(b), we present Eq. ( 7) with dashed lines using the values extracted from the tensors above, and nicely fitting the angular distributions. We deliberately omitted the evolution of ⟨ℓ⟩(Θ ′ ) for η = 1 since the corresponding values for a ′ ℓ are negligible.

Figures 12(a) and 12(b) gather the results for the branch anisotropies a ′ ℓ and preferred orientations Θ ′ ℓ as a function of the loading orientation θ and the levels of inherent anisotropy η. They show that the branch length anisotropy increases with η from a ′ ℓ ≃ 0 for η = 1, up to a ′ ℓ ≃ 0.7 for η = 6. In all of the cases, these anisotropies present only minor variations with the loading orientation θ. We can then say that the variation of the preferred orientation for branch lengths evolves roughly as π/2θ, similarly to P c (Θ ′ ).

The large variations of geometrical anisotropies that we observed are

(a) (b)
Figure 12: Evolution of the normal branch length anisotropy (a) and its preferential angular orientation (b) as a function of the loading orientation and values of inherent anisotropy η. In the inset: angular distribution branch lengths for some values of inherent anisotropy η for loading orientations 0 • and 90 • . We also present Eq. ( 7), fitting the angular branch length distributions with dashed lines.

induced by the inherent anisotropy of the cells -and, furthermore, η is likely also inducing force transmission heterogeneities within the samples. In order to investigate this, we focus next on the interaction forces between cells.

Force orientations

Similarly to the analysis undertaken for branches, we can analyze the interaction forces between cells by using the angular distribution of radial and ortho-radial forces f ′ n and f ′ t , respectively. These angular distributions can be computed as

⟨f ′ n ⟩(Θ ′ ) = 1 N c (Θ ′ ) ∑ c∈δΘ ′ f ′ n , and (9)
⟨f ′ t ⟩(Θ ′ ) = 1 N c (Θ ′ ) ∑ c∈δΘ ′ f ′ t . (10)
In the insets of Fig. 13(b), we present the angular distribution ⟨f ′ n ⟩(Θ ′ ) for loading orientation θ = 45 • which shows the misalignment of the largest forces with respect to the vertical (i.e., the loading orientation). In the insets of Fig. 13(d), we present the distributions of ortho-radial forces for θ = 0 • and θ = 90 • , highlighting how widely these distributions vary as the assembly rotates. Regardless, these angular distributions remain periodic and smooth enough to fit Fourier series for their characterization. We can thus describe the angular variation of forces as

⟨f ′ n ⟩(Θ ′ ) = ⟨f ′ n ⟩ { 1 + a ′ fn cos 2 ( Θ ′ -Θ ′ fn )} , and (11)
⟨f ′ t ⟩(Θ ′ ) = ⟨f ′ n ⟩ { -a ′ ft sin 2 ( Θ ′ -Θ ′ ft )} , (12)
with a ′ fn and a ′ ft being the level of anisotropy for each distribution, and Θ ′ fn and Θ ′ ft the respective preferential orientations. For convenience, we build force tensors that allow us to easily compute the anisotropies and main orientation of each distribution as

H f ′ n ij = ∫ π 0 ⟨f ′ n ⟩(Θ ′ )n ′ i (Θ ′ )n j (Θ ′ )dΘ ′ , and (13)
H f ′ t ij = ∫ π 0 ⟨f ′ t ⟩(Θ ′ )n i (Θ ′ )t j (Θ ′ )dΘ ′ . (14)
This lets us compute the levels of force anisotropy as a ′ fn = 2

( H f ′ n 1 -H f ′ n 2 ) / ( H f ′ n 1 + H f ′ n 2
) for the radial forces and a ′ ft = 2

( H f ′ t 1 -H f ′ t 2 ) / ( H f ′ n 1 + H f ′ n 2 )
for the ortho-radial forces, where H α 1 and H α 2 are the eigenvalues of each one of the tensors. Note that H α 1 > H α 2 , and α stands for either the radial or ortho-radial components of the forces. It is worth mentioning that tr

( H f ′ n ) = ⟨f ′ n ⟩, i.e.
, the average radial force, and tr

( H f ′ t )
= 0 by equilibrium of force moments over the cells.

Figures 13(a In this case, we observe a continuous increase of a ′ ft as a function of η from loading orientations θ = 0 • to θ ≃ 45 • . Beyond that loading orientation, the ortho-radial force anisotropy reaches a plateau and barely varies with θ. For the orientations Θ ′ ft , we can see a variation that becomes more important as the level of inherent anisotropy increases. Although we might have expected a joint evolution of the preferential radial and ortho-radial forces' orientations, these figures show that a non-evident trade-off of force anisotropies occurs for highly anisotropic materials. This fact is, of course, emphasized by the point loading configuration which signs the force transmission at bonds.

Scaling up the strength from the microstructure

Microstructural contributions to the stress tensor

The previous microstructural parameters -concerning bonds, branches, and forces -must act together to produce the macroscopic failure strength we initially measured. This mapping between the micro and macro scales is especially challenging because of the varying shapes and sizes of the cells and the fact that the different microstructural tensors are strongly misaligned and evolving with θ.

In order to reconcile the micro and macro scales, let us consider the granular stress tensor as [START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Nicot | On the defi-835 nition of the stress tensor in granular media[END_REF]

σ ij = 1 V ∑ ∀c f c i ℓ c j , ( 15
)
where V is the volume of the sample, and the sum includes the dyadic product of the force f and branch ℓ vectors for all interactions c. Supposing that the distribution of forces and branches is uncorrelated (which is verified in our simulations), we can rewrite the stress tensor in terms of angular distributions on the frame {n ′ , t ′ } as [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular material[END_REF]

σ ij = n c ∫ π 0 {⟨f ′ n ⟩(Θ ′ )n ′ i (Θ ′ ) -⟨f ′ t ⟩(Θ ′ )t ′ i (Θ ′ )} ⟨ℓ⟩(Θ ′ )n ′ j (Θ ′ )P c (Θ ′ )dΘ ′ , (16)
with n c being the bond number density defined as N c /V = Z/(2⟨V cl ⟩), where ⟨V cl ⟩ is the average volume per cell. Note that we can also write ⟨V cl ⟩ = (π/4)⟨d cl ⟩ 2 , with ⟨d cl ⟩ being the equivalent average diameter of the cells. When we replace Eqs. ( 4), [START_REF] Zhang | Estimating the intensity of rock discontinuities[END_REF], and [START_REF] Pouragha | Micromechanical correlation between elasticity and strength characteristics of anisotropic rocks[END_REF] in the previous expression, focus only on the vertical component of the tensor (i.e., σ yy ), and integrate over the interval [0, π], we find a microstructural definition of the vertical stress at the onset of failure as

σ th yy = Z⟨f ′ n ⟩⟨ℓ⟩ π⟨d cl ⟩ 2 { 1 - 1 2 ∑ a ′ k cos 2 (Θ ′ k ) + 1 2 ∑ a ′ l a ′ m cos 2 (Θ ′ m -Θ ′ l ) + O } . ( 17
)
The term in brackets shows the contributions of the different anisotropies to the strength. In that term, the first sum runs in the set a ′ k ∈ {a ′ c , a ′ ℓ , a ′ fn , a ′ ft }, and in the respective values for Θ ′ k . The second sum is a product of anisotropies in which the combinations of indices l and m belong, respectively, to the set {a ′ fn a ′ ℓ , a ′ fn a ′ c , a ′ ℓ a ′ c }, with the respective angles for Θ ′ m and Θ ′ l . The higherorder term O involves triple products of anisotropies and is purposely neglected for the sake of simplicity. Also note that we added the superscript 'th' to emphasize that this value of strength results from the theoretical decomposition of the stress tensor. For simplicity, the term related to the anisotropies is henceforth written as A.

Equation [START_REF] Marinelli | Anisotropic breakage mechanics: From 38 stored energy to yielding in transversely isotropic granular rocks[END_REF] illuminates the fact that non-trivial microstructural compensations occur between 1) the different anisotropy levels, 2) the preferred orientations of angular distributions, and 3) geometrical and mechanical features. In addition, the choice of the branch frame {n ′ , t ′ } instead of the bond frame {n, t} is deliberate because it allowed us to reduce the number of anisotropies and the number of terms involved in A [START_REF] Azéma | Stress-strain behavior and geometrical properties of packings of elongated particles[END_REF].

In Fig. 14, we summarize the evolution of the different parameters involved in A for single and double anisotropies. On the one hand, we see that the geometrical anisotropies related to the branch orientation a ′ c and branch lengths a ′ ℓ smoothly decrease as a function of the loading orientation θ. On the other, the anisotropies related to the force transmission a ′ fn and a ′ ft have a highly non-linear evolution with θ. For angles between θ = 0 • and θ = 30 • , the radial force anisotropy increases but then finds a relatively steady value for larger loading orientations. For angles greater than θ ≃ 45 • , the ortho-radial anisotropy increases strongly as a function of θ. For the terms involving the product of anisotropies, the trends are all quite similar and not negligible in contribution.

At the bottom of Fig. 14, we see how all of these anisotropies add up.

Given the strong variation of all the anisotropies and preferential orientations, it is notable that the term A ends up fluctuating around the case η = 1.

This is clearly a mechanism involving direct compensations between geometrical microstructural characteristics and the force transmission at bonds. This phenomenon -in which the term A lies close to one -shows that the strong variation of the macroscopic failure strength must lie on the parameters Z⟨f ′ n ⟩⟨ℓ⟩ of the microstructural decomposition of stresses. Note that a version of Eq. ( 17) that neglects the term A has been used many times before for conglomerates or granular assemblies in which particles are of similar size and shape [START_REF] Rumpf | Zur theorie der zugfestigkeit von agglomeraten bei kraftübertragung an kontaktpunkten[END_REF][START_REF] Gröger | Modelling and measuring of cohesion in wet granular materials[END_REF][START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF]. Nonetheless, as we just observed, A cannot be neglected for ellongated bodies.

In Fig. 15(a), we see the evolution of the average branch length at the onset of failure as a function of the inherent anisotropy and the loading orientation, which shows a gradual drop as the loading becomes perpendicular to the layering. Such a variation is accentuated as η grows. In Fig. 15(b), we present the evolution of the average radial force ⟨f ′ n ⟩, which is normalized by the internal cohesion and the average cell equivalent diameter. This curve varies widely and, indeed, carries most of the shape of the macroscopic failure strength.

These observations allow us to conclude that the microstructural mechanisms producing the increase of failure strength with θ are related to the rise of radial forces, the drop of average branch length and coordination number Z, and the complex compensations occurring within the term A.

Finally, in Fig. 16, we plot σ th yy nicely reproducing the macroscopic vertical failure stress measured in Sec. 3. It is remarkable how the approach of decomposition of the granular stress tensor in terms of angular distributions, in spite of the large variability of anisotropies and preferred orientations for each distribution, is capable of providing these set of satisfactorily good predictions for the failure strength. The small differences between the measure and the decomposition are linked to the higher-order terms that were neglected in Eq. [START_REF] Marinelli | Anisotropic breakage mechanics: From 38 stored energy to yielding in transversely isotropic granular rocks[END_REF]. Thus, based on a fine description of the phenomena at bonds and mineral organization in space, our micromechanical description proves capable of describing and scaling up the macroscopic behavior we observe in laboratory.

Summary

545

We developed a series of numerical tests to study the failure strength of brittle materials reminiscent of schists, slates, shales, etc, whose components have a preferential orientation, i.e., an inherent anisotropy. Using a bi-dimensional discrete-element method, we built samples in which we could control the degree of inherent anisotropy by using a modified Voronoi tessellation. This approach allowed us to generate a set of subdivisions (the tessellation) of adjacent irregular polygons that we called cells. The common edges between cells interacted via cohesive bonds, enabling us to control both the failure strength and the cumulated surface energy necessary to produce fissuring. We then measured the macroscopic failure strength by applying a diametrical point load onto circular samples up to breakage. The failure strength turned out to be strongly affected by the layering orientation with respect to the loading direction θ. As observed in many experimental tests, the failure strength in our numerical tests evolved in a parabolic 'U ' shape, with a minimum value around a loading orientation of θ ≃ 25 • .

We also analyzed the variability of the failure strength using the Weibull survival probability, concluding that mixing several anisotropic configurations and loading orientations may lead to misleading conclusions upon the average failure strength and data scatter. This means that experimental test-ing must consider the microstructure of samples to avoid a misinterpretation of the strength of anisotropic brittle materials.

Finally, we performed a thorough characterization of geometrical properties of the cells' assemblies and force transmission mechanisms by means of the fabric, branch, and force tensors, as well as an approximation of their angular distributions using Fourier series. Exploiting the definition of the granular stress tensor in terms of angular contributions, we were able to find the microstructural elements that explain the variability of strength at the macroscopic scale. This was not a straightforward task. The strong geometrical and force anisotropies we found -as well as the misalignment of the different tensors -prompted us to undertake a full description of the contributions of anisotropies involving high-order terms seldom seen when analyzing rocks or granular media. Instead of simplifying particles' shape and size variability, we modelled the complexity of these materials in order to identify the microstructural elements responsible for the macroscopic phenomena. We found that geometrical and mechanical anisotropies present complex compensations, which means they are not the main source of the failure strength variations. Rather, it is the cell coordination, the average branch length, and the average radial forces that present the larger fluctuations -making them the key microstructural elements at the origin of the macroscopic failure strength.

Real materials are tremendously complex. Through this work, we sought to explore this complexity with the most detailed parameters we could gather linked to the granular stress tensor. Note that the circular shape we used for the samples was simply a choice of configuration that allowed us to compare our results to those obtained in rock testing. Our approach, however, is general and can be extended to any sample shape, assemblies of many crushable bodies, and diverse bonding behavior other than pure cohesive. Another very interesting perspective deals with the mechanics of porous rocks or grains for which our numerical model can be modified to not only reproduce adjacent cell configurations, but structures with holes in it. Many questions remain unresolved concerning the compaction or shear properties (rheology) of assemblies composed of several crushable inherently anisotropic bodies, which would benefit from future research.

Figure 1 :

 1 Figure 1: Scheme of a diametrical loading in which the loading orientation varies relative to the orientation of the internal structure.

Figure 2 :

 2 Figure 2: Interaction law for cohesive bonds between cells for the normal (left) and the tangential components (right) relative to the local framework coordinates.

Figure 3 :

 3 Figure 3: Samples presenting increasing average cell aspect ratio η. We varied η from

(a) 1

 1 up to (f) 6 in steps of 1.

Figure 4 :

 4 Figure 4: Different interaction types between convex polygons: (left) vertex-vertex, (center) vertex-edge, and (right) edge-edge.

Figure 5 (Figure 5 :

 55 Figure 5(b) summarizes the average values for σ yy found in our tests as a function of θ and the inherent anisotropy η. Under the same loading rates,

Figure 6

 6 Figure 6 presents the survival probability distribution P s as a function of the applied stress and the different values of η, combining all results by θ. As expected, the stress range σ within which we can expect failure considerably increases with η. The dashed lines correspond to the fitting of Eq. (2) by finding m and σ 0 with a least-squares minimization. The inset on the same figure presents the values found for the Weibull modulus m as a function of η.

Figure 7

 7 Figure7presents the fissuring paths for some of the samples with η = 2 and varying loading orientation θ. When θ = 0 • , the failure is roughly vertical, matching the loading orientation. The failure mechanism is similar when θ = 90 • , although the zones in contact with the platens show more damage. In these two cases, we can infer that tensile stresses are the source of fissuring because many interactions are debonded orthogonally to the loading

Figure 6 :Figure 7 :

 67 Figure 6: Probability of survival P s as a function of the applied stress σ for each value of η and all values of θ combined. The dashed lines correspond to the fit of Eq. (2) to the data. In the inset, we present the corresponding value of the Weibull modulus for each fit of the Weibull distribution. θ = 0 • θ = 30 • θ = 60 • θ = 90 •

Figure 8 :

 8 Figure 8: Schematic representation of the local frames created between two irregular cells i and j. Unitary vectors n and t are linked to the bond, while n ′ and t ′ are defined upon the branch vector. Note that angles Θ and Θ ′ are the orientation of n and n ′ , respectively, measured counterclockwise from the horizontal.

Figure 9 :

 9 Figure 9: (Top) Networks created by bonded adjacent cells for η = 1, 3 and 6 displayed with lines between the center of mass of the corresponding cells. (Bottom) Force networks between bonded cells. The thickness of the lines is proportional to the force intensity.

Figure 10 :

 10 Figure 10: Evolution of the coordination number at the onset of failure Z * as a function of the loading orientation θ and the inherent anisotropy η.

  the eigenvalues of F . The major principal direction of the fabric tensor is Θ ′ c = 1/2 arctan{2F xy /(F xx -F yy )}, with F xx and F yy being the components in the diagonal of F , and F xy the component off the diagonal. Figure 11 presents the corresponding values of preferential branch orientation and anisotropy as a function of η and θ.

Figure 11 :

 11 Figure 11: Bond orientation anisotropy a ′ c (a) and the preferential orientation Θ ′ c (b) for different inherent anisotropy levels η and loading orientations θ. In the inset: probability of branch orientations P c (Θ ′ ) shown with symbols for three values of inherent anisotropy η and loading orientation θ = 0 • and θ = 90 • . The dashed lines are the fitting curves using Eq. (4).

  (b) presents these angular distributions for loading orientations θ = 0 • and θ = 90 • , and three values 385 of η. These angular distributions for branch lengths closely follow the trends previously seen for P c (Θ ′ ).

Figure 13 :

 13 Figures 13(a) and 13(b) display the radial force anisotropies and their preferential orientations for the different levels of inherent anisotropy η and

Figure 14 :Figure 15 :

 1415 Figure 14: Evolution of the terms in Eq. (17) related to single anisotropies (first row) and double anisotropies (second row). We also present the sum of these different parameters in the term A (bottom).

Figure 16 :

 16 Figure 16: Failure strength measure through the wall forces just as in Fig. 5(b) (solid lines), and the same strength found using the microstructural decomposition of the stress tensor using Eq. (1) (dashed lines).

  Typical values for parameter m for silicate materials are found in the range[1.5, 4] [19, 75]. We observe that for microstructures with η = 1, m reaches a value of ≃ 2.5 and then smoothly decreases as values of inherent anisotropy grow. It is remarkable that our experiments satisfactorily reproduce the Weibull modulus for brittle silicate materials despite the strong variation of the cell's configuration, the simplicity of our model, and the fact that we are not able to store elastic energy/deformation in the bulk of the cells.Typical values for parameter m for silicate-rich materials have been found, for instance, in the range m ∈ [1.16, 1.93] for silica sand[START_REF] Mcdowell | The application of weibull statistics to the fracture of soil particles[END_REF], m ∈ [1.23, 3.04] [77] for different silica sands from Japan, m ∈ [3.34, 3.44] [78] for limestone sand, and m ∈ [1.26, 2.52] for calcareous and quartzite shale aggregate rocks [21],

agreeing with our results. Although our modeling strategy was developed targeting the brittle behavior of soil grains and rocks, higher values of m can be obtained for different materials by adjusting parameters γ, C n , and C t .

internal material strength (C n ) to the macroscale, and 2) it can be correctly assumed that for a non-anisotropic brittle structure, the tensile stresses are

Acknowledgments

This research work benefited from the financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) [Ref. RGPIN-2019-06118], the Fonds de recherche du Québec -Nature et technologies (FRQNT) through the Programme de recherche en partenariat sur le développement durable du secteur minier-II [Ref. 2020-MN-281267] and the industrial partners of the Research Institute on Mines and the Environment (RIME) UQAT-Polytechnique (irme.ca/en). This research was enabled in part by support provided by Compute Canada (www.computecanada.ca) under the Resources for Research Groups 2021 program (Project ID 3604).

We also thank Franny McGill for her careful editing of the manuscript.