# Modeling Simultaneous Multiple Fracturing Using the Combined Finite-Discrete Element Method

Authors: Quansheng Liu; Lei Sun; Pingli Liu; Lei Chen
Journal: Geofluids
 (2018)

---

## Abstract

Simultaneous multiple fracturing is a key technology to facilitate the production of shale oil/gas. When multiple hydraulic fractures propagate simultaneously, there is an interaction effect among these propagating hydraulic fractures, known as the stress-shadow effect, which has a significant impact on the fracture geometry. Understanding and controlling the propagation of simultaneous multiple hydraulic fractures and the interaction effects between multiple fractures are critical to optimizing oil/gas production. In this paper, the FDEM simulator and a fluid simulator are linked, named FDEM-Fluid, to handle hydromechanical-fracture coupling problems and investigate the simultaneous multiple hydraulic fracturing mechanism. The fractures propagation and the deformation of solid phase are solved by FDEM; meanwhile the fluid flow in the fractures is modeled using the principle of parallel-plate flow model. Several tests are carried out to validate the application of FDEM-Fluid in hydraulic fracturing simulation. Then, this FDEM-Fluid is used to investigate simultaneous multiple fractures treatment. Fractures repel each other when multiple fractures propagate from a single horizontal well, while the nearby fractures in different horizontal wells attract each other when multiple fractures propagate from multiple parallel horizontal wells. The in situ stress also has a significant impact on the fracture geometry.


---

## Body

## 1. Introduction

In the past two decades, a great number of researchers have been attracted to improving hydraulic fracturing techniques, which provide large surface area contact with formations and facilitate the production of oil and gas, as evident from the experience in North America [1]. Simultaneous multiple fracturing (simul-frac) is a key hydraulic fracturing technology that fractures two or more parallel horizontal wells simultaneously, which can significantly increase the volume and surface area of the flow path to improve the productivity and recovery efficiency [2–4]. A great number of attempts have been made to simultaneously fracture multiple adjacent horizontal wells to generate complex fracture networks. Although field practices have yielded a significantly improvement from simul-frac over stand-along wells [5], microseismic data [6], and numerical simulations [3, 7–12] also show a complex fracturing network produced by simul-frac, the reasons for simul-frac success are still not clear.Multiple hydraulic fractures treatment is a complex process, as it needs to consider not only the hydraulic fracturing process but also fracture interaction between multiple fractures. The hydraulic fracturing process is a typical hydro-mechanical-fracture (H-M-F) coupling problem, in which the following three basic processes require discussion [13]: (i) the rock deformation induced by the fluid pressure acting on the fracture surface, (ii) the fluid flow in the fracture, and (iii) the fracture growth. The significant consequences of multiple fractures interaction are stress-shadow effects that lead to stress field and fracture geometry changes.With the development of computers, numerical techniques have become reliable and convenient tools for studying hydraulic fracturing processes. The finite element method (FEM) [14] is a well-established method for rock engineering problems, and it has been utilized to simulate hydraulic fracture propagation for several decades [15]. However, meshes in the FEM and the fracture surfaces need to be coincident with each other. Therefore, advanced remeshing technologies have been developed for fracturing simulation [16, 17]. FEM with advanced remeshing technologies is still difficult to implement for field-scale reservoirs because of the meshing burden. The extended finite element method (XFEM) [18–20] and the generalized finite element method (GFEM) [21, 22] have been proposed to overcome the meshing burden associated with FEM and applied to model hydraulic fracturing. However, XFEM and GFEM still have difficulties in representing multiple complex fracture networks. The boundary element method (BEM) is also used for hydraulic fracturing problems [23–25]. However, it is difficult for BEM to address interfaces when they are in close proximity [26]. The discrete element method (DEM) [27–29] is well suited for the simulation of hydraulic fracture propagation in rock mass because it inherently incorporates fabric features. However, fracture propagation cannot be simulated directly but is often presented as a function of bond breakage between DEM particles. Furthermore, there is initially interspace between particles regardless of the particle arrangement. The combined finite-discrete element method (FDEM) [30] can well model the transition from continuous to discontinuous behavior in rock mass. FDEM inherits the advantages of FEM in describing elastic deformations and the capabilities of DEM in capturing discontinuities. FDEM is widely used for model collision, fracturing, and fragmentation where both discontinuum and continuum are involved [31, 32]. By coupling it with a fluid flow solver, FDEM has also been applied in handling hydromechanical (H-M) coupling problems [33–36].Although FDEM has been well-developed and applied by many researchers to study the process of hydraulic fracturing, very little literature can be found about the study of stress-shadow effect with FDEM. In this paper, FDEM coupled with a fluid algorithm, referred to as FDEM-Fluid, is developed to handle simultaneous multiple hydraulic fractures propagation. The organization of this paper is as follows. In Section2, the fundamental principles of FDEM are briefly presented. In Section 3, the fluid algorithm in fracture networks is presented. In Section 4, the coupling process and a flowchart of FDEM-Fluid are proposed. In Section 5, verification examples are presented, with comparisons with analytical or experimental results. In Section 6, simultaneous multiple fractures propagation in single and multiple horizontal wells under the stress-shadow effect are investigated. Finally, conclusions are drawn in Section 7.


## 2. FDEM for Solid Phase

### 2.1. Principles of FDEM

As shown in Figure1, the computational model is discretized with 3-node triangular elements and 4-node joint elements, which are inserted between the triangular elements. The stress state and deformation of the bulk material are captured by the linear-elastic triangular elements and the nonelastic joint elements. Fracture initiation and propagation within the continuum are simulated by the breakage of the joint elements on the basis of the combined single and smeared crack model [39] when the slip or opening distance is significantly large. The failure criterion will be explained in great detail in Section 2.4.Figure 1
Discretized intact body with 3-node triangular elements and 4-node joint elements.


### 2.2. Governing Equation

The movement of discrete bodies and the nodal coordinates are updated by the following governing equation:(1)M∂2x∂t2+C∂x∂t+Fintx-Fextx-Fcx=0,where M and C are the lumped mass and damping diagonal matrices of the system, respectively. x is the vector of the nodal coordinates. Fint, Fext, and Fc represent the internal resisting forces, the external loads, and the contact forces, as shown in Figure 2. The internal resisting forces include the contribution from the elastic reaction forces and the interelement bonding forces. The external nodal forces are applied boundary condition. To model quasistatic phenomena by dynamic relaxation or nonlinear material behavior, numerical dampingC is used in the governing equation to account for energy dissipation. The matrixC is calculated as(2)C=ξI,where ξ represents the constant damping coefficient andI is an identity matrix.Figure 2
The schematic diagram of the internal resisting forces, the external loads, and the contact forces (the contact force happens when two blocks contacted).


### 2.3. Contact Detection and Interaction

