A Numerical Study on the Propagation Mechanisms of Hydraulic Fractures in Fracture-Cavity Carbonate Reservoirs

Fang Shi; Daobing Wang; Xiaogang Chen

doi:10.32604/cmes.2021.015384

1Jiangsu Key Laboratory of Advanced Manufacturing Technology, Huaiyin Institute of Technology, Huai’an, 223003, China
2School of Mechanical Engineering, Beijing Key Laboratory of Pipeline Critical Technology and Equipment for Deepwater Oil & Gas Development, Beijing Institute of Petrochemical Technology, Beijing, 102617, China
*Corresponding Author: Fang Shi. Email: shifang@hyit.edu.cn
Received: 15 December 2020; Accepted: 09 February 2021

Carbonate reservoirs, which contain naturally-formed fractures and caves, widely exist in the world and have huge potential for exploitation [1,2]. These caves of different shapes and sizes provide storage space for oil and gas, and fractures are the potential flow paths for hydrocarbon recovery. As a successful stimulation method, the aim of the hydraulic fracturing technique is to connect wellbores and caves through hydraulic fractures as well as natural fractures. In order to achieve maximum oil or gas extraction, it is of great significance to perform the study on the propagation mechanisms of hydraulic fractures in fracture-cavity reservoirs.

During the past decades, based on different numerical approaches like the finite element method (FEM) [3,4], the displacement discontinues method (DDM) [5], the phase field method [6,7], the distinct element method (DEM) [8], the extended finite element method (XFEM) [9–12], and the proper generalized decomposition method (PGD) [13], a variety of numerical models have been established and applied by researchers to study hydraulic fracture propagation in consideration of different kinds of influence factors. Among these studies, Schrefler et al. [3] proposed an adaptive refinement technique to simulate hydraulic fracturing problems based on the generalized finite element formulation. Song et al. [4] performed a series of numerical simulations and examined the influence of some key factors on hydraulic fracturing propagation using RFPA2D-Flow. Zhou et al. [7] proposed a phase field model for fluid-driven dynamic crack propagation simulation in the poroelastic media. Wang et al. [13] numerically studied the hydraulic fracturing problems via the PGD method. Recently, the interaction mechanisms between hydraulic fractures (HFs) and pre-existing natural fractures (NFs) in reservoirs have increased markedly [11,14,15]. When an HF meets an NF, it can be deflected or arrested by the NF, active the NF and branch along one or two sides of the NF, or even bypass the NF [14], leading to very complicated patterns of HF-NF interactions. Readers are referred to [14,15] for a comprehensive review of this topic. However, despite these efforts, a systematic study of mechanisms and factors that govern the propagation behavior of hydraulic fractures around natural cavities has been seldomly reported in the literature. Hence, an attempt is made in this study to present a comprehensive investigation of the interaction mechanism of hydraulic fracture and naturally-formed caves and fractures by establishing a reliable numerical model.

The difference in the propagation mechanisms of hydraulic fractures between fracture-cavity reservoirs and conventional reservoirs is mainly caused by the stress concentration around cavities. As shown in Fig. 1, the analytical solution of stress in the tangential direction, σθ, surrounding the cavity of radius r can be written as [16]

where λ=σh/σH denotes the lateral stress coefficient. From Eq. (1), it can be easily concluded that if λ<1/3, tensile stress zones exist at the left and right sides near the cavity, thus promoting the propagation of a nearby hydraulic fracture towards the cavity located ahead of the fracture tip. If 1/3<λ<1/2, although the tensile stress zone disappears, σθ<σH can still be satisfied around the left and right sides of the cavity, which is also favorable for hydraulic fracturing. However, if λ>1/2, the hydraulic fracture will be arrested or deflected from its original propagation path due to the compressive stress concentration phenomenon. With this in mind, therefore, it can be conjectured that once the cavities are considered in reservoirs, quite different physical mechanisms and propagation-path scenarios of hydraulic fracturing can be obtained. For instance, Liu et al. [6] investigated the effects of full-filled cavities on hydraulic fracturing based on the phase field method. They found that the propagation paths are sensitive to Young’s modulus ratio between rock formation and cavity, spatial location of the cavity as well as differential in-situ stress. Likewise, Luo et al. [17] established a coupled seepage-free flow-mechanical model and studied the influence of the fluid pressure inside a cavity as well as the perforation direction on the hydraulic fracture propagation paths. It was found that higher fluid pressure in the cavity is more likely to reduce fracture path deflection and the perforation direction has a profound influence on the fracture propagation path. Cheng et al. [11] built a numerical model and studied the influence of cavity location as well as the difference of in-situ stress on the propagation path of a hydraulic fracture. He et al. [18] proposed a coupled hydromechanical model and studied the effects of existing pressurized cavities on the propagation of fluid-driven fracture. Besides, a numerical simulator called FEMM-FracFlow was used by Wang et al. [19] to investigate the behaviors of hydraulic fractures and found that fracture paths are closely related to the cavities, pumping rate of fracturing fluid, and the in-situ stress level. In addition to hydraulic fracturing technology, acid fracturing is one of the major measures of reservoir stimulation, especially for carbonate reservoirs. Zhao et al. [20] established a numerical tool for acid fracturing simulation and studied the mechanism of hydraulic fracture propagation in cases involving a single cavity and a pre-existing natural fracture. Simulation results [21] suggest that a smaller approaching angle of the hydraulic fracture towards the natural fracture and longer natural fracture benefit the penetration of hydraulic fracture into the cavity. As another alternative, the explosive pulsed fracturing technique is receiving more attention from researchers [21]. Wang et al. [22] developed a numerical model to investigate the pulsed fracturing process in fracture-cavity reservoirs and found that natural fracture can potentially guide the pulsed fracturing path under the action of pulse pressure.

