Hydromechanical modeling of unrestricted crack propagation in fractured formations using intrinsic cohesive zone model

Journal Pre-proofs Hydromechanical modeling of unrestricted crack propagation in fractured formations using intrinsic cohesive zone model Julio Rueda, Cristian Mejia, Deane Roehl PII: DOI: Reference:

S0013-7944(19)30997-X https://doi.org/10.1016/j.engfracmech.2019.106655 EFM 106655

To appear in:

Engineering Fracture Mechanics

Received Date: Accepted Date:

6 August 2019 31 August 2019

Please cite this article as: Rueda, J., Mejia, C., Roehl, D., Hydromechanical modeling of unrestricted crack propagation in fractured formations using intrinsic cohesive zone model, Engineering Fracture Mechanics (2019), doi: https://doi.org/10.1016/j.engfracmech.2019.106655

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

1

Hydromechanical modeling of unrestricted crack propagation in fractured formations using intrinsic cohesive zone model Julio Ruedaa,b, Cristian Mejiaa, Deane Roehla,b aTecgraf

Institute, Pontifical Catholic University of Rio de Janeiro, 22451-000, RJ, Brazil, www.tecgraf.puc-rio.br bDepartment of Civil and Environmental Engineering, Pontifical Catholic University of Rio de Janeiro, Rua Marquês de São Vicente, 225, Gávea, 22453-900, RJ, Brasil, www.civ.puc-rio.br

Abstract Production from unconventional reservoirs has revolutionized the oil and gas industry strongly through the development of innovative applications, such as horizontal drilling and shale stimulation by hydraulic fracturing. However, the presence of geological discontinuities such as natural fractures (NFs) complicates the hydraulic fracturing treatment and affects the general geometry of the hydraulic fracture (HF). For that reason, it is essential to develop and implement robust numerical methods capable of simulating the phenomena present during hydraulic fracturing in a naturally fractured formation. This paper presents a novel mesh fragmentation technique to simulate unrestricted hydraulic fracture propagation in fractured media. This method is based on the hydromechanical zero thickness interface element combined with the cohesive zone model (CZM). An algorithm for inserting the especial triple-noded interface elements into conventional 2D finite element meshes is described. We introduce a simple but effective topological data structure that considers each finite element independent of others. The proposed topological structure ensures the generation of new nodes and faces which are necessary for the insertion of interface elements. The implementation is validated by simulating some cases of fracturing experiments in laboratory. The numerical procedure provides very good agreement with the laboratory tests and shows the robustness of the proposed computational approach. Some advantages and limitations of the proposed methodology are discussed. 1. Introduction Discontinuities in the rocks present complex behavior and they are a matter of fundamental importance in rock engineering. Several numerical approaches are available to represent the behavior of discontinuities in fractured formations such as the equivalent continuum [1–3], dual porosity/dual permeability [4–7] and interface elements [8–10], among others. The interface element finds a broad range of engineering applications where the crack paths are known as delamination in composite materials [11,12], fault activation [13,14] and grain boundary cracks in brittle materials[15]. Extrinsic or Intrinsic interface element fragmentation is necessary to simulate complex crack patterns during crack propagation. Crack branching [16], mesoscale and microscale arbitrary crack propagation in fiberreinforced material[17] or heterogeneous material characterized by inclusions [18– 20] has been simulated with those approaches. In the extrinsic approach, based on

2

local criterion of stress or strain, an adaptive insertion of cohesive interface elements in the mesh (in space and time) is required. In the intrinsic approach, cohesive interface elements are inserted in the discretization of the structure before the simulation starts. Both approaches have their own advantages and drawbacks. In the extrinsic fragmentation process, new elements are inserted at the edges of regular elements during the simulation. Therefore, this process requires an elaborated data structure to update the modified mesh once new nodes and interface elements are inserted. Pandolfi [21,22] presented a data structure for an extrinsic fragmentation limited to triangle and tetrahedral finite elements. Paulino [23] proposed a more generalized extrinsic fragmentation using several element types. Computational parallelization becomes challenging due to the complexity associated with the mesh topological changes during the simulation process. In order to overcome this problem, Dooley [24] based their implementation on inserting the nodes of the interface elements from the beginning and activating the interface element when the fracture occurs. Alhadeff [25] presented the first proposal of a complete adaptive finite element analysis running on the GPU. Their data structure is specialized in triangular elements. In general, the discrete discontinuity is investigated using the cohesive zone model (CZM) that defines the constitutive response directly in terms of traction versus separation. This model avoids the singularity of the stress field at crack tip, which is present in the classical linear elastic fracture mechanics. The extrinsic cohesive zone model considers that the separation only occurs when interfacial traction meets certain damage initiation criteria 𝑇(𝛿) in terms of stress or strain. Subsequently, the material stiffness is progressively degraded, and finally is null at the critical separation δ𝑐. This behavior is shown in Fig. 1a. On the other hand, the intrinsic cohesive model considers an initial interfacial traction-separation until it reaches the maximal interfacial strength 𝑇(𝛿). Consequently, the material stiffness is progressively degraded, and finally vanishes at the critical separation δ𝑐, as illustrated in Fig. 1b. In general, intrinsic models consider that all cohesive elements are embedded in the mesh since the beginning of the simulation. Therefore, the topological mesh connectivity remains constant during the simulation process. Intrinsic models are easier to implement in terms of mesh and computational parallelization [26]. However, they introduce an artificial compliance that depends on element surface and cohesive parameters related to bulk properties. The phenomenon “artificial compliance” can be reduced increasing the initial stiffness of the cohesive elements. Extrinsic models avoid this artificial effect [27]. Nevertheless, this issue also arises in the extrinsic approach if the inserted fracture closes under compressive load.

Fig. 1. Comparison of two cohesive zone models (CZM): (a) extrinsic model; (b) intrinsic model.

3

It is well known that the presence of discontinuities plays an important role in the hydraulic behavior of the geological formation. The propagation of hydraulic fractures is a relevant problem in rock mechanics once they may alter the hydromechanical properties of the rock. The coupled HM double-noded interface element has been formulated to simulate pre-existing discontinuities [28–30]. However, these formulations do not consider normal flow. Baldoni [31] and Guiducci [32] proposed the hydromechanical triple-noded zero thickness interface element (three coincident nodes) to simulate transversal permeability. Fig. 2 shows a schematic representation of the triple-noded interface element. The corner nodes have displacement and pore-pressure degrees of freedom. The middle node has a pore-pressure degree of freedom to represent the fluid pressure inside the discontinuity [33,34]. This special element has been combined with a cohesive zone model to simulate planar hydraulic fracture propagation [35–37]. In these simulations, hydraulic fractures are modeled by non-intersecting pre-defined cohesive paths.

Fig. 2. Schematic representation of the triple-noded interface element

Recently, the triple-noded zero thickness interface elements have been used to simulate the interaction between hydraulic and natural fractures [34,38–43]. In these cases, the intersected elements needed a special treatment to ensure the continuity in pressure along the hydraulic fracture and natural fractures. Gonzalez [39] defined four triangular elements in the intersection. Haddad et al. [34] and Taleghani et al. [42] applied constraints through additional governing equations to couple the middle-edge pore-pressure degrees of freedom at all inner nodes at the intersection. Guo [38] and Rueda [41] merge the inner nodes on the surfaces of HF and NFs to ensure that the injected fluid flows normal to or along the fracture paths. The implementation of this last configuration, “merge the inner nodes”, seems easier than those proposed by Gonzalez [39], Haddad [34] and Taleghani [42]. Regardless of the choice of those configurations, the hydromechanical CZM has the inability to predict the path and propagation direction of the hydraulic fracture. Accordingly, other numerical techniques have been presented as viable alternatives to simulate nonplanar hydraulic fracture propagation in fractured media. The extended finite element method (XFEM) has been gaining popularity because this method can simulate unrestricted hydraulic fracture propagation [44–48]. However, this method is complex in the treatment of fracture interaction even in 2D cases and its implementation in existing finite element codes is not straightforward. Additionally, parallel computation is not trivial. An intrinsic cohesive approach may be a solution to overcome these limitations. Although the implementation of intrinsic cohesive elements is straightforward in the conventional finite element solvers, the preprocessing step to insert a cohesive interface element at every edge of general finite