A highly efficient No Binary Search (NBS) algorithm [40] is used to detect all the contact couples that are in contact. As soon as contacting couples are detected, contact interaction is performed to calculate forces between discrete bodies by a potential function method. The contact function is given as follows:(3)φc=∫Γc12pert-rcTrt-rcdΓ,where pe is the penalty term, Γc is the boundary of the contactor element, and rt and rc are position vectors of the points falling on the target and contactor elements, respectively.A distributed contact force is generated according to the shape and size of the overlap between the contactor and target elements. As the contactor penetrates an areadA into the target (Figure 3), a contact force (df) is generated as follows:(4)df=gradφcPc-gradφtPtdA,where Pc and Pt are overlapping points of the contactor and target, respectively; φ is the corresponding potential function. The contact force only depends on the endpoints and is independent of the path taken.Figure 3
Contact force due to an infinitesimal overlap around pointsPc and Pt based on the potential function method ([24]).


### 2.4. Material Failure

Failure of rock material is simulated by the combined single and smeared crack model [39], which is developed from nonlinear elastic fracture mechanics (NLEFM) [41, 42]. According to this theory, when the strength of the material is reached, a fracture process zone (FPZ) characterized by nonlinear behavior starts to form ahead of the fracture tip (as shown in Figure 4(a)). The material in the FPZ, albeit damaged, is still able to transfer load across the fracture surfaces. In FDEM, fracture surfaces are assumed to coincide with the edges of the triangular elements, and the relative fracture displacements are evaluated by joint elements, as illustrated in Figure 4(b).Figure 4
(a) Softening part of stress-stain at the FPZ zone; (b) representation in FDEM by using triangular elements and joint elements.
(a)
(b)Based on the local stress and deformation field, fractures initiate and propagate according to three different failure modes, which are often denoted as Mode I (opening failure), Mode II (sliding failure), and Mixed Mode I-II. The constitutive behavior of the joint element is illustrated in detail in the following.(
1
) Mode I: as shown in Figure 5(a), before the fracture initiation, the stress at the fracture tip is linear with the fracture opening. Then, the fracture initiates when the value of the fracture opening, O, reaches a critical value, Op, which corresponds to the tensile strength ft. After the fracture propagates, the normal stress does not disappear immediately but decreases with increasing fracture opening. A stress-free surface (a real fracture) is generated when the opening reaches a certain residual value, Or.Figure 5
Constitutive behavior of the joint element. (a) Mode I. (b) Mode II. (c) Mixed Mode I-II.
(a)
(b)
(c)(
2
) Mode II: as shown in Figure 5(b), the tangential bonding stress between the two fracture surfaces is a function of the slip distance and the normal stress on the fracture surfaces. The critical slip (Sp) corresponds to the intrinsic shear strength (fs), which is calculated according to the Mohr-Coulomb criterion as follows:(5)fs=σn·tan⁡φ+c,where σn is the normal stress acting on the fracture surfaces; φ and c represent the intrinsic intact material friction angle and the internal cohesion, respectively.Upon undergoing critical slip,Sp, the tangential stress is gradually reduced to a residual value, fr, which corresponds to a pure friction resistance:(6)fr=σn·tan⁡φr,where φr is the residual friction angle.(
3
) Mixed Mode I-II is defined by a coupling criterion between Mode I and Mode II (see Figure 5(c)). This mode describes a combination of shear and tensile failure and is defined by the following:(7)O-OpOr-Op2+S-SpSr-Sp2≥1,where O is the opening distance and S is the slip distance.The values of the residual opening and slip depend on Mode I and Mode II fracture energy release rates,GI and GII, which are defined as follows:(8)GI=∫oporσodoGII=∫spsrτs-frds.


## 2.1. Principles of FDEM

As shown in Figure1, the computational model is discretized with 3-node triangular elements and 4-node joint elements, which are inserted between the triangular elements. The stress state and deformation of the bulk material are captured by the linear-elastic triangular elements and the nonelastic joint elements. Fracture initiation and propagation within the continuum are simulated by the breakage of the joint elements on the basis of the combined single and smeared crack model [39] when the slip or opening distance is significantly large. The failure criterion will be explained in great detail in Section 2.4.Figure 1
Discretized intact body with 3-node triangular elements and 4-node joint elements.


## 2.2. Governing Equation

The movement of discrete bodies and the nodal coordinates are updated by the following governing equation:(1)M∂2x∂t2+C∂x∂t+Fintx-Fextx-Fcx=0,where M and C are the lumped mass and damping diagonal matrices of the system, respectively. x is the vector of the nodal coordinates. Fint, Fext, and Fc represent the internal resisting forces, the external loads, and the contact forces, as shown in Figure 2. The internal resisting forces include the contribution from the elastic reaction forces and the interelement bonding forces. The external nodal forces are applied boundary condition. To model quasistatic phenomena by dynamic relaxation or nonlinear material behavior, numerical dampingC is used in the governing equation to account for energy dissipation. The matrixC is calculated as(2)C=ξI,where ξ represents the constant damping coefficient andI is an identity matrix.Figure 2
The schematic diagram of the internal resisting forces, the external loads, and the contact forces (the contact force happens when two blocks contacted).


## 2.3. Contact Detection and Interaction

A highly efficient No Binary Search (NBS) algorithm [40] is used to detect all the contact couples that are in contact. As soon as contacting couples are detected, contact interaction is performed to calculate forces between discrete bodies by a potential function method. The contact function is given as follows:(3)φc=∫Γc12pert-rcTrt-rcdΓ,where pe is the penalty term, Γc is the boundary of the contactor element, and rt and rc are position vectors of the points falling on the target and contactor elements, respectively.A distributed contact force is generated according to the shape and size of the overlap between the contactor and target elements. As the contactor penetrates an areadA into the target (Figure 3), a contact force (df) is generated as follows:(4)df=gradφcPc-gradφtPtdA,where Pc and Pt are overlapping points of the contactor and target, respectively; φ is the corresponding potential function. The contact force only depends on the endpoints and is independent of the path taken.Figure 3
Contact force due to an infinitesimal overlap around pointsPc and Pt based on the potential function method ([24]).


## 2.4. Material Failure

