Geofluids
Volume 2018 (2018), Article ID 4252904, 20 pages
https://doi.org/10.1155/2018/4252904

Research Article

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

Quansheng Liu,¹ Lei Sun,¹ Pingli Liu,² and Lei Chen³

¹School of Civil Engineering, Wuhan University, Wuhan 430072, China
²State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu 610500, China
³Department of Mechanical Engineering, Mississippi State University, Mississippi State, MS 39762, USA

Correspondence should be addressed to Pingli Liu; liupingli@swpu.edu.cn

Received 23 August 2017; Accepted 8 January 2018; Published 18 February 2018

Academic Editor: Wei Wu

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

Abstract.

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.

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 Section 2, 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 Figure 1, 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:

where and are the lumped mass and damping diagonal matrices of the system, respectively. is the vector of the nodal coordinates. , , and 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 damping C is used in the governing equation to account for energy dissipation. The matrix C is calculated as

where represents the constant damping coefficient and I 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:

where is the penalty term, is the boundary of the contactor element, and and 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 area into the target (Figure 3), a contact force () is generated as follows:

where and 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 points and 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).

(a)

(b)

Figure 4: (a) Softening part of stress-stain at the FPZ zone; (b) representation in FDEM by using triangular elements and joint elements.

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.

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, , reaches a critical value, , which corresponds to the tensile strength . 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, .

(a)

(b)

(c)

Figure 5: Constitutive behavior of the joint element. (a) Mode I. (b) Mode II. (c) Mixed Mode I-II.

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 () corresponds to the intrinsic shear strength (), which is calculated according to the Mohr-Coulomb criterion as follows:

where is the normal stress acting on the fracture surfaces; and represent the intrinsic intact material friction angle and the internal cohesion, respectively.

Upon undergoing critical slip, , the tangential stress is gradually reduced to a residual value, , which corresponds to a pure friction resistance:

where is the residual friction angle.

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:

where is the opening distance and is the slip distance.

The values of the residual opening and slip depend on Mode I and Mode II fracture energy release rates, and , which are defined as follows:

3. Fluid Flow in Fracture Networks

As shown in Figure 6, 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 node .

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: 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; the fluid is assumed to be laminar, viscous, and incompressible; and ) 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

where is the equivalent aperture, 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 nodes as an example to illustrate the fluid flow in fracture network, the influence domain of node is the union of all hydraulic elements connected to node (see Figure 6). is a hydraulic element connected to node . and are fluid pressures at nodes and , respectively. The pressure gradient between and can be given by

where and are the vertical coordinates of hydraulic nodes and , respectively. is the length of hydraulic element .

The flow rate from to can be calculated as

Considering the fluid flow saturation, the actual flow rate from to will decrease when the hydraulic node is unsaturated. Equation (11) can be modified as

where is a function of saturation, and ; is the saturation of the node. It should be noted that if , ; if , saturation should be updated as follows:

where is the saturation of node at the previous time step. and are the volumes of influence domain of this node at the current time step and previous time step, respectively. Then, , and. 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 node can be calculated as the summation of the flow rates in the influence domain of node . The fluid pressure of hydraulic node can be updated as follows [44]:

where is the pressure at the previous time step; is the bulk modulus of the fluid; is the total flow rate; and 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:

where is the node volume and is the permeability factor of the th flow path. For all nodes, the minimum value of 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 Figure 7. 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 Figure 8.

Figure 8: Interaction between the solid solver and fluid solver in FDEM-Fluid.

A typical time step in FDEM-Fluid is illustrated in Figure 9. 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 Figure 10, 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.

and are the hydraulic pressure on hydraulic nodes and , respectively. The total force on this surface can be calculated as

where and represent the normal force and shear force, respectively, 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

The shear force on each node is opposite to the fluid flow and can be calculated as

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 Figure 10. and 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:

