# Fracture Propagation Behavior of Jointed Rocks in Hydraulic Fracturing

Authors: Xiaoxi Men; Jiren Li; Zhihui Han
Journal: Advances in Materials Science and Engineering
 (2018)

---

## Abstract

Jointed rocks are typical examples of heterogeneous materials with joints. The existence of joints influences the physical properties of rock mass and propagation of fractures, which can affect production operations in engineering. A series of simulations is performed to understand the failure patterns and fracture propagation behavior of jointed rocks in hydraulic fracturing. Three points, that is, dip-angle joint, joint strength, and complex joints, are considered in the simulations. Results demonstrate three basic kinds of hydraulic fractures on jointed rock, namely, along the joint, across the joint, and partly along the joint and partly across the joint. The maximum principal stress is the control factor of fracture propagation in global scale, and the joint plane is the control factor of fracture propagation in local scale. In the propagation path, when the dip angle is small or the joint is weak, the fracture propagates along the joint; otherwise, the fracture propagates across the joint. In the multijoint interconnection models, hydraulic fractures propagate along joints in the tensile stress regions near the propagating fracture tip without dip angle limitation. Subsequently, the fractures connect with one another to form a complex fracture network based on the joint morphology.


---

## Body

## 1. Introduction

Joint is a structural plane of discontinuities formed by sedimentation. Jointed rock masses are commonly encountered in civil, mining, and petroleum engineering because they are universally distributed over the earth’s surface. As a control factor of rock mass stability [1], joints influence rock failure patterns and fracture propagation behavior. Hydraulic fracturing is commonly used in petroleum engineering, mining, and geotechnical industries, and it is carried out as an important technique to further improve production efficiency. The first experimental hydraulic fracturing treatment in the United States was conducted in the Hugoton Gas Field in Grant County, Kansas [2], and this technology has been widely used worldwide since 1947. In recent years, hydraulic fracturing was also used for jointed rock masses of sandstones, mudstones, shales, and coal seams to form complex fracture networks and improve the permeability of unconventional reservoirs [3–5]. Hydraulic fracture propagation on jointed rocks is complex due to the different morphologies and composition of joints. Extensive experimental and numerical studies on the actual kinematics of hydraulic fracture propagation of jointed rocks were performed by researchers all over the world [6–12]. Discontinuities dominate the geometry, deformation modulus, strength, failure behavior, and permeability of rocks, and even the local magnitude and direction of in situ stress fields [1]. A study on the mechanical influence of joints and the growth process of hydraulic fractures of jointed rocks is essential, and it is even helpful in determining efficient fracturing scenarios or production designs in field operations.The main purpose of this study is to numerically investigate how joints affect the initiation and propagation of hydraulic fractures in different conditions. Numerical simulations have been proven advantageous in modeling. Thus, the hydraulic fracturing process of parallel and complex joint models are simulated by using the Rock Failure Process Analysis (RFPA)2D2.0-Flow. The RFPA is a code to simulate fracture propagation following the continuum damage principle [13, 14] of heterogeneous rock materials. The RFPA has been used extensively, and the creditability of its simulations has been proven many times since its development [15–23].


## 2. Brief Introduction of the RFPA Code

Rocks are typical examples of heterogeneous natural geological materials [24]. These solid media are often described as homogeneous and continuous materials, and their failure processes are explained by elasticity-plasticity principles in classical mechanics. However, recent studies on the mechanical and unique characteristics of rock materials have shown the limitations of the classical method. Subsequently, rock failure is assumed to be a nonlinear process from mesoscopic weakening to macroscopic damaging [25–27].


### 2.1. Nonlinear Macroscopic Model Described by Linear Mesoscopic Elements