Failure of rock material is simulated by the combined single and smeared crack model [39], which is developed from nonlinear elastic fracture mechanics (NLEFM) [41, 42]. According to this theory, when the strength of the material is reached, a fracture process zone (FPZ) characterized by nonlinear behavior starts to form ahead of the fracture tip (as shown in Figure 4(a)). The material in the FPZ, albeit damaged, is still able to transfer load across the fracture surfaces. In FDEM, fracture surfaces are assumed to coincide with the edges of the triangular elements, and the relative fracture displacements are evaluated by joint elements, as illustrated in Figure 4(b).Figure 4
(a) Softening part of stress-stain at the FPZ zone; (b) representation in FDEM by using triangular elements and joint elements.
(a)
(b)Based on the local stress and deformation field, fractures initiate and propagate according to three different failure modes, which are often denoted as Mode I (opening failure), Mode II (sliding failure), and Mixed Mode I-II. The constitutive behavior of the joint element is illustrated in detail in the following.(
1
) Mode I: as shown in Figure 5(a), before the fracture initiation, the stress at the fracture tip is linear with the fracture opening. Then, the fracture initiates when the value of the fracture opening, O, reaches a critical value, Op, which corresponds to the tensile strength ft. After the fracture propagates, the normal stress does not disappear immediately but decreases with increasing fracture opening. A stress-free surface (a real fracture) is generated when the opening reaches a certain residual value, Or.Figure 5
Constitutive behavior of the joint element. (a) Mode I. (b) Mode II. (c) Mixed Mode I-II.
(a)
(b)
(c)(
2
) Mode II: as shown in Figure 5(b), the tangential bonding stress between the two fracture surfaces is a function of the slip distance and the normal stress on the fracture surfaces. The critical slip (Sp) corresponds to the intrinsic shear strength (fs), which is calculated according to the Mohr-Coulomb criterion as follows:(5)fs=σn·tan⁡φ+c,where σn is the normal stress acting on the fracture surfaces; φ and c represent the intrinsic intact material friction angle and the internal cohesion, respectively.Upon undergoing critical slip,Sp, the tangential stress is gradually reduced to a residual value, fr, which corresponds to a pure friction resistance:(6)fr=σn·tan⁡φr,where φr is the residual friction angle.(
3
) Mixed Mode I-II is defined by a coupling criterion between Mode I and Mode II (see Figure 5(c)). This mode describes a combination of shear and tensile failure and is defined by the following:(7)O-OpOr-Op2+S-SpSr-Sp2≥1,where O is the opening distance and S is the slip distance.The values of the residual opening and slip depend on Mode I and Mode II fracture energy release rates,GI and GII, which are defined as follows:(8)GI=∫oporσodoGII=∫spsrτs-frds.


## 3. Fluid Flow in Fracture Networks

As shown in Figure6, two sets of grid systems, namely, the solid grid and the hydraulic grid, are introduced to represent the solid phase and the fluid phase, respectively. The intact rock mass is discretized with solid grids by using 3-node triangular elements (in gray) and 4-node joint elements (in red). The fluid domain is discretized with hydraulic grids by using 2-node hydraulic elements (in blue), which are inserted into the broken joint elements.Figure 6
Fluid flow in the hydraulic fracture network and the influence domain of nodei.


### 3.1. Basic Assumptions of Fluid Flow

The hydraulic fracturing process is complicated, especially in regard to H-M-F coupling. Therefore, the following hypotheses are introduced:(1) the rock mass is impermeable and the fluid flows only in fractures. In the FDEM model, the flow only happens in the broken joint elements; (2) the fluid is assumed to be laminar, viscous, and incompressible; and (3) the fracture apertures can be obtained directly through the displacement of broken joint element surfaces.


### 3.2. Governing Equations

A cubic law based on parallel-plate flow model [43] is used to simulate the fluid flow in hydraulic elements, in which the flow rate can be expressed as (9)q=δ3J12μ,where δ is the equivalent aperture, J is the fluid pressure gradient, and μ is the viscosity of the fluid.


### 3.3. Update of Fluid Pressure and Saturation in the Fracture Network

Hydraulic elements connected with the hydraulic boundary are defined as the hydraulic fracture network, which is the only flow path for the fluid flow. Taking hydraulic nodesi as an example to illustrate the fluid flow in fracture network, the influence domain of node i is the union of all hydraulic elements connected to node i (see Figure 6). Hi-j is a hydraulic element connected to node i. pi and pj are fluid pressures at nodes i and j, respectively. The pressure gradient between i and j can be given by(10)Jj-i=pj-pi+ρgyj-yilj-i,where yj and yi are the vertical coordinates of hydraulic nodes j and i, respectively. lj-i is the length of hydraulic element Hi-j.The flow rate fromj to i can be calculated as(11)qj-i=δj-i3Jj-i12μ.Considering the fluid flow saturation, the actual flow rate fromj to i will decrease when the hydraulic node is unsaturated. Equation (11) can be modified as(12)qj-i=fsi+fsj2·δj-i3Jj-i12μ,where fs is a function of saturation, and fs=s2(3-2s); s is the saturation of the node. It should be noted that if s>1, s=1; if s<1, saturation should be updated as follows: (13)s=s0+QΔtV-ΔVVm,where s0 is the saturation of node i at the previous time step. V and V0 are the volumes of influence domain of this node at the current time step and previous time step, respectively. Then, ΔV=V-V0, andVm=V+V0/2. The volume of a hydraulic node is defined as half of the total volume of all hydraulic elements connected to the node.The total flow rate of nodei can be calculated as the summation of the flow rates in the influence domain of node i. The fluid pressure of hydraulic node i can be updated as follows [44]:(14)p=p0+KwQΔtV-KwΔVVm,where p0 is the pressure at the previous time step; Kw is the bulk modulus of the fluid; Q is the total flow rate; and Δt is the time increase.Note that there is a convergence criterion herein; that is, the hydraulic time step should be less than a critical value:(15)Δtf≤min⁡VKw∑kj,where V is the node volume and kj is the permeability factor of the jth flow path. For all nodes, the minimum value of Δtf is used in the algorithm.


## 3.1. Basic Assumptions of Fluid Flow

The hydraulic fracturing process is complicated, especially in regard to H-M-F coupling. Therefore, the following hypotheses are introduced:(1) the rock mass is impermeable and the fluid flows only in fractures. In the FDEM model, the flow only happens in the broken joint elements; (2) the fluid is assumed to be laminar, viscous, and incompressible; and (3) the fracture apertures can be obtained directly through the displacement of broken joint element surfaces.


## 3.2. Governing Equations

A cubic law based on parallel-plate flow model [43] is used to simulate the fluid flow in hydraulic elements, in which the flow rate can be expressed as (9)q=δ3J12μ,where δ is the equivalent aperture, J is the fluid pressure gradient, and μ is the viscosity of the fluid.


## 3.3. Update of Fluid Pressure and Saturation in the Fracture Network

