Theoretical and computational investigation of the fracturing behavior of anisotropic geomaterials

The fracturing process in geomaterials is studied to characterize a potential host rock for radioactive waste, such as the kaolinite-rich Opalinus Clay formation. Because of its sedimentary genesis, this rock can be considered as a transversely isotropic geomaterial. A semi-circular bending test is here modeled based on the eXtended Finite Element Method (XFEM), to check for the formation and propagation of cracks in the rock, with a particular focus on the effect of notch dimensions and scale effects on the fracturing response of the specimen in terms of peak load. Starting with the XFEM-based results, a novel analytical formulation is also proposed to approximate the response of the material in terms of load-crack mouth opening displacement. The proposed formulation is also capable to provide a reliable estimate of the peak value and time history response, compared to some experimental predictions from literature, starting from a predefined value of initial notch depth, which could represent a useful theoretical tool for design purposes.

underground repositories [15][16][17]. At the same time, it is possible to study the structure and physical properties nearby a 1.5-5 m thick fault zone that cuts through the Opalinus claystone with an offset of about 10m [18]. From a geomechanical standpoint, OPA presents a complex elastic-plastic constitutive behavior, with a marked nonlinearity, anisotropic elasticity, and brittle failure [19,20]. Due to its spatial variability, it is also difficult to assign a unique set of geomechanical parameters to the material, but at the same time such variability must be considered for a proper study of the large-scale response of an underground. Different studies in literature have focused on the layered composition nature of soil and rock. In particular, in Ref. [21][22][23][24][25][26][27], an analytical solution was proposed to obtain the elastic properties of transversely isotropic layered materials from layers properties. Further works from literature analyzed experimentally [28] or numerically [29][30][31] the sensitivity of the mechanical response of geomaterials in terms of strength and stiffness, for varying layers from soft to hard layers. More recently, some empirical correlations were derived to define the elastic properties of transversely isotropic geomaterials [32]. In addition, based on the experimental results performed in MT, the geometry and extent of the excavation-disturbed zone around tunnels has revealed to depend on the orientation of bedding planes and excavation axes [33][34][35][36][37].
From a computational standpoint, despite the large use of finite element method (FEM) to study fracture mechanics problems, this technique is almost unable to follow automatically the crack creation and growth. To this end, the XFEM developed in [38][39][40][41] represents a valid computational way to solve cracking problems with a minimal remeshing, as efficiently verified in [42][43][44]. The method is based on the partition of unity theorem as conventional FEM, while introducing a novel enriched function to model discontinuities in the displacement field due to cracks. This is possible with the introduction of an additional term related to the enriched degrees of freedom [44,45], without any remeshing process during a crack propagation. The crack creation also depends on the mesh quality and density. Moës et al. [40] and Dolbow [46] for the first time introduced the possibility to represent the entire crack independently of the mesh, where the interaction between crack geometry and mesh discretization represents the base for an enriched approximation. A similar approach was implemented by Dolbow et al. [47][48][49], Daux et al. [50] and Sukumar et al. [51] for three-dimensional and arbitrarily branched cracks. In the last years, different works have successfully applied the XFEM for the study of anisotropic rocks and porous media, see Refs. [52][53][54][55]. In all these works, XFEM has revealed to overestimate the number of cracks within materials (along with the corresponding stress and strain field), and to overpredict by 20% the load-deflection curves compared to experimental measurements. Later on, a multiscale modeling approach was applied by Liu and Zheng [56] and Unger and Eckardt [57] for 2D problems, with reliable results compared to the experimental predictions from literature. More recently, Elruby et al. [58] have proposed a novel approach for an accurate prediction of cracking initiation and propagation, along with the corresponding failure load. In finite elements, there is the necessity of embedding an initial crack a priori, within an element mesh to trigger the crack propagation, which requires a certain experience of users to identify the crack propagation. Otherwise, to define XFEM nodal enrichments the cracking regions must be properly determined. An inexperienced user, indeed, could improperly select the entire domain nodes of a finite element mesh as enriching nodes, which could make the runtime longer with a lower computational efficiency. In such a context, a key aspect in the fracturing behavior of heterogeneous materials and structures relies on the possible presence of defects, notches and holes with weaken effects due to high stress concentrations.
In the last decades, different studies from the literature have increased their attention to the influence of notch shapes on the global fracturing response of composite specimens [31,[59][60][61][62][63][64][65][66][67][68][69][70][71][72], while applying incremental methods and advanced computational methods to follow multi-directional fractures [31,59,[61][62][63][64][65]68,71]. As far as OPA-based materials is concerned, the recent work by Dimitri et al. [31] applied the XFEM as computational approach to model semi-circular bending (SCB) tests, compared to the experimental evidences by Barpi et al. [36] and Valente et al. [37]. Starting with such studies, in the present work we reconsider the SCB test and explore the problem both theoretically and computationally, always compared to the existing literature [36,37]. From a computational standpoint, the problem is modeled via the XFEM, at least in the domain portion where the fracture is most likely to occur, in order to reduce the computational effort. The specimen is precracked with a rectangular notch, and it is subjected to a bending loading condition for different notch depths. Based on the preliminary study [31], which focused on the influence of the notch width and shape, in this work we investigate, for the first time, the sensitivity of the load-displacement response and maximum load for different notch depths. To approximate the results obtained for each XFEM simulation, a systematic formulation has been proposed to describe the relative curve shapes also for different scales. The main novelty of the work relies on the analytical formulation based on a best fitting procedure, to predict the fracturing behavior and load-displacement curves of specimens with different notch depths, starting from only one experimental or computational test, with a clear advantage in time and cost. More in detail, the present

Theoretical formulation
In this section, we briefly recall the main basics of XFEM, adopted herein as computational tool for the study of the fracturing process in anisotropic solids as OPA. Based on the proposed computational approach, generalized Heaviside enrichment functions are introduced as extra few degrees of freedom (DOFs).
The general equilibrium equation associated with the problem is defined in strong form as: and the associated boundary conditions are defined as: in which u , t , c stand for the displacement, traction, and crack boundaries, respectively, σ refers to the stress tensor, n is the normal vector, f b is the body force vector and f t is the external traction vector (see Fig. 1). The boundary value problem can be defined in variational form as: where δε and δu, refer to the virtual strain tensor, and virtual displacement field, respectively. Starting with such a variational formulation, the problem can be discretized by approximate functions, and the entire domain can be subdivided into geometrical subdomains [40].
As visible in Fig. 2, let us consider a two-dimensional cracked body, with a propagating crack represented by a dashed line, a finite element mesh and nodes represented by the solid lines and black filled circles (set A), respectively. In XFEM, two types of nodes surround the crack. The first type refers to the "Heaviside" enrichment nodes (blue circles-set B); the second type refers to the "Crack-tip" nodes (red squares-set C). The XFEM problem is based on a standard Galerkin method, as conventional FEM, due to the same global use of shape functions for all regular nodes A, with the introduction of some "enrichment nodes" (sets B and C) and extra terms to the global shape function, in order to account for discontinuities in the displacement field in addition to the crack-tip singularities.
In such a context, the generalized displacement field can be described in discrete form as: where n is the number of total nodes, m represents the number of nodes with a crack face in their support domain, m t1 and m t2 refer to nodes associated with the crack tips 1 and 2, respectively, in their influence area, u j stands for the nodal displacement field, a h is the vector of additional DOFs used to model crack faces, b l1 k and b l2 k refer to the additional DOFs associated with the two crack tips, F i l (x), i = 1, 2, m f , denote the crack tip enrichment functions, and H (ξ(x)) represents the generalized Heaviside function expressed in terms of the mapping function ξ(x). Note that the first term in Eq. (3) refers to the classical FEM-based approximation of the displacement field, while the other terms refer to the enrichment approximation related to any possible discontinuity, such as cracks.
The discretization of Eq. (2) by means of XFEM (3), leads to the following system of linear equilibrium equations: being u h the classical and enriched vector of nodal degrees of freedom, K the stiffness matrix, and f the external force vector. An extended version of the maximum circumferential stress theory [73] is, then, applied in order to treat the mixed-mode crack propagation in homogeneous anisotropic solids, like OPA. For such materials, two different values of mode-I stress intensity factors, K 1 I c and K 2 I c , must be defined ahead of the crack tip to define the brittle behavior of the crack propagation along the principal planes (1 and 2) of the elastic symmetry. As mentioned in Ref. [73], it is possible to determine separately two fracture toughness values, starting from the computation of K 1 I c , thus postulating that K 2 where E 1 and E 2 are the longitudinal elastic constants along the directions 1 and 2, respectively (Fig. 1). Considering the arbitrary orientation of a crack by an angle θ , (see Fig. 1), K θ I c can be defined as a function of K 1 I c and K 2 I c , i.e., According to a maximum circumferential tensile stress model for anisotropic bodies, we can follow the crack propagation in mixed-mode conditions [73], based on the following definition of a maximum tangential stress σ θ : Re where t 2 1 = (s 2 sin θ + cos θ) 3 , t 2 2 = (s 1 sin θ + cos θ) 3 , and K I I = √ 2πr σ xy refers to the stress intensity factor in Mode-II, for r → 0. By maximizing the σ θ σ max θ ratio, the crack propagation angle θ is obtained as [73]: Thus, in XFEM, the enrichment functions at a crack tip can be described analytically in terms of the polar coordinates r and θ , as: where whereas x 1 , x 2 are the Cartesian coordinates in the local coordinate system at each crack tip (Fig. 2), and C i j stand for the elasticity constants of the material. At a micro-scale level, the crack initiation refers to the degradation onset of the cohesive response at an enriched element (Fig. 3). The degradation process begins when the stress field (or the strain field, depending on the selected criterion) satisfies a fixed crack initiation criterion. In such a context, a maximum principal stress criterion is here selected such that where σ 0 max is the maximum principal stress, symbol · represents the Macaulay bracket, which means that a purely compressive stress state does not initiate damage. The damage state is, thus, assumed to initiate when the maximum principal stress ratio reaches a unit value. Then, an additional crack is introduced, or the crack length of an existing crack is extended after an equilibrium increment, when ϕ reaches the unit value with a certain tolerance ϕ tol , i.e., 1.0 ≤ ϕ ≤ 1.0 + ϕ tol .
In the following, a XFEM-based cohesive model is implemented to follow the damage propagation within OPA, similar to formulation based on particle systems [74]. This means that the damage evolution follows the degradation rate of the cohesive stiffness when the initiation criterion is reached. A cohesive traction-separation law is, thus, defined both in normal and tangential directions as: T n , T s and T t being the normal and shear stress components expected by the elastic traction-separation behavior for an undamaged strain field. For a compressive traction, T n < 0, no damage occurs and the material keeps its stiffness constant. Moreover,D is a scalar damage parameter, which represents the averaged overall damage at the intersection between the crack surfaces and the cracked elements edges. Such parameter starts with a null value in undamaged conditions, and evolves monotonically from 0 to 1 during the damage evolution, reaching the unit value, for a complete damage state. In the present work, we select an exponential damage evolution, defined as:D where δ 0 m is the effective displacement at damage initiation, δ max m is the maximum value of the effective displacement achieved during a loading process, δ f m is the effective displacement at complete failure, and α is a dimensionless quantity defining the rate of damage evolution, here set as α = 5.

Computational problem
In what follows, we model the OPA-based SCB test (Fig. 4), as studied experimentally in Refs. [36,37] within a Mont Terri Underground Rock Laboratory project, and analyzed computationally in Ref. [31] based on a classical cohesive zone modeling (CZM) and XFEM. The layered geomaterial is modeled as a transversally isotropic material, with bedding planes rotated by 45°with respect to the x-axis (Fig. 4), and elastic properties E x = E y = 10 GPa, E z = 4 GPa, ν xy = 0.33, ν xz = ν yz = 0.24, G xz =G yz = 1.2 GPa, being x y the isotropic plane. In line with Ref. [31], the specimen is characterized by a diameter d = 80 mm, span l = 62 mm between supports, thickness t = 1 mm, and a varying notch depth a from 1.4 up to 11.75 mm in  order to analyze the sensitivity of its mechanical response, both in terms of load-CMOD curve and peak-load. As plotted in Fig. 5, it can be noticed that for higher values of notch depth, the elastic stiffness of the specimen decreases, along with a decreased peak load. Such maximum value seems, also, to be attained for an increased CMOD. As also shown in Fig. 6, the peak load seems to vary linearly with the notch depth of the specimen, with a clear reduction of 45% of P max within the selected range of a.
The analysis continues considering the possible scale effect on the mechanical response, as plotted in the time-history response of Fig. 7 for specimens modeled with 3D elements (i.e., 10-node quadratic tetrahedron), different diameters d = 80 − 1280 mm, and a proportional increase in thickness and notch dimensions, starting from a notch width of 2 mm and a notch depth of 6 mm for d = 80 mm. Based on the plots in Fig. 7, it is worth observing that bigger specimens enable an increased peak load, and a larger softening branch, except for the highest dimensions here considered for the analyses, which reveal an increased brittleness in the fracturing response. Due to the high computational effort required for 3D XFEM-based discretizations, especially for specimens with higher dimensions, it is possible to limit the analyses in 2D, due to the prevalence of the inplane geometrical dimensions compared to the out-of-plane ones, with a clear computational benefit. For 2D analyses, we adopt 4-node bilinear plane stress quadrilateral elements, while assuming an equivalent isotropic material in lieu of the out-of-plane transversally isotropic 3D approximation. Figure 8 shows the comparative force-CMOD response for a specimen with d = 80 mm and t = 1 mm, as provided by a 2D and 3D XFEM  Comparative results for a specimen modeled with 3D and 2D elements modeling, with a reasonable agreement among results, at least in the ascending branch of the curves and peak load. Such 2D analyses have been, then, repeated systematically for SCBs with higher dimensions, focusing the attention on the scale effect on the maximum load P max and the related mode-I fracture toughness K I c , determined as: [75,76] In Fig. 9, we plot the fracture toughness for an increased core diameter, as provided by a 3D XFEM-based formulation compared to a more simplified 2D XFEM in the double plane-stress/plane-strain version. Such results are also compared to predictions based on a classical 3D FEM cohesive modeling, as proposed in the previous work by Dimitri et al. [31], with a satisfactory correspondence among them especially for lower scales, despite the different computational effort stemming from each approach. A clear increase in K I c is noticed for higher diameters of the specimen, in all the plots of Fig. 9, where the results seem to be conditioned by the selected method for higher geometrical scales. In these cases, indeed, the material becomes even more brittle, and the response features some instability effects during the fracturing process and, in particular, for the finite stress computed at the crack tip and K I c , which is very sensitive to the mesh itself. In such cases, a convergence analysis would be necessary because of the high strain gradient localized around the crack, Fig. 9 Fracture toughness K I c for the specimen modeled with: 3D and cohesive elements; 2D elements using XFEM especially in absence of any penalization within the formulation, see, e.g., Refs. [77][78][79][80][81]. It is worth observing, also, that FEM-based CZM is much more stable than XFEM, although it provides more conservative results than what expected by XFEM.

Numerical investigation
Based on the smooth elasto-softening behavior observed computationally in the load-displacement response of OPA specimens, we want now to determine an analytical approximation that can predict well such behavior, and could serve to designers as useful theoretical tool.
To this end, we start considering an objective function of the type where f represents the applied load, the variable x corresponds to the CMOD, α 1 , α 2 and α 3 are parameters that must be properly calibrated to best fit the computational results. Among such parameters, α 3 can be easily determined considering that the CMOD is null in the absence of loads. This corresponds to the enforcement of f (x = 0) = 0, which yields α 3 = 0. The other two parameters, α 1 and α 2 , must be determined through a fitting procedure and a minimization of the error among theoretical and computational predictions. A least squares method is thus applied such that y i and x i correspond to the load and CMOD computational values, whereas f (x i ) stands for the objective function value calculated at the same CMOD. Let consider now the i-th error defined as: which depends on parameters α 1 and α 2 , for a given y i , x i and f (x i ). According to the least square method, it is necessary to minimize the sum of squares of all these errors, namely [82] where α T = (α 1 , α 2 ) is the unknown parameter vector, and H(α) is the vector function defined as: The sum of squared errors is multiplied by 1 2 to simplify subsequent differentiations.
Since the equations that minimize the sum of squared errors are transcendental, the problem becomes nonlinear, which corresponds to a classic nonlinear least squares data fitting problem. The solution of this problem is approximate and relies on the computation of direction p using the formula for the Newton's method The components of gradient ∇h(α) are obtained by the following chain rule: whereas the Hessian ∇ 2 h(α) is calculated by its further differentiation to get Once computing both terms in Eq. (21), it is possible to solve the initial problem. In order to reduce the computational cost, the previous Hessian's expression can be simplified, as in the case of Gauss-Newton or Levenberg-Marquardt methods [83,84]. Given the small scale of the problem, the calculation effort is limited such that a full Newton method has been finally applied. An automatic procedure has been implemented to approximate all the curves acquired from the XFEM simulations, whereby the values α 1 and α 2 have been estimated for each initial notch depth. As a result, the P-CMOD graphs obtained through the proposed model have been plotted, with a very good agreement against the computational predictions both in the ascending and descending branches of the curves (see Fig. 10 for a specimen with a = 9.45 mm). This is due to the fact that the fitting procedure has been applied to points belonging to the range between the origin and a value just after the peak load, where the computational curves have shown a satisfactory agreement with the experimental data from Ref. [31].
For each curve obtained from the data fitting, the peak load has been determined as a function of the initial notch depth, whose numerical value has been compared to the computational predictions from XFEM (see Fig.  11). Similarly, the measured values of CMOD were extrapolated from both the computational and numerical results at the peak load, whose comparative evaluation is successfully reported in Fig. 12 for different values of a.

Previsional approach
We want now to provide a theoretical prediction of the maximum point of the curves (which corresponds to the maximum load P max supported by the specimen). Thus, the differentiation of the objective function has allowed to derive a relation between the peak load P max and the values of α 1 , α 2 , i.e., Based on a simple mathematical manipulation, the value of α 2 can be related to the CMOD at the peak load as follows: This result suggests that it is possible to know the value of α 1 once the peak load and the corresponding CMOD are known. In particular, α 1 can be expressed as: From the combination of such partial results, it is possible to obtain the general definition of the load-CMOD constitutive law as follows: This expression is then defined once the peak load and the corresponding CMOD value are obtained by an experimental test. From the various results carried out, it is also known how the peak load and parameters Fig. 13 Comparative results in terms of α 1 for different values of a α 1 and α 2 vary with the initial notch depth. Since the relationship between such three parameters has been developed, it is sufficient to derive the variation laws for two of these quantities in order to obtain a previsional tool for any dimension of the notch depth. In particular, it seems that the peak load and α 1 represent the best key parameters that enable easier fitting procedures.
The peak load values follow an almost linear pattern within the selected range. If a simple linear data fitting is performed, the trendline assumes the following expression: . This result would be physically inconsistent, since a peak load of about 8.33 N was found for such notch depth. Nevertheless, the approximation seems to be valid, as the most realistic values of the initial dimension of the notch oscillate around 6mm, which corresponds to the central value of the range considered in this discussion (see Fig. 11). A linear data fitting has been performed, once again, to define the analytical relation between α 1 and the notch depth a, with the only difference of using a natural logarithm function as objective function. As plotted in Fig. 13, the α 1 − a diagram can be defined by the following relation: where k = −2631.272 [N/mm] is obtained by the least squares method, and a 0 corresponds to the notch size value for which α 1 (a 0 ) = 0, in line with the previous findings. Once the variation laws are obtained for the peak load and α 1 , it is possible to derive the variation of α 2 with a, as plotted in Fig. 14, with a reasonable agreement among results from a computational, numerical and theoretical perspectives. Such constitutive law allows the entire curve to be derived from the peak load values and the corresponding CMOD, as finally plotted in Fig. 15. In this figure, it is possible to compare the analytical predictions in terms of force vs. CMOD curves with the experimental results from [37] referring to a specimen with thickness t = 30 mm. The very good matching among the analytical and experimental curves in Fig. 15, both in the ascending and softening branches, confirms the reliability of the proposed theoretical approach to estimate the peak values and the entire constitutive curve, with the only knowledge of the initial notch depth, and favorable advantages for design purposes. The proposed formulation does not account any plastic phenomena, due to the brittle nature of the OPA material, and will be extended in a further work to include such possible phenomena [85].

Conclusions
The present work has focused on the fracturing behavior of an OPA-based three-point bending (SCB) test, from a XFEM-based computational and theoretical perspectives, while focusing on the sensitivity of the response to Comparative evaluation of the analytical results against the experimental predictions from Ref. [36] the notch depth and scale effect. Based on a systematic investigation, it is found that higher precrack depths yield a decreased elastic stiffness and lower P max values, which are reached for higher CMODs. The scale effect has been also investigated by increasing gradually both the in-plane and out-of-plane dimensions of the specimen. This enables a monotonic increase of P max , a larger softening branch, while keeping constant the initial elastic stiffness of all the curves. XFEM tests for larger-scale specimens also show a monotonic increase in the fracture toughness K I c together with an increased attitude of the geomaterial to become quasi-brittle (as a concrete-like material), a severe dependence of the response to the mesh size, ahead of the crack tip, together with an increased numerical instability in the overall response. Smaller dimensions of the specimen, instead, exhibit a smooth load-CMOD response with continuous derivatives both in the ascending and softening branches. Starting with the computational findings based on XFEM, the work has originally described such fracturing behavior of OPA semi-circular specimens, from a theoretical standpoint based on the definition of a best fitting function. According to the proposed procedure, two unknown parameters have been determined theoretically by means of a nonlinear least squares method, whose results have been successfully compared to the computational predictions, especially in the first part of the plots, making possible to determine quite accurately both the peak load and relative CMOD, compared to computational predictions. Simple mathematical manipulations have also allowed to obtain an accurate relation describing the two unknown parameters as a function of P max and its corresponding CMOD at the onset of the softening stage. Based on the same procedure, it has been also possible to predict the overall fracturing response of a SCB test for a wide range of notch depth values, just referring to one experimental or computational test, with a clear advantage in time and costs, due to the reduced number of experimental and/or computational tests required for validation purposes. A further extension of the present study will include the arbitrary distribution of the bedding properties within OPA, together with the influence of the notch inclination within SCB specimens with respect to the bedding planes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open access funding provided by Università del Salento within the CRUI-CARE Agreement.