The RFPA code is based on finite element and statistical damage theory, and it can simulate the failure process of heterogeneous and permeable geomaterials. The heterogeneity of materials and the random distribution of defects of these materials are also considered by the RFPA code. The mesoscopic elements in the numerical model are assumed to be isotropic and homogeneous, and their mechanical properties (e.g., Young’s modulus, Poisson’s ratio, and strength properties, among others) are also assumed to be linear. The mesoscopic elements are statistically distributed (e.g., normal, Poisson, and Weibull distributions) to describe the mechanical properties of the nonlinear macromodel (Figure1). In RFPA, four-node isoparametric elements are used as the basic elements (Figure 2). In the simulation, the mesoscopic elements are defined as damaged if they satisfy the strength criterion. As stress increases, numerous elements fail. Moreover, when these elements connect with one another, the heterogeneous materials reach a failure process.Figure 1
Heterogeneous distribution of mechanical properties from linear mesoscopic mode to nonlinear macroscopic mode.Figure 2
Mesh of RFPA model.The statistical distribution densityfα of the mechanical parameter of the mesoscopic elements is defined by (1). The heterogeneity index m [12, 28] in the RFPA code is used to describe the heterogeneity of solid materials. A high m value indicates the presence of highly homogeneous materials, whereas a small m value denotes the existence of highly inhomogeneous materials (Figure 3). In the equation [29], α represents the mechanical parameters (e.g., Young’s modulus, Poisson’s ratio, tensional strength, and compressive strength) of the elements and α0 denotes the scale parameters related to the average values of the mechanical parameters.(1)fα=mα0αα0m−1exp−αα0m.Figure 3
Weibull’s distribution of the mechanical parameters of materials with different heterogeneitym value indices.


### 2.2. Primitive Phase Transition

The mesoscopic elements of an RFPA-based model are called primitives. Three primitive types (matrix, air, and contact) are considered for describing the deformation and failure processes of solid materials (Figure4). When the strain is located between the maximum tensile and compressive strains, the element is called a matrix primitive. A matrix primitive contains elastic and residual phases. The matrix primitive changes into air primitive when the deformation exceeds the maximum tensile strain, and it changes into contact primitive when the deformation exceeds the maximum compressive strain. Changes occur in the mechanical properties of elements during primitive transformation, and this change is called phase transition.Figure 4
Primitive phase transition.


### 2.3. Main Governing Equations of the RFPA-Flow Code

The geomaterial (rock) in the RFPA-Flow code is assumed to be fully saturated with fluid flow in accordance with Darcy’s law. In addition, the coupled process between stress/strain and fluid flow in the deforming rock mass is governed by Biot’s consolidation theory. The rock medium is an elastic-brittle material with residual strength, and the acoustic emission (AE) of rock also proves it. Thus, the mechanical behavior of rock material in the loading and unloading process conforms with elastic damage theory. The maximum tensile strength criterion and Mohr–Coulomb criterion are used as the damage judgment of the mesoscopic element in the RFPA code.The fundamental equations [30] of fluid-solid coupling in the RFPA code are as follows:(2)equilibrium  equation:∂σij∂xij+ρXj=0i,j=1,2,3,geometrical  equation:εij=12Ui,j+Uj,i,εv=ε11+ε22+ε33,constitutive  equation:σij′=σij−αpδij=λδijεv+2Gεij,seepage  equation:K∇2p=1Q∂p∂t−α∂εv∂t,where ρ is the volume density; σij is the total stress; εij is the total strain and εv is the volume strain; Ui is the displacement; p is pore pressure; λ is the Lame coefficient; δij is the Kronecher constant; G is the modulus of shear deformation; Q is Biot’s constant; K is the permeability coefficient; and α is the coefficient of pore pressure which is determined by the experiment.Under the same loading conditions, fluid flow behavior differs between the propagating fractures and the preexisting fractures. In the failure process of rock mass, the permeability of a rock changes all the time. Thus, the propagation of preexisting fractures, initiation of new fractures, coupling relationship of flow-stress-damage in these two kinds of fractures, and variations of rock permeability in the failure process are considered in the RFPA code.In the elastic state, the relationship of the stress and permeability coefficient is described by(3)Kf=K0e−bσ′,where K0 is the initial permeability coefficient; σ′ is the effective stress; and b is the coupling parameter.In describing the dependency of permeability on the stress and the damage, the following coupling equation [31] is introduced:(4)Kσ,p=ξK0e−βσii/3−αp,where α is the coefficient of pore pressure; ξξ≥1.0 is a mutation coefficient of permeability accounting for the increase in permeability once the element reaches the damage state; and β is the coupling parameter reflecting the influence of stress on the permeability coefficient. The value of the coefficients ξ, α, and β are determined experimentally, and they vary with the stress state and damage evolution of the rock. In the RFPA code, 0≤α<1 and ξ=1 are assumed for the mesoscopic element in the elastic stage (as indicated by D=0 in (7)). Once damage occurs, the permeability of the mesoscopic element increases significantly (i.e., see (8) and (11)). Furthermore, ξ=5 is assumed for the damaged element (denoted by 0<D<1), and α=1, ξ=100 is assumed for the fully damaged element (denoted by D=1).As illustrated in Figure5, the elastic-brittle damage constitutive relation of an element under uniaxial tensile and compressive stress is used to govern the failure of this mesoscopic element. When the stress in an element satisfies the strength criteria, the element begins to accumulate damage and the elastic modulus of this mesoscopic element degrades gradually as damage progresses. The elastic modulus of the damaged element is defined as(5)E=1−DE0,where D represents the damage variable, and E and E0 are the elastic modulus of the damaged and undamaged elements. The element and its damage are assumed to be isotropic and elastic. Therefore, as the stress increases, an increasing number of elements fail. Moreover, as the failed elements connect with one other, the whole failure process of the heterogeneous materials is reached.Figure 5
Elastic-brittle damage constitutive law of mesoscopic elements subjected to uniaxial compressive and tensile stress.When the tensile stress in the mesoscopic element reaches its tensile strengthft,(6)σ3≤−ft,and the damage variable D of the element under uniaxial tension can be expressed as(7)D=0,εt0≤ε,1−ftrE0ε,εtu<ε≤εt0,1,ε≤εtu,where ftr is the residual tensile strength; εt0 is the maximum tensile strain at the elastic limit; and εtu is the ultimate tensile strain of the element at which this element is damaged completely (Figure 5). An ultimate strain coefficient η is employed to describe the relationship of εt0 and εtu, such that εtu=ηεt0.As the damage variableD of the element changes in the damage process, the permeability coefficient of the mesoscopic element is expressed as(8)K=K0e−βσ3/3−αp,D=0,ξK0e−βσ3/3−αp,0<D<1,ξ′K0&e−βσ3/3−p,D=1,.Similarly, the damage variableD and permeability coefficient K of the mesoscopic element in the compressive state (Figure 5) in the damage process of geomaterials can be obtained on the basis of the Mohr–Coulomb strength criterion.When the shear stress in the element satisfies the Mohr–Coulomb failure criterion,(9)σ1−σ31+sinϕ1−sinϕ≥fc,and the damage variable D of the element under uniaxial compression can be expressed as(10)D=0,ε<εc0,1−fcrE0ε,εc0≤ε,where ϕ is the friction angle; fc is the compression strength; fcr is the residual compression strength; and εc0 is the maximum compression strain.As the damage variableD of the element changes in the damage process, the permeability coefficient of the mesoscopic element is expressed as(11)K=K0e−βσ1/3−αp,D=0,ξK0e−βσ1/3−αp,D>0.


