Research Article  Open Access
Quansheng Liu, Lei Sun, Pingli Liu, Lei Chen, "Modeling Simultaneous Multiple Fracturing Using the Combined FiniteDiscrete Element Method", Geofluids, vol. 2018, Article ID 4252904, 20 pages, 2018. https://doi.org/10.1155/2018/4252904
Modeling Simultaneous Multiple Fracturing Using the Combined FiniteDiscrete Element Method
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 stressshadow 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 FDEMFluid, to handle hydromechanicalfracture 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 parallelplate flow model. Several tests are carried out to validate the application of FDEMFluid in hydraulic fracturing simulation. Then, this FDEMFluid 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 (simulfrac) 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 simulfrac over standalong wells [5], microseismic data [6], and numerical simulations [3, 7–12] also show a complex fracturing network produced by simulfrac, the reasons for simulfrac 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 hydromechanicalfracture (HMF) 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 stressshadow 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 wellestablished 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 fieldscale 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 finitediscrete 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 (HM) coupling problems [33–36].
Although FDEM has been welldeveloped and applied by many researchers to study the process of hydraulic fracturing, very little literature can be found about the study of stressshadow effect with FDEM. In this paper, FDEM coupled with a fluid algorithm, referred to as FDEMFluid, 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 FDEMFluid 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 stressshadow 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 3node triangular elements and 4node joint elements, which are inserted between the triangular elements. The stress state and deformation of the bulk material are captured by the linearelastic 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.
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 aswhere represents the constant damping coefficient and I is an identity matrix.
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.
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)
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 III. 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 stressfree surface (a real fracture) is generated when the opening reaches a certain residual value, .
(a)
(b)
(c)
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 MohrCoulomb 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 III 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 3node triangular elements (in gray) and 4node joint elements (in red). The fluid domain is discretized with hydraulic grids by using 2node hydraulic elements (in blue), which are inserted into the broken joint elements.
3.1. Basic Assumptions of Fluid Flow
The hydraulic fracturing process is complicated, especially in regard to HMF 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 parallelplate 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 bywhere 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 aswhere 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. HydroMechanicalFracture (HMF) Coupling Procedure
The hydromechanicalfracture (HMF) 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.
4.1. Coupling Procedure of FDEMFluid
This FDEMFluid 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.
A typical time step in FDEMFluid 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.
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.
and are the hydraulic pressure on hydraulic nodes and , respectively. The total force on this surface can be calculated aswhere 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 HM 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 FDEMFluid
Two numerical tests for the present FDEMFluid 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 FDEMFluid. The mesh sensitivity of fluid flow in FDEMFluid 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
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 FDEMFluid 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.
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 FDEMFluid for simulating fluid flow in fractures are validated.
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
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 FDEMFluid in simulating hydraulic fracturing. A cube plaster sample with preset fractures and a central injection hole is prepared, which is then fractured by injecting highpressure 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
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.

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)
(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 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].
(a) Top view
(b) Sectional view
These two examples validate the ability of FDEMFluid 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 multiplefracture propagation, fracture interaction (i.e., the stressshadow 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^{3}/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].