Hydraulic elements connected with the hydraulic boundary are defined as the hydraulic fracture network, which is the only flow path for the fluid flow. Taking hydraulic nodesi as an example to illustrate the fluid flow in fracture network, the influence domain of node i is the union of all hydraulic elements connected to node i (see Figure 6). Hi-j is a hydraulic element connected to node i. pi and pj are fluid pressures at nodes i and j, respectively. The pressure gradient between i and j can be given by(10)Jj-i=pj-pi+ρgyj-yilj-i,where yj and yi are the vertical coordinates of hydraulic nodes j and i, respectively. lj-i is the length of hydraulic element Hi-j.The flow rate fromj to i can be calculated as(11)qj-i=δj-i3Jj-i12μ.Considering the fluid flow saturation, the actual flow rate fromj to i will decrease when the hydraulic node is unsaturated. Equation (11) can be modified as(12)qj-i=fsi+fsj2·δj-i3Jj-i12μ,where fs is a function of saturation, and fs=s2(3-2s); s is the saturation of the node. It should be noted that if s>1, s=1; if s<1, saturation should be updated as follows: (13)s=s0+QΔtV-ΔVVm,where s0 is the saturation of node i at the previous time step. V and V0 are the volumes of influence domain of this node at the current time step and previous time step, respectively. Then, ΔV=V-V0, andVm=V+V0/2. The volume of a hydraulic node is defined as half of the total volume of all hydraulic elements connected to the node.The total flow rate of nodei can be calculated as the summation of the flow rates in the influence domain of node i. The fluid pressure of hydraulic node i can be updated as follows [44]:(14)p=p0+KwQΔtV-KwΔVVm,where p0 is the pressure at the previous time step; Kw is the bulk modulus of the fluid; Q is the total flow rate; and Δt is the time increase.Note that there is a convergence criterion herein; that is, the hydraulic time step should be less than a critical value:(15)Δtf≤min⁡VKw∑kj,where V is the node volume and kj is the permeability factor of the jth flow path. For all nodes, the minimum value of Δtf is used in the algorithm.


## 4. Hydro-Mechanical-Fracture (H-M-F) Coupling Procedure

The hydro-mechanical-fracture (H-M-F) coupling process is illustrated in Figure7. Fluid flows in the fracture network, which is the only path of fluid flow. The fluid pressure applied on the fracture surfaces influences the fracture aperture, which in turn affects the fluid flow behavior. Furthermore, the fractures may propagate and connect with each other under the drive of fluid pressure, which improves the conductivity of the rock mass.Figure 7
Interaction in H-M-F coupling process.


### 4.1. Coupling Procedure of FDEM-Fluid

This FDEM-Fluid consists of the following two parts: the solid solver, which calculates the deformation of rock mass and fracture propagation, and the fluid solver, which calculates the fluid diffusion process and fluid pressure in fracture networks. The schematic diagram of the coupling between the two solvers is presented in Figure8.Figure 8
Interaction between the solid solver and fluid solver in FDEM-Fluid.A typical time step in FDEM-Fluid is illustrated in Figure9. First, nodal forces and elemental deformation are calculated based on the deformation of the triangular elements and joint elements. Second, fracture initiation and propagation are detected and marked, and, then, the hydraulic elements are updated when the joint elements are broken. Third, the current apertures of fractures are calculated. Fourth, the fluid pressure of the nodes is computed and applied on the triangular element as nodal forces. Then, a new loop starts with the updating of the nodal force.Figure 9
Flowchart of the calculation process in the FDEM-Fluid method.


### 4.2. Fluid Pressure Acting on the Fracture Surfaces

As shown in Figure10, the hydraulic pressure, simplified as linearly distributed loadings, is acting on the triangular element boundaries along the fracture surfaces. Then, the fluid pressure is converted to the nodal forces of the relevant triangular elements in the context of FDEM.Figure 10
Fluid pressure acting on element boundaries.p
i and pj are the hydraulic pressure on hydraulic nodes i and j, respectively. The total force on this surface can be calculated as(16)fn=pi+pj2lfs=pi-pj2δ,where fn and fs represent the normal force and shear force, respectively, l is the length of the fracture surface, and δ is the equivalent aperture.The normal force on each node is perpendicular to the edge and can be calculated as(17)fn1=fn4=pil3+pjl6fn2=fn3=pil6+pjl3.The shear force on each node is opposite to the fluid flow and can be calculated as(18)fs1=fs2=fs3=fs4=pi-pjδ4.


### 4.3. Equivalent Aperture

The fracture aperture can be obtained directly through the displacement of the broken joint element surfaces. Considering the H-M coupling, the fracture aperture is assumed to be a linear distribution along the fracture surface, as shown in Figure10. δi and δj represent the aperture of the two endpoints of the broken joint element. According to Iwai [45], the equivalent aperture of the fracture is as follows: (19)δ=δm16r21+r41/3,where δm=δi+δj/2, andr=δi/δj. It is obvious that the equivalent aperture is less thanδi+δj/2.


## 4.1. Coupling Procedure of FDEM-Fluid

This FDEM-Fluid consists of the following two parts: the solid solver, which calculates the deformation of rock mass and fracture propagation, and the fluid solver, which calculates the fluid diffusion process and fluid pressure in fracture networks. The schematic diagram of the coupling between the two solvers is presented in Figure8.Figure 8
Interaction between the solid solver and fluid solver in FDEM-Fluid.A typical time step in FDEM-Fluid is illustrated in Figure9. First, nodal forces and elemental deformation are calculated based on the deformation of the triangular elements and joint elements. Second, fracture initiation and propagation are detected and marked, and, then, the hydraulic elements are updated when the joint elements are broken. Third, the current apertures of fractures are calculated. Fourth, the fluid pressure of the nodes is computed and applied on the triangular element as nodal forces. Then, a new loop starts with the updating of the nodal force.Figure 9
Flowchart of the calculation process in the FDEM-Fluid method.


## 4.2. Fluid Pressure Acting on the Fracture Surfaces

As shown in Figure10, the hydraulic pressure, simplified as linearly distributed loadings, is acting on the triangular element boundaries along the fracture surfaces. Then, the fluid pressure is converted to the nodal forces of the relevant triangular elements in the context of FDEM.Figure 10
Fluid pressure acting on element boundaries.p
i and pj are the hydraulic pressure on hydraulic nodes i and j, respectively. The total force on this surface can be calculated as(16)fn=pi+pj2lfs=pi-pj2δ,where fn and fs represent the normal force and shear force, respectively, l is the length of the fracture surface, and δ is the equivalent aperture.The normal force on each node is perpendicular to the edge and can be calculated as(17)fn1=fn4=pil3+pjl6fn2=fn3=pil6+pjl3.The shear force on each node is opposite to the fluid flow and can be calculated as(18)fs1=fs2=fs3=fs4=pi-pjδ4.


## 4.3. Equivalent Aperture

The fracture aperture can be obtained directly through the displacement of the broken joint element surfaces. Considering the H-M coupling, the fracture aperture is assumed to be a linear distribution along the fracture surface, as shown in Figure10. δi and δj represent the aperture of the two endpoints of the broken joint element. According to Iwai [45], the equivalent aperture of the fracture is as follows: (19)δ=δm16r21+r41/3,where δm=δi+δj/2, andr=δi/δj. It is obvious that the equivalent aperture is less thanδi+δj/2.


## 5. Validation of FDEM-Fluid