where , and. It is obvious that the equivalent aperture is less than.

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 Figure 11(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  MPa is applied to the left side. The aperture of the fracture  m is constant during the simulation. The hydraulic parameters include the bulk modulus of water  GPa and the viscosity of water  Pa·s.

(a) Sample geometry

(b) Three sets of mesh with mesh sizes of 0.1 m, 0.05 m, and 0.02 m

Figure 11: Rock sample with a single fracture.

The hydraulic pressure distribution with time and location is investigated. The analytical solution of this example is as follows [46]:

where; is the distance to the left side; and is the fluid pressure.

Three sets of meshes, as shown in Figure 11(b), with mesh sizes = 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 Figure 12. 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 Figure 13. 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 Figure 14. 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.

(a) Error analysis at = 0.6 ms

(b) Error analysis at = 2.4 ms

Figure 14: Error between numerical solution and analytical solution under three sets of mesh sizes at = 0.6 ms and 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 Figure 15. 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).

(a) Sample geometry

(b) Sample cross section

(c) Model mesh

Figure 15: A cubic specimen with preset fractures.

Two models are taken in to account: in model I, the confining stress is ignored; in model II, the compressive confining stresses  MPa and  MPa are applied in the and directions. A constant pressure,  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, (s)

Injection pressure (MPa)
0.8

Rock parameters

    Bulk density, (kg/m³)
1100

    Young’s modulus, (GPa)
4.0

    Poisson’s ratio,
0.3

Joint element parameters

    Tensile strength, (MPa)
0.3

    Cohesion, (MPa)
1.0

    Friction angle, (°)
32

    Energy release rate, , (J/m²)
1.0, 10.0

Fluid parameters

    Bulk modulus of fluid, (GPa)
2.2

    Viscosity, μ (Pa·s)
0.001

    Initial value of aperture, (m)

The fracture geometry evolution and fluid pressure distribution are shown in Figures 16(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.

(a)

(b)

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.

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 toward 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].

(a) Top view

(b) Sectional view

Figure 19: Experimental result of Chen et al. [38].

These 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  m³/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, (s)

Injection rate, (m³/s)
0.106

Stress anisotropy, (MPa)
0.9

Rock parameters

    Young’s modulus, (GPa)
30

    Poisson’s ratio,
0.35

Joint element parameters

    Tensile strength, (MPa)
0.3

    Cohesion, (MPa)
1.0

    Friction angle, (°)
32

    Energy release rate, , (J/m²)
1.0, 10.0

Fluid parameters

    Viscosity, (Pa·s)
0.001

    Bulk modulus of fluid, (GPa)
2.2

(a)

(b)

Figure 20: A single horizontal wellbore with two parallel preset fractures.

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 the direction.

Figures 21 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.

(a) Step  4000

(b) Step  8000

(c) Step  12000

(d) Step  16000

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  16000

Figure 22: Propagation of two initially parallel fractures driven by hydraulic pressure: anisotropic in situ stress field.

As shown in Figures 21 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 Figure 23, the fractures under the anisotropic in situ stress condition propagate longer in the direction than those under the isotropic in situ stress condition, while the influence domain in the -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.

(a) Influence domain of fractures in the -direction

(b) Influence domain of fractures in the -direction

Figure 23: Influence domain of fractures under isotropic and anisotropic in situ stress conditions.

The 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 Figures 24 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.

(a) Step  4000

(b) Step  8000

(c) Step  12000

(d) Step  16000

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  16000

Figure 25: Maximum principal stress contours under an anisotropic in situ stress field (positive denotes tension).

As shown in Figure 26, 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.

(a) Under isotropic in situ stress field condition

(b) Under anisotropic in situ stress field condition

Figure 26: Comparison of propagation paths by FDEM-Fluid and Wu et al. [10].

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 Figure 27, 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 Figure 28. The fluid flows into the preset fractures. Then, the fractures start to propagate, driven by the fluid pressure. It can be observed that fractures and 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 , fracture has a larger deviation because of the stress-shadow effect induced by both fracture and . However, and propagate away from each other.