## 2.1. Nonlinear Macroscopic Model Described by Linear Mesoscopic Elements

The RFPA code is based on finite element and statistical damage theory, and it can simulate the failure process of heterogeneous and permeable geomaterials. The heterogeneity of materials and the random distribution of defects of these materials are also considered by the RFPA code. The mesoscopic elements in the numerical model are assumed to be isotropic and homogeneous, and their mechanical properties (e.g., Young’s modulus, Poisson’s ratio, and strength properties, among others) are also assumed to be linear. The mesoscopic elements are statistically distributed (e.g., normal, Poisson, and Weibull distributions) to describe the mechanical properties of the nonlinear macromodel (Figure1). In RFPA, four-node isoparametric elements are used as the basic elements (Figure 2). In the simulation, the mesoscopic elements are defined as damaged if they satisfy the strength criterion. As stress increases, numerous elements fail. Moreover, when these elements connect with one another, the heterogeneous materials reach a failure process.Figure 1
Heterogeneous distribution of mechanical properties from linear mesoscopic mode to nonlinear macroscopic mode.Figure 2
Mesh of RFPA model.The statistical distribution densityfα of the mechanical parameter of the mesoscopic elements is defined by (1). The heterogeneity index m [12, 28] in the RFPA code is used to describe the heterogeneity of solid materials. A high m value indicates the presence of highly homogeneous materials, whereas a small m value denotes the existence of highly inhomogeneous materials (Figure 3). In the equation [29], α represents the mechanical parameters (e.g., Young’s modulus, Poisson’s ratio, tensional strength, and compressive strength) of the elements and α0 denotes the scale parameters related to the average values of the mechanical parameters.(1)fα=mα0αα0m−1exp−αα0m.Figure 3
Weibull’s distribution of the mechanical parameters of materials with different heterogeneitym value indices.


## 2.2. Primitive Phase Transition

The mesoscopic elements of an RFPA-based model are called primitives. Three primitive types (matrix, air, and contact) are considered for describing the deformation and failure processes of solid materials (Figure4). When the strain is located between the maximum tensile and compressive strains, the element is called a matrix primitive. A matrix primitive contains elastic and residual phases. The matrix primitive changes into air primitive when the deformation exceeds the maximum tensile strain, and it changes into contact primitive when the deformation exceeds the maximum compressive strain. Changes occur in the mechanical properties of elements during primitive transformation, and this change is called phase transition.Figure 4
Primitive phase transition.


