Crack Patterns in Heterogenous Rocks Using a Combined Phase Field-Cohesive Interface Modeling Approach : A Numerical Study

Rock fracture in geo-materials is a complex phenomenon due to its intrinsic characteristics and the potential external loading conditions. As a result, these materials can experience intricate fracture patterns endowing various cracking phenomena such as: branching, coalescence, shielding, and amplification, among many others. In this article, we present a numerical investigation concerning the applicability of an original bulk-interface fracture simulation technique to trigger such phenomena within the context of the phase field approach for fracture. In particular, the prediction of failure patterns in heterogenous rock masses with brittle response is accomplished through the current methodology by combining the phase field approach for intact rock failure and the cohesive interface-like modeling approach for its application in joint fracture. Predictions from the present technique are first validated against Brazilian test results, which were developed using alternative phase field methods, and with respect to specimens subjected to different loading case and whose corresponding definitions are characterized by the presence of single and multiple flaws. Subsequently, the numerical study is extended to the analysis of heterogeneous rock masses including joints that separate different potential lithologies, leading to tortuous crack paths, which are observed in many practical situations.


Introduction
Fracture events in geological materials are phenomena of notable importance for the safety and stability of geological masses.This is a relevant issue to engineering applications including tunneling procedures, hydraulic fracture for energy reservoirs, deep foundations, to name a few.Due to the advent of novel modeling techniques and the advances in the computational capabilities, the usage of numerical methods for reliable predictions of failure in rock masses has suffered a tremendous improvement in the last few years, becoming a plausible alternative for the analysis of geo-materials.
In this setting, the presence of flaws such as defects/discontinuities in rocks leads to a substantial reduction of the strength of rock mass, as compared to the intact rock.In line with linear elastic fracture mechanics (LEFM), the singular character of the stress field around flaws can induce the propagation of such defects when the stress state within the rock mass is altered, for instance due to a tunnel excavation.Therefore, the existence, initiation, and propagation of fracture not only affect the mechanical response of the rock mass, but also the construction of the engineering structures.This dependence of the strength and deformability of the rock mass upon the intact rock and the characteristic discontinuities (fracture and joints) make the deep understanding of such fracture processes a challenging task.
Stemming from the complex nature of the rock mass, according to [1,2], it is possible to distinguish between crack initiation and propagation phenomena.Crack initiation is defined as the process through which pre-existing flaws initiate in a confined region in the neighborhood of a crack tip.Crack propagation events are those processes by which the pre-existing cracks continue to grow after the initiation.Thus, propagation phases can lead to crack coalescence and/or crack branching or bifurcation, among different scenarios, depending on the characteristics of the rock mass and the external loading conditions.With the aim of achieving a deep understanding of these phenomena, a great deal of research has been conducted in literature to characterize the crack propagation in rocks, which are succinctly discussed from the experimental and numerical standpoints in Sections 1.1 and 1.2, respectively.

Experimental Investigations
In a landmark investigation, Lajtai [3] performed numerous experimental tests on gypsum of Paris specimens, with single flaws at different inclinations with respect to the direction of the external load.These tests led to the identification of different fracture stages of the rock mass under uniaxial compression.Subsequently Ingraffea et al., [4] have extended the experiments on Limestone and Granodiorite specimens, by orienting the flaws at a particular angle.In [4], the authors reported that the completed crack growth sequence upon failure, which allowed the identification of a complete variety of failure modes.Some recent studies on rocks with several flaws, which were oriented following different inclinations, manifested the complex nature of fracture events in rocks with evident crack coalescence [5][6][7][8].
Recently, several researchers have studied the response of rock mass on crack arrest, considering: (i) the number of defects, (ii) the relative inclination of the flaws and (iii) the applied external loads.Mentionable studies on experimental methods have thoroughly investigated uniaxial and biaxial compression tests in rocks, through which it is possible to identify that these materials may experience a combination of flaw slippage onset and wing propagation, followed by secondary cracks (shear dominated) leading to subsequent potential coalescence [5][6][7][8][9][10].In particular, Sagong and Bobet [11] distinguished 9 coalescence types of cracks in rock masses.
Nevertheless, it is worth mentioning that the experimental studies on this topic have been profusely exploited in recent decades.Standard methods are often limited due to notable difficulties in the simulation of general loading tests and the preparation of specimens for this purpose.Under these circumstances, traditional experimental procedures are also complemented with the use of digital image correlation (DIC) techniques and 3D printing methods [12], fostering a new perspective for the generation of specimens with desired attributes.