4

elements is not trivial. Recently, Nguyen [49] and Truster [50,51] have developed algorithms and open source codes to insert zero thickness interface elements throughout a finite element mesh. However, theses codes have been applied only to mechanical problems using the conventional double-noded interface element. Actually, there is a total lack of discussion on implementation aspects of the insertion of the special hydromechanical triple-noded interface elements into an existing finite element mesh. Also, to our knowledge, no commercial or open source codes have this possibility. In order to fill this gap, the current work presents an innovative methodology of mesh fragmentation to insert the especial triple-noded interface elements into conforming meshes. Its development is motivated by the need of a robust numerical approach to simulate unrestricted hydraulic crack propagation in naturally fractured media. The proposed methodology is independent of the mesh generator and of the finite element solver. Three fracturing experiments in laboratory are simulated to validate its implementation. The first example, studies the influence of natural fracture on the trajectory of the hydraulic fracture propagation. Mesh dependency and computational time of the presented methodology are also evaluated. The second example demonstrates the unrestricted fracture propagation under normal stress regime. In this case, hydraulic fracture propagation is initiated from an unfavorable perforation angle. Finally, the third example shows the mixed effects of the stress contrast and the intersecting angle on the hydraulic and natural fracture interaction. The numerical procedure provides very good agreement with the laboratory tests and shows the robustness of the proposed computational algorithm. Some advantages and limitations of the proposed methodology are discussed. 2. Pore pressure cohesive zone model The governing equations of the fracture involve the interaction between mechanical deformation and fluid flow along and across the interface. The mass conservation equation combines fluid flow along and across the discontinuity according to: ∂𝛿𝑛 ∂𝑡

+

𝑑𝑞𝑙 𝑑𝑠

+ 𝑣𝑡 + 𝑣𝑏 = 0

(1)

where 𝑞𝑙[m2/s] is the longitudinal fluid flow inside the fracture; 𝑣𝑡, 𝑣𝑏 [m/s] are the transversal fluid flow velocities leaking through the top and bottom fracture surfaces, respectively, into the surrounding formation, and δ𝑛[m] is the fracture opening. The fracturing fluid is considered to be incompressible Newtonian. The longitudinal flow 𝑞𝑙 within the fracture is governed by Reynold’s lubrication theory of flow through narrow smooth parallel plates (i.e., Poiseuille flow). The normal flow represents leakoff from fractures into the surrounding porous formation [52]. The cohesive zone model (CZM) has been combined with the triple-noded interface element to simulate hydraulic fracture propagation [35–37]. This model includes a fracture process zone ahead of a fracture tip (Fig. 3a) through a damage model.

5

Material crack tip Mathematical crack tip

𝑣� 𝑞� =

�� � �� ��� ��

Hydraulic pressure 𝑣�

Exponential softening Linear softening

Traction

T(δ)

Crack opening δ�

Point of damage initiation

𝑇(𝛿)

𝐺� � δ�

�

cohesive crack tip

Fracture process zone Fluid-filled fracture (Broken cohesive zone) (unbroken cohesive zone) (a)

𝛿��

Separation

𝛿 � 𝛿�

(b)

Fig. 3. Schematic representation of the cohesive zone in a hydraulic fracture (a) and traction separation response for the cohesive elements (b).

The cohesive material response is specified directly in terms of traction versus separation law (TSL) on the crack surface. Therefore, rock behavior during hydraulic fracturing is represented by two laws: a strain-stress relationship to describe the elastic response of the rock material and the cohesive law that controls fracture initiation and fracture propagation. The process of degradation begins when a function in terms of stress (maximum nominal or quadratic stress ratio) or strain (maximum nominal or quadratic strain ratio) meets its critical values [53]. After that, damage evolution characterizes the progressive degradation of the material stiffness [54] according to linear or exponential softening (Fig. 3b). 3. Mesh Fragmentation: algorithm _MeshFrag There is a total lack of discussion on implementation aspects of the insertion of the special hydromechanical triple-noded interface elements into an existing finite element mesh. Also, to our knowledge, no commercial or open source codes have this capability. Hence, an intrinsic mesh fragmentation program, MeshFrag, is presented in order to simulate unrestricted hydraulic crack propagation. The algorithm is written in Matlab and it supports triangle and quadrilateral elements for 2D analyses. The workflow follows the conventional procedure with the three main phases: model generation, solution, and result interpretation. The insertion of the triple-noded interface elements is based on the fragmentation of conforming meshes and happens in the pre-processing phase, as shown in Fig. 4. For this purpose, node duplication, update of the topological connectivity and generation of extra fracture faces, are required. Subsequently, the conventional doubled-noded interface element is inserted at every edge of the continuum finite elements. Finally, we add the middle-edge nodes of the triple-noded interface elements, which require special treatment in order to ensure continuity of the fluid field. Fig. 5 outlines the fragmentation steps of a conventional finite element mesh.

6

FEM solving

Preprocessing

Mesh generation: Conventional FEM mesh file.msh

Stress initialization: Equilibrate initial stresses, loads and boundary conditions

Mesh fragmentation: Insert interface elements in the file.msh and generate input.inp file

Fracturing model: Hydraulic fracturing analysis and generation of the output file.pos

Postprocessing

FEM Viewer: Results interpretation

Fig. 4. Workflow to simulate unrestricted fracture propagation in fractured media.

(b)

(a)

(c)

Degrees of Freedom Displacement and Pore pressure Pore pressure

(d) Fig. 5. Mesh fragmentation scheme: (a) conventional mesh; (b) duplicate nodes and insert interface elements at the edges of the continuum elements; (c) insertion of intermediate pore pressure nodes; (d) detail of the intersection of several hydromechanical cohesive interface elements.

This methodology is independent of the mesh generation software. Here Gmsh preprocessor [55] is used for practical examples. We make use of the unstructured Delaunay algorithm to generate the final mesh in order to introduce some generality to the description of a heterogeneous structure, crack propagation and fracture paths. This preprocessor combines geometrical entities into more meaningful groups to define domains, boundary conditions, and material properties. By default, in the MSH file format, if physical groups are defined, the output mesh only contains

7

those elements that belong to at least one physical group. Thus, some rules are defined to regularize the process. The topological entities created by Gmsh are identified by italic typeset. Each physical group (attached to elements) is considered to have a different material. Therefore, a physical group is composed of a prefix which defines some procedures and is followed by the material name. For example, it is usually convenient to keep the original nodes in some specific regions denoted by the prefix gce_. On the other hand, we define the prefix intraface_ to fragment a region of continuum elements. Interfaces, joints or predefined fractures are attached to line elements and referenced by the prefix coh_. Finally, we use the prefix bc_ to reference boundary conditions on the nodes. Fig. 6 shows a scheme of a total and a partial fragmentation according to the adopted prefixes.

Fig. 6. Defining element groups - (a) total fragmentation (b) partial fragmentation.