(a) Step  4000

(b) Step  8000

(c) Step  12000

(d) Step  16000

Figure 28: Evolution of fracture geometries and the distribution of fluid pressure.

Figure 29 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.

(a)

(b)

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.

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.

Conflicts of Interest

The authors declare that there are no conflicts of interest related to this paper.

Acknowledgments

This research is supported by the National Basic Research Program of China (973 Program) (Grant no. 2014CB046900) and the National Natural Science Foundation of China (Grant no. 41602296).

References

United States Department of Energy, Modern shale gas development in the United States: A Primer, 2009.

K. Sato, C. A. Wright, and M. Ichikawa, “Post-frac analyses indicating multiple fractures created in a volcanic formation,” SPE Production and Facilities, vol. 14, no. 4, pp. 277–283, 1999.

J. E. Olson, “Multi-fracture propagation modeling: Applications to hydraulic fracturing in shales and tight gas sands,” in Proceedings of the 42nd U.S. Rock Mechanics - 2nd U.S.-Canada Rock Mechanics Symposium 2008, July 2008.

H. Wang, “Numerical investigation of fracture spacing and sequencing effects on multiple hydraulic fracture interference and coalescence in brittle and ductile reservoir rocks,” Engineering Fracture Mechanics, vol. 157, pp. 107–124, 2016.

P. N. Mutalik and B. Gibson, “Case history of sequential and simultaneous fracturing of the Barnett Shale in Parker County,” in Proceedings of the SPE Annual Technical Conference and Exhibition, ATCE 2008, pp. 3203–3209, September 2008.

G. Waters, B. Dean, R. Downie, K. Kerrihard, L. Austbo, and B. McPherson, “Simultaneous hydraulic fracturing of adjacent horizontal wells in the woodford shale,” in Proceedings of the SPE Hydraulic Fracturing Technology Conference 2009, pp. 694–715, January 2009.

Y. Yongtao, T. Xuhai, Z. Hong, L. Quansheng, and L. Zhijun, “Hydraulic fracturing modeling using the enriched numerical manifold method,” Applied Mathematical Modelling, vol. 53, pp. 462–486, 2018.

A. Dahi-Taleghani and J. E. Olson, “Numerical modeling of multistranded-hydraulic-fracture propagation: accounting for the interaction between induced and natural fractures,” SPE Journal, vol. 16, no. 3, pp. 575–581, 2011.

X. Weng, O. Kresse, C. Cohen, R. Wu, and H. Gu, “Modeling of hydraulic-fracture-network propagation in a naturally fractured formation,” SPE Production and Operations, vol. 26, no. 4, pp. 368–380, 2011.

R. Wu, O. Kresse, X. Weng, C. Cohen, and H. Gu, “Modeling of Interaction of Hydraulic Fractures in Complex Fracture Networks,” in Proceedings of the SPE Hydraulic Fracturing Technology Conference, The Woodlands, Tex, USA, 2012.

Y. Yongtao, T. Xuhai, Z. Hong, L. Quansheng, and H. Lei, “Three-dimensional fracture propagation with numerical manifold method,” Engineering Analysis with Boundary Elements, vol. 72, pp. 65–77, 2016.

N. B. Nagel and M. Sanchez-Nagel, “Stress shadowing and microseismic events: a numerical evaluation,” Society of Petroleum Engineers, 2011.

J. Adachi, E. Siebrits, A. Peirce, and J. Desroches, “Computer simulation of hydraulic fractures,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 5, pp. 739–757, 2007.

O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, McGraw-Hill, 2008.

K. Sato, M. Itaoka, and T. Hashida, “FEM simulation of mixed mode crack propagation induced by hydraulic fracturing,” Processing and Properties of Porous Nickel Titanium, 2013.

