A mathematical model of fracturing fluid leak-off based on dynamic discrete grid system

Abstract

Based on the classical PKN two-dimensional fracture propagation mathematical model, the two-dimensional leak-off model of fracturing fluid of fractured dual-medium reservoir is established by considering the time-varying non-Newtonian fracturing fluid leak-off coefficient in the stretching process of fractures. Using the finite element difference method, a dynamic discrete grid system is established and solved by Newton-Raphson iterative method. At the same time, the effect on fracturing fluid leak-off of the fractured reservoir stress sensitivity coefficient, the pumping rate, and the propagating length of the fractures is analyzed. As it is analyzed, under the combined effect of the formation pressure, the fracture pressure, the edge effect, and the fracture permeability, the greater the stress sensitivity coefficient is, the smaller the leak-off rate and coefficient are. However, the greater pumping rate is, the larger leak-off rate and coefficient are. If both of them increase to a certain value, the leak-off coefficient firstly decreases , and then increases; the longer the fracture is, in the same position, the larger fracturing fluid leak-off coefficient is; and the greater boundary effect is, the larger fracturing fluid leak-off coefficient near the pinch point is.

http://crossmark.crossref.org/dialog/?doi=10.1007/s13202-015-0188-4&domain=pdf

Web End = J Petrol Explor Prod Technol (2016) 6:343349 DOI 10.1007/s13202-015-0188-4

http://crossmark.crossref.org/dialog/?doi=10.1007/s13202-015-0188-4&domain=pdf

Web End = ORIGINAL PAPER - PRODUCTION ENGINEERING

A mathematical model of fracturing uid leak-off based on dynamic discrete grid system

Yaokun Yang1,2 Hongwei Jiang1 Mu Li1,2 Shuai Yang1,2 Gang Chen1,2

Received: 17 November 2014 / Accepted: 31 May 2015 / Published online: 23 June 2015 The Author(s) 2015. This article is published with open access at Springerlink.com

Abstract Based on the classical PKN two-dimensional fracture propagation mathematical model, the two-dimensional leak-off model of fracturing uid of fractured dual-medium reservoir is established by considering the time-varying non-Newtonian fracturing uid leak-off coefcient in the stretching process of fractures. Using the nite element difference method, a dynamic discrete grid system is established and solved by NewtonRaphson iterative method. At the same time, the effect on fracturing uid leak-off of the fractured reservoir stress sensitivity coefcient, the pumping rate, and the propagating length of the fractures is analyzed. As it is analyzed, under the combined effect of the formation pressure, the fracture pressure, the edge effect, and the fracture permeability, the greater the stress sensitivity coefcient is, the smaller the leak-off rate and coefcient are. However, the greater pumping rate is, the larger leak-off rate and coefcient are. If both of them increase to a certain value, the leak-off coefcient rstly decreases , and then increases; the longer the fracture is, in the same position, the larger fracturing uid leak-off coefcient is; and the greater boundary effect is, the larger fracturing uid leak-off coefcient near the pinch point is.

Keywords Low-permeability reservoirs Filtration

coefcient Fracture propagation Stress sensitivity

Dynamic grid

Introduction

Low-permeability reservoirs need to be developed with fracturing technology (Balen et al. 1988; Demarchos and Chomatas 2004; Fan and Economides 1995); due to the presence of natural fractures in the reservoirs, the conventional homogeneous reservoir fracturing uid leak-off model is no longer applicable (Settari 1985; Yi and Penden 1993; Settari 1998). For the double-porosity reservoir uid leak-off calculation model, many scholars have studied it. Considering single permeability and dual permeability (Mayerhofer et al. 1991; Nghiem et al. 1984), some established one-dimensional fracturing uid leak-off models for fractured dual-medium rese rvoirs, and then considering the actual situation that fracturing uid leak-off in the formation is two-dimensional owing uid and the fracturing uid is non-Newtonian uid, they established the two-dimensional model of non-Newtonian fracturing uid leak-off, which improved the pressure fracturing uid leak-off model, rendering the results more in line with the actual situation. However, the above models are established on the basis of the situation that the fractures do not extend after the pump stops and the pressure distributes evenly in the fracture, and combining the actual mineral conditions, the fracturing uid leak-off also exists in the propagation process. For the propagation of the fractures, the classical two-dimensional PKN and KGD models and the pseudo-three dimensional (P3D) model are mainly used to conduct the simulation (Al-Shatri et al. 2009; Ouenes and Hartley 2000), and we can make the assumption that the fracturing uid ltration coefcient is constant. However, the leak-off process does not accord with the classic Carter leak-off model. The actual ltration coefcient changes with time and is associated with the uid owing process among the