## 2.3. Main Governing Equations of the RFPA-Flow Code

The geomaterial (rock) in the RFPA-Flow code is assumed to be fully saturated with fluid flow in accordance with Darcy’s law. In addition, the coupled process between stress/strain and fluid flow in the deforming rock mass is governed by Biot’s consolidation theory. The rock medium is an elastic-brittle material with residual strength, and the acoustic emission (AE) of rock also proves it. Thus, the mechanical behavior of rock material in the loading and unloading process conforms with elastic damage theory. The maximum tensile strength criterion and Mohr–Coulomb criterion are used as the damage judgment of the mesoscopic element in the RFPA code.The fundamental equations [30] of fluid-solid coupling in the RFPA code are as follows:(2)equilibrium  equation:∂σij∂xij+ρXj=0i,j=1,2,3,geometrical  equation:εij=12Ui,j+Uj,i,εv=ε11+ε22+ε33,constitutive  equation:σij′=σij−αpδij=λδijεv+2Gεij,seepage  equation:K∇2p=1Q∂p∂t−α∂εv∂t,where ρ is the volume density; σij is the total stress; εij is the total strain and εv is the volume strain; Ui is the displacement; p is pore pressure; λ is the Lame coefficient; δij is the Kronecher constant; G is the modulus of shear deformation; Q is Biot’s constant; K is the permeability coefficient; and α is the coefficient of pore pressure which is determined by the experiment.Under the same loading conditions, fluid flow behavior differs between the propagating fractures and the preexisting fractures. In the failure process of rock mass, the permeability of a rock changes all the time. Thus, the propagation of preexisting fractures, initiation of new fractures, coupling relationship of flow-stress-damage in these two kinds of fractures, and variations of rock permeability in the failure process are considered in the RFPA code.In the elastic state, the relationship of the stress and permeability coefficient is described by(3)Kf=K0e−bσ′,where K0 is the initial permeability coefficient; σ′ is the effective stress; and b is the coupling parameter.In describing the dependency of permeability on the stress and the damage, the following coupling equation [31] is introduced:(4)Kσ,p=ξK0e−βσii/3−αp,where α is the coefficient of pore pressure; ξξ≥1.0 is a mutation coefficient of permeability accounting for the increase in permeability once the element reaches the damage state; and β is the coupling parameter reflecting the influence of stress on the permeability coefficient. The value of the coefficients ξ, α, and β are determined experimentally, and they vary with the stress state and damage evolution of the rock. In the RFPA code, 0≤α<1 and ξ=1 are assumed for the mesoscopic element in the elastic stage (as indicated by D=0 in (7)). Once damage occurs, the permeability of the mesoscopic element increases significantly (i.e., see (8) and (11)). Furthermore, ξ=5 is assumed for the damaged element (denoted by 0<D<1), and α=1, ξ=100 is assumed for the fully damaged element (denoted by D=1).As illustrated in Figure5, the elastic-brittle damage constitutive relation of an element under uniaxial tensile and compressive stress is used to govern the failure of this mesoscopic element. When the stress in an element satisfies the strength criteria, the element begins to accumulate damage and the elastic modulus of this mesoscopic element degrades gradually as damage progresses. The elastic modulus of the damaged element is defined as(5)E=1−DE0,where D represents the damage variable, and E and E0 are the elastic modulus of the damaged and undamaged elements. The element and its damage are assumed to be isotropic and elastic. Therefore, as the stress increases, an increasing number of elements fail. Moreover, as the failed elements connect with one other, the whole failure process of the heterogeneous materials is reached.Figure 5
Elastic-brittle damage constitutive law of mesoscopic elements subjected to uniaxial compressive and tensile stress.When the tensile stress in the mesoscopic element reaches its tensile strengthft,(6)σ3≤−ft,and the damage variable D of the element under uniaxial tension can be expressed as(7)D=0,εt0≤ε,1−ftrE0ε,εtu<ε≤εt0,1,ε≤εtu,where ftr is the residual tensile strength; εt0 is the maximum tensile strain at the elastic limit; and εtu is the ultimate tensile strain of the element at which this element is damaged completely (Figure 5). An ultimate strain coefficient η is employed to describe the relationship of εt0 and εtu, such that εtu=ηεt0.As the damage variableD of the element changes in the damage process, the permeability coefficient of the mesoscopic element is expressed as(8)K=K0e−βσ3/3−αp,D=0,ξK0e−βσ3/3−αp,0<D<1,ξ′K0&e−βσ3/3−p,D=1,.Similarly, the damage variableD and permeability coefficient K of the mesoscopic element in the compressive state (Figure 5) in the damage process of geomaterials can be obtained on the basis of the Mohr–Coulomb strength criterion.When the shear stress in the element satisfies the Mohr–Coulomb failure criterion,(9)σ1−σ31+sinϕ1−sinϕ≥fc,and the damage variable D of the element under uniaxial compression can be expressed as(10)D=0,ε<εc0,1−fcrE0ε,εc0≤ε,where ϕ is the friction angle; fc is the compression strength; fcr is the residual compression strength; and εc0 is the maximum compression strain.As the damage variableD of the element changes in the damage process, the permeability coefficient of the mesoscopic element is expressed as(11)K=K0e−βσ1/3−αp,D=0,ξK0e−βσ1/3−αp,D>0.