A. Paluszny and R. W. Zimmerman, “Numerical simulation of multiple 3D fracture propagation using arbitrary meshes,” Computer Methods Applied Mechanics and Engineering, vol. 200, no. 9, pp. 953–966, 2011.

A. Paluszny, X. Tang, M. Nejati, and R. W. Zimmerman, “A direct fragmentation method with Weibull function distribution of sizes based on finite- and discrete element simulations,” International Journal of Solids and Structures, vol. 80, pp. 38–51, 2016.

N. Moës, J. Dolbow, and T. Belytschko, “A finite element method for crack growth without remeshing,” International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131–150, 1999.

E. Gordeliy and A. Peirce, “Coupling schemes for modeling hydraulic fracture propagation using the XFEM,” Computer Methods Applied Mechanics and Engineering, vol. 253, no. 1, pp. 305–322, 2013.

B. Lecampion, “An extended finite element method for hydraulic fracture problems,” Communications in Numerical Methods in Engineering, vol. 25, no. 2, pp. 121–133, 2009.

K. Su, X. Zhou, X. Tang, X. Xu, and Q. Liu, “Mechanism of cracking in dams using a hybrid FE-meshfree method,” International Journal of Geomechanics, vol. 17, no. 9, Article ID 04017071, 2017.

P. Gupta and C. A. Duarte, “Simulation of non-planar three-dimensional hydraulic fracture propagation,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 13, pp. 1397–1430, 2014.

T. A. Cruse, “Numerical solutions in three dimensional elastostatics,” International Journal of Solids and Structures, vol. 5, no. 12, pp. 1259–1274, 1969.

D. Zhou, P. Zheng, P. He, and J. Peng, “Hydraulic fracture propagation direction during volume fracturing in unconventional reservoirs,” Journal of Petroleum Science and Engineering, vol. 141, pp. 82–89, 2016.

M. Behnia, K. Goshtasbi, G. Zhang, and S. H. M. Yazdi, “Numerical modeling of hydraulic fracture propagation and reorientation,” European Journal of Environmental and Civil Engineering, vol. 19, no. 2, pp. 152–167, 2015.

N. K. Mukhopadhyay, S. K. Maiti, and A. Kakodkar, “A review of SIF evaluation and modelling of singularities in BEM,” Computational Mechanics, vol. 25, no. 4, pp. 358–375, 2000.

F. K. Wittel, H. A. Carmona, F. Kun, and H. J. Herrmann, “Mechanisms in impact fragmentation,” International Journal of Fracture, vol. 154, no. 1-2, pp. 105–117, 2008.

B. Damjanac, I. Gil, M. Pierce, M. Sanchez, As. AV, and J. Mclennan, “A new approach to hydraulic fracturing modeling in naturally fractured reservoirs,” Technology & Health Care, vol. 18, no. 4-5, pp. 325–334, 2010.

S. Marina, I. Derek, P. Mohamed, S. Yong, and EK. Imo-Imo, “Simulation of the hydraulic fracturing process of fractured rocks by the discrete element method,” Environmental Earth Sciences, vol. 73, no. 12, pp. 8451–8469, 2015.

A. Munjiza, The Combined Finite-Discrete Element Method, John Wiley & Sons, 2004.

O. K. Mahabadi, N. X. Randall, Z. Zong, and G. Grasselli, “A novel approach for micro-scale characterization and modeling of geomaterials incorporating actual material heterogeneity,” Research Letters, vol. 39, no. 1, p. 1303, 2012.

Q. Lei, J.-P. Latham, and J. Xiang, “Implementation of an Empirical Joint Constitutive Model into Finite-Discrete Element Analysis of the Geomechanical Behaviour of Fractured Rocks,” Rock Mechanics and Rock Engineering, vol. 49, no. 12, pp. 4799–4816, 2016.

