Geofluids

Volume 2019, Article ID 3958583, 23 pages

https://doi.org/10.1155/2019/3958583

## Discrete Element Analysis for Hydraulic Fracture Propagations in Laminated Reservoirs with Complex Initial Joint Properties

^{1}State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, China^{2}School of Mechanics and Civil Engineering, China University of Mining and Technology, Xuzhou 221116, China^{3}School of Engineering, University of Tasmania, Hobart, Tasmania 7001, Australia

Correspondence should be addressed to J. G. Wang; moc.oohay@gjwsun

Received 30 April 2019; Accepted 17 September 2019; Published 7 November 2019

Academic Editor: Maurizio Barbieri

Copyright © 2019 Fakai Dou et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Previous studies on hydraulic fracturing mainly focus on the effects of the in-situ stress state, permeability, fracturing fluids, and approach angle in homogeneous rocks, but the impacts of joint mechanical properties in laminated shale reservoirs on the propagation and formation of the fracture network are still unclear. In this study, a coupled fluid-mechanical model was developed to investigate the impacts of joint mechanical properties on hydraulic fracture propagation. Then, this model was validated with Blanton’s criterion and some experimental observations on fracture morphology. Finally, a series of numerical simulations were conducted to comparatively analyze the impacts of joint mechanical properties on the total crack number, the percentage and distribution of each fracture type, the process of crack propagation, and the final fracture morphology. Numerical results show that the cracking behaviors induced by joint mechanical properties vary with the approach angle. Joint strength has a significant influence on the generation of matrix tensile cracks. The tensile-to-shear strength ratio plays an even more important role in the shear slips of bedding planes and is conducive to the formation of complex fracture morphology.

#### 1. Introduction

The huge demand for unconventional natural gas has stimulated the exploration and exploitation of shale gas reservoirs. Because shale gas reservoirs have extremely low permeability and porosity, the gas production from a shale reservoir is usually enhanced by a horizontal well plus staged hydraulic fracturing [1–5]. After long-history environmental and geological impacts, shale gas reservoirs are laminated with many natural flaws such as fissures, faults, and bedding planes. These laminated structures of shale reservoirs make hydraulic fracturing behaviors even more complex [6, 7]. Among these natural fractures or flaws, bedding planes are widely considered as the biggest influence on the process of crack initiation, propagation, and coalescence [8]. In particular, layer orientation and joint mechanical properties are the sources of shale anisotropy and heavily affect the hydraulic fracturing performance in a shale reservoir. Therefore, the effects of joint mechanical properties in a shale reservoir on the formation of the fracture network during hydraulic fracturing should be carefully investigated.

Considerable investigations have been conducted on the interaction between induced and natural fractures based on field and laboratory observations. Field data from the Barnett shale showed that the propagation path of the hydraulic fracture presented multiple segments in different directions because of preexisting fractures. Hydraulic fractures can propagate about thousands of feet along the expected direction and hundreds of feet in the orthogonal direction during a staged hydraulic fracturing treatment in Barnett shale [9]. Laboratory observations [10, 11] indicated that hydraulic fractures propagated in both horizontal and vertical directions and the failure pattern was complex rather than a simple bi-wing fracture in the shale specimen. Norman and Jessen [12] reported that the cementing strength of weak planes played an important role in the orientation of crack propagation. Zhou et al. [13] investigated the influence of shear strength of natural fractures on crack propagation behaviors. They concluded that fracture morphology was mainly controlled by in-situ stress and natural fractures. Based on a series of large true triaxial hydraulic fracturing experiments, Tan et al. [14] studied the effects of the injection rate and fracturing fluid viscosity on the initiation and propagation of hydraulic fractures. Their results showed that fracturing fluid viscosity and injection rate also decided the complexity of the fracture network. Based on experimental observations, Liu et al. [15] concluded that an induced-fluid fracture propagated along with the least resistance or the shortest propagation path. Acoustic emission technology and gaseous tracer have also been applied in laboratory experiments to investigate fracturing behaviors [16, 17]. The aforementioned experiments mainly explored the hydraulic fracturing mechanism from external loading conditions and hydraulic fracturing treatments. However, the interaction between induced and natural fractures of shale is extremely complex; hydraulic fracturing mechanism in laminated shale gas reservoirs is still not clear.