Two numerical tests for the present FDEM-Fluid are carried out, and the results are compared with previous analytical and experimental results.


### 5.1. Validation 1: Unsteady Flow in a Single Fracture

This test investigates the fluid flow in a single fracture by FDEM-Fluid. The mesh sensitivity of fluid flow in FDEM-Fluid is also discussed. The geometry and boundary conditions of fluid flowing through a single fracture are shown in Figure11(a). The rock sample is 1 m in length and 0.2 m in height. The right end of the sample is an impermeable boundary, and a hydraulic pressure p0=10 MPa is applied to the left side. The aperture of the fracture a=3×10-5 m is constant during the simulation. The hydraulic parameters include the bulk modulus of water Kw=2.0 GPa and the viscosity of water μ=1×10-3 Pa·s.Figure 11
Rock sample with a single fracture.
(a)
Sample geometry
(b)
Three sets of mesh with mesh sizes of 0.1 m, 0.05 m, and 0.02 mThe hydraulic pressure distribution with time and location is investigated. The analytical solution of this example is as follows [46]:(20)pp0=1+4π∑n=0∞e-2n+12kwa2t/12μl2π2cos⁡2n+1π2ξ-1n+12n+1,whereξ=L-x/L; x is the distance to the left side; and p is the fluid pressure.Three sets of meshes, as shown in Figure11(b), with mesh sizes h = 0.1, 0.05, and 0.02 m, are utilized to investigate the accuracy of FDEM-Fluid and the mesh sensitivity of fluid flow.The hydraulic pressure distributions at different time with the three sets of meshes are shown in Figure12. Note that only hydraulic nodes have fluid pressure. Fluid flows from the left to right side, and the hydraulic pressure decreases from left to right. The pressure at the same location increases with the passing of time.Figure 12
Hydraulic pressure distribution at different time with three sets of mesh size (Pa).A comparison between the analytical and numerical solution is presented in Figure13. A nonlinear decrease of fluid pressure along the fracture can be observed. Fluid pressure increases with the passing of time. The numerical results of the three sets of meshes all agree well with the analytical solution; therefore, the accuracy and reliability of FDEM-Fluid for simulating fluid flow in fractures are validated.Figure 13
Comparison between analytical and numerical solutions of hydraulic pressure distribution with three sets of meshes.The mesh sensitivity of fluid flow is shown in Figure14. Comparing the error between the analytical solution and the numerical solution for the three sets of mesh sizes, it can be found that the error decreases when the mesh size decreases. In other words, the lower the mesh size is, the more accurate the numerical solution.Figure 14
Error between numerical solution and analytical solution under three sets of mesh sizes att = 0.6 ms and 2.4 ms.
(a)
Error analysis att = 0.6 ms
(b)
Error analysis att = 2.4 ms


### 5.2. Validation 2: Hydraulic Fracturing of a Cubic Sample with Preset Fractures

An experiment designed by Jiao et al. [37] is used to validate this FDEM-Fluid in simulating hydraulic fracturing. A cube plaster sample with preset fractures and a central injection hole is prepared, which is then fractured by injecting high-pressure water into the borehole with no other loads applied. In addition, Chen et al. conducted a similar physical experiment [38] where an anisotropic confining stress is applied.The geometry and cross section of the sample are shown in Figure15. Two preset fractures with a 30° inclination angle are symmetrically set along the diameter direction of the borehole. The length and depth of the preset fractures are 20 and 90 mm, respectively. The numerical mesh with 11502 triangular elements is shown in Figure 15(c).Figure 15
A cubic specimen with preset fractures.
(a)
Sample geometry
(b)
Sample cross section
(c)
Model meshTwo models are taken in to account: in model I, the confining stress is ignored; in model II, the compressive confining stressesσH=0.3 MPa and σh=0.05 MPa are applied in the y and x directions. A constant pressure, p0=0.8 MPa is applied on the wellbore surface uniformly. The mechanical parameters and hydraulic parameters are listed in Table 1.Table 1
Parameters of validation 2.
Parameters
Value
Time step,Δt (s)
1
×
e
-
5
Injection pressure,p0 (MPa)
0.8
Rock parameters
Bulk density,ρ (kg/m3)
1100
Young’s modulus,E (GPa)
4.0
Poisson’s ratio,ν
0.3
Joint element parameters
Tensile strength,ft (MPa)
0.3
Cohesion,c (MPa)
1.0
Friction angle,ψ (°)
32
Energy release rate,GI, GII (J/m2)
1.0, 10.0
Fluid parameters
Bulk modulus of fluid,Kw (GPa)
2.2
Viscosity,μ (Pa·s)
0.001
Initial value of aperture,δ0 (m)
5
e
-
5The fracture geometry evolution and fluid pressure distribution are shown in Figures16(a) and 17(a). The displacement magnitude contours of the two models are shown in Figures 16(b) and 17(b). Before fractures initiation, fluid flows into the preset fractures, and the fluid pressure increases gradually. With the passing of time, the prefractures start to propagate from the fracture tips, and fluid flows into the newly generated fracture. The plaster samples are eventually cut into two symmetrical blocks. Note that the prefractures start to propagate at step  3000 in model I and at step  4000 in model II. The plaster samples are cut into two blocks at step  12000 and step  15000 in model I and model II, respectively. This means that the fracture initiation and propagation in model II are harder than those in model I.Figure 16
(a) Fracture geometry and hydraulic pressure distribution of model I; (b) displacement magnitude contours.
(a)
(b)Figure 17
(a) Fracture geometry and hydraulic pressure distribution of model II; (b) displacement magnitude contours.
(a)
(b)Comparing the two fracture geometries, the fractures propagate strictly along the direction of preset fractures in model I. However, in model II, the fractures obviously turn towardy direction, which is parallel to the direction of the maximum principle stress.The numerical result of model I agree well with the experimental result conducted by Jiao et al. [37], as seen in Figure 18. The result of model II is also similar to the existing results (Figure 19), where the preset fractures gradually reorient to the direction of maximum principal stress which is consistent with the observation by Chen et al. [38], Almaguer et al. [47], and Wang [48].Figure 18
Experimental result of model I test [37].Figure 19
Experimental result of Chen et al. [38].
(a)
Top view
(b)
Sectional viewThese two examples validate the ability of FDEM-Fluid to model hydraulic fracture correctly with and without considering the influence of in situ stress.


## 5.1. Validation 1: Unsteady Flow in a Single Fracture