J. P. Latham, J. Xiang, M. Belayneh, H. M. Nick, C. F. Tsang, and M. J. Blunt, “Modelling stress-dependent permeability in fractured rock including effects of propagating and bending fractures,” International Journal of Rock Mechanics and Mining Sciences, vol. 57, pp. 100–112, 2013.

Q. Lei, JP. Latham, J. Xiang, and P. Lang, “Coupled FEMDEM-DFN model for characterising the stress-dependent permeability of an anisotropic fracture system,” International Conference on Discrete Fracture Network Engineering, 2014.

A. Obeysekara, Q. Lei, P. Salinas et al., “A fluid-solid coupled approach for numerical modeling of near-wellbore hydraulic fracturing and flow dynamics with adaptive mesh refinement,” in Proceedings of the 50th US Rock Mechanics/Geomechanics Symposium 2016, pp. 1688–1699, June 2016.

C. Yan, H. Zheng, G. Sun, and X. Ge, “Combined finite-discrete element method for simulation of hydraulic fracturing,” Rock Mechanics and Rock Engineering, vol. 49, no. 4, pp. 1389–1410, 2016.

Y.-Y. Jiao, H.-Q. Zhang, X.-L. Zhang, H.-B. Li, and Q.-H. Jiang, “A two-dimensional coupled hydromechanical discontinuum model for simulating rock hydraulic fracturing,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 39, no. 5, pp. 457–481, 2015.

M. Chen, H. Jiang, G. Q. Zhang, and Y. Jin, “The experimental investigation of fracture propagation behavior and fracture geometry in hydraulic fracturing through oriented perforations,” Petroleum Science and Technology, vol. 28, no. 13, pp. 1297–1306, 2010.

A. Munjiza, K. R. F. Andrews, and J. K. White, “Combined single and smeared crack model in combined finite-discrete element analysis,” International Journal for Numerical Methods in Engineering, vol. 44, no. 1, pp. 41–57, 1999.

A. Munjiza and K. R. F. Andrews, “NBS contact detection algorithm for bodies of similar size,” International Journal for Numerical Methods in Engineering, vol. 43, no. 1, pp. 131–149, 1998.

D. S. Dugdale, “Yielding of steel sheets containing slits,” Journal of the Mechanics and Physics of Solids, vol. 8, no. 2, pp. 100–104, 1960.

G. I. Barenblatt, “The mathematical theory of equilibrium cracks in brittle fracture,” Advances in Applied Mechanics, vol. 7, no. C, pp. 55–129, 1962.

P. A. Witherspoon, J. S. Y. Wang, K. Iwai, and J. E. Gale, “Validity of cubic law for fluid flow in a deformable rock fracture,” Water Resources Research, vol. 16, no. 6, pp. 1016–1024, 2010.

Itasca, UDEC version 4.0 users manuals, Itasca Consulting Group Inc, 2005.

K. Iwai, Fundamental studies of the fluid flow through a single fracture, 1976.

H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Oxford Science Publications, 2nd edition, 1995.

J. Almaguer, J. Manrique, S. Wickramasuriya et al., “Oriented perforating minimizes flow restrictions and friction pressures during fracturing,” Oilfield Rev, vol. 14, no. 1, pp. 16–31, 2002.

H. Wang, “Numerical modeling of non-planar hydraulic fracture propagation in brittle and ductile rocks using XFEM with cohesive zone method,” Journal of Petroleum Science and Engineering, vol. 135, pp. 127–140, 2015.

J. Olson and D. D. Pollard, “Inferring paleostresses from natural fracture patterns: a new method,” Geology, vol. 17, no. 4, pp. 345–348, 1989.

D. Kumar and A. Ghassemi, “A three-dimensional analysis of simultaneous and sequential fracturing of horizontal wells,” Journal of Petroleum Science and Engineering, vol. 146, pp. 1006–1025, 2016.

K. Wu, “Simultaneous multifracture treatments: fully coupled fluid flow and fracture mechanics for horizontal Wells,” SPE Journal, vol. 20, no. 2, 2015.