Numerical Investigations
As mentioned above, many of the previous experimental studies were focused on uniaxial or biaxial compressive loads, since the general loading conditions are difficult to reproduce in the laboratory.This fact motivated the development of computational methods which can potentially capture complex failure mechanisms in rocks.Moreover, numerical methods generally lead to notable cost savings, as compared to arduous experimental campaigns.
Modeling the physical behavior of fracture phenomena in geo-materials has been traditionally tackled using linear elastic fracture mechanics-based methodologies, due to their relative simplicity and representativeness, provided that the confining pressure and loading rate are sufficiently low.
From the modeling point of view, previous experimental studies have been employed as reference investigations to calibrate the numerical methods and to provide a deeper insight into the physical processes associated with fracture events in rocks.However, several difficulties have been identified in order to capture the wide variety of failure modes of such materials under different loading cases such as bidirectional traction-compression, traction-shear loading among many others, which evolve from the complex crack path when it is unknown a priori.In addition, tracking such crack paths will be further complicated by the development of branching and coalescence evolutions, which are typical scenarios encountered in geo-materials.
For this purpose, numerous computational methods have been proposed in the last three decades, especially within the context of nonlinear finite element method (FEM) and boundary element method (BEM), to simulate the complex failure modes in rocks.Later on, advanced computational approaches such as: cracking particles method (CPM) [13][14][15], Peridynamics [16] and dual-horizon Peridynamics [17,18] are developed, which neither do not require an explicit representation of the discontinuities in the displacement field due to fracture nor crack tracking algorithms for the identification of the propagation path.
Exploiting the BEM-based techniques, Einstein and coauthors [19,20] have developed reliable modeling tools to predict crack growth under the loading conditions leading to fracture in Modes I and II, considering the specimens with one or more internal flaws.However, due to the inherent versatility of FEM-based techniques, they are often preferred for fracture applications over alternative methodologies.Within the current modeling techniques to simulate fracture events in solids using FEM, the following categories can be differentiated:

•
Methods based on continuum damage mechanics (CDM) encompass the stiffness degradation, by defining a set of internal-damage variables [21].In their local forms, CDM models suffer from well-known pathological sensitivity of the corresponding estimations on the spatial discretization (finite element mesh) due to the loss of ellipticity of the corresponding governing equations.
To mitigate such deficiency, integral regularization schemes and gradient-enhanced formulations have been proposed [22][23][24][25][26], apart from the alternative procedures introduced by Areias and coauthors [27] by combining the local damage formulation with local re-meshing techniques.

•
Explicit crack tools can be employed to envisage a strong discontinuity in the displacement field.They rely on either the enrichment of the nodal displacement field at the element level using the partition of unity methods (PUM) [28][29][30], or the enrichment of the strain field (enhanced FEM, E-FEM) [31][32][33][34].
Although the above-mentioned methods have been extensively employed to simulate fracture events in different engineering applications, they present several limitations, particularly to capture the onset of multiple cracks and their evolution, apart from tracking the intricate fracture scenarios including branching and coalescence.This latter aspect is of remarkable complexity in strong discontinuity techniques, which suffer from topological operative problems under such conditions.
Due to the advances in the computational capabilities in the last decade, a renovated interest on diffusive damage models complying with nonlocal representation has been witnessed, which endow a scalar-based evolution equation for the damage field.Especially suitable for brittle fracture events, the so-called phase field (PF) approach of fracture envisages a diffusive crack representation with a nonlocal character [42].This new concept shares some common aspects with respect to CDM formulations via the definition of a scalar PF variable, which varies from 0 to 1, and triggers the stiffness degradation upon failure.
PF methods can be understood as a regularized version of the Griffith fracture approach [43], through the advocation of the Γ-convergence [44][45][46], whereby a direct competition between the elastic and the fracture energies is constructed in an integrated global functional of the system.The corresponding minimization of such functional allows capturing the initiation, growth, and propagation of cracks in solids in a robust and reliable manner.Notable contributions in this direction are due to Francfort and Marigo and their coauthors from a fundamental standpoint, see [42,47] and the associated references.Additional formulations exploiting the PF approach for brittle fracture have been rigorously developed by Bourdin et al. [47,48].Whereas, several modified models have been recently proposed complying with a thermodynamically consistent framework, e.g., see [49][50][51][52].Exploiting this potential, posterior studies, the PF approach has been extended to fracture studies in thin-walled structures [53][54][55], ductile failure [56,57], cohesive fracture [58], dynamic applications [59,60], multi-physics environments [61][62][63].Recently, this methodology has also been applied to study the initiation of fracture in rocks with multiple flaws [64][65][66].

Objectives and Organization
Due to the appealing features of the PF method for fracture, the main objective of this work is focused on the application of the PF method and its combination with the interface-like cohesive method to model fracture events in heterogeneous rocks [67].The advocation of the methodology proposed in [68], which was originally developed by the authors, allows the possible interaction of rock masses with different lithologies, which are separated by a series of predefined discontinuities.The analysis and understanding of interfaces in rocks have always been the subject of a great impact in the research community due to their importance in practice.This stems from the fact that they possess a certain amount of cohesion due to interlocking and roughness, along with frictional responses in case of large sliding modifying the response of the rock mass [69].Therefore, the aim of this research is to address the potential predictive capabilities of the combined PF-CZM for its application in rock fracture.Based on this interest, the proposed technique is first validated by reproducing experimental results concerning fracture of intact rock and subsequently extended to the analysis of heterogeneous media.
The manuscript is organized as follows: Section 2 summarizes the basic concepts of the PF method for bulk fracture, whereas, the corresponding simulations are presented in Section 3. Posteriorly, general details concerning the combined PF-CZM for triggering fracture events in heterogenous rock mass is outlined in Section 4. Some representative applications to illustrate the applicability of the proposed method are presented in Section 5. Finally, the main conclusions of the present study are given in Section 6.