## 3. Analysis of Hydraulic Fracture Propagation after Fracture Connection to Joint

Hydraulic fracture propagation, which changes path as it approaches a natural structural plane, is controlled by in situ stress, joint strength, and dip angle between fracture and joint, among others [32]. With the existence of a structural plane, hydraulic fractures either cease to propagate or change direction to find the path of least resistance and energy damage. As depicted by Figure 6(a), the normal and shear stresses of the joint under in situ stress on the basis of Mohr–Coulomb strength theory are(12)σn=σ1+σ32+σ1−σ32cosπ−2α,σn′=σ1+σ32−σ1−σ32cosπ−2α,τ=σ1−σ32sinπ−2α.Figure 6
Reaction of joint on hydraulic fracture.
(a)
(b)As the hydraulic fracture approaches the joint and connects with the initial stage of the fracture-joint intersection, the fluid is injected and flows into the joint. The pressure in the joint is(13)p=σ3+pnet,where pnet is the net pressure. Three types of hydraulic fracture propagations occur as fluid flows into the joint (Figure 6(b)).


### 3.1. Across the Joint

When the hydraulic fracture propagates across the joint and intersects with the joint, the pressure in the intersection point is(14)p>σn′+T,where T is the tensile strength of the rock matrix. Equations (12) and (13) are integrated into (14), and the net pressure in the fracture is obtained as(15)pnet>σ1−σ321+cos2α+T.


### 3.2. Along the Joint

In this situation, the hydraulic fracture connects with the joint and creates a fracture-joint overlap. Then, as the hydraulic fracture propagates along the joint, the fluid is injected. The pressure in the joint is(16)p>σn+Tb,where Tb is the tensile strength of the joint. Equations (12) and (13) are integrated into (16), and the net pressure in the fracture is obtained as(17)pnet>σ1−σ321−cos2α+Tb.After the fracture propagates along the joint and overlaps the joint plane, it propagates along these overlapped joints. Then, the fracture continues to propagate in the maximum principal stress direction in the joint tip.


### 3.3. Partly along the Joint and Partly across the Joint

In the propagation path along the joint plane (Figure6(b)), the hydraulic fracture across the joint is in a weak point. The pressure in this weak point is(18)p>σn′+Tw,where Tw is the tensile strength of the weak point. Equations (12) and (13) are integrated into (18), and the net pressure in the fracture is obtained as(19)pnet>σ1−σ321+cos2α+Tw.


## 3.1. Across the Joint

When the hydraulic fracture propagates across the joint and intersects with the joint, the pressure in the intersection point is(14)p>σn′+T,where T is the tensile strength of the rock matrix. Equations (12) and (13) are integrated into (14), and the net pressure in the fracture is obtained as(15)pnet>σ1−σ321+cos2α+T.


## 3.2. Along the Joint

In this situation, the hydraulic fracture connects with the joint and creates a fracture-joint overlap. Then, as the hydraulic fracture propagates along the joint, the fluid is injected. The pressure in the joint is(16)p>σn+Tb,where Tb is the tensile strength of the joint. Equations (12) and (13) are integrated into (16), and the net pressure in the fracture is obtained as(17)pnet>σ1−σ321−cos2α+Tb.After the fracture propagates along the joint and overlaps the joint plane, it propagates along these overlapped joints. Then, the fracture continues to propagate in the maximum principal stress direction in the joint tip.


## 3.3. Partly along the Joint and Partly across the Joint