As described above, the first stage is to read the MSH file created by the preprocessor program and recognize the physical groups. The second stage consists of the fragmentation process which is followed by the creation of the input file for the FEM solver. 3.1

Mesh fragmentation for modeling of hydraulic crack propagation

The insertion of triple-noded interface elements in a conforming finite element mesh requires topological changes. In this work, we propose the following procedure. 1. Duplicate nodes: element groups denoted as intraface_ are fragmented. In this case, an original node shared by n elements is duplicated by n nodes. Node duplication is the outcome of mesh fragmentation. 2. Update topological connectivity: a simple topological data structure is employed in the mesh fragmentation process. The proposed data structure is used for triangle and quadrilateral elements with 3 and 4 nodes, respectively. We create new tables when fragmentation occurs. A new table of nodes stores the new node indices (Id), the original indices (OId) of the sequential element connectivity and their corresponding x and y coordinates. A second new element table saves the material id and the new nodal connectivity. Fig. 7 shows an example of data structure used in the mesh fragmentation.

8

Original node table Id

x

y

1

0

0

2

0

1

3

2

2

4

0

2

5

0

1

6

1

1

New node table

3

4

2

6

5

1

1

2

7

8

Id

OId

x

y

1

1

0

0

2

2

0

1

3

6

1

1

4

5

0

1

5

5

0

1

6

6

1

1

7

3

2

2

8

4

0

2

2

Face 2

5 4

6 3

Face 1

1

1

2

New element table

Original element table Id material

v1

v2

v3

v4

Id material

v1

v2

v3

v4

1

1

2

6

5

1

1

2

3

4

1

5

6

3

4

1

5

6

7

8

(a)

(b)

Fig. 7. Updating topological connectivity – (a) original; (b) new

In this data structure, each element is assumed to be independent of the others. Therefore, nodes and elements can be updated based on their own information. This procedure ensures the generation of new nodes and faces which are necessary for the insertion of interface elements. The original node indices (OId) are essential to identify faces shared by two elements (Step 4). For example, face 1 formed by Id nodes 4 and 3 and face 2 formed by Id nodes 5 and 6 share the same OId nodes 5 and 6 (Fig. 7b). 3. Generate extra fracture faces: in order to create the interface elements, a face table saves the material id of its corresponding element and the updated nodal incidence of every face. 4. Create the double-noded interface element: the conventional double-noded interface element is created when two faces share the original nodal incidence. 5. Check predefined interface material: interfaces, joints or predefined fractures are attached to line elements and saved in the table of lines. This table stores the material id and the original nodal incidence of the line element. The material id of the created interface element is updated if it shares the same original nodal incidence with some predefined line element. 6. Create and insert the middle nodes: the middle nodes of the interface elements need a special treatment in order to ensure the continuity of the fluid field. The insertion of the middle nodes of each interface element is based on the verification of the coordinates and nodal incidence of the created middle nodes. Therefore, a new middle node is created provided that its coordinates and nodal

9

incidence are different from those of the previously created middle nodes. Finally, the middle nodes are inserted in the interface elements. Fig. 8 shows the nodal incidence changes to insert the triple-noded interface element in a conventional mesh (Steps 1, 2, 4, and 6). 7 6

4

6

6 21 14

13

4

2

4

2 9 8 1

1

7

4

2 7

13

14

3

8

5

15

16

4

3

9

12

5

2

(a)

1

2

10

(b)

11

1

16 20 12

9 3

1

3

1

3

15

5 18 3

8 17 4

2 19 10

(c)

Fig. 8. Updating nodal incidence to insert the triple-noded interface element: (a) original FEM mesh; (b) mesh fragmentation and insertion of double –noded interface element; (c) insertion of middle nodes.

7. Create the input file for the FEM solver: the final step is to write an input file for the processing in a finite element solver. Thus, user functions can be implemented for this purpose. For such, the triple-noded interface element must be available in the software library. In this work, we have used the commercial software Abaqus®/Standard (2017) for the practical examples. In this software, the triple-noded displacement and pore pressure cohesive interface element (COH2D4P) is used to simulate the fractures while the 3-node displacement linear triangle (CPE3) or 4-node bilinear quadrilateral displacement and pore pressure (CPE4P) plane strain elements are used to model the surrounding medium. More details of cohesive-element size, meshing scheme, and far-field boundary conditions can be found in [56]. 3.2

Stress initialization

In naturally fractured media, the in-situ stresses play an important role in hydraulic fracture propagation. Abaqus/Standard provides a geostatic procedure as the first step of geotechnical problems. In this work, stress initialization corresponding to the geostatic step occurs in two stages. The first stage calculates the loads to equilibrate the in-situ stress field. Therefore, the displacements of all nodes are restricted in the model. Subsequently, the initial normal stress 𝜎𝑛 and shear stress 𝜏 in all cohesive interface elements are defined in the Abaqus user subroutine sigini.for [57] according to the following expressiona: 1

𝜏 = 2(𝜎ℎ,𝑚𝑖𝑛,𝑒𝑓𝑓 ― 𝜎𝐻,𝑚𝑎𝑥,𝑒𝑓𝑓)sin 2𝜃

(2)

𝜎𝑛 = 𝜎𝐻,𝑚𝑎𝑥,𝑒𝑓𝑓cos2 𝜃 + 𝜎ℎ,𝑚𝑖𝑛,𝑒𝑓𝑓sin2 𝜃

(3)

where 𝜃 is the inclination angle of the interface element; 𝜎𝐻,𝑚𝑎𝑥,𝑒𝑓𝑓 and 𝜎ℎ,𝑚𝑖𝑛,𝑒𝑓𝑓 are the effective maximum and minimum horizontal stresses attributed to solid

10

elements. Finally, the nodal reactions in the system can be obtained. The second stage allows to verify that an initial geostatic stress field (assigned to the interface and solid elements) is in equilibrium with the applied boundary conditions and loads (nodal reactions obtained in the first stage). The hydromechanical analysis is started when the in-situ stress field is equilibrated and produce zero deformation. More details about this procedures can be found in [57]. 4

Effect of the interface element stiffness on the elastic response of the material

The intrinsic cohesive model introduces an artificial compliance that alters the elastic response of the material prior to the onset of a hydraulic fracture. This phenomenon is called “artificial compliance” and can be reduced by increasing the initial stiffness of the cohesive elements. The normal stiffness (𝐾𝑛) and shear stiffness (𝐾𝑠) can be defined in terms of the elastic properties of the surrounding rock material, as shown in Eq (4): Elastic stiffness