Fundamentals of the Phase Field Approach for Brittle Fracture
The basic concepts of the PF approach for fracture in rocks are briefly introduced in this section.The current formulation advocates the methodology proposed in [47,49], therefore, detailed derivations are herewith omitted for brevity.

Basic Concepts
Within the context of multi-dimensional analysis, see Figure 1, we consider an arbitrary body Ω ∈ R nd in the Euclidean space of dimension n d , whose material points are denoted by the position vectors x in the Cartesian setting.Body actions are identified by the vector f v : Ω −→ R n d .The delimiting boundary Ω is defined such that ∂Ω ∈ R n d −1 , where ∂Ω u and ∂Ω t denote the boundary portions in which kinematic and static surface actions are prescribed, respectively, complying to: The kinematic field is denoted by u, whereas the second-order Cauchy stress tensor is identified by σ.Boundary conditions in this arbitrary body render: where n is the outward normal to the body.
where ∆d is the Laplacian of d, and ∇ x d stands for the spatial gradient.At this point, relevant contributions on this matter are due to the contributions by Wu and coauthors [70][71][72], who extended the applicability the PF method for cohesive and length-scale insensitive fracture.
To account for the regularization, a crack density functional γ(d, ∇ x d) is defined, which can be estimated using the relation: As discussed in [42,45,73], Equation (3) represents the crack surface Γ c , in the limit l → 0. Furthermore, to prevent crack healing upon unloading, the thermodynamic consistency condition given by Γc (d) ≥ 0, is required.The crack density functional can be then expressed as [49]: The above definition of γ(d, ∇ x d) allows the approximation of the crack surface integral in the original Griffith formulation by a volume integral as follows: Therefore, the total potential of the body can be split into the internal and external counterparts represented, Π int (u, d) and Π ext (u), respectively: where the first term in the integral represents the elastic energy of the body and the second term accounts for the energy dissipation.
To simulate the crack growth under tensile load, we adopt the spectral decomposition of the infinitesimal strain tensor between its corresponding positive and negative counterparts [49]: ) where λ and µ are the Lamé's constants, and ε + and ε − respectively identify the tensile and compressive counterparts of the strain tensor.In Equation ( 7), the symbol tr[•] represents the trace operator, while • ± stands for the so-called Macaulay bracket: Finally, the degradation function g(d), which introduces a stiffness deterioration with the growth of the defect (crack), can be defined as: where K is a residual stiffness parameter to prevent numerical instability issues for fully degraded states, i.e., when d ≈ 1.Also, note that according to Equation (7a), the degradation function (Equation ( 8)) exclusively affects the positive part of the elastic energy of the body.Thus, the corresponding Cauchy stress tensor in a decomposed fashion for the PF formulation renders: In the Equation ( 9), 1 represents the second-order identity matrix and σ ± identify the positive and negative parts of the stress tensor.

Numerical Implementation
In this Section, the numerical implementation of the PF approach for the fracture of intact rock, within the context of nonlinear FEM are briefly addressed.In the sequel, the corresponding to bulk operators will be denoted by the label b in the superscript.
Based on the isoparametric interpolation, the LaGrangian shape functions at the element level N I (ξ), defined in the natural space ξ = {ξ 1 , ξ 2 }, are used for the interpolation of the geometry (x).The displacement field (u), as well as the variation (δu) and the linearization (∆u) of u are defined as: where n indicates the number of nodes at the element level, and x I y d I denotes the nodal coordinates and displacements, respectively, which are collected into the operator x y d.Accordingly, the interpolation functions are arranged into the vector form N.
In the similar lines, the strain field (ε), its variation (δε) and its linearization (∆ε) can be approximated using the compatibility operator B d associated with the kinematic field: Similarly, the interpolation of the PF variable (damage-like variable), its variation and linearization render: where d I are the nodal values of the PF variable, which are arranged into the vector d.The spatial gradient of the PF (∇ x d), its variation (∇ x δd) and its linearization (∇ x ∆d) can be interpolated through the operator B d : Therefore, the discrete version of Equation ( 6) at the element level, identified by the super-index el, can be expressed as: where the residual vector and f b d,ext represent the internal and external forces, respectively, given by: Correspondingly, f b d is the residual vector associated with the PF variable, which can be estimated as: In this study, a staggered solution scheme in the form of Jacobi-type iterative procedure is adopted for the numerical treatment of the governing equations.Thus, the consistent linearization of Equation ( 14) yields the following coupled system [57]: The particular forms of the matrices K b dd , K b dd , K b dd and K b dd are omitted here for brevity [74].The numerical implementation of the PF approach is carried out within the framework of the FE codes FEAP and ABAQUS, by generating the user-defined element routines.