In the propagation path along the joint plane (Figure6(b)), the hydraulic fracture across the joint is in a weak point. The pressure in this weak point is(18)p>σn′+Tw,where Tw is the tensile strength of the weak point. Equations (12) and (13) are integrated into (18), and the net pressure in the fracture is obtained as(19)pnet>σ1−σ321+cos2α+Tw.


## 4. Numerical Modeling

The simulation of hydraulic fracture propagation on jointed rocks is divided into three groups. The first group is the simulation of dip angle models, as shown in Figure7(a). The dip angle α is the angle between the joint plane and the maximum principal stress direction. The α values of the seven models are 0°, 15°, 30°, 45°, 60°, 75°, and 90° (fc = 20 MPa). The second group is the simulation of different joint strength models. In this group, the dip angle of the joint is 60°, and the joint strengths of the five models are 10, 20, 30, 40, and 80 MPa. As shown in Figures 7(b) and 7(c), the complex joint in the third group consists of two multijoint models and two masonry-shaped joint models after the maximum principal stress directions are changed. In the models, a wellbore with a radius of 32 mm is assumed to be located at the center of the 640 mm × 640 mm specimen discretized into 102,400 elements. The process of defining joints: first draw lines on the existing rock material model and then assign the joint material properties to the mesoscopic elements on the lines. The thickness of the joint is the size of one element, 2 mm. In the parallel joint models, the distance between the adjacent planes is 40 mm. In the multijoint models, the two dip angle joints and the two strength joints are positioned horizontally and vertically around the well. In the masonry-shaped joint models, the distances of the horizontally and vertically adjacent parallel joint planes are 40 and 80 mm, respectively. The confining pressures of σH= 15 MPa and σh= 10 MPa (in complex joint models: σH= 15 MPa, σh= 10 MPa; σH= 10 MPa, σh= 15 MPa) are imposed on the outward boundaries of the models. An increasing water pressure, with an initial value of 12 MPa followed by 0.1 MPa increments per step, is considered for the wellbore. The rock and joint parameters employed in the simulations are shown in Table 1.Figure 7
Schematic of the model set. (a) Parallel joint model. (b) Multijoint model. (c). Masonry-shaped joint model.
(a)
(b)
(c)Table 1
Physicomechanical parameters of rock and joint employed in simulations.
Parameters
Unit
Rock
Joint
Heterogeneity index (m)
—
2
2
Young’s modulus (E)
GPa
20
3
Possion’s ratio (ν)
—
0.25
0.25
Friction angle (φ)
°
37
37
Compressive strength (fc)
MPa
120
20 (10, 30, 40, 80)
Permeability coefficient (K)
md−1
0.000864
0.00864


## 5. Fracture Propagation of Parallel Joint Models with Different Dip Angles

As shown in Figure8, it depicts the hydraulic fractures propagation of different dip angle specimens. Hydraulic fracture initiation and propagation are significantly changed when the dip angle is increased.Figure 8
Minimum principal stress field of dip angle models. (a)α = 0°. (b) α = 15°. (c) α = 30°. (d) α = 45°. (e) α = 60°. (f) α = 75°. (g) α = 90°.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure8 shows that when the dip angle is set between 0° and 15°, the fracture initiated in the wellbore propagates in a straight path along the joint to form a well-connected pathway for fluid flow because of the small fracture toughness of the joint material. In the process of hydraulic fracturing, a tensile stress region (highlighted region) emerges around the propagating fracture tip, which causes the fracture to propagate along the weakened joint plane. In this situation, hydraulic fracture propagation is mainly controlled by the joint plane.When the dip angle is slightly increased (30° to 45°), the main fracture and some secondary ones propagate along the parallel joint planes in the tensile stress region to form multiple parallel straight fractures. Within this dip-angle range, the enlarged dip angle develops more secondary fractures. Fracture propagation is partially controlled by the joint plane.When the dip angle is increased to 60°, the fractures propagate along the joint plane. Moreover, a secondary fracture is induced in the maximum principal stress direction and connects with other secondary parallel fractures. When fractures further propagate, a complex fracture network is formed. In this situation, fracture propagation is controlled by the joint plane and the maximum principal stress.For large dip angles (75° to 90°), the hydraulic fracture ceases to propagate along the joint planes. In the propagation process, when the fracture meets discontinuities (joints) of strongly contrasting mechanical properties, it propagates across multiple joints and forms a jagged fracture with multiple branches in the maximum principal stress direction. In the numerical models, the mechanical properties of mesoscopic elements are presented as a Weibull distribution [33] owing to the heterogeneity of jointed rocks. Then, fracture initiation and propagation occur on the weak elements, that is, these weak elements fail under high stress when pressure increases, and they may connect with one another and with other numerous isolated failure elements to form macroscopic fractures [33, 34]. The termination of fracturing at the frequently contacting joint interfaces limits horizontal flow and instead produces complex flow paths. As a result, some secondary tensile fractures propagate along the joint near the terminated main fracture tip. Furthermore, the main fracture length is shorter than those depicted by previous simulations. In this situation, the fracture propagation is mainly controlled by the maximum principal stress.Figure9 presents the linear increase in breakdown pressure as the dip angle is increased. The smaller the dip angle, the easier the specimen is fractured. However, a small dip angle does not necessarily lead to good production operations in engineering, that is, a single fracture unrealistically portrays the efficient economic development of unconventional reservoirs, and thus, multiple fractures are needed for reservoir stimulation [35–37]. The simulation of this study suggests that the dip angles of 30° to 60° for parallel jointed rocks can efficiently propagate a main fracture and induce secondary ones to form a complex fracture network.Figure 9
Breakdown pressures of dip angle specimens.