This test investigates the fluid flow in a single fracture by FDEM-Fluid. The mesh sensitivity of fluid flow in FDEM-Fluid is also discussed. The geometry and boundary conditions of fluid flowing through a single fracture are shown in Figure11(a). The rock sample is 1 m in length and 0.2 m in height. The right end of the sample is an impermeable boundary, and a hydraulic pressure p0=10 MPa is applied to the left side. The aperture of the fracture a=3×10-5 m is constant during the simulation. The hydraulic parameters include the bulk modulus of water Kw=2.0 GPa and the viscosity of water μ=1×10-3 Pa·s.Figure 11
Rock sample with a single fracture.
(a)
Sample geometry
(b)
Three sets of mesh with mesh sizes of 0.1 m, 0.05 m, and 0.02 mThe hydraulic pressure distribution with time and location is investigated. The analytical solution of this example is as follows [46]:(20)pp0=1+4π∑n=0∞e-2n+12kwa2t/12μl2π2cos⁡2n+1π2ξ-1n+12n+1,whereξ=L-x/L; x is the distance to the left side; and p is the fluid pressure.Three sets of meshes, as shown in Figure11(b), with mesh sizes h = 0.1, 0.05, and 0.02 m, are utilized to investigate the accuracy of FDEM-Fluid and the mesh sensitivity of fluid flow.The hydraulic pressure distributions at different time with the three sets of meshes are shown in Figure12. Note that only hydraulic nodes have fluid pressure. Fluid flows from the left to right side, and the hydraulic pressure decreases from left to right. The pressure at the same location increases with the passing of time.Figure 12
Hydraulic pressure distribution at different time with three sets of mesh size (Pa).A comparison between the analytical and numerical solution is presented in Figure13. A nonlinear decrease of fluid pressure along the fracture can be observed. Fluid pressure increases with the passing of time. The numerical results of the three sets of meshes all agree well with the analytical solution; therefore, the accuracy and reliability of FDEM-Fluid for simulating fluid flow in fractures are validated.Figure 13
Comparison between analytical and numerical solutions of hydraulic pressure distribution with three sets of meshes.The mesh sensitivity of fluid flow is shown in Figure14. Comparing the error between the analytical solution and the numerical solution for the three sets of mesh sizes, it can be found that the error decreases when the mesh size decreases. In other words, the lower the mesh size is, the more accurate the numerical solution.Figure 14
Error between numerical solution and analytical solution under three sets of mesh sizes att = 0.6 ms and 2.4 ms.
(a)
Error analysis att = 0.6 ms
(b)
Error analysis att = 2.4 ms


## 5.2. Validation 2: Hydraulic Fracturing of a Cubic Sample with Preset Fractures

An experiment designed by Jiao et al. [37] is used to validate this FDEM-Fluid in simulating hydraulic fracturing. A cube plaster sample with preset fractures and a central injection hole is prepared, which is then fractured by injecting high-pressure water into the borehole with no other loads applied. In addition, Chen et al. conducted a similar physical experiment [38] where an anisotropic confining stress is applied.The geometry and cross section of the sample are shown in Figure15. Two preset fractures with a 30° inclination angle are symmetrically set along the diameter direction of the borehole. The length and depth of the preset fractures are 20 and 90 mm, respectively. The numerical mesh with 11502 triangular elements is shown in Figure 15(c).Figure 15
A cubic specimen with preset fractures.
(a)
Sample geometry
(b)
Sample cross section
(c)
Model meshTwo models are taken in to account: in model I, the confining stress is ignored; in model II, the compressive confining stressesσH=0.3 MPa and σh=0.05 MPa are applied in the y and x directions. A constant pressure, p0=0.8 MPa is applied on the wellbore surface uniformly. The mechanical parameters and hydraulic parameters are listed in Table 1.Table 1
Parameters of validation 2.
Parameters
Value
Time step,Δt (s)
1
×
e
-
5
Injection pressure,p0 (MPa)
0.8
Rock parameters
Bulk density,ρ (kg/m3)
1100
Young’s modulus,E (GPa)
4.0
Poisson’s ratio,ν
0.3
Joint element parameters
Tensile strength,ft (MPa)
0.3
Cohesion,c (MPa)
1.0
Friction angle,ψ (°)
32
Energy release rate,GI, GII (J/m2)
1.0, 10.0
Fluid parameters
Bulk modulus of fluid,Kw (GPa)
2.2
Viscosity,μ (Pa·s)
0.001
Initial value of aperture,δ0 (m)
5
e
-
5The fracture geometry evolution and fluid pressure distribution are shown in Figures16(a) and 17(a). The displacement magnitude contours of the two models are shown in Figures 16(b) and 17(b). Before fractures initiation, fluid flows into the preset fractures, and the fluid pressure increases gradually. With the passing of time, the prefractures start to propagate from the fracture tips, and fluid flows into the newly generated fracture. The plaster samples are eventually cut into two symmetrical blocks. Note that the prefractures start to propagate at step  3000 in model I and at step  4000 in model II. The plaster samples are cut into two blocks at step  12000 and step  15000 in model I and model II, respectively. This means that the fracture initiation and propagation in model II are harder than those in model I.Figure 16
(a) Fracture geometry and hydraulic pressure distribution of model I; (b) displacement magnitude contours.
(a)
(b)Figure 17
(a) Fracture geometry and hydraulic pressure distribution of model II; (b) displacement magnitude contours.
(a)
(b)Comparing the two fracture geometries, the fractures propagate strictly along the direction of preset fractures in model I. However, in model II, the fractures obviously turn towardy direction, which is parallel to the direction of the maximum principle stress.The numerical result of model I agree well with the experimental result conducted by Jiao et al. [37], as seen in Figure 18. The result of model II is also similar to the existing results (Figure 19), where the preset fractures gradually reorient to the direction of maximum principal stress which is consistent with the observation by Chen et al. [38], Almaguer et al. [47], and Wang [48].Figure 18
Experimental result of model I test [37].Figure 19
Experimental result of Chen et al. [38].
(a)
Top view
(b)
Sectional viewThese two examples validate the ability of FDEM-Fluid to model hydraulic fracture correctly with and without considering the influence of in situ stress.


## 6. Simultaneous Multifracture Treatments

Simultaneous multiple hydraulic fracturing in combination with horizontal drilling has been the key technology in the exploitation of unconventional gas/oil. This technology can significantly increase the volume and surface area of the flow path, so as to improve the productivity and recovery efficiency. For simultaneous multiple-fracture propagation, fracture interaction (i.e., the stress-shadow effect) plays an important role in fracture geometry. In this section, two examples are studied to investigate fracture interactions during simultaneous fracturing processes.


### 6.1. Simultaneous Multifracture Propagation from a Horizontal Well