Figure 1: A cylindrical cavity of radius a under the action of biaxial stresses. Plane strain model. Stress σH is acting in the horizontal direction, and stress σh is acting in the vertical direction. r and θ represent polar coordinates

Although some works aiming at explaining the mechanisms of interaction among HF, NF, and cavity have been done by researchers, it still lacks a systematic study and deeper understanding of crucial factors that influence the hydraulic fracturing efficiency in fracture-cavity reservoirs. On the other hand, existing studies mainly focus on simple geometrical configuration, and the study of the influence of cavity on hydraulic fractures in a wellbore is not available in the literature. In this study, the XFEM proposed by Belytschko et al. [23] is adopted to describe the displacement field around fractures and cavities, the intersection between an HF and an NF, and the penetration of a fracture into a cavity. Consequently, the tedious remeshing process can be avoided after the propagation of fractures [10]. The fluid flow within fractures is described by Reynold’s equation [9] which is discretized using the FEM. Afterward, the displacements of all solid nodes including enriched nodes, and the pressure distribution of all fluid nodes along the hydro-fracture are obtained by solving the fully-coupled governing equations [10] using the Newton–Raphson (N–R) method.

This paper is presented as follows. The description of the numerical model is given in Section 2. Implementation and validation of the proposed numerical model are presented in Section 3. Mechanisms and factors that govern the propagation behavior of hydraulic fractures around natural cavities are systematically studied in Section 4. Major conclusions are summarized in Section 5.

A shown in Fig. 2, we consider the domain Ω with a hydraulic fracture ΓHF, a frictional [24] natural fracture ΓNF, and a cavity inside. The boundary of the domain Ω is Γ with a normal vector nΓ. Ω is fixed at boundary Γu and traction t is applied at boundary Γt. s represents the coordinate system defined along the hydraulic fracture. The incompressible Newtonian fluid is pumped at a rate of Qinj into ΓHF. To make the numerical model more tractable, the following assumptions are made. The fracture growth process is treated as quasi-static [25]. The fluid lag phenomenon [14] is ignored. Natural cavities usually exist in three forms, i.e., non-filled, partly-filled, and fully-filled cavities [2]. Non-filled cavities are studied in this article. In addition, although natural cavities are of quite irregular shapes, they are treated as regular circles or ellipses in this article for the sake of modeling convenience.

Figure 2: Illustration of the hydraulic fracturing problem in the fracture-cavity reservoir

According to the linear elastic assumption of rock formation, the governing equations for quasi-static deformation can be written as [10]

where σ denotes the Cauchy stress tensor, tcont represents the traction vector caused by contact behaviors between fracture surfaces.

The fluid flow along the hydraulic fracture can be described by Reynold’s equation [25]

where w, t, p, and μ represent fracture aperture, time, fluid pressure, and fluid viscosity, respectively. The fracture aperture w can be calculated by w=[[u]]⋅nΓHF, where [[u]] denotes the displacement jump across the fracture. The term gL represents the fluid leak-off from the hydraulic fracture into the surrounding rock formation, which can be written as [9,26]

where CL is the fluid leak-off coefficient, and t0(s) represents the time instant when fluid starts to leak for a point s along the hydraulic fracture.

After introducing the trial function u(x, t) and test function δu(x, t), we can obtain the weak form of Eq. (2):

∫Ωδε:σdΩ+∫ΓHF[[δu]]⋅pnΓHFdΓ+∫ΓNF[[δu]]⋅tcontdΓ=∫Γtδu⋅tdΓ.(5)

Besides, by introducing the test functionδp(s, t) and integrating by parts, the weak form of Eq. (3) can be written as:

Spatial and time discretization of Eqs. (5) and (6) can be found in our previous work [10]. Afterward, the solid–fluid coupling problem is governed by Eqs. (2) and (3) can be iteratively solved by using the Newton-Raphson method [10].

In this paper, we use the XFEM to describe the displacement field around fractures and cavities. By simply introducing additional degrees of freedoms (DOFs), i.e., the enriched DOFs, the remeshing and data mapping between old and new meshes can be avoided [27] and the demand on mesh refinement around fracture tip can be greatly alleviated owing to the fracture tip enrichment function, leading to a significant reduction of the problem complexity and computational time. This is a substantial improvement for numerical simulation of problems containing strong or weak discontinuities, especially for simulation cases where a large number of fractures and caves are involved. Since its introduction in 1999 [23,27], the XFEM has been adopted by many researchers to study various kinds of problems related to hydraulic fracturing. For problems where the intersection of an HF and an NF, as well as the intersection of a fracture and a cavity, are considered, the displacement of a point x inside the domain Ω can be written as

u(x)=∑I∈SallNI(x)uI+∑I∈SfracNI(x)H(x)aI+∑I∈ScavityNI(x)C(x)bI+∑I∈StipNI(x)∑l=14Fl(x)cIl+∑I∈SjuncNI(x)J(x)dI+∑I∈SpenNI(x)P(x)eI,(7)