Various criteria for the prediction of failure morphology have been proposed based on field and laboratory experimental results. Blanton [18] proved that an induced hydraulic fracture preferred to cross a natural fracture under high approach angles and deviatoric stresses. Warpinski and Teufel [19] proposed a relationship between deviatoric stress and approach angle by considering the effect of shear strength of natural fractures as well as pore pressure distribution. Based on laboratory experiments, Renshaw and Pollard [20] proposed a criterion for the propagation path of the hydraulic fracture when it perpendicularly penetrated a natural fracture. This Renshaw and Pollard’s criterion was extended to consider the interaction at nonorthogonal approach angles [21] and a more general case of the nonorthogonal cohesive natural interface [22]. Though these criteria did not consider the rock anisotropy, they are still effective tools to validate the feasibility of hydraulic fracturing models.

Various numerical methods have been developed for hydraulic fracturing modeling [23, 24]. For example, the Extended Finite Element Method (X-FEM) can highly refine mesh around the crack tip through the enrichment of freedom degrees in discontinuous behaviors [25, 26]. Meshless methods do not need mesh during the process of crack propagation and can obtain smooth stress around the crack tip [27]. A Boundary Element Method (BEM) approach was developed to solve the brittle anisotropic problems [28]. The Rock Failure Process Analysis (RFPA) was used to account for the interaction between hydraulic fractures and preexisting fractures [29]. Compared to these aforementioned continuum-based methods, the Particle Flow Code (PFC) proposed by Potyondy and Cundall [30] has also been developed to investigate the fracturing behaviors in a naturally fractured reservoir. Contact models are the essence of PFC. The PFC has great advantages in simulating the materials with particle-like behaviors just like rock and soil. After continuous improvements by Itasca Consulting Group [31, 32], the PFC is now efficient in investigating the interaction between fluid-induced fractures and preexisting fractures. Based on PFC^{2D}, Shimizu et al. [33] investigated the influence of rock brittleness on the complexity of fracture morphology. Their simulation results observed a positive correlation between rock brittleness and the complexity of fracture morphology. Yoon et al. [34] studied the effect of stress shadowing in a low permeability reservoir by using PFC^{2D} and proposed a method to optimize the hydraulic fracture network. Their optimization method has been successfully applied to a geothermal reservoir. Zhou et al. [35] studied the feasibility of the smooth joint model to mimic the preexisting fractures. Chong et al. [36] investigated the influence of natural fracture density on the complex fracture network based on the smooth joint model. Zhang et al. [37] studied the influence of deviatoric stress, the strength of the preexisting interface, permeability, fluid injection rate, and the viscosity of fracturing fluid on the interaction between induced and natural fractures. Their numerical model was used to optimize the design of hydraulic fracturing. However, most of the previous investigations only considered one or two main preexisting fractures in their model, while the shale gas reservoir as a sedimentary formation has massive layered structures. In this sense, a particle-based numerical method is a good choice when the interaction between induced and natural fractures is to be captured under various geological and environmental conditions.

In this study, a coupled fluid-mechanical model was established based on PFC^{2D} and validated by Blanton’s criterion. This model was then used to comparatively study the influence of the joint strength and tensile-to-shear strength ratio on the crack propagation as well as hydraulic fracture morphology under different approach angles. This paper is organized as follows. Section 2 presents the numerical model for discrete element analysis. Section 3 validates this numerical model through three scenarios of interaction between induced and preexisting fractures and Blanton’s criterion. Section 4 numerically investigates the impacts of the joint strength and tensile-to-shear strength ratio on fracture propagation and fracture network morphology. Their relative roles are explored, too. The conclusions are drawn in Section 5.