## 6. Fracture Propagation in Joint Strength Models

To understand the hydraulic fracture propagation behavior of the joint strength specimens, a dip angle of 60° is considered on the basis of the first group simulation, that is, the hydraulic fracture propagation of the 60° dip model is simultaneously controlled by the joint plane and the maximum principal stress.Figure10 shows that as joint strength increases, the main fracture shifts from along-the-joint to across-the-joint propagation. The influence of joints on fracture evolution is weakened gradually. In this situation, the initiation and propagation of hydraulic fractures are controlled initially by the joint plane, followed by the maximum principal stress.Figure 10
Pore pressure field of joint strength models. (a)fc = 10 MPa. (b) fc = 20 MPa. (c) fc = 30 MPa. (d) fc = 40 MPa. (e) fc = 80 MPa.
(a)
(b)
(c)
(d)
(e)A joint strength of 10 MPa implies that the joint is much weaker than the rock material. In this situation, fracture propagation is mainly controlled by the weak joint plane. When the joint strength is set to 20 MPa and hydraulic pressure is increased, secondary fractures are induced along and across the joints, and they connect with one another upon intersection. In this situation, fracture evolution is controlled by the joint and the maximum principal stress. When the joint strength is increased to 30 MPa, the joint influence on fracture growth is significantly weakened, and the fracture propagates partly along the joint and partly across the joint to form a horizontally crooked fracture. When the joint strength is set to 40 MPa, the main fracture propagates horizontally and a jagged fracture with multiple branches is formed. A joint strength of 80 MPa is similar to that of a rock matrix. In this situation, the joint influence on fracture evolution is negligible, and the maximum principal stress mainly controls fracture propagation. When the fluid flows into the fracture, the fracture continues to propagate and later attains a relatively smooth interface.


## 7. Numerical Simulation of Complex Joint Models

The previous simulations in this study show that the joint plane and the maximum principal stress both control the propagation of hydraulic fractures in different situations. The multijoint model with different dip angles and strength joints (Figure7(b)), and the masonry-shaped joint model with horizontal and vertical joints (Figure 7(c)), are simulated to capture the hydraulic fracture propagation on jointed rocks.Figure11(a) presents the simulated hydraulic fracture initiation and propagation of the multijoint model in the horizontal direction under the confining pressures of σH= 15 MPa and σh= 10 MPa. The fracture deterministically selects the path of least resistance and energy damage through the rock, and it propagates to the proximal tip of the joint after the hydraulic fracture is initiated in the right side of the well. When fluid flows into the fracture, the fracture propagates along the small-dip (45°) joint. Then, the fracture turns horizontally under maximum principal stress condition in the distal tip of the joint. At the other side of the well, the fracture propagates across the large-dip (75°) joint and forms a jagged fracture with branches. When the confining pressure is changed to σH= 10 MPa and σh= 15 MPa (Figure 11(b)), the hydraulic fracture initiates and propagates vertically and spreads along the weak joint (60°, fc = 10 MPa) and across the high-strength joint (60°, fc = 80 MPa).Figure 11
Pore pressure field of multijoint models. (a)σH= 15 MPa, σh= 10 MPa. (b) σH= 10 MPa, σh= 15 MPa.
(a)
(b)The simulation of the multijoint models indicates that the maximum principal stress controls fracture propagation in global scale, whereas the joint plane controls fracture propagation in local scale. However, all these models consider only independent joints. Subsequently, masonry-shaped joint models are numerically simulated to capture the fracture propagation processes of interconnecting jointed rocks.As shown in Figure12(a), the hydraulic fracture propagates along the straight joint plane under the confining pressures of σH = 15 MPa and σh = 10 MPa. The breakdown pressure and the failure mode are similar to those of the dip angle model of 0° (Figure 8). When the confining pressure is changed to σH = 10 MPa and σh = 15 MPa (Figure 12(b)), the main fracture propagates vertically, while secondary hydraulic fractures initiate and propagate along the horizontal and vertical joints under tensile stress conditions (Figures 12(b) and 12(d)). When the hydraulic pressure is further increased, the propagating fractures connect with one another and merge to form a masonry-shaped fracture network.Figure 12
Pore pressure field (a, c) and minimum principal stress distribution (b, d) of masonry-shaped joint models:σH = 15 MPa, σh = 10 MPa (a, b) and σH = 10 MPa, σ = 15 MPa (c, d).
(a)
(b)
(c)
(d)