in which NI denotes the shape function of node I, Sall denotes the set of all nodes whose displacement vector is denoted by uI. Sfrac, Scavity, Stip, Sjunc, and Spen are the set of enriched nodes for Heaviside enrichment, cavity enrichment, tip enrichment, junction enrichment, and penetration enrichment, respectively, and their corresponding vectors of enriched DOFs are denoted by aI, bI, cI, dI and eI, respectively. Heaviside enrichment function H(x) equals 1 on one side of the fracture and −1 on the other side. The cavity enrichment function C(x) takes the value of 0 if the point x lies inside the cavity and 1 if outside the cavity. The fracture tip enrichment function Fl(x) is composed of four parts [27]

where (r, θ) represents the local polar co-ordinates defined at fracture tips [27]. As depicted in Fig. 3, in addition to the Heaviside enrichment nodes, the junction enrichment nodes are necessary to describe the intersection of two fractures. The T-shaped junction considered in this article, J(x) takes different values of 1, −1, and 0 in different subdomains [10]. Besides the Heaviside enriched nodes and cavity enriched nodes, the penetration enriched nodes are presented to describe the penetration of a fracture into the cavity, just as illustrated in Fig. 4. The penetration enrichment function P(x) takes the following form

which indicates that P(x) equals 0 if x lies inside the cavity, and otherwise equals 1 or −1 in accordance with the Heaviside enrichment function.

Figure 3: Depiction of the enriched nodes for a T-shaped junction formed by two fractures. The fluid nodes are also shown. Nodes whose support domain contains two intersected fractures should be selected as junction enriched nodes

Figure 4: Depiction of the enriched nodes for penetration of a fracture into a cavity. The fluid nodes are also shown. Nodes whose support domain contains both the boundary of the cavity and the fracture should be selected as penetration enriched nodes

Frictional and cemented fractures are two common types of natural fractures in reservoirs [24]. In this article, the frictional natural fractures are considered due to their stronger influence on the creation of the fracture network in comparison with cemented natural fractures [24]. Thus, to avoid embedding of fracture surfaces under the effect of compressive stress, the frictional behaviors of natural fractures must be considered. The no-embedding conditions of fracture surfaces can be expressed as [28]

where [[u]]N and tNcont represent the normal component of [[u]] and tcont, respectively. Based on Coulomb’s friction law, the contact status can be determined according to the following expression [29]

where tTcont represents the tangential component of tcont, ϕf represents the friction angle of the fracture surface, and cf is the cohesive strength between fracture surfaces. If Ff < 0, then the fracture surfaces are in a bonding state. If Ff = 0, then the fracture surfaces are in a slipping state.

In this article, the penalty function method together with the Newton–Raphson method [10] is adopted to consider the frictional force and the contact status in the simulation. For each Newton–Raphson iteration Step i, the following linear system is solved

and the global conventional DOFs vector ũ and enriched DOFs vector ã are updated according to

In Eq. (12), Kuu, Kua, Kau, Kaa, Ψ1 and Ψ2 take the following forms, respectively