Validation of the Developed PF Approach for Fracture Using the Brazilian Splitting Test
The developed PF methodology for fracture of intact rocks is validated by comparing the results with the standard Brazilian splitting test reported in [64].The Brazilian test is extensively employed to characterize the tensile strength of rocks.The geometric definition of the benchmark problem involves a specimen of diameter equal to 52 mm, see Figure 2a.Furthermore, the following material parameters are considered: Young's modulus E = 31.5GPa, Poisson's ratio ν = 0.25 and G b c = 100 J/m 2 .Furthermore, the length scale in the PF approach can be related to the material strength σ c of one-dimensional bar: Therefore, in the current application, we set l = 1 mm.Half of the specimen is discretized using 55,452 elements, the simulations are performed by controlling the displacement.Figure 2b depicts the crack path estimated by implementing the present PF method in the ABAQUS platform, in line with [64].During the simulations, it is observed that the crack initiated at the center of the disc, which agrees with the experimental evidence.Upon further loading, the initiated crack starts to grow towards the upper end of the Brazilian disc.From the qualitative standpoint, Figure 3 shows the experimental-numerical correlation, whereby a very satisfactory agreement can be observed along the loading path, as compared to the PF methods in [64].In particular, the simulation response is characterized by a linear evolution until the abrupt failure and hence the fracture.Deviation of the simulated results from the experimental data is attributed to the absence of contact modeling features in the simulations.

Numerical Predictions of Rock Fracture: Estimation of Damage Patterns Various Load Cases for Specimens with Multiple Flaws
In this section, several examples are presented to illustrate the potential of the developed PF approach to simulate fracture phenomena in rock specimens, with an inclined echelon flaws (with respect to the horizontal direction), subjected to different load cases.
The basic geometry under consideration is a rock specimen with dimensions 200 mm × 200 mm, consisting of four inclined echelon flaws which can be simplified due to symmetry, to a quarter specimen of size 50 mm × 50 mm with a single inclined flaw, as shown in Figure 4.The adopted mechanical properties in the simulation are: E = 30 GPa, ν = 0.3, G c = 3 J/m 2 and l = 0.1 mm [64].
The domain is discretized using first-order iso-parametric elements equipped with PF attributes, and a uniform mesh of size h = 0.2 mm.The specimen is subjected to displacements imposed along the vertical and horizontal directions, see Figure 4.In particular, the maximum vertical displacement is prescribed as ūy = ±3 × 10 −2 mm (tensile/compression), whereas ūx can have different ratios with respect to ūy depending on the specific case.Simulations are performed under controlled displacement, with a maximum displacement increment equal to ∆u = 1 × 10 −5 mm using a monolithic PF formulation.

Uniaxial Compression Test
We first analyze the deformation patterns in a uniaxial compression test, which is achieved by specifying ūy = −3 × 10 −2 mm, in several steps and ūx = 0, see Figure 4. Figure 5 shows the failure maps, highlighting the fracture patterns, when the specified displacement reaches ūy = −2.4× 10 −2 mm (left) and ūy = −2.78× 10 −2 mm (right), respectively.Moreover, according to Figure 5, two cracks are observed to first initiate around the tips of the flaw, which are subsequently propagated parallel to the loading direction.These predictions are in good agreement with experimental observations, see [64].

Uniaxial Tensile Test
In the next step, the fracture patterns in a uniaxial tensile test are captured using the developed PF approach.Therefore, a tensile load along the y direction equal to ūy = 3 × 10 −2 mm, is specified on to the system shown in Figure 4.The simulation is carried out by specifying the displacement load in several steps.Figure 6a.shows the failure maps highlighting the fracture patterns, when the specified displacement reaches ūy = −2.4× 10 −2 mm (left) and ūy = −2.78× 10 −2 mm (right), respectively.
As compared to Figure 5, the predicted crack pattern is perpendicular to the loading direction, which agrees with the results reported in [64].
To summarize, the estimated results from compression and tension tests, using the present PF method are in close agreement with the reported experimental results [64], illustrating the robustness of the proposed methodology.

Compressive-Traction and Traction-Traction Tests
To complement the previous results, the developed methodology is extended to predict the fracture patterns under compression-traction and traction-traction load cases.The particular load cases are achieved by suitably combining the imposed vertical and horizontal displacements.
Figure 6b,c depict the failure maps under uniaxial tensile load, vertical compression and horizontal traction with ūy / ūx = −1 and vertical compression and horizontal traction cases, respectively, setting ūy = −2.5 ×10 −2 mm.Notable incidence of the external loadings with respect to the estimated crack patterns can be observed.Figure 6b,c correspond to the failure maps of the samples subjected to compression-traction loads with ūy / ūx = −1 and ūy / ūx = −2, respectively.Comparing Figure 6b,c, the applied horizontal displacement is observed to have a moderate influence on the final damage pattern, which is principally aligned with the vertical imposed displacement.
The analysis is also applied to samples subjected to tensile-tensile loads.Figure 7a-c correspond to the failure maps estimated by specifying the displacement ratios ūy / ūx = 1, 2, and 10, respectively.We note that the predicted crack patterns evolve from a perfectly aligned path with the flaw orientation for ūy / ūx = 1, to almost horizontal path for ūy / ūx = 10.In line with previous cases, the estimated results agree with the reported results in [64].