In this section, a single horizontal wellbore with two parallel initial fractures is simulated and compared to the results of Wu et al. [10]. As shown in Figure 20, two preset fractures are orthogonal to the wellbore axes, with the offset between the two parallel preset fractures being 10 m. A constant injection rate Q=0.106 m3/s is applied at the middle of the horizontal wellbore. The numerical model with 123466 unstructured meshes is shown in Figure 20(b). The mechanical parameters and hydraulic parameters are listed in Table 2 [10].Table 2
Parameters for simultaneous multifracture propagation from a horizontal well.
Parameters
Value
Time step,Δt (s)
1
×
1
0
-
5
Injection rate,Q (m3/s)
0.106
Stress anisotropy,σy-σx (MPa)
0.9
Rock parameters
Young’s modulus,E (GPa)
30
Poisson’s ratio,ν
0.35
Joint element parameters
Tensile strength,ft (MPa)
0.3
Cohesion,c (MPa)
1.0
Friction angle,ψ (°)
32
Energy release rate,GI, GII (J/m2)
1.0, 10.0
Fluid parameters
Viscosity,μ (Pa·s)
0.001
Bulk modulus of fluid,Kw (GPa)
2.2Figure 20
A single horizontal wellbore with two parallel preset fractures.
(a)
(b)This model is simulated under isotropic and anisotropic in situ stress conditions. The differential stress for the anisotropic stress is 0.9 MPa, with the maximum horizontal stress oriented in they direction.Figures21 and 22 show the fracture geometry and hydraulic pressure distribution of fracturing in the two models. The fluid first flows into the preset fractures from the borehole, and fluid pressure is applied on the fracture surfaces. After step  4000, the two fractures start to propagate from the fracture tips, driven by the fluid pressure. Then, fluid flows into the newly generated fractures.Figure 21
Propagation of two initially parallel fractures driven by hydraulic pressure: isotropic in situ stress field.
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000Figure 22
Propagation of two initially parallel fractures driven by hydraulic pressure: anisotropic in situ stress field.
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000As shown in Figures21 and 22, in both cases, two preset fractures propagate away from each other with a side-by-side configuration because of the mechanical fracture interaction (i.e., stress-shadow effect). Compared to the isotropic in situ stress condition, the fractures under the anisotropic in situ condition have smaller curvature because of the confining stress (see Figure 22).As shown in Figure23, the fractures under the anisotropic in situ stress condition propagate longer in the y direction than those under the isotropic in situ stress condition, while the influence domain in the x-direction is smaller. The complex influence domain is critical to the exploitation of unconventional gas/oil because it influences the surface area contact with the formation or fracture network.Figure 23
Influence domain of fractures under isotropic and anisotropic in situ stress conditions.
(a)
Influence domain of fractures in they-direction
(b)
Influence domain of fractures in thex-directionThe reason fractures turn is because of the mechanical fracture interaction (i.e., the stress-shadow effect), which exerts additional stress on the surrounding rock, resulting in a local change in direction of the maximum horizontal stress and a deviation of the fracture path.As shown in Figures24 and 25, there is a domain surrounding the two adjacent hydraulic fractures, where the stress-shadow effect obviously changes the stress field around the fracture tip and the direction of fracture propagation. Under the anisotropic in situ stress condition, the anisotropic in situ stress can reduce the influence of the stress-shadow effect. In other words, anisotropic far-field stresses tend to suppress fracture curving [49], which indicate that the stress shadow and anisotropic stress compete with each other.Figure 24
Maximum principal stress contours under an isotropic in situ stress field (positive denotes tension).
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000Figure 25
Maximum principal stress contours under an anisotropic in situ stress field (positive denotes tension).
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000As shown in Figure26, the FDEM-Fluid results agree well with previous numerical results [10] in both isotropic and anisotropic in situ stress conditions. Therefore, this FDEM-Fluid method is able to simulate the multiple hydraulic fracturing, with considering the H-M-F coupling, the stress-shadow effect, and in situ stress.Figure 26
Comparison of propagation paths by FDEM-Fluid and Wu et al. [10].
(a)
Under isotropic in situ stress field condition
(b)
Under anisotropic in situ stress field condition


### 6.2. Simultaneous Multifracture Propagation from Multiple Horizontal Wells

In this example, simultaneous propagation of two sets of initial fractures from two horizontal wells is investigated. As shown in Figure27, the horizontal wells are 4 m apart, and two preset fractures are orthogonal to each wellbore. The offset between the two parallel perforations in the same wellbore is 4 m, and the offset between two fractures in different horizontal wells is 1.5 m. The fluid pressure in the horizontal wellbores is constant at 10 MPa, and the fluid viscosity is taken as 0.001 Pa·s. The properties of the rock mass are the same as in Table 2.Figure 27
Two horizontal wells with two sets of parallel preset fractures.The evolution of the preset fracture geometry and the distribution of fluid pressure at different time steps are shown in Figure28. The fluid flows into the preset fractures. Then, the fractures start to propagate, driven by the fluid pressure. It can be observed that fractures f1 and f3 which propagate from each horizontal well turn toward each other because the induced stress-shadow effect changes the maximum principle stress at the fracture tips (see Figure 28(c)). The phenomenon that hydraulic fractures from different horizontal wells attract each other is consistent with the observation by Wang [4], Kumar and Ghassemi [50], and Wu [51]. Compared to fracture f3, fracture f1 has a larger deviation because of the stress-shadow effect induced by both fracture f2 and f3. However, f1 and f2 propagate away from each other.Figure 28
Evolution of fracture geometries and the distribution of fluid pressure.
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000Figure29 shows the distribution of the maximum principle stress and displacement magnitude. Figure 30 shows the direction of the maximum principal stress map (the green line segment denotes the direction) at steps 4000 and 8000. As the fracture is propagating, it alters the local in situ stresses around it, and the direction of the maximum principle stress can be reversed close to the fracture surface. When the hydraulic fracture propagates into the “stress alteration zone” that induced by the other fracture, the fracture reorients itself to propagate along the new direction of the local maximum principal stress.Figure 29
Counter of (a) maximum principle stress and (b) displacement magnitude.
(a)
(b)Figure 30
Direction of maximum principal stress at steps 4000 and 8000.
(a)
(b)


## 6.1. Simultaneous Multifracture Propagation from a Horizontal Well