& Yaokun [email protected]

1 China National Petroleum Corporation Drilling ResearchInstitute, Beijing 102206, China

2 College of Petroleum Engineering, China University of Petroleum, Beijing 102249, China

123

344 J Petrol Explor Prod Technol (2016) 6:343349

Kf lef

oPfoy

formation and the fractures (Economides and Demarchos 2008; Gidley et al. 1989)).

At present, the main point is the stress sensitivity when the formation pressure decreases in production. Therefore, the fracturing uid leak-off makes the reservoir pressure increase, and the permeability of the stress sensitive reservoir may increase. So it is necessary to discuss the increasing volume of leak-off uid. This paper combines the dual-porosity formation ow equations with PKN two-dimensional fracture-stretching model, and considering the ltration coefcient dynamic variation during the fracture-stretching process, and the effect of the reservoir stress sensitivity, we will establish two-dimensional fracturing uid leak-off model. Compared with the actual site, the results of simulation get better adaptability.

Mathematical model

Formation

According to the Darcy law, when fracturing uid ows in low-permeability reservoirs fractures and matrix, the differential equations are derived asFlowing in the fractures:

r vf a

Kmlem Pm Pf /Ctf

p

t sx

Cx

; 0\x\Lt; y 0

oPfoy 0;

Lt\x\xe; y 0

6

oPfoy 0; 0\x\xe; y ye

oPfox 0;

7

x 0; 0\y\ye

oPfox 0; x xe; 0\y\ye: 8

Articial fractures

The classic PKN two-dimensional model is used to simulate the propagation of fractures in the formation, and the fracture height is constantly equal to the effective reservoir thickness. While considering the fracturing uid leaking off into the formation, the ltration coefcient is not the same at different locations and changes with time.

Considering the articial fracture vertical faces strain and combining with the England and Green equation, we can derive the fracture width equation:

Wx; t

oPfot 1

Flowing in the matrix:

r vm a

Kmlem Pm Pf /Ctm

1 vHPFx; t rhG : 9

Considering the articial fracture vertical prole as oval and the fracturing uid as non-Newtonian uid, when the fracturing uid ows in the fracture, we can obtain the pressure drop equation as

oPFox

64 p

oPmot : 2

Wherein, when the fracturing uid is considered as non-Newtonian uid, the equation of motion can be derived as

vi

Kilei rPi; i f ; m 3

n 1: 10

when the fracturing uid ows in the articial fracture, the continuity equation can be expressed as

oq ox

pH 4

q W3H Kn

2n 1

3n

n

6q HW2

oWot

2CH

p 0 11

t s

3n 1

8n

/i 8Ki

1 n

2n

Initial condition:

Wx; t 0; x [ Lt 12 Boundary condition:

Wx; 0 0 13 oW4ox

x0

lei

2kn

rPi

; i f ; m: 4

For the reservoirs of which the depth is relatively shallow and the deformation of the skeleton particles is obvious, such as the coalbed methane reservoirs, when the fracturing uid ltrates into the reservoirs, the formation pressure increases, which may cause fractures or matrix porosity to expand and the permeability to increase. For the low-permeability reservoir characteristics, only the stress sensitivity is considered, namelyInitial condition:

Kf Kfie bPi Pf 5 Pjx; y; 0 Pi; j f ; m

Boundary condition

1n

n 1

n

2561 vQ

pG ; x 0: 14

The model solution

Formation

Combined with Eqs. (1)(5), and considering the quasi steady-state channeling in matrix and the two-dimensional ow of the fracturing uid in the formation, the synthesis

123

J Petrol Explor Prod Technol (2016) 6:343349 345

Fig. 1 Dynamic grid

ow equation in fracture and matrix can be derived asFracture:

r

Kfie bPi Pf

lef rPf

! a

Kmlem Pm Pf /Ctf

oPf ot

15

Matrix:

a Km

oPmot : 16

Equation (15) shows a strong non-linear characteristic, and it is difcult to the analytical solution. Therefore, the nite-difference distribution combined with the Newton Raphson iterative method is used to solve the equation. With the fracturing proceeding, the fractures keep propagating, which makes it difcult to divide meshes. So a dynamic discrete grid is established in Fig. 1, which works as follows: according to the increase of the fractures, the number of meshes is increased and all the time steps are assumed to be equal to T0. During the rst period T0, the fracture length is LF1, and the entire area is divided into the two grids in the fracture crack orientation. The length of the grids are, respectively, LF1 and (Xe - LF1). During the

second period T0, the length of the fracture increases by LF2, and the entire region is divided into three grids in the direction of the fracture, and the length of the grids are, respectively, LF1, LF2, and (Xe - LF1 - LF2). We can nd that the pressure of the grid of LF2 during the last period is the pressure on the grid of (Xe - LF1) during the last period. It can be successively obtained according to this method that during the n period of T0, the length of the fracture increases by LFn, and the whole region is divided into (n ? 1) grids in the direction of the fracture, the length of which are, respectively, LF1, LF2,, LFn, and [Xe - (LF1

? LF2 ? ? LFn)] and the pressure on the grid LFn during

this period is equal to the pressure on the last grid of [Xe - (LF1 ? LF2 ? ? LFn-1) ] during last period.

Dynamic meshes dividing diagram is shown as follows in Fig. 2:

Now the differential discretization is conducted for the formulas (15) and (16), which is the process of simple discretization. We regard the fracture and the matrix of channeling as the same grid, while we can make use of the pressure value of last period to calculate the effective crack permeability and effective viscosity of this period. Then we can get the nal differential equation as follows:

ai;jPN1fi;j 1 bi;jPN1fi 1;j ci;jPN1fi;j di;jPN1fi1;j

ei;jPN1fi;j1 gi;j:

lem Pm Pf /Ctm

17

Among them:

2Kfie bPi P

Nf i;j

ai;j

Dyjlefj 1=2Dyj Dyj 1

bi;j

2Kfie bPi P

Nf i;j

Dxilefi 1=2Dxi Dxi 1

2Kfie bPi P

Nf i;j

Dxilefi1=2Dxi Dxi1

ei;j

di;j

2Kfie bPi P

Nf i;j

Dyjlefj1=2Dyj Dyj1

Fig. 2 Dynamic grid schematic diagram in different time

Lf1

Xe-Lf1 Lf1 Lf2 Xe-Lf-Lf2 Lf1 Lf2 Lfn Xe-Lf1-

Lf2-Lfn

(1) T0 (2)2T0 (3) nT0

123

346 J Petrol Explor Prod Technol (2016) 6:343349

Forma on: matrix and fracture ow equa on

PKN fracture extension model

Fracture length and stress distribu on

Forma on pressure distribu on last me

Fracture ltra on coecient this me

Dynamic grid

Inner boundary dieren al equa on

Forma on Dieren al equa on and coecient matrix

Outer boundary dieren al equa on

Fig. 3 Solving owchart

gi;j

Ct/maKmPNmi;j lemi;jCt/m aT0Km

ci;j ai;j bi;j di;j ei;j

Ct/f

T0

On the basis of the value of the pressure distribution of last period, the value of the leak-off in Eq. (18) can be approximated as

C

t s

Ct/f PNfi;j

T0

Nf i;1PNF PNfi;1lef Dy1 : 19

2Kfie bPi P

Ct/maKm lemi;jCt/m aT0Km

lefi 1=2

3n 1

8n

p

!

1 n

2n

/f 8Kfie bPi P

Nf i;j

Thus, by the t period (which can also be expressed as NT0), we can get the ltration coefcient at the fracture tip as

C

2Kfie bPi P

2kn

1n PNfi 1 PNfi Dxi Dxi 1=2

n 1

n

Nf i;1PNF PNfi;1 lef Dy1

NT0

ci;j ai;j bi;j di;j ei;j

Ct/f

T0

Ct/maKm lemi;jCt/m aT0Km

p : 20

With the combination of the Eqs. (18) and (20), according to Carter method, at the t period (which can be expressed as the rst (N ? 1) one time T0), the full length and pressure distribution of the fracture are derived as

LN 1T0 Q

N 1

N

!

1 n

2n

lemi;j

3n 1

8n

/f 8Kfie bPi P

N mi;j

r

lef Dy1

2Kfie bPi P

N

f i;1PNF PNfi;1

1n PNmi1j PNmi 1;j Dxi 1 2Dxi Dxi1

21

2kn

14

t18

n 1

n

:

PF

4G

1 vH

21 vlQ2

p3GCH

22

Then the fracture is regarded as the boundary condition of formation pressure, which can be expressed as

PNmi;j1 PNmi;j 1

Dyj 1 2Dyj Dyj1

( )

14

rH:

xL sin 1

x L

14

1

x L

2

p

2

x L

Articial fracture

The fracture synthetical owing equation can be derived of the combination from Eqs. (9) to (11) as

G641 vHl

o2W4

ox2

Clef Dy1

2Kfie bPi P

oWot

8Cp

p 0: 18

PN1fi;1 PN1F

Nf i;1

: 23

t s

p

N 1T0

123

J Petrol Explor Prod Technol (2016) 6:343349 347

(a)

(b)

Based on all the equations above, considering the additional boundary condition and initial condition, we can use the strong implicit NewtonRaphson iterative method to derive the pressure value at each point. By using the Eq. (20), we can get fracturing uid ltration coefcient and rate at different locations. The overall calculation owchart is shown as follows: rst, we need to discrete the formation owing equation, then combine it with the formation pressure distribution at last period to obtain the correlation coefcient and the outer boundary condition of differential equations; then we make use of PKN model and combine it with the pressure distribution of grids nearby the fracture on the last moment to derive the dynamic ltration coefcient, and the length of the fracture and pressure distribution in fractures are obtained. Then the pressure distribution of fractures is regarded as a boundary condition of the formation; the differential equations, the outer boundary conditions, and the inner boundary conditions are integrated to obtain the current reservoir pressure by using the NewtonRaphson iterative method and the fracturing uid ltration coefcient is nally obtained. Solving owchart is shown in Fig. 3.

Calculation and analysis

Based on the models above, making use of non-uniform grid mesh which is divided into 20 perpendicular to the fracture, we may analyze the examples and factors. The other basic data are shown in Table 1.

As is shown in Fig. 4, considering stress sensitivityduring the fracture propagation process, we can get the lawof fracturing uid leak-off. It can be concluded from Fig. 4,as the stress sensitivity coefcients increase, the fracturing uid rate and ltration coefcient decrease; as time passes by, the leak-off rate decreases but the ltration coefcient increases. In the beginning of fracturing, both the fracturing uid rate and ltration coefcient appear obvious variation, and both of them remain unchanged later. As is shown in Fig. 4b, when the stress sensitivity coefcient increases to a certain value, the ltration coefcient will rstly decreases and then increases. Based on comprehensive analysis of the reason, it is obtained that, along with the continuous injection of the fracturing uid, the articial fractures pressure continues to increase, and due to the fracturing uid leak-off and the increase of the formation pressure, the stress sensitivity makes the formation fractures open and the permeability increase, which improve the capacity of the increase of the formation pressure. Also as the fracture permeability increases, the comprehensive effect results in the decrease of leak-off. When the pressure wave reaches the boundary, the sealed boundary will weaken the effect; therefore, the fracturing uid ltration coefcient will increase again.

Table 1 related basic data

Parameters Value

Well control area (m2) xe=100 m, ye=100 m

Initial formation pressure (Mpa) 25

Initial fracture permeability, K(mD) 10

Stress sensitivity coefcient, b(MPa-1) 0.5

Minimum horizontal stress (MPa) 20.5

Consistency coefcient, Kn (Pa sn) 0.01

Fracture porosity, Uf 0.01 Fluid index, n 0.85

Matrix porosity, Um 0.01 Matrix permeability, Km(mD) 1

Fracture compressibility coefcient, Ctf (MPa-1) 0.00025

Matrix compressibility coefcient Ctm (MPa-1) 0.00030

Pumping volume, Q(m/min) 103

Fig. 4 The effect of stress sensitivity on fracturing uid leak-off: a leak-off rate, b leak-off coefcient

123

348 J Petrol Explor Prod Technol (2016) 6:343349

(a)

(b)

(a)

(b)

Fig. 5 The effect of pumping rate on fracturing uid leak-off: a leak-off rate, b leak-off coefcient

Fig. 6 The effect of different propagating lengths on fracturing uid leak-off: a leak-off rate, b leak-off coefcient

Figure 5 shows the laws of fracturing uid leak-off considering the different fracturing uid pumping rate during the fracture propagation process. It can concluded from Fig. 5, as the fracturing uid pumping rate increases, the leak-off velocity and ltration coefcient will increase. As time passes by, the leak-off rate decreases but the ltration coefcient continues to increase. In the beginning of the fracturing, both of them appear obvious variation rst and remain unchanged later. As is shown in Fig. 4b, when the stress sensitivity coefcient increases to a certain value, the ltration coefcient will rstly decrease and then increase. Based on comprehensive analysis of the mechanism, the increased volume of uid pumped in increases the uid pressure of the fractures, and then increases the fracturing uid leak-off. However, when the volume increases to a certain value, the formation pressure will increase, the fracture permeability resulting from the stress sensitivity will increase, and the articial fracture pressure will also increase. This comprehensive effect makes the ltration coefcient decrease, and the sealed boundary

effect results in the increase of the ltration coefcient later. Thus, an optimal displacement volume of fracturing uid pumped in exists during the process of fracturing.

Figure 6 shows the laws of fracturing uid leak-off considering the different stretching lengths of the articial fractures. In Fig. 6, the lengths of the fractures are, respectively, 20, 40, and 60 m. It can be concluded that, for a certain length of the fracture, the leak-off rate and the ltration coefcient are different, respectively, at the different locations of the fractures. The closer the distance to the end of the fracture is, the smaller the leak-off rate is and the greater the ltration coefcient is. For the different lengths of the fractures, the ltration coefcients are different at the same position of the fracture, the longer the fractures are, the greater the leak-off rate and the ltration coefcient are. Based on the analysis of the mechanism, the longer the fracture and the time of pumping in fracturing uids are, the larger the articial fracture pressure will be, and the amount of leak-off will also increase. While the

123

J Petrol Explor Prod Technol (2016) 6:343349 349

closer to the end of the fracture, the more signicant the effect of the sealed boundary will be, and the worse the formation pressure supply will be. All of the phenomena above will result in larger fracturing uid leak-off rate and greater ltration coefcient.

Conclusion

1. Based on the classical PKN two-dimensional fracture-stretching mathematical model, the two-dimensional leak-off model of fracturing uid of fractured dual-medium reservoir is established by considering the time-varying non-Newtonian fracturing uid leak-off coefcient in the stretching process of fractures;

2. Using the nite element difference method, a dynamic discrete grid system is established and solved by NewtonRaphson iterative method, and the relevant factors are analyzed;

3. As the stress sensitivity coefcient increases, the fracturing uid ltration coefcient and fracturing uid leak-off rate will decrease, while as the volume of fracturing uid pumped in enlarges, the fracturing uid ltration coefcient and the leak-off rate will increase. When the stress sensitivity coefcient or the displacement volume of fracturing uid pumped in increases to a certain value, under the comprehensive effects of the formation pressure, the fracture pressure, the boundary, and the fracture permeability, the ltration coefcient will rstly decrease and then increase; and

4. Considering the propagation length of the fracture, the longer the length is, the more signicant the leak-off effect will be, and the boundary effect close to the end of the fracture will increase the fracture ltration coefcient and the leak-off rate.

Acknowledgments The authors acknowledge a fund called the engineering theory and risk control of Deepwater drilling and

completion from the National Basic Research Program of China (No. 2015CB251206) and supports from the MOE Key Laboratory of Petroleum Engineering.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

References

Al-Shatri A, Bellah S, Furuya K, Umezawa YT (2009) Successful application of 3D seismic technology on middle cretaceous carbonate reservoir-signicant OIP Increasing. SPE paper 12539

Balen RM, Meng H-Z, Economides MJ (1988) Applications of the net present value (NPV) in the optimization of hydraulic fracture.SPE paper 18541

Demarchos AS, Chomatas AS (2004) Pushing the limits of hydraulic fracture design.SPE paper 86483

Economides MJ, Demarchos AS (2008) Benets of a p-3D Over a 2D model for unied fracture design. SPE paper 112374

Fan Y, Economides MJ (1995) Fracturing uid leak off and net pressure behavior in Frac and Pack stimulation. SPE paper29988

Gidley JL, Holditch SA, Nierode DE, Veatch Jr RW (1989) Recent advances in hydraulic fracturing. SPE Monogr, vol 12

Mayerhofer MJ, Economides MJ, Nolte KG (1991) An experimental and fundamental interpretation of fracturing lter-cake uid loss. SPE paper 22873

Nghiem LX, Forsyth PA Jr, Behie A (1984) A fully implicit hydraulic fracture model. J Petroleum Technol 36(7):11911198

Ouenes ARC, Hartley LJ (2000) AEA Technology. Integrated fractured reservoir modeling using both discrete and continuum approaches. SPE paper 62939

Settari A (1985) A new general model of uid loss in hydraulic fracturing. SPE-JPT pp 491501

Settari A (1998) General model of uid ow (leak off) from fractures induced in injection operations. SPE paper 18197

Yi TC, Penden JM (1993) A comprehensive model of uid loss in hydraulic fracturing. Paper SPE 25493

123