General Aspects
In many engineering applications, rock mechanics in particular, heterogenous media separated by existing interfaces are often observed.The simulation of fracture events in such media is a challenging task.In this context, a plausible alternative is the combined phase field-cohesive zone model (PF-CZM) numerical approach [68], where a consistent coupling between the PF approach for brittle fracture in the matrix and the interface crack model are proposed.This technique can be conceived as a new paradigm for the consistent coupling of damage events of different signature, allowing complex crack evolution with multiple branches.
The point of departure of the PF-CZM consists of the consideration of a system with a predefined crack Γ b and a prescribed interface Γ i , see Figure 8.A generic point of the interface is denoted by the vector x c .Considering the internal potential of the body (omitting the subscript int to alleviate the notation), the governing functional of the system (Equation ( 6)) can be written as: The main attribute of the current method renders the additive decomposition of the dissipative part of the functional G c into two components related to: (i) the fracture energy dissipated in the matrix rock, which is characterized by G b c and is modeled using the PF method, and (ii) the interface energy G i that identifies the crack propagation along the existing interfaces and is modeled using a cohesive zone representation.Correspondingly, the governing functional in Equation ( 21), can be rewritten as: where g is used to identify the displacement gaps between the two flanks of the interface, h is a history parameter [58], and d is the PF variable in the matrix.
In this work, we account for the tension cut-off interface behavior given in [68], whose apparent stiffness is modified by taking into consideration the matrix damage state in the surrounding region via d.Moreover, it is worth mentioning that the interface fracture energy, Equation (22), can be decomposed into the corresponding contributions to the fracture Modes I and II in a 2D setting of analysis, which are represented by G i IC and G i I IC , respectively.According to [68], it is assumed that the critical interface opening (g c ) obeys a linear law with respect to the damage state in the matrix.Therefore, the following relation can be postulated: , where g c,0 = g c (d = 0) and g c,1 = g c (d = 1).Then the traction-separation law at the interface (Figure 9 for mode I), are given by: where σ and τ are the critical traction values associated with the fracture Modes I and II, respectively, g are the displacements whose sub-indices n and t refer to fracture Modes I and II, respectively.Consequently, the apparent interface stiffnesses for both fracture modes as a function of the bulk damage d render where k 0 and g 0 represent for the interface stiffness corresponding to intact matrix rock, i.e., d = 0. We adopt a quadratic-based fracture criterion for triggering propagation events, given by: where G i I and G i I I are the individual energy release rates (ERRs) for fracture Modes I and II, respectively, which take the form: Critical values of the ERRs defined above are compared to their critical values via Equation ( 26), whereby G i IC and G i I IC are the critical values, given by: Finally, the interface formulation to encompass geo-mechanical aspects can be arrived in the similar lines of [75], where we incorporate a modified version of the normal penalty stiffness (to prevent inter-penetration between the interface flanks), which obeys the following hyperbolic expression: where g n is the interface closure, k n,0 is the initial stiffness, g m is the maximum interface closure.The values k n,0 and g m , can be related in rock mechanics to the joint compressive strength (JCS), which is provided in the literature for many different geological masses, using the following correlations (Figure 9): where JCS is the joint wall compression strength, which is the joint roughness coefficient introduced by Barton and Choubey [76], see also a recent overview on the use of this index in [77] and a 0 = JRC 0 /50.These coefficients can be characterized by means of experimental campaigns on specific rock masses.In line with the description of the current model, it is important to note that the independent parameters that characterize the interface response are: (i) the critical normal and tangential displacements corresponding to Mode I and II of fracture g n and g t , (ii) the associated interface fracture energies and the fracture criterion, and (iii) the corresponding geo-mechanical parameters.Additionally, the present model includes an additional capability that couples the apparent interface stiffness with respect to the PF (bulk) damage variable, and therefore two extra parameters are required corresponding to the pairs: g n,0 , g n and g t,0 , g t .

Variational Form and Finite Element Formulation
In this section, we present the particular aspects concerning the interface contribution to the total potential of the system, Equation (22).Following a standard Galerkin procedure, the weak form of the interface contribution can be expressed as: where δu are the test functions of the displacement field V u = δu | u = u on ∂Ω u , u ∈ H 1 , and δd stands for the test functions of the PF variable d The discretization of the present interface formulation is performed using first-order elements with linear interpolation functions.Thus, analogous to the PF approach for bulk fracture presented in Section 2.2, d denotes the nodal displacement vector and d is the phase field nodal vector, both at the element level.Accordingly, the discrete form of Equation ( 32) at the element level The vector of displacement gaps (g) for any arbitrary point on the surface Γ el i is the result of the difference in displacement vector between opposite flanks, which can be obtained from the displacement field d by performing a pre-multiplication by the operator L: where N is a matrix including the element shape functions of the displacement field and Bd = NL stands for the compatibility operator of the displacement field associated with the interface.The transformation of the global gap vector into the local setting referred to the interface [38,41] is conducted via the rotation matrix R from Equation (35) leading to the local gap vector g loc , given by: In a similar manner, the following discretization forms can be deduced for the average PF variable d at the interface, at the element level Γ el i as: where M d is an operator computing the average value of the PF variable between the two interface flanks and Bd = N d M d identifies the corresponding compatibility operator.The specific forms of these matrices are given in [39,41].
Therefore, the discrete weak form renders leading to the residual vectors: Finally, through the directional derivative, the linearization of the residual vectors yields the following tangent matrices: where: and α and β are given by: Please note that in line with Equation ( 19), the resulting coupled equations are solved using the monolithic Newton-Raphson scheme, obeying the form: 5. Fracture Analysis in Heterogeneous Rock Masses Using the PF-CZM Approach

