Numerical simulation of fracture path and nonlinear closure for simultaneous and sequential fracturing in a horizontal well

Computers and Geotechnics 88 (2017) 242–255

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Numerical simulation of fracture path and nonlinear closure for simultaneous and sequential fracturing in a horizontal well Wan Cheng a,⇑, Guosheng Jiang a, Yan Jin b a b

Faculty of Engineering, China University of Geosciences, Wuhan, PR China State Key Laboratory of Petroleum Resources and Prospecting, Beijing, PR China

a r t i c l e

i n f o

Article history: Received 1 June 2016 Received in revised form 28 February 2017 Accepted 30 March 2017

Keywords: Horizontal well Sequential fracturing Hydraulic fracturing Stress interaction Multi-stage fracturing Simultaneous fracturing

a b s t r a c t Reservoir stimulation requires a model to evaluate the fracture path and closure for the simultaneous or sequential propagation of the hydraulic fracture (HF). This paper presents a fluid-solid coupled model to simulate multi-stage HF propagation. A non-linear joint model is proposed to evaluate the fracture closure when the created fractures are elastically propped. HF closure continues until the balance of external stress matches the proppant’s resistance. The reservoir along the horizontal wellbore is not stimulated equally by the multi-stage fracturing. The HFs in the subsequent stage are ‘repelled’ and restrained by the HFs in the previous stage. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Over the past few years, hydraulic fracturing has been proven to be an effective reservoir stimulation technique in developing shale gas resources. One of the main objectives of multi-stage fracturing in a horizontal well is to create transverse fractures or even fracture networks in the subsurface [8,10,28]. To hold the fracture open and facilitate hydraulic conductivity after the hydraulic fracture (HF) closure, proppant is injected in nearly all hydraulic fracturing treatments [24,25]. The opening of a propped fracture causes a reorientation of stress on the surrounding rock, which in turn affects the fracture propagation direction. This phenomenon is often referred to as stress shadow [32,13]. The stress shadow among multiple HFs is often known as stress interaction, which is an important subject in the hydraulic fracturing of horizontal wells [20]. Extensive research has been carried out to evaluate the impact of the stress shadow on the multi-fracture geometry and fracture design in recent years. The impacts of the in situ stress ratio, elastic parameters and net fluid pressure on the fracture spacing were evaluated numerically by a finite element model [13]. Meyer and Bazan [17] used a discrete fracture network model to investigate the numerical correlation between the dimensionless fracture ⇑ Corresponding author at: No. 388 Lumo Road, Wuhan 430074, PR China. E-mail address: [email protected] (W. Cheng). http://dx.doi.org/10.1016/j.compgeo.2017.03.019 0266-352X/Ó 2017 Elsevier Ltd. All rights reserved.

spacing and width. Considering the stress interaction, Olson [21] and Olson and Dahi-Taleghani [22] simulated multi-fracture propagation in a naturally fractured reservoir. The stress interaction around multi-fracture has been identified as a limiting factor in determining the fracture spacing and perforation number [5], or even in choosing different fracturing methods, such as consecutive fracturing, alternative fracturing and zipper fracturing [20,2]. However, these papers did not couple the fluid flow with multi-fracture deformation because they assumed a constant pressure inside the multi-fracture. More recent papers successfully solved the fluid flow during the multi-fracture propagation [29,31,30]. Sesetty and Ghassemi [26,27] provided a series of examples to illustrate the impact of fracture spacing on the expected stimulated reservoir volume (SRV) and the multi-fracture path. The multi-fracture path is also affected by the proppant and fluid pressure in the previously created HF [8]. However, little attention has been given to the HF closure [18,19]. Shiozawa and McClure [24,25] developed a model to describe the fracture closure against proppant after the end of the injection. The residual fracture opening causes the stress interaction in its neighborhood and the nearby HFs after the injection is stopped [20]. The stress interaction and fracture conductivity will decrease drastically if the HF is subjected to the confining stresses that tend to close. Therefore, the residual opening of the HF along with the compressible proppant is an important factor in the success of multi-fracture treatment.

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

243

Nomenclature M Qk

number of HF fluid injection rate of the kth fracture, k = 1, 2, 3, . . ., M, m3/s Nk number of segments in the kth HF N total number of segments of all of the HFs Dsj , Dnj shear and normal DDs, respectively, m Aijss ; Aijsn ; Aijns ; Aijnn coefficients, representing the stresses on the ith segment due to the constant DD on the jth segment rh, rH the minimum and maximum in situ stress, respectively, MPa hi angle of the ith segment, rad pi fluid pressure on the ith segment, MPa Q volumetric flow rate, m3/s n0 and l fluid power-law index and consistency index, respectively H fracture height, m s distance along fracture path, m w fracture width or opening, m p pressure of fracturing fluid flow, MPa Qc total injection rate, m3/s T total injection time, s

Crouch and Starfield [7] proposed a linear joint element (assumed as a linear spring) to describe the subsurface crack with a compressible filling for compression or shear. Similarly, this paper developed a nonlinear joint element (assumed as a nonlinear spring) to evaluate the HF closure after the injection is stopped. To evaluate the influence of propped HF on the multi-fracture geometry and stress interaction in the multi-stage fracturing technology, an HF model coupling the fracturing fluid flow along the fracture and the rock deformation is adopted in this paper.

QT Lk (T) qL CL t0(s) Lf pj

a

½Dsj 0 ; ½Dnj 0 Kn D Dn Ep

qo

½Dsj p ; ½Dnj p Sk E v KIC

rxx, ryy, sxy