## 8. Comparison of Numerical and Experimental Results

To prove the reliability of the RFPA code in the numerical simulation of rocks during hydraulic fracturing, the numerical model of shale was also considered. The dimension of the shale model was 300 mm × 300 mm, and the radius of the well was 12 mm. The model was discretized into 90,000 elements. The confining pressure ofσH = 20 MPa and σh = 17 MPa was imposed on the outward boundaries of the model. Then, the simulation results were compared with experimental results [38] (Figure 13). Plane-strain calculation was used for the simulation. The shale and joint parameters employed in the simulation are shown in Table 2.Figure 13
Fracture propagation in the numerical and experimental [38] results of the shale model. (a) Numerical results. (b) Experimental results.
(a)
(b)Table 2
Physicomechanical parameters of shale and joint employed in simulation.
Parameters
Unit
Shale
Joint
Heterogeneity index (m)
—
2
2
Young’s modulus (E)
GPa
30
3
Possion’s ratio (ν)
—
0.25
0.25
Friction angle (φ)
°
35
35
Compressive strength (fc)
MPa
150
20
Permeability coefficient (K)
md−1
0.000864
0.00864Figure13 shows the main hydraulic fracture initiating and propagating perpendicularly to the joint under the control of the maximum principal stress σH. Given the weak cementation effect and small fracture toughness of the joint material, some secondary fractures propagate along the joint planes when the main fracture propagates. This finding suggests that hydraulic fracture propagation is controlled by the joint plane in local scale. However, because the intensity of fracture toughness is large in the perpendicular direction of the joint, hydraulic fracture propagation is prevented after it comes across multiple joints. Then, a jagged main fracture with turning and intersection points is formed. In the hydraulic fracturing process of the shale, a complex fracture network is formed when the hydraulic fractures and structural planes connected with one another.However, some differences are observed between the 2D numerical simulation and the 3D experiment. For instance, the number of fractured joints and the degree of fracturing are lower in the simulation compared with those in the experiment. Nonetheless, the initiation and propagation of fractures are essentially consistent between the simulated and experimental results.


## 9. Conclusions

In this study, the RFPA code is applied to simulate the hydraulic fracturing of heterogeneous jointed rocks. However, reality is often much more complex than simulations to which numerical models are applied. Nevertheless, the present study provides several interesting insights into fracture initiation and propagation mechanisms associated with hydraulic fracturing. A number of conclusions are derived from the numerical simulations.Fracture propagation is controlled by the maximum principal stress in global scale and the joint plane in local scale. In the maximum principal stress direction, hydraulic fractures propagate along small-dip joints (α ≤ 60° in this paper) or weak joints (fc≤ 20 MPa in this paper) and across large-dip joints (α > 60° in this paper) or strong joints (fc > 20 MPa in this paper).As depicted by the multijoint interconnection models, the main fracture propagates in the maximum principal stress direction. When multiple joints interconnect, secondary hydraulic fractures propagate along the joints within tensile stress regions without dip angle limitation. When hydraulic pressure is increased, the fractures connect with one another to form a complex fracture network based on the joint morphology.The intensity of fracture toughness is large in the perpendicular direction of the joint, and this prevents the propagation of hydraulic fractures. When hydraulic fractures frequently appear across the joints, an apparent decrease in fracture length is observed. During the main fracture propagation process, some joints appear with secondary fractures because of the weak cementation effect and limited fracture toughness of the joint material. As hydraulic fractures connect with structural planes, a complex fracture network can be formed.



---

*Source: 9461284-2018-07-29.xml*