General Considerations
The presence of interfaces in joint propagation is especially relevant in sedimentary rocks composed of siltstone and shale layers.The analysis of cracking events in multilayered materials has been a recurrent research in the last few years, as described in [68].Similar to composite materials, there are three fundamental aspects of crucial importance in fracturing across interfaces, called joints in the rock mechanics literature: (i) the strength of the joint, (ii) the material properties of the layers and (iii) the loading conditions.
As discussed in [78], when an approaching crack impinges on a tough interface this crack generally tends to propagate into the adjacent bulk without experience significant deviations from its initial path.Conversely, weak joints are much more prone to fail, leading to noticeable delamination events.This basic scenario can be reproduced in a straightforward manner by the PF-CZM approach outlined in the previous section.Thus, Figure 10 qualitatively depicts the simulation of a homogeneous system of span L = 50 mm with the mechanical properties given above with a prescribed horizontal interface, which is then subjected to a uniform displacement at the lateral edges ūx .We set the ratio between the fracture toughness of the bulk and the interface as G b c /G i c = 50 for weak joint scenario, whereas, G b c /G i c = 0.05 for tough joints.In this graph, it can be observed that the proposed PF-CZM faithfully reproduces the basic behavior discussed in [78].Specifically, a discontinuous displacement map is shown for the weak joint case illustrating delamination events, whereas a well-defined crack front perpendicular to the loading direction is reproduced for the strong joint response.However, for bi-material systems with weak joints, more complex situations can be encountered according to LEFM studies [68].Particularly, three different propagation scenarios can be found based on the specific Dundurs's parameters of the system when a crack impinges on an interface: (1) single deflection along the interface, (2) double deflection along the interface and (3) penetration.Figure 11 succinctly reproduces the results discussed in [68] using the current PF-CZ approach for a bi-material system subjected to uniform tensile load with a prescribed vertical joint between both media.Observing Figure 11, the numerical predictions led to very satisfactory qualitative agreement with respect to theoretical LEFM results in terms of fracture patterns for different values of the Dundurs's parameter α with respect to Relying on the illustrated potential capabilities, the next Sections addresses the applicability of the PF-CZ method to simulate fracture events in multilayered rock masses.

Application to Rock Salt Fracture with Inclined Interfaces
One of the potential applications of the developed numerical framework is the investigation of the response of bedded rocks salt.Particularly, these rocks perform in a different manner depending upon the corresponding lithology, which can obey either single or multi-layer arrangements.Following [79], mixed lithotypes present horizontal or inclined beds of anhydrite mudstone alternating with halite, leading to anhydrite-halite composite rocks.
From the experimental evidence in [79], it was observed that in such configurations, the presence of interlayers tends to promote crack initiation at this location, this fact being also promoted due to the inclined interface definition.Thus, these authors reported that typical failure sequence can be described by three different phases: (i) near/along the region of the interface, initial micro-defects stemming from the differences between sedimentation and consolidation evolve; (ii) increasing the applied loading, such micro-cracks propagate along the interface and within the rock salt region; (iii) finally, the coalescence of different crack occur, breaking the specimen.
The current numerical method is applied to study the deformation behavior and the fracture characteristics (initiation and coalescence) in bedded rocks salt with inclined interfaces under uniaxial compression in line with [79].For this purpose, uniaxial compression simulations of composite rocks salt with interlayers of 20 mm in thickness are numerically analyzed, which showed a mechanical behavior close to brittle character.Due to the lack of reliable data regarding the geometry of the specimens analyzed in [79], the study is restricted to a qualitative investigation following a central interlayer disposal, see Figure 12.A domain of dimensions 100 × 200 mm is defined with a central interlayer of 20 mm in thickness following a macroscopic representation with two small flaws close to the interface.This domain is dsicretized using 84388 elements combining bulk and interface elements.The mechanical and fracture properties of both materials are listed in Table 1.Please note that for the interface between the rock salt and the anhydrite we choose an intermediate value between the transgranular boundary crack properties reported in [80], so that we set G i c = 1.2715J/m 2 .Figure 12 depicts the qualitative numerical-experimental correlation of the predicted crack path of the current method and that observed in the experiments, showing a good qualitative agreement.In particular, during the simulations, it is observed that the main vertical cracks initiate at the flaws near the interface and propagate through the interlayer and subsequently both cracks propagate along the interface coalescing.The final failure pattern at this location agrees from a qualitative point of view with that reported in [79].However note that the development of wing and shear crack paths are slightly deviated from the experimental evidence, so that alternative PF methods for bulk fracture as proposed in [65,66] can be employed in forthcoming studies.Nevertheless, a more comprehensive investigation falls beyond the scope of the current paper due to the lack of the required material data.