In this section, a single horizontal wellbore with two parallel initial fractures is simulated and compared to the results of Wu et al. [10]. As shown in Figure 20, two preset fractures are orthogonal to the wellbore axes, with the offset between the two parallel preset fractures being 10 m. A constant injection rate Q=0.106 m3/s is applied at the middle of the horizontal wellbore. The numerical model with 123466 unstructured meshes is shown in Figure 20(b). The mechanical parameters and hydraulic parameters are listed in Table 2 [10].Table 2
Parameters for simultaneous multifracture propagation from a horizontal well.
Parameters
Value
Time step,Δt (s)
1
×
1
0
-
5
Injection rate,Q (m3/s)
0.106
Stress anisotropy,σy-σx (MPa)
0.9
Rock parameters
Young’s modulus,E (GPa)
30
Poisson’s ratio,ν
0.35
Joint element parameters
Tensile strength,ft (MPa)
0.3
Cohesion,c (MPa)
1.0
Friction angle,ψ (°)
32
Energy release rate,GI, GII (J/m2)
1.0, 10.0
Fluid parameters
Viscosity,μ (Pa·s)
0.001
Bulk modulus of fluid,Kw (GPa)
2.2Figure 20
A single horizontal wellbore with two parallel preset fractures.
(a)
(b)This model is simulated under isotropic and anisotropic in situ stress conditions. The differential stress for the anisotropic stress is 0.9 MPa, with the maximum horizontal stress oriented in they direction.Figures21 and 22 show the fracture geometry and hydraulic pressure distribution of fracturing in the two models. The fluid first flows into the preset fractures from the borehole, and fluid pressure is applied on the fracture surfaces. After step  4000, the two fractures start to propagate from the fracture tips, driven by the fluid pressure. Then, fluid flows into the newly generated fractures.Figure 21
Propagation of two initially parallel fractures driven by hydraulic pressure: isotropic in situ stress field.
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000Figure 22
Propagation of two initially parallel fractures driven by hydraulic pressure: anisotropic in situ stress field.
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000As shown in Figures21 and 22, in both cases, two preset fractures propagate away from each other with a side-by-side configuration because of the mechanical fracture interaction (i.e., stress-shadow effect). Compared to the isotropic in situ stress condition, the fractures under the anisotropic in situ condition have smaller curvature because of the confining stress (see Figure 22).As shown in Figure23, the fractures under the anisotropic in situ stress condition propagate longer in the y direction than those under the isotropic in situ stress condition, while the influence domain in the x-direction is smaller. The complex influence domain is critical to the exploitation of unconventional gas/oil because it influences the surface area contact with the formation or fracture network.Figure 23
Influence domain of fractures under isotropic and anisotropic in situ stress conditions.
(a)
Influence domain of fractures in they-direction
(b)
Influence domain of fractures in thex-directionThe reason fractures turn is because of the mechanical fracture interaction (i.e., the stress-shadow effect), which exerts additional stress on the surrounding rock, resulting in a local change in direction of the maximum horizontal stress and a deviation of the fracture path.As shown in Figures24 and 25, there is a domain surrounding the two adjacent hydraulic fractures, where the stress-shadow effect obviously changes the stress field around the fracture tip and the direction of fracture propagation. Under the anisotropic in situ stress condition, the anisotropic in situ stress can reduce the influence of the stress-shadow effect. In other words, anisotropic far-field stresses tend to suppress fracture curving [49], which indicate that the stress shadow and anisotropic stress compete with each other.Figure 24
Maximum principal stress contours under an isotropic in situ stress field (positive denotes tension).
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000Figure 25
Maximum principal stress contours under an anisotropic in situ stress field (positive denotes tension).
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000As shown in Figure26, the FDEM-Fluid results agree well with previous numerical results [10] in both isotropic and anisotropic in situ stress conditions. Therefore, this FDEM-Fluid method is able to simulate the multiple hydraulic fracturing, with considering the H-M-F coupling, the stress-shadow effect, and in situ stress.Figure 26
Comparison of propagation paths by FDEM-Fluid and Wu et al. [10].
(a)
Under isotropic in situ stress field condition
(b)
Under anisotropic in situ stress field condition


## 6.2. Simultaneous Multifracture Propagation from Multiple Horizontal Wells

In this example, simultaneous propagation of two sets of initial fractures from two horizontal wells is investigated. As shown in Figure27, the horizontal wells are 4 m apart, and two preset fractures are orthogonal to each wellbore. The offset between the two parallel perforations in the same wellbore is 4 m, and the offset between two fractures in different horizontal wells is 1.5 m. The fluid pressure in the horizontal wellbores is constant at 10 MPa, and the fluid viscosity is taken as 0.001 Pa·s. The properties of the rock mass are the same as in Table 2.Figure 27
Two horizontal wells with two sets of parallel preset fractures.The evolution of the preset fracture geometry and the distribution of fluid pressure at different time steps are shown in Figure28. The fluid flows into the preset fractures. Then, the fractures start to propagate, driven by the fluid pressure. It can be observed that fractures f1 and f3 which propagate from each horizontal well turn toward each other because the induced stress-shadow effect changes the maximum principle stress at the fracture tips (see Figure 28(c)). The phenomenon that hydraulic fractures from different horizontal wells attract each other is consistent with the observation by Wang [4], Kumar and Ghassemi [50], and Wu [51]. Compared to fracture f3, fracture f1 has a larger deviation because of the stress-shadow effect induced by both fracture f2 and f3. However, f1 and f2 propagate away from each other.Figure 28
Evolution of fracture geometries and the distribution of fluid pressure.
(a)
Step  4000
(b)
Step  8000
(c)
Step  12000
(d)
Step  16000Figure29 shows the distribution of the maximum principle stress and displacement magnitude. Figure 30 shows the direction of the maximum principal stress map (the green line segment denotes the direction) at steps 4000 and 8000. As the fracture is propagating, it alters the local in situ stresses around it, and the direction of the maximum principle stress can be reversed close to the fracture surface. When the hydraulic fracture propagates into the “stress alteration zone” that induced by the other fracture, the fracture reorients itself to propagate along the new direction of the local maximum principal stress.Figure 29
Counter of (a) maximum principle stress and (b) displacement magnitude.
(a)
(b)Figure 30
Direction of maximum principal stress at steps 4000 and 8000.
(a)
(b)


## 7. Conclusions

The interaction (i.e., the stress-shadow effect) during simultaneous multiple-fracture propagation plays an important role in fracture geometry, which can significantly affect the productivity and recovery efficiency. Therefore, understanding and controlling the interaction effects between multiple fractures are critical. In this study, a hydraulic fracturing simulator, named FDEM-Fluid, based on FDEM and the parallel-plate flow model is presented to capture the development of complex hydraulic fracturing propagation. This FDEM-Fluid is verified by comparing it with existing analytical and experimental results, and good agreements are observed. Then, the multifracture treatment in the hydraulic fracturing process is investigated with FDEM-Fluid. The main conclusions from the study include the following:(1)
When multiple fractures are orthogonal to the horizontal well and propagate from a single horizontal well simultaneously, fractures repel each other as a result of the stress-shadow effect, resulting in a complex fracture geometry. Furthermore, the anisotropic in situ stress, in which the maximum horizontal stress is orthogonal to the horizontal well, has a great influence on the curvature and length of the fracture geometry because of the competition between the stress-shadow effect and the anisotropic stress.(2)
When multiple fractures are orthogonal to the horizontal wells and propagate from multiple parallel horizontal wells simultaneously, the nearby fractures in different horizontal wells attract each other as a result of the stress-shadow effect, which alters the local in situ stresses and the direction of the maximum principle stress. The fracture tips in these fractures tend to reorient themselves to be perpendicular to the main fractures. This leads to the development of a more complex fracture geometry and prevent fractures from penetrating deeper into the opposite wells.



---

*Source: 4252904-2018-02-18.xml*