{Kuu=∫Ω(Bstd)TDBstddΩKua=∫Ω(Bstd)TDBenrdΩKau=∫Ω(Benr)TDBstddΩKaa=∫Ω(Benr)TDBenrdΩ+∫Γc(Nstd)TDcont(Nstd)TdΓ︸KcontactΨ1=∫Ω(Bstd)TσdΩ−∫Γt(Nstd)TtdΓΨ2=∫Ω(Benr)TσdΩ−∫Γt(Nenr)TtdΓ+∫Γc(Nstd)TtcontdΓ︸Fcontact(15)

In Eq. (15), N is the matrix of shape function, B is the matrix of shape function derivatives, D is the stress-strain matrix, and Dcont is related to contact status of fractures and can be written as [10]

Dcont={kN(nΓNF⊗nΓNF)+kT(mΓNF⊗mΓNF)forbondingstatekN(nΓNF⊗nΓFF)+ tanϕfkT(mΓNF⊗nΓNF),slippingstate(16)

where kN and kT are the penalty parameters in the normal and tangential directions, respectively, and mΓNF is the unit vector in the tangential direction of the natural frictional fracture. In Eq. (15), the frictional force between fracture surfaces is considered by the last term Fcontact of Ψ2.

The widely adopted maximum hoop-stress criterion [30] related to Mode-I (KI) and Mode-II (KII) stress intensity factors (SIFs) is used to check whether and along which direction the fracture will propagate

If Ke exceeds the fracture toughness of reservoirs (KIC), the fracture will propagate and deflect by the angle α compared to its original direction. The SIFs are calculated by using the interaction integral approach [31]. When an HF meets a frictional NF, whether the HF will propagate along the NF, cross the NF, or be arrested by the NF is determined using the extended Renshaw and Pollard criterion [32]. Implementation details of this criterion can be found in our previous article [24].

The proposed numerical model is implemented in an in-house Fortran code called PhiPsi (http://phipsi.top/): A general-purpose XFEM-based program. The level set method [33] is used to describe fractures and cavities, and to track the growth of fractures. The compressed sparse row (CSR) format is used to store global system matrices with optimized memory consumption. For the simulation cases where frictional natural fractures are considered (Section 4.1.6), two types of Newton–Raphson iteration loops are performed [10]. The outer loop solves the solid–fluid coupling equations (Eqs. (2) and (3)), and the inner loop determines the contact status of natural frictional fractures [10]. For each Newton–Raphson iteration step, an iterative solver called Lis [34] is adopted to solve the linear system. The flowchart of the simulation algorithm [10] is shown in Fig. 5.

The validation of the presented numerical model without considering cavities has been sufficiently conducted in our previous articles [9,10,24,35]. In this section, the propagation of an initial fracture towards a cavity in a panel under the action of tensile load is studied to verify the capacity of the proposed model for predicting interactions between fractures and cavities. As illustrated in Fig. 6, the bottom of a plane stress panel is fixed and the top of the panel is subjected to tensile stress σ. The geometric parameters H, h1, h2, l, b1, b2, and a are taken as 3.0, 1.5, 0.1, 1.2, 0.7, and 0.2 m, respectively. Besides, the values of Young’s modulus and Poisson’s ratio are given as 30 GPa and 0.2, respectively. The tensile stress σ is 5 KPa. The discretized model consists of 3,400 quadrilateral elements in total. The fracture propagates by an increment of 0.045 m and penetrates the cavity after 15 propagation steps. The resulting vertical displacement field of the panel and the enriched nodes are presented in Fig. 7 where six penetration enriched nodes can be seen. This problem has been studied by Nguyen et al. [36] using the extended meshfree Galerkin method. The comparison of the calculated propagation path with [36] is visualized in Fig. 8. As expected, the propagation path predicted by the proposed method is well consistent with that predicted by other numerical method, indicating that the interactions between fracture and cavity can be well captured using the proposed model.

Figure 6: Illustration of a panel that contains a fracture of initial length l and a cavity of radius a. The bottom of the panel is fixed and the top of the panel is subjected to tensile stress σ. Plane stress condition is assumed

Figure 7: Vertical displacement field and enriched nodes of the last propagation step when the fracture penetrates into the cavity. The deformation scale factor is 200

Figure 8: Comparison of the propagation path with the existing simulation result in the literature. The black dots represent locations of the fracture tip during the propagation process. The black circle represents the cavity

The simulation of a frictional natural fracture of length 10 m in a plane strain plate under the action of compressive stress will be performed to verify the frictional model presented in Section 2.3. As shown in Fig. 9a fracture at 50∘ with respect to x-direction is placed at the center of the model. The bottom edge of the plate is fixed in vertical direction and the left bottom corner is fixed in both vertical and horizontal directions. The material parameters E, υ of the plate are 20 GPa and 0.2, respectively. The friction angle of the fracture surface and the cohesive strength of the frictional fracture are π/10 and 0 MPa, respectively. kN and kT are taken as 104 GPa/m. The applied stress σy on the top edge is 5 MPa. The discretized model consists of 2,184 quadrilateral elements in total. As a comparison, this problem is also simulated using the FEM (ANSYS) with 6,806 quadrilateral elements. The contours of the maximum principal stress obtained from both the XFEM and the FEM solutions are shown in Fig. 10, from which good agreements can be seen. In order to make a more accurate comparison, the relative shear displacement (displacement in the tangential direction) between fracture surfaces obtained from both the XFEM and the FEM solutions are presented in Fig. 11. Again, good agreements can be observed, which means that the proposed scheme to simulate contact behaviors [37] of frictional fractures is accurate and reliable.

Figure 9: Illustration of a frictional natural fracture in a plane strain plate under the action of compressive stress

Figure 10: Contours of maximum principal stress obtained from: (a) The XFEM solution and (b) The FEM solution

In this section, interactions between hydraulic fracture and cavity will be studied via a simple model shown in Fig. 12. The influence of factors including in-situ stress, viscosity of the injected fluid, pumping rate of the fluid, size of the cavity, shape of the cavity, and natural fracture will be thoroughly investigated. The size of the model is 5m×10m. The radius of the wellbore is 0.1 m. The length of the initial hydraulic fracture is 0.5 m. The radius of the cavity denoted by a is set to 0.5 m. The vertical distance between the cavity and the fracture denoted by d is taken as 0.25 m (d = 0.5 a). The in-situ stresses σh and σH are 40 and 50 MPa, respectively. Thus, the lateral stress coefficient λ equals 0.8, and the in-situ stress difference σD equals 10 MPa. The Young’s modulus of the rock formation is E = 30 GPa, Poisson’s ratio is ν=0.2, and fracture toughness is KIC=1.0MPa⋅m1/2. In addition, the viscosity of the injected fluid is μ=0.01 Pa⋅s, the fluid pumping rate is Qinj = 0.001 m2/s, and the fluid leak-off coefficient CL is taken as 10−5 m/s1/2. The model is discretized into 5,618 quadrilateral elements. The hydraulic fracture propagates by an increment of 0.15 m.

Figure 12: Geometry and boundary conditions of the model studied in Section 4.1. Plane strain condition is assumed

The strongly deflected fracture propagation path and the stress distribution in y-direction are shown in Fig. 13. It can be conjectured that if there is no cavity in this example, the hydraulic fracture will propagate horizontally, i.e., in the direction orthogonal to the minimum principal stress σh. However, as can be seen in Fig. 13, because of the strong stress concertation surrounding the cavity, the hydraulic fracture propagates and bypasses the cavity. In the meantime, the injection pressure at the wellbore required to drive the hydraulic fracturing process changes in a significantly different manner from the scenario without the cavity, as shown in Fig. 14. Specifically, higher fluid pressure is required at the beginning of the simulation because of the intense compressive stress (in y-direction) at the left side of the cavity, and the required fluid pressure drops quickly once the fracture grows into the area of low compressive stress above the cavity.

Figure 13: Propagation path of the hydraulic fracture and the stress distribution in y-direction. A positive value indicates compressive stress

Figure 14: Comparison of the injection pressure at the wellbore between the scenarios with and without the cavity

The level of in-situ stress increases with the increase of the reservoir depth. On the other hand, as discussed in the introduction section, the stress field surrounding the cavity is strongly affected by the lateral stress coefficient λ. Hence, three different cases of in-situ stress are investigated in this section, as shown in Tab. 1. In Case 1, we choose a higher level of in-situ stress: σh is 90 MPa and σH is 100 MPa. In Case 2, we choose a lower lateral stress coefficient: σh is 10 MPa and σH is 50 MPa. In Case 3, we choose the same lateral stress coefficient as Case 2 but low in-situ stress difference: σh is 1 MPa and σH is 5 MPa. The simulated fracture propagation paths are presented in Fig. 15. It can be found that both the lateral stress coefficient and the level of confining stress (or in-situ stress difference) have a strong influence on propagation paths. Specifically speaking, a higher lateral stress coefficient can enhance the influence of the natural cavity, causing a more curved fracture path. However, lower confining stress or smaller in-situ stress difference can reduce the influence, and thus contributes to the penetration of the hydraulic fracture towards the cavity.

Figure 15: Comparison of fracture propagation paths for three different cases of in-situ stress

Fluid viscosity is an important factor in hydraulic fracturing treatments. In this section, a fluid of higher viscosity, 0.1 Pa⋅s, is injected into the wellbore. The simulation results are presented in Fig. 16, from which it can be concluded that increasing the fluid viscosity can reduce the impact of the cavity on the propagation path. The main reason for this is that the fluid pressure inside the hydraulic fracture is higher for higher viscosity fluid [38], and the higher fluid pressure can counteract the effect of the stress concertation caused by the cavity.

Figure 16: Comparison of fracture propagation paths for higher fluid viscosity and higher fluid pumping rate with the base case

In this section, we investigate the effect of another key factor, the fluid pumping rate. Different from the base simulation case, a higher pumping rate, 0.01 m2/s, is applied. The simulated propagation path is also shown in Fig. 16 in which a smoother path can be seen compared to the base case. Therefore, it can be concluded that a high fluid pumping rate is able to attenuate the effect of the cavity, and thus benefits the propagation of the hydraulic fracture along its original direction. This phenomenon can be explained with the aid of the analytical solution of the KGD model [38] which reveals that a higher pumping rate results in increased fluid pressure.

Cavities with greatly different scales are widely distributed in fracture-cavity carbonate reservoirs [2]. In this section, we consider a smaller cavity of radius 0.3 m to study the effect of cavity size. The simulated path is illustrated in Fig. 17. Just as expected, for a smaller cavity, its effect is significantly weakened in comparison with the base case. Besides, it is worth noting that the hydraulic fracture does not start to deflect until its tip comes into the zone near enough to the cavity.

Figure 17: Comparison of fracture propagation paths for different sizes of cavities. The red dotted circle represents the smaller cavity

Natural cavities are of quite irregular shapes in fracture-cavity reservoirs [2]. In this section, two elliptical cavities oriented in different directions (a horizontal elliptical cavity and a vertical elliptical cavity) are studied to investigate the effect of cavity shape. The elliptical and circular cavities have the same area. The eccentricity of both ellipses takes a value of 3/2. From Fig. 18 we can conclude that cavities of irregular shape have a stronger influence on the propagation path than regular cavity and the horizontal elliptical cavity leads to a more twisty propagation path. The stronger influence of the irregular cavity is caused by its higher stress concentrations compared to the regular circle cavity.

Figure 18: Comparison of fracture propagation paths for cavities of different shapes

As shown in Fig. 19, the influence of natural fractures will be investigated through four different cases. All of the natural fractures have the same length of 0.6 m and are positioned along the normal direction of the cavity. The friction angle of the fracture surface and the cohesive strength of the natural fracture are π/6 and 0.5 MPa, respectively. kN and kT are taken as 104 GPa/m. In Figs. 19a and 19c, natural fractures are connected to the cavity; In Figs. 19b and 19d, however, natural fractures are disconnected from the cavity with a small distance of 0.1 m. The angle of natural fractures in Cases (a) and (b) is set to π/4 with respect to the x-axes, and in Cases (c) and (d) the angle is taken as -π/4. Simulation results are depicted in Figs. 20 and 21. It can be seen that the connectivity between existing natural fracture and cavity is a very important factor that influences the propagation behavior of the hydraulic fracture. The connected natural fracture dramatically deviates the hydraulic fracture from its original propagation direction. However, for the disconnected natural fracture, the influence is limited since no sliding occurs between the surfaces of the frictional fracture. In Case (d), a T-shaped junction is formed since the hydraulic fracture is arrested by the natural fracture.

Figure 19: Illustration of different cases of natural fracture around the cavity: (a) The natural fracture with an angle of π/4 is disconnected from the cavity with a small distance; (b) The natural fracture with an angle of π/4 is connected to the cavity; (c) The natural fracture with an angle of -π/4 is disconnected from the cavity with a small distance; (d) The natural fracture with an angle of -π/4 is connected to the cavity

Figure 20: Comparison of fracture propagation paths for cavities with connected and disconnected natural fractures: Case (a) and Case (b)

Figure 21: Comparison of fracture propagation paths for cavities with connected and disconnected natural fractures: Case (c) and Case (d). In Case (d), the hydraulic fracture is arrested by the natural fracture

This section is aimed to study the influence of cavity on hydraulic fracturing in a wellbore. As illustrated in Fig. 22, for the sake of simplicity, only the first two stages (Stages 1 and 2) are considered [9]. The size of the model is 100m×100m. The spacing between stages is set to 10 m. The length of the initial hydraulic fracture is 2.5 m. The radius of the cavity denoted by a is set to 2.0 m. The vertical distance between the cavity and the fracture denoted by d is taken as 3 m. The in-situ stresses σh and σH are 20 and 30 MPa (the lateral stress coefficient λ equals 2/3), respectively. The Young’s modulus E is 30 GPa, Poisson’s ratio ν is 0.2, and fracture toughness KIC of the rock formation equals 1.0 MPa⋅m1/2. The viscosity μ and pumping rate Qinj of the injected fluid is 0.01 Pa⋅s and 0.001 m2/s, respectively. The fluid leak-off coefficient CL is taken as 10−5 m/s1/2. The proppant volume concentration of the injected fluid [39] is set to 0.3. The distribution of proppant along the hydraulic fracture is obtained according to the method proposed in our previous papers [9,35]. Then, the width of the propped fracture (wp) can be determined by [9,35] wp=woc/η, where wo denotes the fracture opening before significant pressure drops, c is the proppant volume concentration, and η=0.74 is the proppant packing density [40]. Afterwards, the penalty function method [9] is adopted to simulate the propped fracture, namely, the Stage-1 hydro-fracture in this example. The model is discretized into 10,281 quadrilateral elements. The fractures of both stages propagate by an increment of 0.85 m for each propagation step. It should be noted that simulation without a cavity is also performed for comparison purposes.

The propagation paths of Stages 1 and 2 fractures for both cases without and with cavity are illustrated in Fig. 23. The stress contours in x-direction for cases without and with cavity are shown in Figs. 24a and 24b, respectively. These simulation results reveal that the existing cavity between two sequential fracturing stages significantly influences the propagation paths. For the case without the cavity, the Stage-2 hydraulic fracture curves away from the straight Stage-1 fracture because of the stress shadow effects [41] caused by the propped Stage-1 fracture, and the final paths are symmetrical about the wellbore. For the case with cavity, the Stage-1 fracture deviates leftward for both the upper and lower fracture tips, and the fracture path below the wellbore (y < 0) is more curved than that above the wellbore (y > 0). For the subsequent Stage-2 fracture, its fracture tip above the wellbore (y > 0) grows away from the previous Stage-1 fracture; However, the fracture tip below the wellbore (y < 0) firstly grows towards and then away from Stage-1 fracture. Besides the propagation paths, the stress fields are also quite distinct. The maximum stress values in x-direction for the cases without and with cavity are 36.1 and 43.2 MPa, respectively. For the case shown in Fig. 24b, both the maximum stress and the minimum stress in x-direction occur around the cavity. The maximum stress occurs on the upper and lower sides of the cavity, and the minimum stress occurs on the left and right sides of the cavity. The distinct propagation paths shown in Fig. 23 are a direct result of the stress concentration caused by the cavity.

Figure 22: Illustration of two-stages sequential hydraulic fracturing in a wellbore. The cavity is placed between Stage-1 hydro-fracture and Stage-2 hydro-fracture with a vertical distance d to the wellbore

Figure 23: Propagation paths of Stage-1 fracture and subsequent Stage-2 fracture within a zone of size 50m×50m for cases with and without cavity. The red solid lines represent propagation paths for the case with cavity, and the dotted black lines represent propagation paths for the case without the cavity

Figure 24: Stress contours in x-direction within a zone of size 50m×50m: (a) For the case without cavity; and (b) For the case with cavity. A positive value indicates compressive stress

In this paper, we established a fully-coupled numerical model to investigate the mechanisms of hydraulic fractures in fracture-cavity reservoirs using the XFEM. The Heaviside, cavity, fracture tip, T-shaped junction, and penetration enrichment functions are proposed to describe the displacement jump across the fracture surface, displacement discontinuity over the cavity boundary, singular displacement field near the fracture tip, intersection between the hydro-fracture and the natural fracture, and penetration of a fracture into the cavity, respectively. Hence, tedious remeshing can be avoided. The fluid flow within fractures is described by Reynold’s equation which is discretized using the FEM. Afterwards, the fully-coupled governing equations are solved iteratively using the Newton–Raphson method. After the validation of the proposed model in Section 3, several cases are simulated to investigate the effects of factors such as in-situ stress, fluid viscosity, fluid pumping rate, cavity size and shape, and natural fractures in Section 4.1. Besides, the effects of a cavity on the sequential hydraulic fracturing in a wellbore are studied in Section 4.2. According to the cases studied in this paper, the major conclusions can be reached as follows:

1. Both the lateral stress coefficient and the level of confining stress (or in-situ stress difference) have a strong influence on propagation paths of hydraulic fractures near cavities. A higher lateral stress coefficient can enhance the influence of the natural cavity, causing a more curved fracture path. However, lower confining stress or smaller in-situ stress difference can reduce this influence, and thus contributes to the penetration of the hydraulic fracture towards the cavity.

2. The fluid viscosity and fluid pumping rate are two dominant factors on the propagation path in hydraulic fracturing treatments when natural cavities are considered. Higher fluid viscosity and high fluid pumping rate are both able to attenuate the effect of the cavity, and thus benefit the propagation of the hydraulic fracture along its original direction.

3. The influence of a cavity depends not only on its size but also on its shape. Cavities of irregular shape (ellipse, for example) have a stronger influence on the propagation path of hydraulic fracture than regular circle cavity.

4. The frictional natural fracture disconnected from the cavity, even with a very small distance between the fracture tip and the cavity, has limited influence on the stress field around the cavity. Nevertheless, a frictional natural fracture connected to the cavity can significantly change the stress distribution around the cavity, thus dramatically deviates the hydraulic fracture from its original propagation direction.

5. Natural cavity existing between two adjacent fracturing stages will significantly influence the stress distribution between fractures and is more likely to result in irregular propagation paths compared to the case without cavity.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author on reasonable request.

Funding Statement: This research was jointly funded by the National Natural Science Foundation of China (No. 51904111), the Natural Science Foundation of Jiangsu Province (No. BK20170457), the Open Fund for Jiangsu Key Laboratory of Advanced Manufacturing Technology (No. HGAMTL-1712), and the Natural Science Research of Institution of Higher Education of Jiangsu Province (No. 17KJA460003).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

Khvatova, I. E., Renaud, A., Malyutina, G., Sansiev, G., Kuzilov, I. et al. (2012). Simulation of complex carbonate field: Double media vs. single media Kharyaga field case (Russian). SPE Russian Oil and Gas Exploration and Production Technical Conference and Exhibition, Society of Petroleum Engineers. Moscow, Russia.
Li, Y., Hou, J., & Li, Y. (2016). Features and classified hierarchical modeling of carbonate fracture-cavity reservoirs. Petroleum Exploration and Development, 43(4), 655-662. [Google Scholar] [CrossRef]
Schrefler, B. A., Secchi, S., & Simoni, L. (2006). On adaptive refinement techniques in multi-field problems including cohesive fracture. Computer Methods in Applied Mechanics and Engineering, 195(4–6), 444-461. [Google Scholar] [CrossRef]
Song, C., Chen, Y., & Wang, J. (2019). Experiment and simulation for controlling propagation direction of hydrofracture by multi-boreholes hydraulic fracturing. Computer Modeling in Engineering & Sciences, 120(3), 779-797. [Google Scholar] [CrossRef]
Tang, H., Wang, S., Zhang, R., Li, S., & Zhang, L. (2019). Analysis of stress interference among multiple hydraulic fractures using a fully three-dimensional displacement discontinuity method. Journal of Petroleum Science and Engineering, 179(3), 378-393. [Google Scholar] [CrossRef]
Liu, Z., Lu, Q., Sun, Y., Tang, X., & Shao, Z. (2019). Investigation of the influence of natural cavities on hydraulic fracturing using phase field method. Arabian Journal for Science and Engineering, 44(12), 10481-10501. [Google Scholar] [CrossRef]
Zhou, S., Zhuang, X., & Rabczuk, T. (2019). Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 350(6), 169-198. [Google Scholar] [CrossRef]
Shimizu, H., Murata, S., & Ishida, T. (2011). The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution. International Journal of Rock Mechanics and Mining Sciences, 48(5), 712-727. [Google Scholar] [CrossRef]
Shi, F., Wang, X., Liu, C., Liu, H., & Wu, H. (2016). A coupled extended finite element approach for modeling hydraulic fracturing in consideration of proppant. Journal of Natural Gas Science and Engineering, 33(5), 885-897. [Google Scholar] [CrossRef]
Shi, F., Wang, X., Liu, C., Liu, H., & Wu, H. (2017). An XFEM-based method with reduction technique for modeling hydraulic fracture propagation in formations containing frictional natural fractures. Engineering Fracture Mechanics, 173, 64-90. [Google Scholar] [CrossRef]
Cheng, L., Luo, Z., Yu, Y., Zhao, L., & Zhou, C. (2019). Study on the interaction mechanism between hydraulic fracture and natural karst cave with the extended finite element method. Engineering Fracture Mechanics, 222(9), 106680. [Google Scholar] [CrossRef]
He, B. (2019). Hydromechanical model for hydraulic fractures using XFEM. Frontiers of Structural and Civil Engineering, 13(1), 240-249. [Google Scholar] [CrossRef]
Wang, D., Zlotnik, S., Diez, P., Ge, H., & Zhou, F. (2020). A numerical study on hydraulic fracturing problems via the proper generalized decomposition method. Computer Modeling in Engineering & Sciences, 122(2), 703-720. [Google Scholar] [CrossRef]
Kolawole, O., & Ispas, I. (2020). Interaction between hydraulic fractures and natural fractures: Current status and prospective directions. Journal of Petroleum Exploration and Production Technology, 10(4), 1613-1634. [Google Scholar] [CrossRef]
Taleghani, A. D., Gonzalez, M., & Shojaei, A. (2016). Overview of numerical models for interactions between hydraulic fractures and natural fractures: Challenges and limitations. Computers and Geotechnics, 71(3), 361-368. [Google Scholar] [CrossRef]
Yu, H. S. (2020). Cavity expansion methods in geomechanics. Springer, Dordrecht: Kluwer.
Luo, Z., Zhang, N., Zhao, L., Zeng, J., & Liu, P. (2020). Interaction of a hydraulic fracture with a hole in poroelasticity medium based on extended finite element method. Engineering Analysis with Boundary Elements, 115(4), 108-119. [Google Scholar] [CrossRef]
He, B., & Zhuang, X. (2018). Modeling hydraulic cracks and inclusion interaction using XFEM. Underground Space, 3(3), 218-228. [Google Scholar] [CrossRef]
Wang, H., Tang, X., Luo, Z., Liu, P. (2018). Investigation of the fracture propagation in fractured-vuggy reservoirs. 52nd U.S. Rock Mechanics/Geomechanics Symposium, American Rock Mechanics Association. Seattle, Washington, USA.
Zhao, H., Xie, Y., Zhao, L., Liu, Z., & Li, Y. (2020). Simulation of mechanism of hydraulic fracture propagation in fracture-cavity reservoirs. Chemistry and Technology of Fuels and Oils, 55(6), 814-827. [Google Scholar] [CrossRef]
Al-Nakhli, A. R. (2015). Chemically-Induced pressure pulse: A new fracturing technology for unconventional reservoirs. SPE Middle East Oil & Gas Show and Conference, Society of Petroleum Engineers. Manama, Bahrain.
Wang, Y., Li, X., Zhao, B., & Zhang, Z. (2020). 3D numerical simulation of pulsed fracture in complex fracture-cavitied reservoir. Computers and Geotechnics, 125(6), 103665. [Google Scholar]
Belytschko, T., & Black, T. (1999). Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 45, 601-620. [Google Scholar]
Wang, X., Shi, F., Liu, C., Lu, D., & Liu, H. (2018). Extended finite element simulation of fracture network propagation in formation containing frictional and cemented natural fractures. Journal of Natural Gas Science and Engineering, 50(5), 309-324. [Google Scholar] [CrossRef]
Adachi, J., Siebrits, E., Peirce, A., & Desroches, J. (2007). Computer simulation of hydraulic fractures. International Journal of Rock Mechanics and Mining Sciences, 44(5), 739-757. [Google Scholar] [CrossRef]
Carter, R. D. (1957). Derivation of the general equation for estimating the extent of the fractured area. Drilling and Production Practices, pp. 261–270. Dallas: American Petroleum Institute.
Moës, N., Dolbow, J., & Belytschko, T. (1999). A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46(1), 131-150. [Google Scholar] [CrossRef]
Wrigger, P. (2006). Computational contact mechanics. Second edition. Berlin, Heidelberg: Springer.
Khoei, A. R., & Nikbakht, M. (2007). An enriched finite element algorithm for numerical computation of contact friction problems. International Journal of Mechanical Sciences, 49(2), 183-199. [Google Scholar] [CrossRef]
Erdogan, F., & Sih, G. C. (1963). On the crack extension in plates under plane loading and transverse shear. Journal of Biomechanical Engineering-Transactions, ASME, 85(4), 519-525. [Google Scholar] [CrossRef]
Moran, B., & Shih, C. F. (1987). Crack tip and associated domain integrals from momentum and energy balance. Engineering Fracture Mechanics, 27(6), 615-642. [Google Scholar] [CrossRef]
Gu, H., Weng, X., Lund, J. B., Mack, M. G., & Ganguly, U. (2012). Hydraulic fracture crossing natural fracture at nonorthogonal angles: A criterion and its validation. SPE Production & Operation, 27(1), 20-26. [Google Scholar] [CrossRef]
Chessa, J., & Belytschko, T. (2004). Arbitrary discontinuities in space-time finite elements by level sets and X-FEM. International Journal for Numerical Methods in Engineering, 61(15), 2595-2614. [Google Scholar] [CrossRef]
Nishida, A. (2010). Experience in developing an open source scalable software infrastructure in Japan. Computational Science and Its Applications, vol. 6017, pp. 448–462. Berlin, Heidelberg: Springer.
Shi, F., Wang, X., Liu, C., Liu, H., & Wu, H. (2018). An XFEM-based numerical model to calculate conductivity of propped fracture considering proppant transport, embedment and crushing. Journal of Petroleum Science and Engineering, 167(5), 615-626. [Google Scholar] [CrossRef]
Nguyen, N. T., Bui, T. Q., Zhang, C., & Truong, T. T. (2014). Crack growth modeling in elastic solids by the extended meshfree Galerkin radial point interpolation method. Engineering Analysis with Boundary Elements, 44, 87-97. [Google Scholar] [CrossRef]
Peng, S., Zhang, Z., Mou, J., Zhao, B., & Liu, Z. (2018). Hydraulic fracture simulation with hydro-mechanical coupled discretized virtual internal bond. Journal of Petroleum Science and Engineering, 169(3), 504-517. [Google Scholar] [CrossRef]
Hu, J., & Garagash, D. I. (2010). Plane-strain propagation of a fluid-driven crack in a permeable rock with fracture toughness. Journal of Engineering Mechanics, 136(9), 1152-1166. [Google Scholar] [CrossRef]
Jeon, J., Bashir, M. O., Liu, J., Wu, X. (2016). Fracturing carbonate reservoirs: Acidising fracturing or fracturing with proppants? SPE Asia Pacific Hydraulic Fracturing Conference, Society of Petroleum Engineers. Beijing, China.
Rivers, M., Zhu, D., Hill, A. D. (2012). Proppant fracture conductivity with high proppant loading and high closure stress. SPE Hydraulic Fracturing Technology Conference, Society of Petroleum Engineers. The Woodlands, Texas, USA.
Roussel, N. P., Sharma, M. M. (2011). Strategies to minimize frac spacing and stimulate natural fractures in horizontal completions. SPE Annual Technical Conference and Exhibition, Denver, Colorado, USA.