Numerical Investigation: Specimen with Multiple Flaws and Presence of Interfaces
Continuing with the assessment of the predictive capabilities of the current PF-CZM approach for heterogenous rock masses, we adopt the specimen definition described in Section 3, including the prescribed weak joints and multiple flaws with different lengths and orientations forming 45 • with respect to the horizontal direction.Material properties of the specimens also replicate those given in Section 3.For the sake of conciseness in the results, we assume that the different lithologies are from the same material.Also note that bi-and tri-material systems can be directly analyzed using the current numerical methodology in a straightforward manner.
Additionally, with respect to the joint definitions that separate different bulk domains, we consider two cases: (i) a single horizontal joint, (ii) a single vertical joint.Joint mechanical properties are given in Table 2, whereas, contact stiffness properties relying on geo-mechanical parameters are listed in Table 3.
It is worth mentioning that due to the lack of reliable experimental data, after validating the proposed method in the previous Sections, the present results are conducted to analyze the potential fracture patterns for different geometric characteristics and loading conditions.The specific description of the flaw definitions within the domain, the joint positions and the loading conditions are shown in Figure 13a.The position of the joint coincides with the medial horizontal edge of the specimen.A compression-tension load case, where ūy = −0.2ūx , with compression along the vertical direction is considered (we set ūy = 5 × 10 −2 mm in the following analysis here presented) Figure 13b,c depict the predicted contour plots of the vertical displacement and the damage pattern at the end of the simulation.Analyzing the displacement field, Figure 13b, a clear discontinuity at the position of the joint can be appreciated.This fact indicates the occurrence of joint failure and therefore delamination.However, almost no increase in the value of PF variable within the specimen is predicted, see Figure 13c.These results clearly indicate that delamination along the existing joint can be identified as the predominant failure phenomenon for the given loading conditions.To examine the clear influence of the applied loading on the failure characteristics, we analyze the specimen by prescribing the traction-traction loading conditions, i.e., through the definition of ūy = 0.5 ūx (both positive values along the corresponding edges based on the Cartesian setting given in the graph).The computed results clearly pinpoint a substantial change in the fracture conditions within the specimen, see Figure 14.In this case, the development of nucleation and propagation phases around the flaw tips are noticed.Upon load increasing, both crack events emanating from the flaws coalesce and subsequently reach the joint.A third phase in the failure path evolution is manifested by the crack propagations along the existing joint and subsequently progressing into the bottommost substrate, which concomitantly occurs with the crack propagation in the topmost substrate.These phenomena can be identified by analyzing the discontinuities in the contour map associated with the vertical displacement field and the PF evolution.This concomitant simulation of joint and bulk damage can be identified as one of the most appealing attributes of the present numerical methodology, since this allows interaction without user intervention throughout the computations.The simultaneous evolution of failure events from different signature during the analysis can be identified in Figure 14 regarding the load-displacement curve, whereby a first linear elastic evolution is followed by crack growth from the two flaws with incipient delamination along the interface in a later stage.Failure initiation is around u y / ūy =0.53.
Finally, for a most severe traction-traction case we define ūy = ūx .Analyzing the contour map corresponding to the vertical displacement, Figure 15a, it is observed that for this loading configuration the crack pattern tends to progress with gradual change in its path.In particular, bulk fracture is reoriented along the direction of the flaws, and then reduce the extent of the delamination front, as seen in Figure 15b.This effect is also noticed above for the homogeneous crack example with a single inclined flaw, in Section 5.3.This behavior can be also noticed in the load-displacement evolution curve, Figure 15c, which in comparison with the previous case, the failure initiation is predicted to occur at around u y / ūy =0.5.