{

𝐾𝑛 = 𝐾𝑠 =

𝛼 𝑡 𝛼 𝑡

𝐸

(4)

𝐺

where 𝐸 and 𝐺 are the Young’s and the shear modulus of the adjacent material; 𝑡 is an artificial mechanical thickness; and 𝛼 is a hardening factor much larger than 1 (𝛼 ≫ 1). This factor should be high enough to provide a reasonable stiffness but small enough to avoid numerical problems such as spurious oscillations of the interface element [12,58,59]. Values of 𝛼 𝑡 = 0.1, 1, 10 and 100 were considered in a compression test in order to study the influence of the stiffness parameter on the elastic response of the material. A rectangle of 0.4 x 1 m is compressed by applying a distributed pressure of 1 MPa at the model top surface. The base is constrained in the y-direction while left and right sides are constrained in the x-direction. Fig. 9 shows these conditions. In the numerical model, it is assumed a Young’s Modulus equal to 1*104 MPa and Poisson’s ratio equal to 0.20. Pressure 0.2 m

1m Path

y

x

0.4 m

Fig. 9. Compression test and finite element mesh of the numerical model

11

Fig. 10 shows the numerical results associated with the effect of 𝛼 𝑡 on the elastic response of the model and their comparisons against a numerical solution considering only continuum elements. Fig. 11 shows the vertical stress distribution for 𝛼 𝑡 = 0.1 and 𝛼 𝑡 = 10. As one would expect, the general response of the model is sensitive to the relation 𝛼 𝑡 < 1, but tends to describe the continuum response as 𝛼 𝑡 > 1. In this test, 𝛼 𝑡 ≥ 10 ensures the appropriate representation of the elastic response of the material prior to the onset of hydraulic fracture (Fig. 10 and Fig. 11).

1 0.9 0.8

Path(m)

0.7 0.6 0.5 0.4 Continuum CZM (α/t) = 0.1 CZM (α/t) = 1 CZM (α/t) = 10 CZM (α/t) = 100

0.3 0.2 0.1 0 0 Fig. 10. Effect of

-0.05 -0.1 Displacement y(mm)

𝛼 𝑡

-0.15

of the interface element on the elastic response of the model

12

σy (MPa)

Deformation scale = 1500

0.252 -0.113 -0.479 -0.845 -1.210

𝛼 = 0.1 𝑡

Deformation scale = 1500 σy (MPa) 0.177 -0.154 -0.487 -0.819 -1.151

𝛼 = 10 𝑡

𝛼

Fig. 11. Vertical stress distribution for 𝑡 = 0.1 and

5

𝛼 𝑡

= 10.

KGD fracture problem

The presented numerical approach is verified against the analytical KGD solution and FEM solution presented by Carrier and Granet [36]. A Newtonian incompressible fluid is injected at a constant rate 𝑄𝑜 = 5 *10-4 m2/s. The 60 x 45 m domain is subjected to a compressive initial stress state 𝜎𝑜 = 3.7 MPa. Fig. 12 shows a schematic representation of the model. The comparison results between intrinsic Cohesive Zone Model (CZM) against the numerical and analytical results reported by Carrier and Granet [36] are displayed in Fig. 13 and Fig. 14. Table. 1 lists the rock, hydraulic fracture and pumping parameters used in the numerical model.

13

Injection point

σ0

60 m

45 m

σ0

Fig. 12. Schematic model and FEM mesh. Table. 1. Hydraulic and mechanical properties for the numerical tests Categories

Rock

Hydraulic fracture Pumping parameters

Variables Young's Modulus, E Poisson's ratio, υ Permeability, k Biot Coefficient Biot Modulus Porosity, ϕ Tensile strength, 𝜎0𝑛 Fracture energy, 𝐺𝑐 Leak-off, cb=ct Fluid viscosity, μ Injection rate, 𝑄𝑜

Units (kPa) (-) (m2) (-) (kPa) (%) (kPa) (kN/m) (m*(kPa.s)-1) (centipoise) (m2/s)

Values 1.70E+07 0.2 10-15,10-16 0.75 68700 20 1250 0.120 1.00E-6 0.1 5.00E-04

Fig. 13. Injection net pressure predicted by the intrinsic Cohesive Zone Model (CZM) against the numerical and analytical results reported by Carrier and Granet [36] (a) and pore pressure distribution (b) at t = 10 s.

For the net fluid injection pressure, 𝑝𝑓 + 𝜎𝑜 ( Fig. 13a), a small deviation from the analytical solution is observed in both numerical results. These differences can be related with the so-called back stress effect [60,61]. The pore pressure 𝑝𝑓 is

14

uniform in the fracture, as illustrated in Fig. 13b. This effect is directly related to the low viscosity of the injected fluid.

Fig. 14. Comparison results between the intrinsic Cohesive Zone Model (CZM) against the numerical and analytical results reported by Carrier and Granet [36]: fracture length (a) and fracture profile at t = 10 s (b)

The predicted fracture length (Fig. 14a) and the lips displacement (Fig. 14b) at t = 10 s, perfectly match the numerical and analytical solutions reported by Carrier and Granet [36]. 6

Methodology validation against laboratory tests

Several experimental examples have been simulated in order to validate and demonstrate the capabilities of the proposed numerical methodology. The first example studies the hydraulic fracture propagation in a prefractured rock. Mesh dependency and computational time of the presented methodology are also evaluated. The second example demonstrates the unrestricted fracture propagation under normal stress regime. In this case, hydraulic fracture propagation is initiated from an unfavorable perforation angle. Finally, the effect of the stress contrast and intercepted angle on the hydraulic and natural fracture interaction is simulated. 6.1 Hydraulic fracture propagation in a fractured rock sample This example is presented to give a better understanding of the hydraulic fracture behavior in a naturally fractured formation. Khoei [62] performed a hydraulic fracturing experimental test in two naturally fractured rocks under plane strain conditions. Moreover, the authors [62] compared the measured experimental data with those of numerical results obtained from XFEM method. Here, we compare our numerical results with CZM elements and those of Khoei [62] from two laboratory tests. A schematic representation of the geometry, boundary conditions of the simulated test are illustrated in Fig. 15. Table. 2 presents the dimensions of the models as well as the coordinates of the initial notch and the natural fracture, and material properties of the rock specimen. The fluid is injected with a constant pressure 𝑃𝑖𝑛𝑗 = 39.3 MPa through the hydro-fracture mouth. The fracturing fluid is assumed to be water with viscosity 𝜇 = 1e ― 03 Pa s.

15

Fig. 15. Schematic view of the geometry and boundary conditions of hydraulic fracturing experimental tests by Khoei [62] Table. 2. Geometry and mechanical properties of the rock sample

Specimen 1 Specimen 2

Specimen 1 Specimen 2

Width w(mm)

Height H(mm)

(x1,y1) (mm)

(x2,y2) (mm)

yo (mm)

111 110

45 54

(30.5;4) (12.13;8.82)

(55.4;42.16) (98.12;46.07)

27.9 26.9

Young´s Modulus E (GPa)

Poisson ratio ν

36.5 32.5

0.25 0.25

Tensile strength ft(mm) 29.2 22.3

Fracture energy Gc(J/m2) 330 330

We consider that the natural fracture stiffness and strength are equivalent to 10% of the corresponding rock properties, hence natural fractures are more sensitive to stress changes than the rock matrix. As referred by Khoei [62], the rock is considered impermeable in the first test. Accordingly, four different triangular meshes are evaluated in order to study the mesh dependency and the computational time of the proposed approach. We use unstructured meshes based on Delaunay algorithm in order to reduce certain bias on the crack propagation. Only the left region, where the crack will appear, was refined and fragmented with constant element size 𝑙𝑖 (red elements in Fig. 16). Interface element length 𝑙𝑖 = 2 mm, 1.2 mm, 0.9 mm and 0.7 mm were considered. Table. 3 shows the total number of elements, the nodes and the processing time on a desktop computer equipped with an Intel® Xeon® CPU X5650 @ 2.67 GHz (2 processors) for the adopted meshes. In this work, no attempt was made to optimize the CPU time in the different analyzes.

Fig. 16. Finite element triangular mesh. Interface element length 𝑙𝑖 = 2 mm

16 Table. 3. Mesh specifications and processing time

Mesh 1 Mesh 2 Mesh 3 Mesh 4

li (mm)

Nodes

Elements

Time (s)

2 1.2 0.9 0.7

5277 14405 25227 43535

3963 10877 19098 32918

390 751 1572 3673 Mesh 2

Mesh 1 σx (MPa)

σx (MPa)

8.385

8.348

0.749

0.051

-6.888

-8.247

-14.52

-16.54

-22.16

-24.84

σx (MPa)

Mesh 3 σx (MPa)

14.31

15.27

5.790

7.047

-2.730

-1.170

-11.25

-9.388

-19.77

-17.61

Mesh 4

Fig. 17. Stress distribution in the 𝑥 direction for meshes 1, 2, 3 and 4 when the HF meets the NF. Displacements are magnified 20 time for demonstration purposes.

The influence of mesh size on the evolution of the hydraulic fracture propagation pattern is depicted in Fig. 17. It is clear that the fracture patterns are influenced by mesh refinement. Fig. 18 shows the comparison of hydraulic fracture trajectories between laboratory test 1 [62] and the numerical solutions. For Mesh 1 (𝑙𝑖 = 2 mm), the fracture is deviated from its correct path. This effect is reduced for more refined meshes. However, higher refinement results in higher computational cost (Table 3). We notice that for Mesh 3 (𝑙𝑖 = 0.9 mm) and Mesh 4 (𝑙𝑖 = 0.7 mm), the hydraulic fracture paths differ very little between the experimental test 1 [62], Mesh 3 (𝑙𝑖 = 0.9 mm) and Mesh 4 (𝑙𝑖 = 0.7 mm). On the other hand, the processing time for the model with Mesh 4 is 2.34 times higher than for Mesh 3 (Table. 3). Therefore, the initial mesh should be well balanced in terms of refinement in order to obtain a reliable solution and optimize the computational cost.

Experimental (Khoei,2015) Mesh 1 Mesh 2 Mesh 3 Mesh 4

Fig. 18. Comparison of hydraulic fracture trajectory between laboratory test 1 [62] and intrinsic CZM fragmentation.

In the second test, we have considered rock permeability κ = 0.1 mD and porosity ϕ = 20% in order to demonstrate the capability of the presented approach to simulate hydraulic fracture propagation in a porous media. 4-node bilinear quadrilateral displacement and pore pressure plane strain elements (CPE4P) are

17

used to model the porous rock. Fig. 19 shows the pore pressure distribution and the minimum principal stress for this test. In Fig. 19 a, it is observed that the fluid flow inside the hydraulic fracture is transferred to the natural fracture. Fig. 20 shows the hydraulic fracture trajectory between laboratory test 2 and numerical simulations with XFEM [62] against intrinsic CZM fragmentation. POR(M Pa) 20.98 13.33 9.333 5.300 0.000

(a) Deformation scale: 10 σmin (M Pa) 8.376 -7.198 22.77 -38.34 -53.92

(b) Fig. 19. Hydraulic fracture and natural fracture interaction: pore pressure distribution (a); minimum principal stress (b)..

Experimental (Khoei,2015) XFEM (Khoei,2015) Intrinsic CZM

Fig. 20. Comparison of hydraulic fracture trajectory between experimental test 2, numerical simulations with XFEM [62] and intrinsic CZM fragmentation.

As shown in Fig. 20, the numerical simulations predict the hydraulic fracture trajectory accurately. Experimental conditions, rock specimen heterogeneity, permeability, fracture propagation criteria, among others may explain the slight discrepancies. It is evident that the natural fracture affects the hydraulic fracture trajectory, which tends to curve when approaching the NF. Therefore, unrestricted fracture propagation is essential to predict complex fracture patterns in the simulation of hydraulic fracturing in fractured formations.

18

6.2

Hydraulic fracture propagation from oriented perforation

Abbas [63,64] studied the effect of the perforation orientation on hydraulic fracturing. A series of laboratory fracturing experiments were performed using rectangular blocks of hydrostone (gypsum cement) with dimensions of 150 x 150 x 250 mm. The wellbore was drilled in the center of the block aligned with the vertical stress. A series of perforation orientations was considered: 𝜃 = 0, 15, 30, 45, 60, 75, and 90 degrees from the maximum horizontal stress. All samples were confined in a triaxial loading vessel subject to the principal stresses of 20.7 MPa in vertical direction, and 17.2 and 9.65 MPa minimum and maximum horizontal stresses, respectively. No pore fluid was present in the sample blocks [63,64]. Fig. 21 shows a schematic representation of the core sample and the quadrilateral mesh. Table. 4 presents the physical and mechanical properties. The numerical simulation is performed considering 𝜃 = 45 and 60 degrees. This choice is made since the hydraulic fracture initiates and propagates from perforations only. In addition, photographs are available to compare with the numerical models. The dissipated fracture energy was not specified in the experimental test. Therefore, we assumed a fracture energy of 0.1 N/mm. Homogenous and isotropic rock material were assumed and the casing or cement was not considered in the numerical models. 𝜎�

𝜎�

𝜎�

250 𝑚𝑚

Wellbore

Wellbore 𝜎�

θ°

𝜎�

Perforation 𝜎� z

y x

𝜎� 150 𝑚𝑚

Fig. 21. Schematic view of the geometry and finite element mesh of the numerical model Table. 4. Geometry and mechanical properties of the rock sample Rock sample properties Dimension well radius Perforation length Young's Modulus Poisson ratio Min Horizontal Stress Max Horizontal Stress vertical Stress Permeability Porosity Tensile strength (Brazilian) Fluid viscosity Injection rate

Units (mm) (mm) (mm) (MPa) (-) (MPa) (MPa) (MPa) (mD) (%) (MPa) (cp) (cm3/s)

Values 150 x 150 x 250 9.5 6.35 1.43*104 0.21 9.650 17.240 20.680 9.5 25 5.564 14 0.5

19

Fig. 22 shows the pore pressure distribution of the numerical simulation and some photographs from available experimental tests. θ = 45º

Deformation scale: 5

θ = 45º

POR(MPa) 22.30 8.330

𝜎�

𝜎�

5.833 2.500

𝜎�

0.000

𝜎�

θ = 60º

θ = 60º Deformation scale: 5 POR(MPa) 21.40 16.00

𝜎�

𝜎�

10.61 5.210 -0.18

𝜎�

Numerical model

𝜎� Experimental test

Fig. 22. Comparison of the hydraulic fracture trajectory between: intrinsic CZM fragmentation approach and laboratory tests [63]; for the perforation orientations at 𝜃 = 45° and 60°.

There is a good agreement between the numerical and the experimental results in terms of hydraulic fracture path (Fig. 22). It can be noted that the fracture initially propagates along the orientation of the perforations and gradually turns towards the most favorable direction. 6.3

Interaction between hydraulic and natural fractures

Fracture interaction is one of the most relevant issues to consider in the management of naturally fractured reservoirs. This includes the interaction between hydraulic fractures and natural fractures as well as the interaction among themselves. In order to study the effect of the orientation of the natural and horizontal stress contrast on hydraulic fracture propagation, Blanton [65,66] performed a group of hydraulic fracturing experiments in laboratory. Hydrostone blocks of 30 cm x 30 cm were fabricated with a prefracture oriented at different angles (30°, 45°, 60°, 90°) and subject to different triaxial compressive stresses. After that, fluid was injected to propagate the hydraulic fracture. Fig. 23 shows the test schematic model and the quadrilateral mesh. Table 5 presents the experimental conditions and interaction types identified in the tests. The prefractured surface is

20

frictional and non-cohesive with friction coefficient of 0.75. The tensile strength of the hydrostone is 3.1 MPa; Young’s Modulus is 17 GPa; Poisson’s ratio is 0.22 and the fracture energy is 0.1 kN/m. The fracturing fluid is injected at a constant flow rate of 8.20e-07 m3/s. The fluid viscosity used in these experiments is not specified in the reference. In this work, the fluid is assumed to be water. The natural fracture stiffness is equivalent to 10% of the rock stiffness. The behavior of the NF is governed by the Mohr-Coulomb criterion with zero tensile cutoff. The numerical simulations were carried out for the experimental tests CT - 7, 9, 11, 8, 4 (Table. 5). These tests have been chosen because photographs are available for comparison with the numerical models. Fig. 24-Fig. 28 show the fracture opening for these numerical cases.

Fig. 23. Test schematic model [65] and finite element mesh. Table. 5. Experimental condition and results for hydraulic fracture experiments on hydrostone [65] Test #

NF orientation (°)

CT- 7 CT- 9 CT-11 CT-12 CT-13 CT-14 CT-22 CT- 8 CT- 21 CT- 4 CT- 20

30 30 45 45 45 45 45 60 60 60 90

Horizontal stresses (MPa) 𝜎𝐻𝑚𝑎𝑥 𝜎ℎ𝑚𝑖𝑛 19 20 20 18 16 14 10 20 14 12 14

10 5 5 5 5 5 5 5 5 10 5

interaction type Opening Arrest Arrest Arrest Arrest Arrest Opening Crossing Arrest Opening Crossing

21 Deformation s cale: 100

CT-7

HF propa gation Opening (μm) 54.01

CT-7

HF opens the NF

HF opens the NF

40.25 27.17

HF propa gation

13.87 0.000

Fig. 24. Hydraulic fracture - natural fracture interaction for intrinsic CZM fragmentation approach and test CT-7 [65]. Deformation s cale: 100

Opening (μm) 73.80

HF propa gation

HF i s a rres ted

HF i s a rres ted

56.40 39.00

HF propa gation

21.60

CT-9

CT-9

0.000

Fig. 25. Hydraulic fracture- natural fracture interaction for intrinsic CZM fragmentation approach and test CT-9 [65]. Deformation s cale: 100

Opening (μm) 53.01

HF propa gation

HF i s a rres ted

HF propa gation

HF i s a rres ted

41.27 29.55

CT-11

CT-11

17.82 0.000

Fig. 26. Hydraulic fracture- natural fracture interaction for intrinsic CZM fragmentation approach and test CT-11 [65].

22 Deformation s cale: 100 NF s ta ys cl osed

NF s ta ys cl osed Opening (μm) 95.97

HF propa gation

HF propagation

HF cros s the NF

HF cros ses the NF

72.52 49.08

CT-8

CT-8 NF s ta ys cl osed

25.63 NF s ta ys cl os ed

0.000

Fig. 27. Hydraulic fracture- natural fracture interaction for intrinsic CZM fragmentation approach and test CT-8 [65]. Deformation s cale: 100 NF s ta y cl osed Opening (μm) 39.07 31.12 23.20

HF propa gation

NF s ta y cl osed HF propagation

HF open the NF

CT-4

HF open the NF

CT-4

15.23 0.000

Fig. 28. Hydraulic fracture- natural fracture interaction for intrinsic CZM fragmentation approach and test CT-4 [65].

For test CT-7 (Fig. 24), the differential stress is 9 MPa and the intersecting angle is 30 degrees. In this case, the natural fracture is dilated because the fluid pressure is larger than the normal compressive stress on the natural fracture. In tests CT-9 and CT-11 (Figs. Fig. 25 and Fig. 26), the differential stress is 15 MPa and the intersecting angle is 30° and 45°, respectively. In those scenarios, the hydraulic fracture is arrested by the natural fracture once the NF dilates by slippage increasing fluid leakoff. Test CT-8 (Fig. 27) has an intersecting angle of 60 ° and the differential stress is 15 MPa. Although the differential stress is the same as that in test CT-8, the interaction between HF and NF is different. In this test, the HF crosses the NF due to a higher NF angle of approach and the dilatation of the natural fracture is restricted by the maximal horizontal stress. Finally, test CT-4 test (Fig. 28) has a differential stress of 2 MPa and the intersecting angle is 60°. In this case, the natural fracture opens once the differential stress is low and the fluid pressure is larger than the normal compressive stress on the natural fracture. In general, these results show that the angle of approach and the horizontal differential stresses have a first-order effect on the interaction between hydraulic and natural fractures (cross, open and arrest). Other parameters such as mechanical properties of the rock and of the fractures, injection flow, fluid viscosity, NF length and opening, distance between wellbore and NF, geometry of the natural fracture

23

network, among others can also influence on the final geometry of the hydraulic fracture. 7

Discussion and conclusion

This paper describes a novel methodology for the implementation of the preprocessing algorithm which is used to insert the triple-noded interface element into conforming meshes. Its development is motivated by the need for robust numerical approaches capable of simulating nonlinear coupling processes which are present in hydraulic fracturing treatments. We introduce a simple but effective topological data structure that considers each finite element independent from others. The proposed topological structure ensures the generation of new nodes and faces which are necessary to create the conventional double-nodded interface element. Subsequently, a special insertion of middle nodes ensures the continuity of the fluid field. Finally, the workflow for unrestricted hydraulic fracture simulation is presented where three systems have been linked in a non-intrusive and versatile manner. The developed scheme is validated against three cases of fracturing experimental studies performed in laboratory. The first example studies the influence of the natural fracture on the trajectory of the hydraulic fracture propagation. Mesh dependency and computational time of the proposed methodology are also evaluated. The second example demonstrates the unrestricted fracture propagation under normal stress regime. In this case, hydraulic fracture propagation is initiated from an unfavorable perforation angle. Finally, the third example shows the effects of the stress contrast and the intersecting angle on the hydraulic and natural fracture interaction. The results show that the final trajectory of the hydraulic fracture propagation is influenced by changes in the perforation orientation, horizontal differential in situ stresses and the presence and interaction with natural fractures. Overall, it may be stated that the numerical procedure provides very good agreement with the laboratory tests and shows the robustness of the proposed computational algorithm. Some advantages and limitations of the proposed methodology are discussed briefly. Simplicity in terms of mesh data structure and successful representation of a localized failure and complex hydraulic crack patterns are the most important features of the developed intrinsic CZM. Furthermore, parallel implementation is straightforward once the topological mesh connectivity remains constant during the simulation process. With the proposed methodology, it is possible to simulate a reinitiated fracture from the crack tip, crossing with an offset, branching and multiple cracks, not worrying about multiple crack interactions. The adopted cohesive crack approach avoids the singularity issues of the stress at crack tip, which is present in linear elastic fracture mechanics. The main drawback in the developed method is mesh dependence because the fracture can only propagate along general element edges. The numerical results reveal that coarse meshes can affect the correct fracture patterns. This effect is reduced for more refined meshes. However, higher refinement results in higher computational cost. Therefore, the initial mesh should be well balanced in terms of refinement aiming to obtain reliable solution and optimize the computational cost. In addition, we use unstructured meshes based on Delaunay algorithm in order to reduce certain bias on the crack propagation. Another well-known problem of the

24

intrinsic CZM is that it introduces an artificial compliance that alters the elastic response of the material prior to the onset of hydraulic fracture. This phenomenon “artificial compliance” was reduced increasing the initial stiffness of the cohesive elements in terms of a hardening factor. The proposed mesh fragmentation is extensible to three dimensions. We leave the 3D simulation of unrestricted hydraulic crack propagation for future work. However, the current state of the implementation can be applied to study the effect of geomechanical properties of the rock and natural fractures, hydraulic properties of the fluid and injection rate among others - on the hydraulic fracture propagation. In addition, the methodology can simulate multi-staged fracturing into multilayered fractured rock formations. That is an important issue in order to reduce the environment associated risks and increase the oil recovery in unconventional reservoirs. Acknowledgement The authors gratefully acknowledge support from Petrobras and from the Carlos Chagas Filho Foundation for Research Support of Rio de Janeiro State (FAPERJ). References

[1] [2] [3]

[4] [5] [6] [7] [8] [9] [10]

Sitharam TG, Sridevi J, Shimizu N. Practical equivalent continuum characterization of jointed rock masses. Int J Rock Mech Min Sci 2001;38:437– 48. doi:https://doi.org/10.1016/S1365-1609(01)00010-7. Zhao J, Zhu W. Stability Analysis and Modelling of Underground Excavations in Fractured Rocks. 1st ed. 2004. Mejia Sanchez EC, Pezo Zegarra E, Oliveira MFF, Roehl D. Application of a 2D Equivalent Continuum Approach To the Assesment of Geological Fault Reactivation in Reservoirs. Proc XXXVI Iber Lat Am Congr Comput Methods Eng 2016. doi:10.20906/cps/cilamce2015-0880. Barenblatt GI, Zheltov IP, Kochina IN. BASIC CONCEPTS IN THE THEORY OF SEEPAGE OF HOMOGENEOUS LIQUIDS IN FISSURED ROCKS [STRATA]. J Appl Math Mech 1960;24:852–64. Kazemi H, Merrill LS, Porterfield KL, Zeman PR. Numerical Simulation of Water-Oil Flow in Naturally Fractured Reservoirs. Soc Pet Eng J 1976;16:317– 26. doi:10.2118/5719-PA. Hill AC, Thomas GW. A New Approach for Simulating Complex Fractured Reservoirs. Middle East Oil Tech Conf Exhib 1985:429–36. doi:10.2118/13537MS. Rueda Cordero JA, Mejia Sanchez EC, Roehl D. Integrated discrete fracture and dual porosity - Dual permeability models for fluid flow in deformable fractured media. J Pet Sci Eng 2019;175:644–53. doi:10.1016/j.petrol.2018.12.053. Goodman RE, Taylor RL, Brekke TL. A MODEL FOR THE MECHANICS OF JOINTED ROCKS. J Soil Mech Found Div 1968. Day RA, Potts DM. Zero thickness interface elements—numerical stability and application. Int J Numer Anal Methods Geomech 1994;18:689–708. doi:10.1002/nag.1610181003. Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434. doi:10.1016/00225096(94)90003-5.

25

[11]

[12] [13] [14]

[15] [16] [17]

[18] [19]

[20]

[21] [22] [23] [24] [25] [26]

Ortiz M, Pandolfi A. Finite-Deformation Irreversible Cohesive Elements for Three-Dimensional Crack-Propagation Analysis. Int J Numer METHODS Eng Int J Numer Meth Engng 1999;1282:1267–82. doi:10.1002/(SICI)10970207(19990330)44:9<1267::AID-NME486>3.0.CO;2-7. Turon A, Dávila CG, Camanho PP, Costa J. An Engineering Solution for using Coarse Meshes in the Simulation of Delamination With Cohesive Zone Models. Nasa/Tm-2005-213547 2005:1–26. doi:10.1007/978-3-540-79056-3_2. Rueda JC, Noreña N V, Oliveira MFF, Roehl D. Numerical Models for Detection of Fault Reactivation in Oil and Gas Fields. 48th US Rock Mech Symp 2014:8. Pereira FLG, Roehl D, Paulo Laquini J, Oliveira MFF, Costa AM. Fault reactivation case study for probabilistic assessment of carbon dioxide sequestration. Int J Rock Mech Min Sci 2014;71:310–9. doi:10.1016/j.ijrmms.2014.08.003. Zavattieri PD, Espinosa HD. Grain Level Analysis of Crack Initiation and 2001;4:4291–311. doi:10.1016/S1359-6454(01)00292-0. Geubelle PH, Baylor JS. Impact-induced delamination of composites: A 2D simulation. Compos Part B Eng 1998;29:589–602. doi:10.1016/S13598368(98)00013-4. Wu L, Tjahjanto D, Becker G, Makradi A, Jérusalem A, Noels L. A micro–mesomodel of intra-laminar fracture in fiber-reinforced composites based on a discontinuous Galerkin/cohesive zone method. Eng Fract Mech 2013;104:162– 83. doi:https://doi.org/10.1016/j.engfracmech.2013.03.018. Dai Q, Ng K. 2D cohesive zone modeling of crack development in cementitious digital samples with microstructure characterization. Constr Build Mater 2014;54:584–95. doi:10.1016/j.conbuildmat.2013.12.095. Wang XF, Yang ZJ, Yates JR, Jivkov AP, Zhang C. Monte Carlo simulations of mesoscale fracture modelling of concrete with random aggregates and pores. Constr Build Mater 2015;75:35–45. doi:https://doi.org/10.1016/j.conbuildmat.2014.09.069. Benedetto MF, Caggiano A, Etse G. Virtual elements and zero thickness interface-based approach for fracture analysis of heterogeneous materials. Comput Methods Appl Mech Eng 2018;338:41–67. doi:https://doi.org/10.1016/j.cma.2018.04.001. Pandolfi A, Ortiz M. Solid Modeling Aspects of Three-Dimensional Fragmentation. Eng Comput 1998;14:287–308. doi:10.1007/BF01201761. Pandolfi A, Ortiz M. An Efficient Adaptive Procedure for Three-Dimensional Fragmentation Simulations. Eng Comput 2002;18:148–59. doi:10.1007/s003660200013. Paulino GH, Celes W, Espinha R, Zhang ZJ. A general topology-based framework for adaptive insertion of cohesive elements in finite element meshes. Eng Comput 2008;24:59–78. doi:10.1007/s00366-007-0069-7. Dooley I, Mangala S, Kale L, Geubelle P. Parallel simulations of dynamic fracture using extrinsic cohesive elements. J Sci Comput 2009;39:144–65. doi:10.1007/s10915-008-9254-0. Alhadeff A, Celes W, Paulino GH. Mapping Cohesive Fracture and Fragmentation Simulations to Graphics Processor Units. Int J Numer Methods Eng 2015;103:859–93. doi:10.1002/nme.4842. Radovitzky R, Seagraves A, Tupek M, Noels L. A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method. Comput Methods Appl Mech Eng 2011;200:326–44.

26

[27] [28] [29] [30] [31]

[32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]

[43]

doi:10.1016/j.cma.2010.08.014. Zhang Z (Jenny), Paulino GH, Celes W. Extrinsic cohesive modelling of dynamic fracture and microbranching instability in brittle materials. Int J Numer Methods Eng 2007;72:893–923. doi:10.1002/nme.2030. Segura JM, Carol I. Coupled HM analysis using zero-thickness interface elements with double nodes—Part II: Verification and application. Int J Numer Anal Methods Geomech 2008;32:2103–23. doi:10.1002/nag.730. Segura JM, Carol I. Coupled HM analysis using zero-thickness interface elements with double nodes. Part I: Theoretical model. Int J Numer Anal Methods Geomech 2008;32:2083–101. doi:10.1002/nag.735. Ng KLA, Small JC. Behavior of joints and interfaces subjected to water pressure. Comput Geotech 1997;20:71–93. doi:https://doi.org/10.1016/S0266352X(96)00015-8. Baldoni F, Millard A. A finite element formulation for the coupled hydromechanical behaviour of porous rock joints. Poromechanics A Tribut. to Maurice A Biot. Proc. 1st Biot Conf., Taylor & Francis: Louvain-la-Neuve, Belgium; 1998, p. 339–44. Guiducci C, Pellegrino A, Radu J-P, Collin F, Charlier R. Numerical modeling of hydro-mechanical fracture behavior. Numer Model Geomech 2002:293–9. doi:10.1111/j.1752-8062.2012.00338.x. Segura JM, Carol I. On zero-thickness interface elements for diffusion problems. Int J Numer Anal Methods Geomech 2004;28:947–62. doi:10.1002/nag.358. Haddad M, Du J, Vidal-Gilbert S. Integration of Dynamic Microseismic Data With a True 3D Modeling of Hydraulic-Fracture Propagation in the Vaca Muerta Shale. SPE J 2017;22:1714–38. doi:10.2118/179164-PA. Chen Z, Bunger AP, Zhang X, Jeffrey RG. Cohesive zone finite element-based modeling of hydraulic fractures. Acta Mech Solida Sin 2009;22:443–52. doi:10.1016/S0894-9166(09)60295-0. Carrier B, Granet S. Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model. Eng Fract Mech 2012;79:312–28. doi:10.1016/j.engfracmech.2011.11.012. Zielonka MG, Searles KH, Ning J, Buechler SR. Development and Validation of Fully-Coupled Hydraulic Fracturing Simulation Capabilities. SIMULIA Community Conf SCC2014 2014:1–31. Guo J, Zhao X, Zhu H, Zhang X, Pan R. Numerical simulation of interaction of hydraulic fracture and natural fracture based on the cohesive zone finite element method. J Nat Gas Sci Eng 2015;25:180–8. doi:10.1016/j.jngse.2015.05.008. Gonzalez M, Taleghani AD, Olsen JE. SPE-173384-MS A Cohesive Model for Modeling Hydraulic Fractures in Naturally Fractured Formations 2015;d:3–5. doi:10.2118/173384-MS. Nikam A, Awoleke OO, Ahmadi M, Fairbanks A. Modeling the Interaction Between Natural and Hydraulic Fractures Using. Soc Pet Eng 2016. doi:10.2118/180364-MS. Rueda J, Mejia C, Roehl D. Numerical simulation of hydraulic and natural fracture interaction and propagation. Cilamce 2017 2017. Dahi Taleghani A, Gonzalez-Chavez M, Yu H, Asala H. Numerical simulation of hydraulic fracture propagation in naturally fractured formations using the cohesive zone model. J Pet Sci Eng 2018;165:42–57. doi:10.1016/j.petrol.2018.01.063. Rueda J, Mejia C, Roehl D. Hydro-mechanical modeling of hydraulic fracture

27

[44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61]

propagation and its interactions with frictional natural fractures. Comput Geotech n.d. Dahi-Taleghani A, Olson JE. Numerical Modeling of Multistranded-HydraulicFracture Propagation: Accounting for the Interaction Between Induced and Natural Fractures. SPE J 2011;16:575–81. doi:10.2118/124884-PA. Keshavarzi R, Jahanbakhshi R. Investigation of Hydraulic and Natural Fracture Interaction: Numerical Modeling or Artificial Intelligence? Eff Sustain Hydraul Fract 2013:1039–58. doi:10.5772/56382. Khoei AR, Vahab M, Hirmand M. Modeling the interaction between fluid-driven fracture and natural fault using an enriched-FEM technique. Int J Fract 2016;197:1–24. doi:10.1007/s10704-015-0051-0. Cruz F, Roehl D, Vargas E do A. An XFEM element to model intersections between hydraulic and natural fractures in porous rocks. Int J Rock Mech Min Sci 2018;112:385–97. doi:10.1016/j.ijrmms.2018.10.001. Gutierrez ER, Mejía SEC, Roehl D, Romanel C. XFEM Modeling Of Stress Shadowing In Multiple Hydraulic Fractures In Multi-Layered Formations. J Nat Gas Sci Eng 2019. Nguyen VP, Lian H, Rabczuk T, Bordas S. Modelling hydraulic fractures in porous media using flow cohesive interface elements. Eng Geol 2017;225:68–82. doi:10.1016/j.enggeo.2017.04.010. Truster TJ. On interface element insertion into three-dimensional meshes. Eng Fract Mech 2016;153:171–4. doi:10.1016/j.engfracmech.2015.12.019. Truster TJ. DEIP, discontinuous element insertion Program — Mesh generation for interfacial finite element modeling. SoftwareX 2018;7:162–70. doi:10.1016/j.softx.2018.05.002. Settari A, Price H.S. Simulation of Hydraulic Fracturing In Low-Penneability Reservoirs. Soc Pet Eng J 1984;24:141–52. doi:10.2118/8939-PA. Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. doi:https://doi.org/10.1016/0020-7683(95)00255-3. Camanho PP, Dávila CG. Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials. Nasa/Tm-2002-211737 2002:1–37. doi:10.1177/002199803034505. Geuzaine C, Remacle J-F. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Int J Numer Methods Eng 2009;79:1309–31. doi:10.1002/nme.2579. Chen Z. Finite element modelling of viscosity-dominated hydraulic fractures. J Pet Sci Eng 2012;88–89:136–44. doi:https://doi.org/10.1016/j.petrol.2011.12.021. Smith M. ABAQUS/Standard User’s Manual, Version 16 2016. Jin W, Arson C, Busetti S. Simulation of Mode II Unconstrained Fracture Path Formation Coupled with Continuum Anisotropic Damage Propagation in Shale. 50th US Rock Mech Symp 2016:10. Jin W, Xu H, Arson C, Busetti S. Computational model coupling mode II discrete fracture propagation with continuum damage zone evolution. Int J Numer Anal Methods Geomech 2017;41:223–50. doi:10.1002/nag.2553. Vandamme LM, Roegiers J-C. Poroelasticity in Hydraulic Fracturing Simulators. J Pet Technol 1990;42:1199–203. doi:10.2118/16911-PA. Kovalyshen Y. Fluid-Driven Fracture in Poroelastic Medium. University of Minnesota, 2010.

28

[62]

[63] [64] [65] [66]

Khoei AR, Hirmand M, Vahab M, Bazargan M. An enriched FEM technique for modeling hydraulically driven cohesive fracture propagation in impermeable media with frictional natural faults: Numerical and experimental investigations. Int J Numer Methods Eng 2015;104:439–68. doi:10.1002/nme.4944. Abass HH, Brumley JL, Venditto JJ. Oriented Perforations - A Rock Mechanics View. SPE Annu Tech Conf Exhib 1994. doi:10.2118/28555-MS. Abass HH, Hedayati S, Meadows DL. Nonplanar Fracture Propagation From a Horizontal Wellbore: Experimental Study. SPE Prod Facil 1996;11:133–7. doi:10.2118/24823-pa. Blanton TL. An Experimental Study of Interaction Between Hydraulically Induced and Pre-Existing Fractures. SPE Unconv Gas Recover Symp 1982. doi:10.2118/10847-MS. Blanton TL. Propagation of Hydraulically and Dynamically Induced Fractures in Naturally Fractured Reservoirs. SPE Unconv Gas Technol Symp 1986. doi:10.2118/15261-MS.

29

Highlights  A novel mesh fragmentation approach to simulate unrestricted hydraulic fracture propagation.  Effective topological data structure to insert the especial triple-noded interface elements  Special treatment of the middle node of the interface to ensure the continuity of the fluid field at fracture intersections.  Workflow for simulation of unrestricted hydraulic fracturing in naturally fractured formations.