(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 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
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
As shown in Figures 21 and 22, in both cases, two preset fractures propagate away from each other with a sidebyside configuration because of the mechanical fracture interaction (i.e., stressshadow 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
The reason fractures turn is because of the mechanical fracture interaction (i.e., the stressshadow 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 stressshadow 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 stressshadow effect. In other words, anisotropic farfield 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
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
As shown in Figure 26, the FDEMFluid results agree well with previous numerical results [10] in both isotropic and anisotropic in situ stress conditions. Therefore, this FDEMFluid method is able to simulate the multiple hydraulic fracturing, with considering the HMF coupling, the stressshadow effect, and in situ stress.
(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 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.
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 stressshadow 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 stressshadow 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 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)
(a)
(b)
7. Conclusions
The interaction (i.e., the stressshadow effect) during simultaneous multiplefracture 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 FDEMFluid, based on FDEM and the parallelplate flow model is presented to capture the development of complex hydraulic fracturing propagation. This FDEMFluid 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 FDEMFluid. 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 stressshadow 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 stressshadow 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 stressshadow 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. View at: Publisher Site
 K. Sato, C. A. Wright, and M. Ichikawa, “Postfrac analyses indicating multiple fractures created in a volcanic formation,” SPE Production and Facilities, vol. 14, no. 4, pp. 277–283, 1999. View at: Publisher Site  Google Scholar
 J. E. Olson, “Multifracture 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. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Google Scholar
 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. View at: Google Scholar
 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. View at: Google Scholar
 A. DahiTaleghani and J. E. Olson, “Numerical modeling of multistrandedhydraulicfracture propagation: accounting for the interaction between induced and natural fractures,” SPE Journal, vol. 16, no. 3, pp. 575–581, 2011. View at: Publisher Site  Google Scholar
 X. Weng, O. Kresse, C. Cohen, R. Wu, and H. Gu, “Modeling of hydraulicfracturenetwork propagation in a naturally fractured formation,” SPE Production and Operations, vol. 26, no. 4, pp. 368–380, 2011. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 Y. Yongtao, T. Xuhai, Z. Hong, L. Quansheng, and H. Lei, “Threedimensional fracture propagation with numerical manifold method,” Engineering Analysis with Boundary Elements, vol. 72, pp. 65–77, 2016. View at: Google Scholar
 N. B. Nagel and M. SanchezNagel, “Stress shadowing and microseismic events: a numerical evaluation,” Society of Petroleum Engineers, 2011. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar
 O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, McGrawHill, 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. View at: Google Scholar
 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. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Publisher Site  Google Scholar
 K. Su, X. Zhou, X. Tang, X. Xu, and Q. Liu, “Mechanism of cracking in dams using a hybrid FEmeshfree method,” International Journal of Geomechanics, vol. 17, no. 9, Article ID 04017071, 2017. View at: Publisher Site  Google Scholar
 P. Gupta and C. A. Duarte, “Simulation of nonplanar threedimensional hydraulic fracture propagation,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 38, no. 13, pp. 1397–1430, 2014. View at: Publisher Site  Google Scholar
 T. A. Cruse, “Numerical solutions in three dimensional elastostatics,” International Journal of Solids and Structures, vol. 5, no. 12, pp. 1259–1274, 1969. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 F. K. Wittel, H. A. Carmona, F. Kun, and H. J. Herrmann, “Mechanisms in impact fragmentation,” International Journal of Fracture, vol. 154, no. 12, pp. 105–117, 2008. View at: Publisher Site  Google Scholar
 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. 45, pp. 325–334, 2010. View at: Google Scholar
 S. Marina, I. Derek, P. Mohamed, S. Yong, and EK. ImoImo, “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. View at: Google Scholar
 A. Munjiza, The Combined FiniteDiscrete Element Method, John Wiley & Sons, 2004. View at: Publisher Site
 O. K. Mahabadi, N. X. Randall, Z. Zong, and G. Grasselli, “A novel approach for microscale characterization and modeling of geomaterials incorporating actual material heterogeneity,” Research Letters, vol. 39, no. 1, p. 1303, 2012. View at: Google Scholar
 Q. Lei, J.P. Latham, and J. Xiang, “Implementation of an Empirical Joint Constitutive Model into FiniteDiscrete Element Analysis of the Geomechanical Behaviour of Fractured Rocks,” Rock Mechanics and Rock Engineering, vol. 49, no. 12, pp. 4799–4816, 2016. View at: Publisher Site  Google Scholar
 J. P. Latham, J. Xiang, M. Belayneh, H. M. Nick, C. F. Tsang, and M. J. Blunt, “Modelling stressdependent 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. View at: Google Scholar
 Q. Lei, JP. Latham, J. Xiang, and P. Lang, “Coupled FEMDEMDFN model for characterising the stressdependent permeability of an anisotropic fracture system,” International Conference on Discrete Fracture Network Engineering, 2014. View at: Google Scholar
 A. Obeysekara, Q. Lei, P. Salinas et al., “A fluidsolid coupled approach for numerical modeling of nearwellbore 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. View at: Google Scholar
 C. Yan, H. Zheng, G. Sun, and X. Ge, “Combined finitediscrete element method for simulation of hydraulic fracturing,” Rock Mechanics and Rock Engineering, vol. 49, no. 4, pp. 1389–1410, 2016. View at: Publisher Site  Google Scholar
 Y.Y. Jiao, H.Q. Zhang, X.L. Zhang, H.B. Li, and Q.H. Jiang, “A twodimensional 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 A. Munjiza, K. R. F. Andrews, and J. K. White, “Combined single and smeared crack model in combined finitediscrete element analysis,” International Journal for Numerical Methods in Engineering, vol. 44, no. 1, pp. 41–57, 1999. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 G. I. Barenblatt, “The mathematical theory of equilibrium cracks in brittle fracture,” Advances in Applied Mechanics, vol. 7, no. C, pp. 55–129, 1962. View at: Publisher Site  Google Scholar
 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. View at: Publisher Site  Google Scholar
 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. View at: MathSciNet
 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. View at: Google Scholar
 H. Wang, “Numerical modeling of nonplanar 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. View at: Publisher Site  Google Scholar
 J. Olson and D. D. Pollard, “Inferring paleostresses from natural fracture patterns: a new method,” Geology, vol. 17, no. 4, pp. 345–348, 1989. View at: Publisher Site  Google Scholar
 D. Kumar and A. Ghassemi, “A threedimensional analysis of simultaneous and sequential fracturing of horizontal wells,” Journal of Petroleum Science and Engineering, vol. 146, pp. 1006–1025, 2016. View at: Publisher Site  Google Scholar
 K. Wu, “Simultaneous multifracture treatments: fully coupled fluid flow and fracture mechanics for horizontal Wells,” SPE Journal, vol. 20, no. 2, 2015. View at: Publisher Site  Google Scholar
Copyright
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.