#### 2. Model of Discrete Element Analysis

##### 2.1. Fundamental Algorithm of the Bonded Particle Model

PFC is based on one of discontinuum methods, whose constitutive models are described by different bonds between particles or blocks. In two-dimensional problems, shale can be simulated by bonded circular discs. This study uses the parallel bond model to represent the shale matrix and uses the smooth joint model to mimic the mechanical properties of bedding planes. The fundamental algorithm for these contact models is briefed below.

###### 2.1.1. Parallel Bond Model

The contact algorithm of the parallel bond model is firstly proposed by Potyondy and Cundall [30]. This parallel bond model updates bond force and moment with the following algorithms: where the parallel bond force is divided into normal force and shear force , represents the contact area, and represent the normal and shear stiffness in per unit area, respectively, and represent the increments of normal and shear displacements, respectively, and represents the bending moment at the contact which is the product of bending stiffness and the increments of bending rotation .

Further, the tensile stress and shear stress in the parallel bond model are obtained as follows:

The moment-contribution factor is used to reflect the influence of bending moment on the normal stress at the parallel-bond periphery [38]. It ranges from 0.0 to 1.0. This paper takes the value of as 0.2, which is also the recommended value in PFC^{2D} user’s manual (version 5.0) [32]. is the bond radius, and is the inertia moment of the bond. The crack type and number are counted as follows: if the tensile stress exceeds the bond tensile strength , a tensile crack will be formed. If the shear stress exceeds the bond shear strength , a shear crack will be generated at the contact surface. The shear stress on the contact surface is calculated according to the Mohr-Coulomb strength criterion when the contact surface slips. Under this failure criterion, micro cracks will be generated in the shale matrix.

###### 2.1.2. Smooth Joint Model

The smooth joint model was originally proposed by Pierce et al. [39] and used to simulate the shear behavior of rock materials. Recently, it was used to model hydraulic fracturing [34]. Its ability in modeling the interaction between the induced and single preexisting fractures was validated by Zhou et al. [35]. A smooth joint model assumes that the joint orientation is perpendicular to the slide surface rather than aligned with the contact normal direction: where the smooth joint force is resolved into normal force and shear force . The normal force is oriented to the joint normal direction, and the shear force is along the shear direction. The updating algorithm of the smooth joint model is where and represent normal and shear stiffness in per unit area, respectively, and and represent the elastic deformation of the normal and shear displacement increments, respectively.

Simultaneously, the shear force is computed by

In our numerical simulations, the bedding planes are bonded before breakage. The normal and shear stresses acting on the smooth joint are calculated by

If the tensile stress exceeds the tensile strength , the bond will be broken by tension. After tensile breakage, the normal and shear forces of the smooth joint drop to zero. If the shear stress exceeds the shear strength , the bond will be broken in shear, but the contact force will change to trial shear force. If the sliding occurs, the forces are updated as

Here, is the dilation angle on the bedding plane. Under this failure criterion, the cracks of bedding planes will be judged and counted.

##### 2.2. Coupled Fluid-Mechanical Algorithm in PFC^{2D}

Incompressible fluid flow in the laminated reservoir is simulated according to the following fluid-mechanical algorithm. Al-Busaidi et al. [40] documented the development history of this fluid-mechanical algorithm and indicated that this algorithm was modified and introduced into a discrete element method by Cundall. This algorithm has two concepts: flow channel and reservoir. The fluid-mechanical algorithm assumes that contacts in the assembled particles are regarded as the “flow channels” and the pores enclosed by the adjacent particles are the “reservoirs” as shown in Figure 1. The formed reservoirs have been connected by the adjacent flow channels. The volume of each reservoir is calculated by surrounding short lines which are connecting the center of surrounded particles. All the reservoirs in the numerical model are connected by flow channels. Reservoirs and flow channels are combined to form the whole flow network. Fluid is stored in these reservoirs, and the differential pressure between adjacent domains is the driving force of fluid flow.