Crack patterns in heterogenous rocks using a combined phase field-cohesive interface modeling approach: A numerical study

energies Article Crack Patterns in Heterogenous Rocks Using a Combined Phase Field-Cohesive Interface Modeling Approach: A Numerical Study José Reinoso 1,* , Percy Durand 2, Pattabhi Ramaiah Budarapu 3and Marco Paggi 4 1Elasticity and Strength of Materials Group, School of Engineering, Universidad de Sevilla, 41092 Seville, Spain 2Department of Building Structures and Geotechnical Engineering, Universidad de Sevilla, 41012 Seville, Spain; [email protected] 3School of Mechanical Sciences, Indian Institute of Technology, Bhubaneswar 752050, India; [email protected] 4IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy; [email protected] *Correspondence: jr[email protected] Received: 26 November 2018; Accepted: 4 March 2019; Published: 13 March 2019   Abstract: 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. Keywords: rock mechanics; phase field approach to fracture; fracture of geo-materials; cohesive zone model; interface modeling 1. 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 Energies 2019,12, 965; doi:10.3390/en12060965 www.mdpi.com/journal/energies Energies 2019,12, 965 2 of 28 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. 1.1. 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–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 – 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. 1.2. 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. Energies 2019,12, 965 3 of 28 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 – 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 – 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 – 30 ], or the enrichment of the strain field (enhanced FEM, E-FEM) [31–34]. • Adaptive insertion of interface-like cohesive elements within the delimiting edges of the existing mesh, using extrinsic or intrinsic representations of the cohesive law [35–41]. 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 – 46 ], whereby a direct competition between the elastic and the fracture energies is constructed in an integrated global functional of the system. The corresponding Energies 2019,12, 965 4 of 28 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 – 52 ]. Exploiting this potential, posterior studies, the PF approach has been extended to fracture studies in thin-walled structures [ 53 – 55 ], ductile failure [ 56 , 57 ], cohesive fracture [ 58 ], dynamic applications [ 59 , 60 ], multi-physics environments [ 61 – 63 ]. Recently, this methodology has also been applied to study the initiation of fracture in rocks with multiple flaws [64–66]. 1.3. 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 2summarizes 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. 2. 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. 2.1. Basic Concepts Within the context of multi-dimensional analysis, see Figure 1, we consider an arbitrary body Ω∈Rnd in the Euclidean space of dimension nd , whose material points are denoted by the position vectors x in the Cartesian setting. Body actions are identified by the vector fv:Ω−→ Rnd . The delimiting boundary Ω is defined such that ∂Ω∈Rnd−1 , where ∂Ωu and ∂Ωt denote the boundary portions in which kinematic and static surface actions are prescribed, respectively, complying to: ∂Ωt∪∂Ωu=∂Ω and ∂Ωt∩∂Ωu=∅ . The kinematic field is denoted by u , whereas the second-order Cauchy stress tensor is identified by σ. Boundary conditions in this arbitrary body render: u=uon ∂Ωu, and t=σ·non ∂Ωt, (1) where nis the outward normal to the body. Energies 2019,12, 965 5 of 28 1 3 2 1 3 2  1 3  2 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: d−l2∆d=0 in Ωand ∇Xd·n=0 in ∂Ω(2) where ∆d is the Laplacian of d , and ∇xd stands for the spatial gradient. At this point, relevant contributions on this matter are due to the contributions by Wu and coauthors [ 70 – 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 , ∇xd) is defined, which can be estimated using the relation: Γc(d):=ZΩγ(d,∇xd)dΩ. (3) 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]: γ(d,∇xd) = 1 2ld2+l 2|∇xd|2. (4) The above definition of γ(d , ∇xd) allows the approximation of the crack surface integral in the original Griffith formulation by a volume integral as follows: ZΓc Gcd∂Ω≈ZΩGcγ(d,∇xd)dΩ. (5) Therefore, the total potential of the body can be split into the internal and external counterparts represented, Πint(u,d)and Πext(u), respectively: Π(u,d) = Πint(u,d) + Πext(u) = ZΩψ(ε,d)dΩ+ZΩGcγ(d,∇xd)dΩ+Πext(u), (6) where the first term in the integral represents the elastic energy of the body and the second term accounts for the energy dissipation. Energies 2019,12, 965 6 of 28 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]: ψ(ε,d) = g(d)ψe +(ε) + ψe −(ε), (7a) ψe +(ε) = λ 2(htr[ε]i+)2+µtr[ε2 +], (7b) ψe −(ε) = λ 2(htr[ε]i−)2+µtr[ε2 −], (7c) 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 h•i±stands for the so-called Macaulay bracket: h•i±= (•±|•|)/2. Finally, the degradation function g(d) , which introduces a stiffness deterioration with the growth of the defect (crack), can be defined as: g(d) = (1−d)2+K, (8) 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: σ:=∂ˆ ψ ∂ε=g(d)σ++σ−; where σ±=λ(htr[ε]i±)1+2µε±. (9) In the Equation (9) , 1 represents the second-order identity matrix and σ± identify the positive and negative parts of the stress tensor. 2.2. 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 bin the superscript. Based on the isoparametric interpolation, the LaGrangian shape functions at the element level NI(ξ) , 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: x∼ =∑n I=1NIexI=Nex;u∼ =∑n I=1NIdI=Nd;δu∼ =∑n I=1NIδdI=Nδd;∆u∼ =∑n I=1NI∆dI=Nδd,(10) where n indicates the number of nodes at the element level, and xI y dI denotes the nodal coordinates and displacements, respectively, which are collected into the operator ex 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 Bdassociated with the kinematic field: ε∼ =Bdd;δε∼ =Bdδd;∆ε∼ =Bd∆d. (11) Similarly, the interpolation of the PF variable (damage-like variable), its variation and linearization render: d∼ = n ∑ I=1 NIdI=Nd;δd∼ = n ∑ I=1 NIδdI=Nδd;∆d∼ = n ∑ I=1 NI∆dI=Nδd, (12) Energies 2019,12, 965 7 of 28 where dI are the nodal values of the PF variable, which are arranged into the vector d . The spatial gradient of the PF ( ∇xd ), its variation ( ∇xδd ) and its linearization ( ∇x∆d ) can be interpolated through the operator Bd: ∇xd∼ =Bdd;∇x(δd)∼ =Bdδd;∇x(∆d)∼ =Bd∆d. (13) Therefore, the discrete version of Equation (6) at the element level, identified by the super-index el, can be expressed as: δ˜ Πel b(d,δd,d,δd)(14) =δdTZΩel h(1−d)2+KBT dσ++BT dσ−idΩ−Z∂Ωel NTtd∂Ω−ZΩel NTfvdΩ +δdTZΩel −2(1−d)NTψe +(ε)dΩ+ZΩel Gb clBT d∇xd+1 l2NTddΩ =δdTfb d+δdTfb d(15) where the residual vector fb d=fb d,ext −fb d,int , in which, fb d,int and fb d,ext represent the internal and external forces, respectively, given by: fb d,int =ZΩel h(1−d)2+KBT dσ++BT dσ−idΩ, (16) fb d,ext =Z∂Ωel NTtd∂Ω+ZΩNTfvdΩ. (17) Correspondingly, fb d is the residual vector associated with the PF variable, which can be estimated as: fb d=ZΩel −2(1−d)NTψe +(ε)dΩ+ZΩel Gb clBT d∇xd1 l2NTddΩ. (18) 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]: "Kb dd 0 0 Kb dd#"∆d ∆d#="fb d,ext 0#−"fb d,int fb d#. (19) The particular forms of the matrices Kb dd , Kb dd , Kb dd and Kb 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. 3. Bulk Failure Predictions 3.1. 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.5 GPa, Poisson’s ratio ν = 0.25 and Gb c= 100 J/m2. Furthermore, the length scale in the PF approach can be related to the material strength σcof one-dimensional bar: l=27EGb c 256σ2 c . (20) Energies 2019,12, 965 8 of 28 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 3shows 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. Figure 2. (a) Brazilian splitting test: loading condition. (b) Damage pattern of the specimen. Figure 3. Brazilian splitting test: Load-displacement evolution curve. 3.2. 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, Gc= 3 J/m2 and l= 0.1 mm [ 64 ]. Energies 2019,12, 965 9 of 28 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 ¯ uy=± 3 × 10 −2 mm (tensile/compression), whereas ¯ ux can have different ratios with respect to ¯ uy 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. Figure 4. Specimen with inclined flaws: geometry and loading conditions. 3.2.1. Uniaxial Compression Test We first analyze the deformation patterns in a uniaxial compression test, which is achieved by specifying ¯ uy=− 3 × 10 −2 mm, in several steps and ¯ ux= 0, see Figure 4. Figure 5shows the failure maps, highlighting the fracture patterns, when the specified displacement reaches ¯ uy=−2.4 ×10−2mm (left) and ¯ uy=− 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]. Figure 5. Uniaxial compression case: Failure pattern at an initial stage (a) and at advanced stage (b). 3.2.2. 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 ¯ uy= 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 ¯ uy=− 2.4 × 10 −2 mm (left) and ¯ uy=− 2.78 × 10 −2 mm (right), respectively. Energies 2019,12, 965 16 of 28 Please note that in line with Equation (19) , the resulting coupled equations are solved using the monolithic Newton-Raphson scheme, obeying the form: "Ki dd Ki dd Ki ddKi dd #"∆d ∆d#="fi d fi d#. (43) 5. Fracture Analysis in Heterogeneous Rock Masses Using the PF-CZM Approach 5.1. 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 ¯ ux . We set the ratio between the fracture toughness of the bulk and the interface as Gb c/Gi c = 50 for weak joint scenario, whereas, Gb c/Gi 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. 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. Energies 2019,12, 965 17 of 28 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 Π1=Gb c/Gi c. 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). 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. 5.2. 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 Energies 2019,12, 965 18 of 28 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 Gi c=1.2715 J/m2. Figure 12. Salt rock composite specimen under uniaxial compressive loading. Geometry and boundary conditions, experimental failure pattern and numerical predictions. Table 1. Mechanical and fracture properties: halite and anhydrite. Material E (GPa) νGb c(Jm−2) Halite 36.87 0.254 1.155 Anhydrite 74.36 0.269 1.805 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. Energies 2019,12, 965 19 of 28 5.3. 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. Table 2. Interface fracture properties. Critical traction Mode I interface σc075 MPa Critical traction Mode II interface τc090 MPa Fracture toughness Mode I 0.002 N/mm Fracture toughness Mode II 0.008 N/mm Table 3. Penalty stiffness properties of the joint based on geo-mechanical aspects. gm10 gn0 a00.194 mm JCS0169 MPa JRC09.7 Friction angle 32o L00.2 m 5.3.1. Specimen with Two Inclined Flaws and a Horizontal Joint 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 ¯ uy=− 0.2 ¯ ux , with compression along the vertical direction is considered (we set ¯ uy= 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. Energies 2019,12, 965 20 of 28 Figure 13. ( a ) Specimen with 2 inclined flaws and a horizontal joint. Case ¯ uy=− 0.2 ¯ ux : ( b ) horizontal displacement field; (c) damage pattern. 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 ¯ uy= 0.5 ¯ ux (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 uy/¯ uy=0.53. Finally, for a most severe traction-traction case we define ¯ uy=¯ ux . 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 uy/¯ uy=0.5. Energies 2019,12, 965 21 of 28 Discontinuity in the phase field variable Figure 14. Specimen with 2 inclined flaws and a horizontal joint. Case ¯ uy= 0.5 ¯ ux traction-traction: ( a ) vertical displacement; ( b ) damage pattern; ( c ) Load-displacement evolution curve including snapshots of the phase field variable. Discontinuity in the phase field variable Figure 15. Specimen with 2 inclined flaws and a horizontal joint. Case ¯ uy=¯ ux traction-traction: ( a ) vertical displacement; ( b ) damage pattern; ( c ) Load-displacement evolution curve including snapshots of the phase field variable. Energies 2019,12, 965 22 of 28 5.3.2. 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. Figure 16. ( a ) Specimen with 3 inclined flaws and a horizontal joint. ( b ) Phase field damage pattern at the end of the simulation. With reference to the case with the horizontal interface under traction-traction loading conditions given by ¯ uy=¯ ux , 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 ¯ uy=− 0.1 ¯ ux and ¯ uy=−¯ ux , respectively, showing completely different failure mechanisms in the predictions. Clear bulk failure emanating from the two larger flaws can be observed for the case ¯ uy=−¯ ux , whereas the scenario setting ¯ uy=−0.1 ¯ ux 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 ¯ uy=−0.1 ¯ ux 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 ¯ uy=−¯ ux , whereby the failure mechanism is due to the bulk failure without any interaction with the interface decohesion. Energies 2019,12, 965 23 of 28 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. Figure 18. Specimen with 3 inclined flaws and a vertical interface: ( a ) Loading and geometric description. (b) Phase field pattern for ¯ uy=−0.1 ¯ ux. (c) Phase field pattern for ¯ uy=−¯ ux. Energies 2019,12, 965 24 of 28 Figure 19. Specimen with 3 inclined flaws and a vertical interface: ( a ) Load-displacement curve for ¯ uy=−0.1 ¯ ux. (b) Load-displacement curve for ¯ uy=−¯ ux. 6. 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. Author Contributions: Conceptualization, J.R., P.D.; Methodology, J.R., M.P.; Investigation, P.R. Budarapu, J.R.; Resources, J.R.; Writing-Original Draft Preparation, J.R., P.R.B., M.P.; Writing-Review and Editing, J.R., P.R.B., M.P., P.D.; Funding Acquisition, J.R. Funding: J. Reinoso is grateful to the support of the Spanish Ministry of Economy and Competitiveness/FEDER (Project MAT2015-71036-P). Conflicts of Interest: The authors declare no conflict of interest. References 1. Hoek, E.; Bieniawski, Z.T. Brittle Fracture Propagation in Rock under Compression. Int. J. Fract. 1984 ,26, 276–294. [CrossRef] 2. Hoek, E. Brittle failure of rocks in rock mechanics in engineering practice. In Rock Mechanics in Engineering Practice; Stagg, K.G., Zienkiewicz, O.C., Eds.; Wiley: New York, NY, USA, 1968; pp. 99–124. 3. Lajtai, E.Z. Brittle Fracture in Compression. Int. J. Fract. 1974,10, 525–536. [CrossRef] Energies 2019,12, 965 25 of 28 4. Ingraffea, A.R.; Heuze, F.E. An Analysis of Discrete Fracture Propagation in Rock Loaded in Compression. In Proceedings of the 18th U.S. Symposium on Rock Mechanics, Keystone, Colorado, CO, USA, 22–24 June 1977. 5. Suits, L.D.; Sheahan, T.C.; Wong, L.N.Y.; Einstein, H.H. Using high speed video imaging in the study of cracking processes in rock. Geotech. Test. J. 2009,32, 164–180. [CrossRef] 6. Bobet, A. Numerical simulation of initiation of tensile and shear cracks. In Proceedings of the DC Rocks (2001) the 38th US Symposium on Rock Mechanics (USRMS), Washington, DC, USA, 7–10 July 2011. 7. Yang, S.Q.; Yang, D.S.; Jing, H.W.; Li, Y.H.; Wang, S.Y. An experimental study of the fracture coalescence behaviour of brittle sandstone specimens containing three fissures. Rock Mech. Rock Eng. 2012 ,45, 563–582. [CrossRef] 8. Reyes, O.; Einstein, H.H. Failure mechanisms of fractured rock—A fracture coalescence model. In Proceedings of the 7th International Conference on Rock Mechanics, International Society for Rock Mechanics, Aachen, Germany, 16–20 September 1991. 9. Wong, L.; Einstein, H. Crack coalescence in molded gypsum and carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech. Rock Eng. 2009,42, 475–511. [CrossRef] 10. Wong, R.; Chau, K.; Tang, C.; Lin, P. Analysis of crack coalescence in rock-like materials containing three flaws—Part I: Experimental approach. Int. J. Rock Mech. 2001,38, 909–924. [CrossRef] 11. Sagong, M.; Bobet, A. Coalescence of multiple flaws in a rock-model material in uniaxial compression. Int. J. Rock Mech. Min. Sci. 2002,39, 229–241. [CrossRef] 12. Ma, G.W.; Dong, Q.Q.; Fan, L.F.; .Gao, J. An investigation of non-straight fissures cracking under uniaxial compression. Eng. Fracture Mechanics 2018,15, 300–310. [CrossRef] 13. Rabczuk, T.; Belytschko, T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. Int. J. Numer. Meth. Eng. 2004,61, 2316–2343. [CrossRef] 14. Rabczuk, T.; Belytschko, T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput. Methods Appl. Mech. Eng. 2007,196, 2777–2799. [CrossRef] 15. Rabczuk, T.; Zi, G.; Bordas, S.; Nguyen-Xuan, H. A simple and robust three-dimensional cracking-particle method without enrichment. Comput. Methods Appl. Mech. 2010,199, 2437–2455. [CrossRef] 16. Rabczuk, T.; Ren, H. A peridynamics formulation for quasi-static fracture and contact in rock. Eng. Geol. 2017,225, 42–48. [CrossRef] 17. Ren, H.; Zhuang, X.; Cai, Y.; Rabczuk, T. Dual-horizon peridynamics. Int. J. Numer. Meth. Eng. 2016 ,108, 1451–1476. [CrossRef] 18. Ren, H.; Zhuang, X.; Rabczuk, T. Dual-horizon peridynamics: A stable solution to varying horizons. Comput. Methods Appl. Mech. Eng. 2017,318, 762–782. [CrossRef] 19. Bobet, A.; Einstein, H.H. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int. J. Rock Mech. Min. Sci. 1998,35, 863–888. [CrossRef] 20. Gonçalves da Silva, B.; Einstein, H.H. Modeling of crack initiation, propagation and coalescence in rocks. Int. J. Fract. 2013,182, 167–186. [CrossRef] 21. Lemaitre, J.; Desmorat, R. Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures; Springer-Verlag: Berlin, Germany, 2005. 22. Bažant, Z.P.; Pijaudier-Cabot, T.G.P. Nonlocal continuum damage, localization instability and convergence. J. Appl. Mech. 1988,55, 287–293. [CrossRef] 23. Jirásek, M. Nonlocal models for damage and fracture: Comparison of approaches. Int. J. Solids Struct. 1988,35, 4133–4145. [CrossRef] 24. Forest, S. Micromorphic approach for gradient elasticity, viscoplasticity, and damage. J. Engnr. Mech. 2009 ,35, 117–131. [CrossRef] 25. Dimitrijevic, B.J.; Hackl, K. A regularization framework for damage-plasticity models via gradient enhancement of the free energy. Int. J. Numer. Methods Biomed. Energy 2011,27, 1199–1210. [CrossRef] 26. Peerlings, R.; Geers, M.; de Borst, R.; Brekelmans, W. A critical comparison of non local and gradient-enhanced softening continua. Int. J. Solids Struct. 2001,38, 7723–7746. [CrossRef] 27. Areias, P.; Msekh, M.A.; Rabczuk, T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Eng. Fract. Mech. 2016,158, 116–143. [CrossRef] 28. Moës, N.; Dolbow, J.; Belytschko, T. A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 1999,46, 131–150. [CrossRef]