Specimen with Three Inclined Flaws
The previous analysis can be further extended to specimen definitions with 3 parallel flaws, 2 in the topmost substrate and 1 in the bottommost substrate.Moreover, regarding the joint position, we analyze the cases given by: (1) a single horizontal joint whose position is coincident with that described in the previous sections (Figure 16), and (2) a single vertical joint, placed at a distance of 20 mm with respect to the rightmost vertical edge, see Figure 18a.With reference to the case with the horizontal interface under traction-traction loading conditions given by ūy = ūx , it can be observed that the inclusion of the third flaw in the bottommost substrate of the model notably modifies the damage pattern.It is also noticed that the smallest flaw is not activated by the applied loading, and no damage events occur at the corresponding tips (Figure 16).Thus, crack phenomena are associated with the remaining larger flaws, featuring clear propagation paths which tend to coalesce.However, such two main cracks interact with the existing interface, and correspondingly the final damage pattern is a combination of bulk and interface failure.This sequence of damage can be identified analyzing the load-displacement curve (Figure 17) where different snapshots of the PF variable and the horizontal displacement are given.According to this graph, as usual, an initial linear evolution is depicted up to the fracture propagation of the main two flaws, which coincide with the first peak of the curve.Current results predict that these two main cracks initially propagate without involving any decohesion along the interface (note the continuity of the horizontal displacement field).However, when these two cracks are very close to the interface, they also induce the development of decohesion events, which are responsible for the final rupture of the specimen.This fact is characterized by the final abrupt failure in the load-displacement curve due to the tension cut-off interface model herein envisaged.
Regarding the vertical joint configuration, Figure 18b,c depict the final PF patterns for the cases ūy = −0.1 ūx and ūy = − ūx , respectively, showing completely different failure mechanisms in the predictions.Clear bulk failure emanating from the two larger flaws can be observed for the case ūy = − ūx , whereas the scenario setting ūy = −0.1 ūx shows that the main failure mechanism in this later specimen is due to joint decohesion.These two failure patterns can be confirmed by the analysis of the corresponding load-displacement curves, see Figure 19.Thus, the case ūy = −0.1 ūx features a standard evolution up to the first peak.At this point, the decohesion failure starts in unstable manner, followed by a second linear evolution up to the complete decohesion without any bulk failure, see the snapshot of the vertical displacements in Figure 19b.This response contrasts with that corresponding to the case ūy = − ūx , whereby the failure mechanism is due to the bulk failure without any interaction with the interface decohesion.

Conclusions
In this paper, the application of the combined PF-CZ approach for fracture heterogeneous in rock masses has been outlined with the incorporation of geo-mechanical parameters into the interface definition.The present methodology has been first validated for matrix rock damage, exhibiting an excellent accuracy with respect to experimental data.Subsequently, the developed method has been applied to the simulation of complex damage events in heterogeneous rock masses with weak joints and multiple flaws.
Results from the current model illustrate that the PF-CZ approach permits the prediction of crack initiation, propagation, coalescence, and branching in geo-materials in an autonomous manner without the necessity of user intervention.Moreover, differing from alternative PF formulations, its coupling with interface-like modeling techniques allows the characterization of tortuous crack patterns along predefined weak joints in such materials.
Future directions in this concern will regard the extension of the proposed method for anisotropic rocks, hydraulic fracture, porous media, and modified formulations for triggering shear failure.

Figure 1 .
Figure 1.Regularization of the sharp crack topology within the spirit of the phase field approach for fracture.The PF approach for fracture envisages the definition of a scalar-based phase field variable d, where d : Ω × [0, t] −→ [0, 1] [48], which quantifies the stiffness degradation at the material point level, in order to smear out the sharp crack topology.The intact and fully deteriorated states are identified based on the value of the PF parameter, i.e., d = 0 and d = 1, respectively.Therefore, d adopts continuous values between 0 and 1.The one-dimensional PF approximation is given by the expression: d = e −|x|/l , where l indicates the characteristic material length related to the material strength and governs the regularization width.According to [49], the PF problem is governed by:

Figure 4 .
Figure 4. Specimen with inclined flaws: geometry and loading conditions.

Figure 5 .
Figure 5. Uniaxial compression case: Failure pattern at an initial stage (a) and at advanced stage (b).

Figure 8 .
Figure 8. Coexistence between brittle fracture in the bulk and cohesive debonding of an interface within the context of the phase field approach of fracture.

Figure 9 .
Figure 9. Modification of interface penalty stiffness under normal compression due to geo-mechanical aspects: schematic representation under negative normal gap g n .

Figure 10 .
Figure 10.String system subjected to uniaxial displacement.Predictions using PF-CZM for two types of interfaces between similar materials: weak joints with delamination events along the interface and strong joints with continuous crack path.

Figure 11 .
Figure 11.Geometry and boundary conditions for the bi-material problem.Crack pattern scenarios for brittle interface.The contour plots of the dimensionless vertical displacement field correspond to three different cases labeled A (double deflection), B (single deflection), C (penetration).

Figure 12 .
Figure 12.Salt rock composite specimen under uniaxial compressive loading.Geometry and boundary conditions, experimental failure pattern and numerical predictions.

Figure 15 .
Figure 15.Specimen with 2 inclined flaws and a horizontal joint.Case ūy = ūx traction-traction: (a) vertical displacement; (b) damage pattern; (c) Load-displacement evolution curve including snapshots of the phase field variable.

Figure 16 .
Figure 16.(a) Specimen with 3 inclined flaws and a horizontal joint.(b) Phase field damage pattern at the end of the simulation.

Figure 17 .
Figure 17.Specimen with 3 inclined flaws and a horizontal joint.Load-displacement evolution curve including snapshots of the bulk (phase field) failure and the horizontal displacement fields.

Table 1 .
Mechanical and fracture properties: halite and anhydrite.

Table 3 .
Penalty stiffness properties of the joint based on geo-mechanical aspects.