accumulate injection volume on time T, m3 length of the kth HF on time T fluid loss volume, m3 fluid loss coefficient, m2.s0.5 time at which HF first reaches the portion s, s half-length of HF, m pressure of the jth segment, MPa weighting coefficient, dimensionless shear and normal DDs before the pump’s shutoff, respectively, m normal stiffness of fracture, MPa/m normal displacement increment, m elastic modulus of the propped HF, MPa porosity of proppant filled HF shear and normal DDs after the pump’s shutoff, respectively, m number of segments in the kth stage Young’s modulus of the formation, MPa Poisson’s ratio of the formation, dimensionless fracture toughness, MPa m0.5 stress components, MPa

tion is an isotropic, linearly elastic medium; (2) the HF is allowed to propagate if the equivalent stress intensity factor exceeds the fracture toughness of rock (e.g., [14,4]; (3) the HF is a vertical fracture with a constant height, and the horizontal plane satisfies the plane strain condition; (4) the fluid flow inside an HF is equivalent to Poiseuille flow (e.g., [33]; (5) the fluid leakoff in the direction normal to the HF surface is given by Carter’s model [3]; (6) the pressure loss along the borehole and perforation is ignored; and (7) the fracture closure [38,39] after pump shutoff is determined by the assumption that the elastic proppant is a nonlinear spring.

2. Numerical model 2.1. Model assumption

2.2. Fracture deformation

A horizontal well is generally fractured in multiple stages, starting from the toe to the heel. We assume that each stage includes several perforation clusters, but each cluster creates only one fracture during the fracturing process (Fig. 1). To capture the characteristics of simultaneous fracture propagation and the stress perturbation in a horizontal well, an HF model is proposed under the following assumptions: (1) the target forma-

Under two-dimensional plane strain conditions, we adopt the well-established displacement discontinuity (DD) method [6] to calculate the fracture deformation. The multiple fractures in the formation are divided into N line segments in total (Fig. 2). Each fracture has Nk (k = 1, 2, 3, . . ., M) segments. Each line segment represents a boundary element. According to the superposition principle, the stress perturbation created by the fracture deformation on

Fig. 1. Simultaneous fracturing in the horizontal well. Qi represents the injection rate of fracturing fluid. M is the number of HFs. M fractures are produced simultaneously.

244

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

Fig. 2. Elements in the simultaneous fractures in a horizontal well. Only half of the entire fractures are shown because of the symmetry along the borehole axis.

each element can be found [6]. The stress equilibrium equation on each segment is given by

8 N X > > > ðAij D j þ Aijsn Dnj Þ ¼  12 ðrh  rH Þ sin 2hi > > < j¼1 ss s N > X > 2 > > ðAijns Dsj þ Aijnn Dnj Þ ¼ pi  rh sin hi  rH cos2 hi > :

the multiple HFs and the volume leaking into formation (e.g., [33,35].

ði; j ¼ 1; 2; . . . ; NÞ; N ¼

M X Nk

ð1Þ

k¼1

j¼1

where Dsj and Dnj represent the shear and normal displacement Aijss ; Aijsn ; Aijns ; Aijnn

discontinuities (DDs), respectively. are the coefficients [7], representing the stresses on the ith segment caused by the constant DD on the jth segment. Note that the physical meaning of Dnj represents the opening of segment j. rh and rH represent the minimum and maximum in-situ stress, respectively; hi is the angle starting anticlockwise from the rh direction to the ith segment. pi is the fluid pressure on the ith segment. Nk (k = 1, 2, 3, . . . M) is the number of segments in the kth HF. 2.3. Flow of the fracturing fluid The power-law fluid flow between two parallel plates is (e.g., [33])


n0 0 0 @p 1 þ 2n0 Qn ¼ 2n þ1 l 0 @s n0 Hn w2n0 þ1

Z

Z

T

Q c ðtÞdt ¼ 0

ð3Þ

T

ðQ 1 ðtÞ þ Q 2 ðtÞ þ . . . þ Q M ðtÞÞdt

k¼1

Lk ðTÞ

Hwds þ

0

k¼M Z X k¼1

0

Lk ðTÞ

Z

T

t 0 ðsÞ

qL ðs; tÞdtds ði ¼ 1; 2; . . . ; NM Þ ð5Þ

Carter’s [3] leak-off model is given as

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qL ¼ 2HC L = t  t 0 ðsÞ

ð6Þ

where CL is the fluid loss coefficient; t0(s) is the time at which HF first reaches the portion s. Basically, the HF width equals zero before the fracture initiation. If the injection rate at the wellbore remains constant and no flow exists at the fracture tip, the boundary conditions and the initial conditions are

Q T ð0; tÞ ¼ Q c ðtÞ wðLf ; tÞ ¼ 0 wðs; 0Þ ¼ 0

The total injection volume during the entire treatment is given by the integration of the injection rate over the injection time.

QT ¼

k¼M Z X

ð7Þ

ð2Þ

where Q is the volumetric flow rate; n’ and l are the fluid powerlaw index and consistency index, respectively; H is the fracture height; s is the distance along the fracture path; w is the fracture width; @p=@s is the pressure gradient of the fracturing fluid flow. We assume that multiple HFs are initially perpendicular to the horizontal wellbore at the perforation-cluster location. Only one fracture is created from each perforation cluster. The injection rate Qc is dynamically divided into injection rates of each fracture (Fig. 1). Ignoring the compressibility of fluid under pressure, the fracturing fluid satisfies

Q c ðtÞ ¼ Q 1 ðtÞ þ Q 2 ðtÞ þ . . . þ Q M ðtÞ

QT ¼

ð4Þ

0

where T is the total injection time. M is the number of HFs (Fig. 1). The global conservation of fluid volume over the entire multifracture is required to be satisfied in order to obtain a unique solution. The global mass conservation requires that the total volume of fluid injected during treatment equals the volume preserved in

where Lf is the half-length of HF; Qc is the half injection rate of the pump. 2.4. Fluid-solid coupling Dnj in Eq. (1) represents the aperture w in Eq. (2). The fluid pressure solved from Eq. (2) must satisfy Eq. (1). The flow equation determines the fluid pressure, while the fluid pressure distribution decides the fracture aperture. Then, the fracture aperture affects the flow equation in turn. Hence, the flow equation and deformation equation cannot be solved separately. This fluid-solid coupled problem gives a system of nonlinear equations. They are solved numerically by the Picard iteration method [34,12,30]. The flow and pressure are solved with a finite difference, while the fracture width has the form of Eq. (1). At each time step, an initial guess of pressure and injection rate for each HF is set as the initial solution vector. The solution vector is given by

 ½X ¼ pj

Qk

T

k ¼ 1; 2; . . . ; M; j ¼ 1; 2; . . . ; N

ð8Þ

where N is the total number of elements in the fracture model; M is the fracture number. In this article, pressure loss along the borehole is neglected, and thus, the pressure is equal along the borehole.

245

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

There are M Qk unknowns and N pj unknowns in the solution vector. Hence, we need (M + N) equations to solve all of the variables. Fracture width w is given by Eq. (1). A residual function vector is given by the combination of Eqs. (2), (3)

½F ¼ ½ F 1

F 2 T

ð9Þ

where

8 0 0 Q nk n0 þ1 1þ2n0 n > @p > 0 > 0 < F 1 ¼ @s  2 l n0 Hn w2n þ1

where [J] is the Jacobian matrix of residual function vector [F]. The finite difference is conducted on Eq. (10) in order to obtain the Jacobian matrix. The initial and new solution vectors are compared to check the convergence. If convergence is not achieved, the pressure distribution is modified as a weighted average by the present and former pressures. ðkþ1Þ

pj ð10Þ

X > > Qi > : F2 ¼ Q c  i¼M i¼1

There are N F1 equations in Eq. (10), while there is only one F2 equation. Since the fluid pressures at the M fracture centers are equal, there are (M  1) equations. If we put these equations together, we have (M + N) equations. The number of equations is equal to the number of unknowns. The fracture width from the rock deformation is obtained by using Eq. (1). A new solution vector is obtained by

Xðkþ1Þ ¼ XðkÞ  ½J1 ½FðXðkÞ Þ k ¼ 0; 1; 2; . . .

ð11Þ

¼ ð1  aÞpj

ðkþ1Þ

þ apj

ðkÞ

j ¼ 1; 2; . . . ; N; 0 < a < 1

ð12Þ

where a is the weighting coefficient. The weighting coefficient is suggested to be from 0.1 to 0.3 [34]. The weighting coefficient mainly decides the convergent velocity of iteration. After many attempts, this model reaches the fastest convergence when a is approximately 0.2. In order to improve the convergence, the analytical or asymptotic solution to a similar problem is suggested as the initial guess of the solution. In addition, after the convergence is achieved, the solution at the current propagation-step is suggested as the initial guess solution of the next propagation-step. After obtaining a new pressure, the process of Eqs. (11), (12) is repeated until the convergence is achieved. The time increment is calculated by the adoption of the global mass balance equation. A

Fig. 3. Workflow of the fluid-solid coupled procedure based on Picard iteration [34,30].

246

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

(a) Pump on

(b) Pump off

Fig. 4. The influence of pump shutoff on HF width.

series of new fracture increments (one for each fracture tip) are added to the fracture segments. Hence, the hydraulic fracture path is built step-by-step. At every time step, the HF length increases by the elementary step added to the end of the crack [1,23]. The fluid-solid coupling is required for every propagation step. When the accumulated time increment reaches the treatment time, the hydraulic fracturing is finished. A workflow of fluid-solid coupling is demonstrated in Fig. 3. The iteration method in Fig. 3 does not locally enforce mass balance (solve mass balance within each element) because F1 in Eq. (10) is a constitutive equation for fluid flow. It only globally enforces mass balance, which means that theoretically, fluid could move instantaneously from one place to another in the model. McClure [16] outlines a method of solving these equations that locally conserves mass at each element. The singularity of boundary integral equation usually makes the iteration difficult to be convergent. However, the iteration method in Fig. 3 is more efficient and convergent because it is reduced by one unknown needed at each element.

normal closure. In this study, we only consider the normal closure of the propped fracture. The DDs of the HF before the pump’s shutoff are labelled as ½Dsj 0 ; ½Dnj 0 (j = 1, 2, 3, . . ., N), which can be easily solved by the fluid-solid coupling method mentioned in Section 2.4. The joint deformation [7] is simplified as a linear relationship with normal stress acting on the joint, and reads

rn ¼ K n DDn

ð13Þ

where Kn represents the normal stiffness [7]; DDn is the normal increment of DD. Physically, DDn represents the closure increment of the joint. However, the joint closure DDn cannot exceed the normal opening Dn. Obviously, the DDn in Eq. (13) has the possibility of being greater than Dn if the compressive stress acting on the fracture surface is large. To avoid the possibility of the closure becoming greater than the aperture, we proposed a tangent relation between the compressive stress and relative closure (Fig. 5). The compressive stress satisfies




2.5. Nonlinear closure of elastically propped HF When the fracturing treatment reaches the target treatment time, the injection pump will shut off. The HF width before the pump’s shutoff is larger than that after the pump’s shutoff (Fig. 4a and b). Before the pump is shut off, the HFs are propped by the fluid pressure. After the pump is shut off and the fluid pressure is released, the HFs are propped by the proppant. The proppant holds the fracture opening and increases the hydraulic conductivity after closure [25]. In low permeability formations, the mechanical closure of HF may take hours or even days until the balance with proppant is achieved [24]. The proppant transport depends on a number of factors such as the proppant density, flow rate of fluid, fluid properties and fracture width. The proppant particles are deposited when the injection is stopped. Since this model is a 2D model based on the plane strain condition, it is difficult to consider the vertical distribution of proppant. Hence, we simplify the proppant transport process. In addition, the ultra-low density, small size proppant is developed for the oil and gas development. The most ideal situation of proppant is uniformly distributed along the vertical profile. As a result, the HFs’ width will shrink until the balance with the stressed proppant is achieved. Compared to the normal DD, the shear DD causes a much smaller influence on the



2 p DDn rn ¼  Ep ðqo Þ tan ; 0 6 DDn < ½Dn 0 p 2 ½Dn 0

ð14Þ

where DDn represents the normal closure of the boundary element induced by the pump’s shutoff. Ep is the elastic modulus of the propped HF. qo represents the porosity of the proppant-filled HF. Obviously, the compressive stress will go to infinity when the relative closure goes to 1 [11,26,27]. If the relative closure goes to zero, Eq. (14) will be simplified as a linear relation,

Fig. 5. Nonlinear joint closure versus the dimensionless compression.

247

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

rn ¼ Ep

DD n ½Dn 0

ð15Þ

(

½Dsj p ¼ ½Dsj 0 þ DDsj ½Dnj p ¼ ½Dnj 0 þ DDnj

We are making a simplified assumption that the fracture has zero frictional strength when the HF walls are in contact. The stress on the proppant consists of the in situ stresses, stress interaction induced by the HF itself and the resistance of the proppant agents. Therefore, the propped HF satisfies

ðj ¼ 1; 2; . . . ; NÞ

ð19Þ

where ½Dsj p ; ½Dnj p represent shear and normal DDs after the pump’s shutoff, respectively.

8 N X > > > fAij ð½D j  þ DDsj Þ þ Aijsn ð½Dnj 0 þ DDnj Þg þ 12 ðrh  rH Þ sin 2hi ¼ 0 > > < j¼1 ss s 0 N > X > 2 > > fAijns ð½Dsj 0 þ DDsj Þ þ Aijnn ð½Dnj 0 þ DDnj Þg þ rh sin hi þ rH cos2 hi ¼  p2 Ep tan p2 > : j¼1

DDin ½Din 0



ð16Þ

ði ¼ 1; 2; . . . ; NÞ where DDs is the shear DD increment. Since ½Dsj 0 ; ½Dnj 0 (j = 1, 2, . . ., N) are solved by the fluid-solid coupling method in Section 2.4,

8 N X > > > fAij ½D j  þ Aijsn ½Dnj 0 g ¼  12 ðrh  rH Þ sin 2hi > > < j¼1 ss s 0 N > X > 2 > > fAijns ½Dsj 0 þ Aijnn ½Dnj 0 g ¼ pi  rh sin hi  rH cos2 hi > :

ði ¼ 1; 2; . . . ; NÞ

ð17Þ

j¼1

Substituting Eq. (17) into Eq. (16) yields

If the relative closure is small, Eq. (18) is simplified as

8 N X > > > fAij DD j þ Aijsn DDnj g ¼ 0 > > < j¼1 ss s

N > X > ij ij j j > 2 p > fA D D þ A D D g ¼ p  E tan > p i s n ns nn p 2 : j¼1

DDin ½Din 0



ði ¼ 1; 2; . . . ; NÞ

8 N X > > > fAij DD j þ Aijsn DDnj g ¼ 0 > > < j¼1 ss s

N > X i > > > fAijns DDsj þ Aijnn DDnj g ¼ pi  Ep ½DDDi n > : n 0

ði ¼ 1; 2; . . . ; NÞ

ð20Þ

j¼1

ð18Þ Eq. (18) is a system of non-linear equations with 2N equations and 2N unknowns. The Newton-Raphson iteration method is used to solve Eq. (18). After the stress balances the proppant with the HF surface, the new DD is given by the superposition of DD before the pump’s shutoff and the DD increment after the pump’s shutoff,

Eq. (20) is a linear equation system with 2N unknowns and 2N equations. 2.6. Multi-stage fracturing A horizontal well is generally fractured in multiple stages to enhance the hydraulic conductivity for oil and gas flowing from the reservoir to the wellbore. One stage is fractured at a time, starting from the toe to the heel (Fig. 6).

Fig. 6. Multi-stage fracturing in a horizontal well (starting from the left to the right).

248

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

Multi-fracturing propagates simultaneously in each stage, but the multi-stage HFs propagate sequentially. The fractures’ propagation at the first stage is affected by the in situ stresses, fluid pressure and stress interaction. The fractures’ propagation in the subsequent stage is affected by the in situ stresses, fluid pressure, re-distributed stress in the previous stage and the stress interaction. Hence, we still use the fluid-solid coupling procedure in Section 2.4 to simulate the simultaneous fracture propagation in the 2nd stage but consider the induced stress from the fractures in the 1st stage. Each stage has Sk (k = 1, 2, 3, . . .) elements in total. The boundary element equation of the 2nd stage fracture is given by

j¼S1 þ1

(

½Dsj p ¼ ½Dsj 0 þ DDsj ½Dnj p ¼ ½Dnj 0 þ DDnj

S 1 þ 1 6 j 6 S1 þ S 2

ð25Þ

Similarly, all of the fractures in the previous stages should be considered in the propagation of fractures in the kth stage, 8 S þS þ...þS S1 þS2X þ...þSk1 1 X 2 k > > > ðAijss Dsj þ Aijsn Dnj Þ ¼ ris  fAijss ½Dsj p þ Aijsn ½Dnj p g > > < j¼S þS þ...S þ1 j¼1 1

> > > > > :

8 SX S1 1 þS2 X > > ij ij j j i > fA D þ A D g ¼ r  fAijss ½Dsj p þ Aijsn ½Dnj p g > s ss s sn n > > < j¼S1 þ1 j¼1 > SX S1 1 þS2 > X > > > fAijns Dsj þ Aijnn Dnj g ¼ rin  fAijns ½Dsj p þ Aijnn ½Dnj p g > :

After the stress balances the proppant with the HF surface, the new DDs of the 2nd stage are given by using the superposition of the DDs before the pump’s shutoff and the increment of the DDs after the pump’s shutoff,

2

k1

S1 þSX 2 þ...þSk

ðAijns Dsj þ Aijnn Dnj Þ ¼ rin 

S1 þS2X þ...þSk1

j¼S1 þS2 þ...Sk1 þ1

fAijns ½Dsj p þ Aijnn ½Dnj p g

j¼1

ðS1 þ S2 þ . . . þ Sk1 þ 1 6 i 6 S1 þ S2 þ . . . þ Sk Þ ð26Þ

ð21Þ

Coupling Eq. (26) with fluid flow equations in Section 2.3, the DDs in the kth stage can be solved numerically and are labelled as

j¼1

ðS1 þ 1 6 i 6 S1 þ S2 Þ Coupling Eq. (21) with the fluid flow equation in Section 2.3, the DDs in the 2nd stage can be solved and are labelled as ½Dsj 0 ; ½Dnj 0 (j = S1 + 1, S1 + 2, . . ., S1 + S2). When the stimulation of the 2nd stage fracturing reaches the target time, the injection is stopped. The stress on the proppant consists of in-situ stresses, stress interaction induced by the HF itself and the resistance of the proppant agents. Therefore, the propped HFs in the 2nd stage satisfy

When ½Dsj 0 ; ½Dnj 0 ðS1 þ S2 þ . . . þ Sk1 þ 1 6 j 6 S1 þ S2 þ . . . þ Sk Þ. stimulation of the kth stage fracturing reaches the target time, the injection will be stopped. The closure of the kth stage can be solved by

8 > > > > > < > > > > > :

S1 þSX 2 þ...þSk

ðAijss DDsj þ Aijsn DDnj Þ ¼ 0

j¼S1 þS2 þ...þSk1 þ1 S1 þSX 2 þ...þSk j¼S1 þS2 þ...þSk1 þ1



i DDin ðAijns DDsj þ Aijnn DDnj Þ ¼ pi  p2 Ep tan p2 ½D i 

ð27Þ

n 0

ðS1 þ S2 þ .. . þ Sk1 þ 1 6 i 6 S1 þ S2 þ . .. þ Sk Þ

8 SX 1 þS2 > > > fAij ð½D j  þ DDsj Þ þ Aijsn ð½Dnj 0 þ DDnj Þg þ 12 ðrh  rH Þ sin 2hi ¼ 0 > > < S þ1 ss s 0 1

SX > 1 þS2 > 2 > > fAijns ð½Dsj 0 þ DDsj Þ þ Aijnn ð½Dnj 0 þ DDnj Þg þ rh sin hi þ rH cos2 hi ¼  p2 Ep tan p2 > : S1 þ1

DDin ½Din 0



ð22Þ

ðS1 þ 1 6 i 6 S1 þ S2 Þ where DDsj ; DDnj (j = S1 + 1, S1 + 2, . . ., S1 + S2) are the normal and shear DDs. The physical meaning of DDn represents the normal closure of a propped HF induced by the pump’s shutoff. Since ½Dsj 0 ; ½Dnj 0 (j = S1 + 1, S1 + 2, . . ., S1 + S2) are solved by the fluid-solid coupling method,

8 S þS 1 2 X > > > fAij ½D j  þ Aijsn ½Dnj 0 g ¼  12 ðrh  rH Þ sin 2hi > > < j¼S þ1 ss s 0 1

SX > 1 þS2 > 2 > > fAijns ½Dsj 0 þ Aijnn ½Dnj 0 g ¼ pi  rh sin hi  rH cos2 hi > :

After the stress balances the proppant with the HF surface, the new DDs of the kth stage are given using the superposition of the DDs before the pump’s shutoff and the DD increment after the pump’s shutoff,

ðS1 þ 1 6 i 6 S1 þ S2 Þ

ð23Þ

j¼S1 þ1

Substituting Eq. (23) into Eq. (22) yields

8 S þS 1 2 X > > > fAij DD j þ Aijsn DDnj g ¼ 0 > > < j¼S þ1 ss s 1

SX > 1 þS2 > > > fAijns DDsj þ Aijnn DDnj g ¼ pi  p2 Ep tan p2 > : j¼S1 þ1

DDin ½Din 0



ðS1 þ 1 6 i 6 S1 þ S2 Þ

ð24Þ

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

(

½Dsj p ¼ ½Dsj 0 þ DDsj ½Dnj p ¼ ½Dnj 0 þ DDnj

249

ðS1 þ S2 þ ... þ Sk1 þ 1 6 j 6 S1 þ S2 þ . .. þ Sk Þ ð28Þ

2.7. Model verification 2.7.1. Single fracture The KGD model [15,9] requires the horizontal plane to satisfy the plane strain condition. This model also assumes that the horizontal plane satisfies the plane stain condition. Hence, we verify this model with the analytical solution of the KGD model using the parameters measured from the shale gas field in the Longmaxi

Table 1 Input parameters measured from the shale gas field in the Longmaxi formation. Parameter

Value

Parameter

Value

rh rH

48.97 MPa 54.19 MPa 51.3 GPa 0.26 1

H

100 m 10 mPa s 2 MPa m0.5 15 m3/min 0

E v n0

l KIC Qc CL

Fig. 8. Comparisons of leakoff (CL = 0.00005 m/s0.5) and no leakoff (CL = 0) in this model.

formation, Sichuan basin, China (in Table 1). Using the given parameters, both the fracture length and width obtained from this model match well with the results from the KGD model (Fig. 7).

(a) Width comparison (Element length is 0.1m)

(b) Fracture-length comparison (Element length is 0.1m) Fig. 7. Comparisons of this model and the KGD model (Element number increases during the propagation).

250

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

(a) Two fractures in one horizontal well

(b) Two horizontal wells with one fracture each

Fig. 9. Two initiation fractures in a horizontal well [36,37]

(a) Two fractures in one horizontal well

(b) Two horizontal wells with one fracture each

Fig. 10. Comparison of fracture path [36,37]

Since the KGD model did not consider the leakoff normal to the HF surface, this model indicates that the leakoff is important if the rock is permeable (Fig. 8). Leakoff of fracturing fluid into the formation results in a lower fracture propagation velocity and a narrower fracture width. However, neglecting the leakoff will result in a simpler form of Eq. (5) and a stable iteration procedure. Note that neglecting the leakoff is still acceptable when the reservoir has ultra-low permeability [26,27,36,37].

3. Simulation example In this part, we still use the input parameters in Table 1 (the loss coefficient excluded). The loss coefficient is 5.0  1010 m/s0.5. Simultaneous fracturing for five fractures and sequential fracturing for multi-stage fractures are given as two examples in Sections 3.1 and 3.2, respectively. 3.1. Simultaneous propagation of five fractures

2.7.2. Double fracture Since this model is designed to simulate the multiple fracture propagation, the verification of a single fracture propagation with the KGD model is technically insufficient. Therefore, we further compared this model with Zhang’s model [36,37] for the fracture geometry of a horizontal well with two initial fractures (Fig. 9a) and two horizontal wells with one initial fracture each (Fig. 9b). Numerical simulations from this model present similar fracture paths as Zhang’s model for both scenarios (Fig. 10).

Five fractures propagate simultaneously with fracture spacing of 30 m. The fluid pressure distribution in each fracture and the fracture path before the pump’s shutoff are given in Fig. 11. The fluid pressure decreases sharply on the segment of fracture tip. The geometry of five fractures has two symmetrical lines: one is the wellbore axis and the other is along the central HF. The outside fractures (Frac 1 and Frac 5) are repelled but the inside fractures (Frac 2 and Frac 4) are attracted from the central fracture (Frac 3). In addition, the outside fractures are much longer than the

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

251

(a) Fluid pressure distribution along the fracture path

(b) Fluid pressure distribution along the y-axis Fig. 11. Fluid pressure distribution in the hydraulic fracture before pump shutoff. The colored lines represent the fracture paths, and the color on the path represents the pressure value with respect to the color bar in the right side. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

inside fractures. The inside fractures stop propagating because of the stress interaction. The geometry of five fractures indicates that the inside area of the reservoir is not stimulated as well as the outside area. The outside HFs (Fracs 1 and 5) are much wider than the inside ones (Fig. 12a and b). HF width distribution implies the complexity of proppant design in simultaneous fracture propagation, which is totally different from single fracture propagation. The width fluctuation in the range from y = 10 to y = 30 is a numerical instability during the Picard iteration. Since the path and width of the HFs are both distributed symmetrically, the stress field is also symmetrical (Fig. 13). Stress interaction among multiple HFs determined the propagation direction of the HF, and the HF propagation releases the stress interaction in turn. The dynamic compatibility in each propagation step between the stress interaction and HF path results in the symmetrical geometry and stress field. After the injection is stopped, the fractures will close due to the fluid leakoff. Further closure will occur when the pressure is drawn down during production. The HF will shrink until a balance with the resistance of proppant agents is achieved. With respect to the model in Section 2.5, the closure of each fracture segment can be solved. The relative closure (the closure over the fracture width before the pump’s shutoff) in each HF is given in Fig. 14.

3.2. Multi-stage fracturing A horizontal well is separated into 5 stages with each stage comprising 3 perforation clusters. The fracture spacing (equal to perforation spacing) is 30 m and the stage spacing is 40 m. The elastic modulus of the proppant is set as equal to the elastic modulus of formation. The fluid pressure and fracture geometry in all of the stages are shown in Fig. 15. The geometry of fracture and fluid pressure distribution are symmetrical along the wellbore axis. The fracture geometry in stage 2 is quite different from that in stage 1 but similar to the other three stages. Frac 3 in the subsequent stage is much longer than Fracs 1 and 2. This geometry implies that the reservoir is not stimulated equally. The HFs in the subsequent stage are ‘repelled’ and restrained by HFs in the previous stage. The fracture width distribution is symmetrical along the wellbore axis (Fig. 16). Frac 3 in the subsequent stage is much wider than the other fractures. Since the induced stress of HF on the neighboring HF decreases drastically as the distance between them increases, Frac 3 in the subsequent stage experiences a smaller compressive stress from the previous stage than Fracs 1 and 2. As the injection continues, Fracs 1 and 2 in the subsequent stage are restrained by Frac 3 and the HFs in the previous stage. This is why Frac 3 in the subsequent stage is much wider and longer than Fracs 1 and 2. It is believed that Frac 3 in the subsequent stage

252

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

(a) Fracture width along the HF path

(b) Fracture width along y-axis Fig. 12. Fracture width distribution in the hydraulic fracture before pump shutoff. Colored lines represent the fracture paths, and the color on the path represents the width value according to the color bar in the right side. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

absorbs most of the fluid volume. The HF width distribution implies the complexity of the proppant design in a multi-stage fracturing of horizontal well. Since the path and width of the HFs are both distributed symmetrically along the wellbore axis, the stress field is also symmetrical (Fig. 17). The rxx among the multiple HFs is much greater than that in other areas. The stress field among the multi-stage fractures demonstrates the stress interaction among them. This HF model can simulate the multi-stage fracturing and find the stress interaction strength at each time step. 4. Conclusions This paper presents an HF model that couples the flow of fracturing fluid and rock deformation through the Picard iteration method. A displacement discontinuity (DD) model with a nonlinear joint element is used in simulating the closure of the propped HF in multi-stage fracturing. The verification study against the published results indicates that this model reliably

tracks the evolution of fracture width and length during single fracture propagation and predicts the fracture propagation path during double fracture propagation. This model can simulate multi-stage sequential fracturing in horizontal wells and can determine the stress interaction among multiple HFs. Some interesting findings include the following: (1) Leakoff of fracturing fluid into the formation results in a lower fracture propagation velocity and a narrower fracture width, which is not negligible if the rock is relatively permeable. (2) The coupled DD model can analyse the complex behavior of the fluid transport along the HF and non-linear HF closure using the joint model. The HF geometry in simultaneous fracturing is very complicated and different from the single fracture model because of the stress interaction among the HFs. Dynamic compatibility in each propagation step between the stress interaction and HF path leads to the symmetrical geometry and stress field.

253

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

(1)

xx

(MPa)

(2)

yy

(MPa)

(3)

xy (MPa)

Fig. 13. Stress field around five fractures before pump shutoff.

Fig. 14. Relative closure along the HF path due to the pump shutoff. (Ep = 0.1 * E).

Fig. 15. Fluid pressure distribution before the pump’s shutoff in multi-stage fracturing (Fracturing process is from left to right, the same order as in Fig. 6).

254

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255

Fig. 16. Final width distribution in a multi-stage fracturing.

1

xx

(MPa)

2

yy

(MPa)

Fig. 17. Stress field and fracture geometry of multi-stage fracturing in a horizontal well.

(3) The reservoir area along the horizontal well is not stimulated equally by the multi-stage fracturing. The fracture in the right side of the subsequent stage is much longer and wider than the other fractures. The HFs in the subsequent stage are repelled and restrained by the HFs in the previous stage. The HF width distribution implies the complexity of the proppant design in multi-stage fracturing of horizontal wells.

Acknowledgements This research was supported by the China National High Technology Research and Development Program 863 (No. 2013AA064503). The authors wish to thank Dr. Xiaojing Wang for English text editing of this article and Dr. Cheng Zhu for providing so many excellent suggestions.

References [1] Andreev AA, Galybin AN, Izvekov OY. Application of complex SIE method for the prediction of hydrofracture path. Eng Anal Boundary Elem 2015;50:133–40. [2] Bunger AP, Zhang X, Jeffrey RG. Parameters affecting the interaction among closely spaced hydraulic fractures. SPE J 2012;17(01):292–306. [3] Carter RD. Derivation of the general equation for estimating the extent of the fractured area. Appendix I of ‘‘optimum fluid characteristics for fracture extension”. In: Howard GC, Fast CR, editors. Drilling and production practice. New York (New York, USA): American Petroleum Institute; 1957. p. 261–9. [4] Cheng W, Jin Y, Chen M. Reactivation mechanism of natural fractures by hydraulic fracturing in naturally fractured shale reservoirs. J Natural Gas Sci Eng 2015;23:431–9. [5] Cheng Y. Impacts of the number of perforation clusters and cluster spacing on production performance of horizontal shale gas wells. SPE Reservoir Eval Eng 2012;15(01):31–40. [6] Crouch SL. Solution of plane elasticity problems by the displacement discontinuity method. Int J Numer Meth Eng 1976;10:301–43. [7] Crouch SL, Starfield AM. Boundary element methods in solid mechanics. London: Goerge Allen and Unwin Publishers; 1983.

W. Cheng et al. / Computers and Geotechnics 88 (2017) 242–255 [8] Damjanac B, Cundall P. Application of distinct element methods to simulation of hydraulic fracturing in naturally fractured reservoirs. Comput Geotech 2015;71:283–94. [9] Geertsma J, De Klerk F. A rapid method of predicting width and extent of hydraulically induced fractures. J Petrol Technol 1969;21(12):1571–81. [10] Gerolymatou E, Galindo-Torres SA, Triantafyllidis T. Numerical investigation of the effect of preexisting discontinuities on hydraulic stimulation. Comput Geotech 2015;69:320–8. [11] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. J Soil Mech Found Div 1968;94:637–59. [12] Gordeliy E, Peirce A. Coupling schemes for modeling hydraulic fracture propagation using the XFEM. Comput Methods Appl Mech Eng 2013;253:305–22. [13] Morrill JC, Miskimins JL. Optimizing hydraulic fracture spacing in unconventional shales. In: SPE hydraulic fracturing technology conference, 6–8 February, The Woodlands, Texas, USA. [14] Kim JH, Paulino GH. On fracture criteria for mixed-mode crack propagation in functionally graded materials. Mech Adv Mater Struct 2007;14:227–44. [15] Khristianovic SA, Zheltov YP. Formation of vertical fractures by means of highly viscous liquid. In: Proceedings of the fourth world petroleum congress, Rome. [16] McClure MW. Modeling and characterization of hydraulic stimulation and induced seismicity in geothermal and shale gas reservoirs PhD Thesis. Stanford (California): Stanford University; 2012. [17] Meyer BR, Bazan LW. A discrete fracture network model for hydraulically induced fractures theory, parametric and case studies. In: SPE hydraulic fracturing technology conference and exhibition, 24–26, January, the Woodlands, Texas, USA.. [18] Neto LB, Kotousov A. Residual opening of hydraulically stimulated fractures filled with granular particles. J Petrol Sci Eng 2012;100:24–9. [19] Neto LB, Kotousov A. Residual opening of hydraulic fractures filled with compressible proppant. Int J Rock Mech Min Sci 2013;61:223–30. [20] Roussel NP, Sharma MM. Strategies to minimize frac spacing and stimulate natural fractures in horizontal completions. In: SPE annual technical conference and exhibition, 30 October-2 November, Denver, Colorado, USA. [21] Olson JE. Multi-fracture propagation modeling: applications to hydraulic fracturing in shales and tight gas sands. In: ARMA, The 42th US rock mechanics symposium, June 29–July 2, San Francisco, USA. [22] Olson JE, Dahi-Taleghani A. Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures.. In: SPE hydraulic fracturing technology conference and exhibition, 19–21 January, The Woodlands, Texas. [23] Shi J, Shen B, Stephansson O, Rinne M. A three-dimensional crack growth simulator with displacement discontinuity method. Eng Anal Boundary Elem 2014;48:73–86.

255

[24] Shiozawa S, McClure MW. Simulation of proppant transport with gravitational settling and fracture closure in a three-dimensional hydraulic fracturing simulator. J Petrol Sci Eng 2016;138:298–314. [25] Shiozawa S, McClure MW. Comparison of pseudo-3d and fully-3d simulations of proppant transport in hydraulic fractures, including gravitational settling, formation of proppant banks, tip-screen out, and fracture closure. In: SPE hydraulic fracturing technology conference, 9–11 February, The Woodlands, Texas, USA. [26] Sesetty V, Ghassemi A. A numerical study of sequential and simultaneous hydraulic fracturing in single and multi-lateral horizontal wells. J Petrol Sci Eng 2015;132:65–76. [27] Sesetty V, Ghassemi A. Simulation of Simultaneous and Zipper Fractures in Shale Formations. In: 49th US rock mechanics/geomechanics symposium. american rock mechanics association. 28 June–1 July, San Francisco, California, USA. [28] Dahi-Taleghani A, Gonzalez M, Shojaei A. Overview of numerical models for interactions between hydraulic fractures and natural fractures: challenges and limitations. Comput Geotech 2015;71:361–8. [29] Wu K, Olson JE. Investigation of critical in situ and injection fractors in multifrac treatments: guidelines for controlling fracture complexity. In: SPE hydraulic fracturing technology conference, 4–6 February, Woodlands, Texas, USA.. [30] Wu K. Numerical modeling of complex hydraulic fracture development in unconventional reservoirs PhD thesis. Texas (USA): University of Texas at Austin; 2014. [31] Wu K, Olson JE. Simultaneous multifracture treatments: fully coupled fluid flow and fracture mechanics for horizontal wells. SPE J 2014;20(2):334–46. [32] Wu R, Kresse O, Weng X, et al. Modeling of interaction of hydraulic fractures in complex fracture networks. In: SPE hydraulic fracturing technology conference. 6–8 February, The Woodlands, Texas, USA. [33] Valko P, Economides MJ. Hydraulic fracture mechanics. New York: Wiley; 1995. [34] Yew CH. Mechanics of hydraulic fracturing. Houston (Texas): Gulf Publishing Company; 2014. [35] Zhang X, Jeffrey RG, Bunger AP, et al. Initiation and growth of a hydraulic fracture from a circular wellbore. Int J Rock Mech Min Sci 2011;48:984–95. [36] Zhang Z, Li X, He J, et al. Numerical analysis on the stability of hydraulic fracture propagation. Energies 2015;8:9680–877. [37] Zhang Z, Li X, Yuan W, et al. Numerical analysis on the optimization of hydraulic fracture networks. Energies 2015;8:12061–79. [38] Zhu C, Arson C. A thermo-mechanical damage model for rock stiffness during anisotropic crack opening and closure. Acta Geotech 2014;9:847–67. [39] Zhu C, Pouya A, Arson C. Prediction of viscous cracking and cyclic fatigue of salt polycrystals using a joint-enriched finite element model. Mech Mater 2016;103:28–43.