- Email: [email protected]
Horizontal well transient rate decline analysis in low permeability gas reservoirs employing an orthogonal transformation method
Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
Horizontal well transient rate decline analysis in low permeability gas reservoirs employing an orthogonal transformation method Li-Na Cao a, *, Xiao-Ping Li a, Cheng Luo b, Lin Yuan c, Ji-Qiang Zhang d, Xiao-Hua Tan a, ** a
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, No. 8 Xindu Road, Xindu Street, Chengdu, Sichuan 610500, China b Geological Exploration and Development Research Institute of CNPC Chuanqing Drilling Engineering Ltd, Chengdu, Sichuan 610056, China c Northeastern Sichuan Production Gas Plant of Sinopec Southwest Oil and Gas Field Co., Langzhong, Sichuan 637402, China d CNOOC Zhanjiang Co., Ltd., Zhanjiang, Guangdong 524057, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 December 2015 Received in revised form 26 May 2016 Accepted 1 June 2016 Available online 6 June 2016
Transient rate decline analysis (TRDA) is an important analysis in modern reservoir engineering and applications. Most studies have focused on the threshold pressure gradient (TPG) of low permeability reservoirs and have placed an emphasis on the transient pressure response or steady productivity, whereas transient production decline has not been a focus. This paper presents a transient production decline analysis model of a horizontal well in low permeability gas reservoirs considering the threshold pressure gradient. An accurate solution is determined by employing the Laplace transformation, SturmLiouville eigenvalue method, and orthogonal transformation. The bi-logarithmic type curves of dimensionless production and the derivative are plotted by the Stehfest numerical inversion method. Five different flow regimes were recognized, and the effects of the influencing factors, such as the threshold pressure gradient, wellbore storage coefficient, skin factor, horizontal section length, and formation thickness, are discussed. This research contributes to the understanding of transient flow behavior and provides the theoretical basis for efficiently exploiting low permeability reservoirs. © 2016 Elsevier B.V. All rights reserved.
Keywords: Transient rate decline analysis Horizontal well Orthogonal transformation Low permeability
1. Introduction As a hot research topic in modern reservoir engineering and applications, transient rate decline analysis (TRDA) has become the focus of study in recent years. TRDA can quantitatively calculate the physical parameters of gas wells and reservoirs and can also be used to monitor reservoir dynamics (Aminian and Ameri, 1989). Horizontal wells are suitable for the development of low permeability reservoirs. TRDA of horizontal wells in low permeability gas reservoirs has both theoretical and practical significance to understand fluid transport behavior characteristics and correctly guide reservoir production. Numerous studies on the production rate transient analysis of horizontal wells have been performed. Ehlig-Economides and Ramey Jr. (1981) studied the transient rate response of a well that
* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (L.-N. Cao), [email protected] (X.-H. Tan). http://dx.doi.org/10.1016/j.jngse.2016.06.001 1875-5100/© 2016 Elsevier B.V. All rights reserved.
produced at a constant pressure according to three outer boundary conditions and several values of the wellbore skin factor by using the numerical Laplace transformation inversion algorithm. Fraim and Wattenbarger (1987) introduced the concept of normalized time for a gas reservoir producing against a constant wellbore pressure during (external) boundary-dominated flow. Therefore, type curve matching can be applicable to a reservoir with any shape. Aminian and Ameri (1989) discussed the application of the production type curves and issues concerning the production performance of horizontal wells in low permeability gas reservoirs. They obtained consistent results between the pre-stimulation data and type curve. Poon (1991) developed production decline curves for horizontal wells under various reservoir boundary conditions using Green’s and source functions. He noted that the wellbore stimulation, well spacing, and well length are probably the most important factors determining the effectiveness of a horizontal well. Doublet et al. (1994) presented a type of method to perform a decline curve analysis for both the variable rate and variable bottom hole pressure cases, without regard to the structure of the reservoir (shape and size) or the reservoir drive mechanisms. The
704
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
parameters determined using this method include the skin factor for near-well damage or stimulation, formation permeability, reservoir drainage area, and original oil-in-place. Accounting for the relative permeability and capillary pressure, Kewen and Horne (2003) proposed a decline analysis model based on fluid flow mechanisms. The model reveals a linear relationship between the oil production rate and reciprocal of the oil recovery or accumulated oil production. Medeiros et al. (2007) employed a semi-analytical model that incorporates the key features of reservoir heterogeneity as well as the details of hydraulic fracture and the wellbore flow to compute the production decline. They proposed a transient productivity index to represent the production decline characteristics. Ozkan et al. (2009) presented a discussion of the fractured horizontal-well performance in conventional (milli-Darcy permeability) and unconventional (micro-to nanoDarcy permeability) reservoirs and provided an interpretation of the different production performances of fracturing horizontal wells in both types of formations. Cai et al. (2010, 2012) (Cai and Yu, 2011)analyzed the capillary-driven flow in gas-saturated porous media based on fractal theory and modified the classical LucasWashburn equation. The factors influencing the imbibition process upon approaching the equilibrium weight were also analyzed. Later, the hydraulic conductivity and effective porosity in porous media were investigated. Nie et al. (2011) proposed the negative skin problem in a transient well test and rate decline analysis; they established a new mathematical model for a horizontal well in a multiple-zone composite reservoir. The contribution of their work was the removal of the oscillation in type curves, thus solving the negative skin problem. Based on the widely familiar Arp’s equation, Bahadori (2012) proposed a simple-to-use method, which was formulated to arrive at an appropriate estimation of the nominal (initial) decline rate and Arp’s decline-curve exponent. This method is quite simple to apply and accurately generates the coefficients of the equations instead of using already generated coefficients with uncertainty. Without applying the empirical concepts from Arps decline models and the explicit calculations from pseudofunctions, Ayala and Ye (2013) suggested a new generation analytical decline equation for boundary-dominated flow in a gas well producing at constant pressure for both full and partial production potential. Tan et al. (2013) created a dual fractal reservoir transient flow model by embedding a fracture system simulated by a tree-shaped fractal network into a matrix system simulated by fractal porous media. They plotted the bilogarithmic type curves of the dual fractal reservoirs and discussed the influence of different fractal factors on the pressure transient responses. Xie et al. (2014, 2015) focused on the two-phase flow caused by flow back after hydrofracturing shale gas reservoirs and established a two-phase pressure transient analysis model of a multi-stage fractured horizontal well. They studied not only the transient pressure but also the transient production rate decline features. Li et al. (2014) established a characteristic value method of the well test analysis for horizontal gas wells. Following the characteristic lines that are manifested by the pseudo-pressure derivative curve of each flow period, they developed formulas to determine the horizontal permeability, vertical permeability, skin factor, reservoir pressure, and pore volume of the gas reservoir. However, flow in low permeability reservoirs does not obey Darcy’s law due to the threshold pressure gradient (TPG), which makes the flow mechanism of horizontal wells in low permeability reservoirs more complicated. The concept of a threshold pressure gradient was suggested by Florin in 1951 (Florin, 1965). In the middle of the 20th century, numerous scholars concluded that the non-Darcy percolation that exists in low permeability reservoirs is
caused by the threshold pressure gradient (TPG), and they investigated the flow problems related with TPG (Miller and Low, 1963; Thomas et al., 1968; Das, 1997; Boukadi et al., 1998; Bennion et al., 2000; Dey, 2007). Guo et al. (2012) presented a horizontal gas well model accounting for the threshold pressure gradient in tight gas reservoirs. They note that the effect of the threshold pressure gradient on the pressure behavior for horizontal gas wells mainly occurs during the mid-late time periods. Song et al. (2014) studied the steady productivity equations of fractured wells in waterbearing tight gas reservoirs under low-velocity non-Darcy flow. The influence of the threshold pressure gradient on the volumetric flow rate of a gas well is significant and should not be ignored. Zhao et al. (2014) established a nonlinear steady flow mathematical model of a horizontal well for oil-water two-phase flow, accounting for the permeability stress sensitivity and threshold pressure gradient. Xu et al. (2015) conducted TPG and stress sensitivity experiments and investigated the existing condition of the TPG and change rule of permeability modulus. Then, they presented a transient pressure analysis mathematical model for a horizontal well in low-permeability offshore reservoirs. Liu and Yao (2015) constructed a moving boundary model of the radial flow in lowpermeability reservoirs, considering the influence of the threshold pressure gradient and quadratic pressure gradient term. The model was solved using numerical methods because it is nonlinear. Cai (2014) and Tan et al. (2015) investigated the low velocity non-Darcy flow, which starts the pressure gradient based on fractal theory. Based on the above discussion, most studies on the TPG of low permeability reservoirs were mainly focused on the transient pressure or steady productivity, and solutions of an unsteady model were determined using point source functions theory following Ozkan and Raghavan (1991) and Ozkan et al. (2010). This paper presents a transient production decline analysis model of horizontal wells in low permeability gas reservoirs accounting for the threshold pressure gradient. An accurate solution is determined by applying the Laplace transformation, SturmLiouville eigenvalue method, and orthogonal transformation. The bilogarithmic type curves of the dimensionless production and derivative are plotted using the Stehfest numerical inversion method. The characteristics of the transient rate decline behavior are thoroughly analyzed, and different flow regimes are observed in each type of curve. In addition, the effects of the relevant parameters are discussed, which improves the understanding of transient flow behavior and efficient exploitation of low permeability reservoirs.
2. Mathematical model A physical model of a low permeability gas reservoir, as shown in Fig. 1, is based on fluid flow in a horizontal well. To make the mathematical model more tractable and easier to understand, the following hypothesis and descriptions are outlined: (1) The reservoir is homogeneous with impermeable top and bottom boundaries, and the lateral boundary is infinite. The initial formation pressure is pi and is equal everywhere; (2) The reservoir is anisotropic, and the horizontal and vertical permeability are Kh and Kv, respectively; (3) The horizontal section of the well is parallel to the closed top and bottom boundaries, and its location is unlimited and represented by zw; (4) The single-phase gas flow obeys a non-Darcy law that accounts for the threshold pressure gradient;
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
705
Fig. 1. Sketch map of a horizontal gas well in a formation (revised from Xiao-ping Li et al. (2014)).
(5) A slightly compressible gas has a constant viscosity and compressibility coefficient, and the gas flow is under isothermal conditions; (6) The influence of gravity and the capillary force are ignored, but both the wellbore storage effect and the skin effect are taken into consideration.
obeyed only if the pressure gradient reaches a certain limit value (threshold pressure gradient), lB (Fig. 2). There is no longer a linear relationship between the flow velocity and pressure gradient because of the threshold pressure gradient. The flow velocity can be expressed by (Zhao et al., 2014):
.
v ¼
2.1. Basic differential equation of flow 2.1.1. Mass conservation equation Based on the mass conservation principle, the generalized form of the continuity equation for single-phase fluid flow can be written as:
vðrfÞ . ¼0 V$ r v þ vt
(1)
In cylindrical coordinates, the partial differential form of Eq. (1) is:
! ! ! v ¼ vr r þ vq q þ vz z
(2)
vðrvr Þ rvr 1 vðrvq Þ vðrvz Þ vðrfÞ þ ¼ þ þ vr r vq vz vt r
(3)
.
8 > < K ðVp lB Þ; > :
m
0;
jVpj > lB
(6)
jVpj lB
where lB is the threshold pressure gradient, MPa/m. In cylindrical coordinates, the components of Eq. (6) is:
8 Kh vp > > ¼ 3:6 l v > r B < m vr > Kv vp > > : vz ¼ 3:6 lB m vz
(7)
Gas flow only exists along the radial r direction and vertical z direction around the horizontal well in the reservoir. There is no rotational motion along the q axis, which is represented by vq ¼ 0. Thus, the flow velocity vector (2) and the continuity equation (3) can be simplified to:
! ! v ¼ vr r þ vz z
(4)
vðrvr Þ rvr vðrvz Þ vðrfÞ þ ¼ þ vr vz vt r
(5)
.
2.1.2. Non-Darcy equation The pores of a low permeability reservoir are very small. The thickness of the fluid adsorbed layer on the wall of pores is as large as the size of the pore, and the flow resistance of the adsorbed layer is larger still, so fluid flow is difficult when the pressure gradient is not large enough. In a low permeability reservoir, Darcy’s law is
Fig. 2. Non-Darcy flow relationship with the threshold pressure gradient in a low permeability gas reservoir.
706
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Table 1 Definitions of the dimensionless variables. Dimensionless pseudo pressure
hh mD ¼ 78:489K ðmi mÞ Tqsc
Dimensionless pseudo threshold pressure gradient
h hL lmBD ¼ 78:489K lmB Tqsc
Dimensionless time
3:6Kh t tD ¼ 4m 2 Ct rw CD ¼ 2p4CC hL2 t qffiffiffiffi ffi LD ¼ hL KKhv qffiffiffiffiffi hD ¼ rhw KKhv
Dimensionless wellbore storage coefficient Dimensionless horizontal section length Dimensionless gas reservoir thickness Dimensionless radial distance Dimensionless outer boundary distance
rw L
Dimensionless horizontal section position
zwD ¼ zhw
f ¼ f0 1 þ Cf ðp p0 Þ
Mp RTZ
2p
mZ
lB
(18)
Eq. (17) can then be rewritten as:
(8)
(9)
All of the dimensionless definitions are shown in Table 1. The dimensionless flow differential equation, converted from Eq. (19), is:
The density of natural gas can be expressed in another form according to the equation of state for real gas:
r¼
lmB ¼
1 v vm l Kv v2 m fmCt vm r mB þ ¼ r vr vr r Kh vz2 3:6Kh vt
The gas state equation:
r ¼ r0 1 þ Cr ðp p0 Þ
(17)
The pseudo threshold pressure gradient is introduced as follows:
2.1.3. State equation Porous rock and gas are both elastic and compressible, so both the rock and gas state equations should be considered. The rock state equation:
The value of lB is commonly on the order of magnitude of 103. p v As a quadratic gradient term, the product term of lB and vr mZ or p v vz mZ can be ignored for the sake of simplicity. Therefore, Eq (16)
1 v vm 1 2p K v2 m fmCt vm r lB þ v 2 ¼ r vr vr r mZ Kh vz 3:6Kh vt
reD ¼ rLe rwD ¼ zD ¼ hz
(16)
can be rearranged into:
rD ¼ Lr
Dimensionless vertical distance
Dimensionless wellbore radius
1 v vm 2p Kv v vm 2p fmCt vm r lB lB ¼ þ r vr vr mZ Kh vz vz mZ 3:6Kh vt
1 v vm 1 v2 m vm rD D þ lmBD þ L2D 2D ¼ ðhD LD Þ2 D rD vrD rD vrD vtD vzD
(19)
(20)
(10) 2.2. Definite conditions
2.1.4. Dimensionless differential equation of flow Combining Eq. (5) and Eq. (7)e(10) and simplifying leads to the real gas flow differential equation that accounts for the TPG:
Kh v p vp v p vp fmCt p vp r lB lB þ Kv ¼ vz mZ vz r vr mZ vr 3:6 mZ vt (11) where the total compressibility coefficientCt ¼ Cr þ Cf , MPa1. The pseudo pressure is introduced in the following mathematical model:
Zp mðpÞ ¼ p0
2p dp mZ
Initial condition:
mD ðrD ; zD ; 0Þ ¼ 0
(21)
Inner boundary condition:
2
ε
lim 4 lim
εD /0
rD /0
εD 2
D zwD Zþ2
rD
ε zwD 2D
vmD þ lmBD vrD
3 1 dzwD 5 ¼ ; jzD zwD j 2
(22) Infinite lateral boundary condition:
(12)
mD ðrD /∞; zD ; tD Þ ¼ 0
(23)
Closed top and bottom boundaries: where p0 is the reference pressure and can be variable. Therefore:
vm 2p vp ¼ vr mZ vr
(13)
vm 2p vp ¼ vz mZ vz
(14)
vm 2p vp ¼ vt mZ vt
(15)
by
Combining Eq. (13)e(15) into Eq. (11) and multiplying both sides 2 Kh, we obtain:
vmD
vmD
¼ 0; ¼0 vzD zD ¼0 vzD zD ¼1
(24)
2.3. Problem for definite solutions By simultaneously solving Eq. (20)e(24), the mathematical flow model of a horizontal well in a low permeability gas reservoir can be written as follows:
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
707
8 1 v vm 1 v2 m vm > rD D þ lmBD þ L2D 2D ¼ ðhD LD Þ2 D > > rD vrD vtD > rD vrD vz > D > > > > > mD ðrD ; zD ; 0Þ ¼ 0 > > > > εD > 2 3 > zwD > Zþ2 < vmD 1 ε 4 þ lmBD rD dzwD 5 ¼ ; jzD zwD j D lim lim > 2 vr 2 ε /0 r /0 D D > D > ε > zwD 2D > > > > > > > mD ðrD /∞; zD ; tD Þ ¼ 0 > > >
> > vmD
> : vmD
¼ 0; ¼0 vzD zD ¼0 vzD zD ¼1
(25)
3.2. Sturm-Liouville eigenvalue problem
3. Solutions of the mathematical model Importantly, Eq. (20) is so nonlinear that it is hard to calculate an exact solution by conventional methods. For the initial boundary value problem (25), the first step is to apply the Laplace transform to the time variable. The second step is to apply the orthogonal transform to the spatial variable, thus reducing the nonlinear degree of the problem to help calculate the final solution.
The differential equation of model (28) can be rewritten as:
v2 mD 1 vmD 1 RmD þ þ f ðrD ; zD Þ ¼ 2 rD vrD gðzD Þ vrD
where R is the differential operator, and pðzD Þ, qðzD Þ, and gðzD Þ are all real functions on [0,1]. The expressions are:
3.1. Laplace transformation The Laplace transformation on mD ðrD ; zD ; tD Þ from tD to s defines mD as the Laplace transformation function of mD , that is:
Zþ∞ mD ðrD ; zD ; sÞ ¼ L½mD ðrD ; zD ; tD Þ ¼
mD ðrD ; zD ; tD ÞestD dtD
0
(26) Applying the initial conditions, the derivative of mD with respect to tD in the Laplace domain can be written as:
dmD ðrD ; zD ; tD Þ ¼ s$L½mD ðrD ; zD ; tD Þ mD ðrD ; zD ; tD ÞjtD ¼0 L dtD
(29)
R¼
v v pðzD Þ qðzD Þ vzD vzD
(30)
pðzD Þ ¼ L2D
(31)
qðzD Þ ¼ 0
(32)
gðzD Þ ¼ 1
(33)
f ðrD ; zD Þ ¼
1 l ðhD LD Þ2 smD srD mBD
(34)
The integral transformation on the space variable zD is as follows:
¼ smD (27) The dimensionless mathematical model (25) in the Laplace domain can be written as:
MðrD ; mÞ ¼ F½mD ðrD ; zD Þ ¼
Z1
gðzD ÞZðzD ; mÞmD ðrD ; zD ÞdzD
(35)
0
where gðzD ÞZðzD ; mÞ is known as the integral transformation kernel.
8 1 v vm 1 v2 m > rD D þ lmBD þ L2D 2D ¼ ðhD LD Þ2 smD > > r vr sr vr > D D vzD D D > > > > εD > 2 3 > zwD > Zþ2 > > vmD lmBD 1 ε > < lim 4 lim dzwD 5 ¼ ; jzD zwD j D rD þ 2s vrD s 2 εD /0 rD /0 ε > zwD 2D > > > > > > > mD ðrD /∞; zD ; sÞ ¼ 0 > > >
> > > vmD
: vmD
¼ 0; ¼0
vzD zD ¼0 vzD zD ¼1
(28)
708
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
The Sturm-Liouville theorem solves an eigenvalue problem that is composed of a second-order variable coefficient differential equation and corresponding boundary conditions. To determine the transformation kernel, we need to solve a Sturm-Liouville eigenvalue problem relating to ZðzD ; mÞ, which is:
lmBD
F
Z1 ¼
s
lmBD s
0
8 < lmBD ; n ¼ 0 s cos npzD dzD ¼ : 0 ; ns0
(41)
"
8 lgZ; 0 < zD < 1 > < RZ ¼
dZ
dZ
> ¼ ¼0 :
dzD zD ¼0 dzD zD ¼1
(36)
# Z1 v2 mD v2 mD F cos npzD dzD ¼ n2 p2 mD ¼ vz2D vz2D
(42)
1 v vmD 4 rD þ n n2 p2 L2D mD ¼ ðhD LD Þ2 smD rD vrD vrD rD
(43)
8 < lmBD ; n ¼ 0 s , the differential equation of Defining 4n ¼ : 0 ; ns0 model (28),can be orthogonally transformed to: 0
Merging the constant coefficient, Eq. (36) can be simplified to:
8 00 > < Z þ lZ ¼ 0; 0 < zD < 1
dZ
> dZ
¼ ¼0 : dzD zD ¼0 dzD zD ¼1
(37)
where l is the corresponding eigenvalue when Eq. (37) has a nonzero solution. Characteristic equation: r 2 þ l ¼ 0 ①When l ¼ 0, r1 ¼ r2 ¼ 0, then ZðzD Þ ¼ c1 þ c2 zD . Z’ð0Þ ¼ 00c2 ¼ 0, naturally, meet the requirement of Z’ð1Þ ¼ 0. So ZðzD Þ ¼ c1 ðc1 s0Þ. pffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffi lzD . ②When l < 0, r1;2 ¼ ± l, then ZðzD Þp¼ffiffiffiffiffic1 e lzDpþ ffiffiffiffiffic2 e Z’ð0Þ ¼ 00c1 c2 ¼ 0; Z’ð1Þ ¼ 00c1 e l c2 e l ¼ 0. Solve simultaneously and obtainc1 ¼ c2 ¼ 0, such that ZðzD Þ has a non-zero solution. pffiffiffi ③When pffiffiffi l > 0,pffiffiffi r1;2 ¼ ± li, then ZðzD Þ ¼ c1 cos lzD þ c2 sin lzD . Z’ð0Þ ¼ 00c2p¼ffiffiffi 0; pffiffiffi Z’ð1Þ ¼ 00c1 sin l ¼ 00 l ¼ np; n ¼ 1; 2; 3:::. So ln ¼ ðnpÞ2 ðn ¼ 1; 2; :::Þ, Zn ðzD Þ ¼ cn cos npzD . According to the above discussion on the S-L problem (37), the corresponding eigenvalues are ln ¼ ðnpÞ2 , and the corresponding eigenfunctions are Zn ðzD Þ ¼ cos npzD (where n ¼ 0; 1; 2:::), which are composed of a weighted function gðzD Þ into a complete orthogonal system in the [0,1] domain. Thus, the integral transformation kernel is:
gðzD ÞZðzD ; mÞ ¼ 1$Zn ðzD Þ ¼ cos npzD ; n ¼ 0; 1; 2:::; m ¼ n
(38)
The integral transformation whose kernel is a complete orthogonal function system can collectively be known as the orthogonal transformation.
The next step is to take the orthogonal transformation based on the inner boundary conditions. The right side of the inner boundary 1 , which is independent of ε . Therecondition of model (28) is 2s D fore, the following expression is established:
1 1 ¼ lim 2s εD /0 2s
(44)
Then, the inner boundary condition can be written as:
2
ε
D zwD Zþ2
lim 4 lim
εD /0
rD
rD /0
ε
¼ lim
εD /0
zwD 2D
vmD lmBD þ vrD s
3 dzwD 5
1 ε ; jzD zwD j D 2s 2
(45)
Simplifying Eq. (45), we obtain:
lim
lim εD rD
εD /0 rD /0
vmD 1 ε ; jzD zwD j D ¼ lim 2s vrD 2 εD /0
Because the integral calculation of the variable ZD is independent of variable εD, the orders of the orthogonal transformation and the limit calculation of εD can be exchanged, leading to:
ℱ
lim
εD /0
lim εD rD
rD /0
vmD vrD
¼ℱ
lim
εD /0
1 2s
i
vm 1 lim ℱ lim εD rD D ¼ lim ℱ 2s vrD εD /0 rD /0 εD /0
3.3. Orthogonal transformation
(46)
(47)
Because: The orthogonal transformation of variable zD , with eigenfunctions cos npzD as the kernel and mD as the orthogonal transformation function of mD , can be written as:
Z1 mD ðrD ; n; sÞ ¼ F½mD ðrD ; zD ; sÞ ¼
Based on the differentiability properties, we obtain:
Z1 vmD v vmD ¼ mD ðrD ; zD ; sÞcos npzD dzD ¼ vrD vrD vrD 0
rD /0
1 ℱ 2s (39)
lim εD rD
mD ðrD ; zD ; sÞcos npzD dzD
0
F
ℱ
(40)
vmD vrD
¼ lim εD rD rD /0
vmD vrD
(48)
ε
D zwD Zþ2
¼
ε
cos npzD dzD 2s
zwD 2D
ε ε sin np zwD D sin np zwD þ D 2 2 ¼ 2snp 3 2 cos npzwD sin np D 2 ¼ 2snp Thus:
(49)
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
lim εD rD
rD /0
cos npzwD sin np ε2D vmD ¼ vrD snp
(50)
Because variable εD is independent of rD, moving εD to the right side of Eq. (50), we obtain:
lim rD
rD /0
cos npzwD sin np ε2D vmD ¼ vrD snpεD
εD /0
lim rD
rD /0
vmD vrD
¼ lim
εD /0
cos npzwD sin np ε2D snpεD
(51)
rD /0
cos npzwD sin np ε2D vmD ¼ lim vrD snpεD εD /0
(52)
lim
εD /0
cos npzwD sin np snpεD
εD 2
ε D d cos npz wD sin np dðsnpεD Þ 2 ¼ lim dεD dεD εD /0 εD np cos npzwD cos np 2 ¼ cos npzwD ¼ lim 2s 2snp εD /0
(54)
lim rD
vmD cos npzwD ¼ vrD 2s
v2 mD 2 vrD
¼
v vrD
vmD vx $ vx vrD
(55)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v vmD $ n2 p2 L2D þ h2D L2D s vrD vx qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 v mD vx $ ¼ n2 p2 L2D þ h2D L2D s$ vx2 vrD v2 m D ¼ n2 p2 L2D þ h2D L2D s vx2
¼
(61) Combining Eq. (60) and Eq. (61) into Eq. (59), we have:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2 m vmD D 2 n2 p2 L2D þ h2D L2D s rD þ r n2 p2 L2D þ h2D L2D s D vx vx2 2 2 2 2 2 2 rD n p LD þ hD LD s mD ¼0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mD ¼ AI0 rD n2 p2 L2D þ h2D L2D s þ BK0 rD n2 p2 L2D þ h2D L2D s (62)
(56)
Simultaneously solving Eqs. (43), (55), and (56) for a mathematical model on mD can be written as:
8 > > > > > > > > <
4 vmD 1 v vmD þ n n2 p2 L2D mD ¼ ðhD LD Þ2 smD lim rD rD > rD vrD rD /0 r vr vr > D D D > > > > > > :
where I0 is the first zero order modified Bessel function, K0 is the second zero order modified Bessel function, and A and B are constants. If n ¼ 0, 4n ¼ lmBD s s0, the basic differential equation of model (57) is inhomogeneous because its free term is not equal to zero. The particular solution can be expressed by the Green functions. The general solution can be written as:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mD ¼ AI0 rD n2 p2 L2D þ h2D L2D s þ BK0 rD n2 p2 L2D þ h2D L2D s Zþ∞ þ
GðrD ; tÞdt
0
cos npzwD ¼ mD ðrD /∞; n; sÞ ¼ 0 2s
(63) (57)
If ns0 and 4n ¼ 0, the basic differential equation of model (57) can be converted into a zero order modified Bessel equation. The derivation process is shown in the following steps:
v2 mD 1 vmD 2 2 2 þ n p LD þ h2D L2D s mD ¼ 0 2 rD vrD vrD
(60)
2
Taking the orthogonal transformation on the outer boundary conditions of model (28):
mD ðrD /∞; n; sÞ ¼ 0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vmD vmD vx vmD $ $ n2 p2 L2D þ h2D L2D s ¼ ¼ vrD vx vrD vx
2 2 D Namely: x2 vvxm2D þ x vm vx ðx þ y ÞmD ¼ 0, where y ¼ 0 and y represents the order number of the modified Bessel equation. The general solution of the equation is written as:
Therefore, the inner boundary condition is simplified to:
rD /0
n2 p2 LD þhD LD s
(53)
Applying the L’Hospital Rule to the right side of Eq. (53) is “00” uncertain limit type, that is:
(59)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x ffi. Following the Let x ¼ rD n2 p2 L2D þ h2D L2D s, so rD ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2
In addition, because variable εD is independent of rD, Eq. (52) can be rewritten as:
lim rD
v2 mD vmD 2 2 2 2 2 2 þ r n p L þ h L s rD mD ¼ 0 D D D D 2 vrD vrD
compound function derivation rule, we obtain:
Taking the limit calculation of εD on both sides of Eq. (51):
lim
2 rD
709
(58)
2 , and Eq. (58) can be Both sides of Eq. (58) are multiplied byrD rewritten as:
where the Green function
R∞ 0
GðrD ; tÞ is:
qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi 8 lmBD 2 2 2 2 > > K ð0 < t < rD Þ > 0 rD hD LD s I0 t hD LD s < s GðrD ; tÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi > lmBD qffiffiffiffiffiffiffiffiffiffiffiffiffi > > : K0 t h2D L2D s I0 rD h2D L2D s ðrD < t < þ ∞Þ s (64) Combining Eq. (56) with an infinite outer boundary condition and I0 ð∞Þ ¼ ∞, we obtain A ¼ 0. So, Eq. (63) is rewritten as:
710
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
mD ¼ BK0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Zþ∞ rD n2 p2 L2D þ h2D L2D s þ GðrD ; tÞdt
(65)
0
Combining the inner boundary conditions of Eq. (55) and the properties of the Bessel functions, if z/0, then I1 ðzÞy0, so npzwD B ¼ cos 2s . Replacing B into Eq. (65) and taking the inverse orthogonal transformation, we determine the pressure solution in the Laplace domain:
qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 mD ¼ K0 rD h2D L2D s 2s qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∞ 1X K0 rD n2 p2 L2D þ h2D L2D s cos npzwD cos npzD þ s n¼1 Zþ∞ þ
GðrD ; tÞdt
0
(66) Following the source function idea method, let yD ¼ 0, and use Eq. (66) to mark the point-source solution. Then, compute the integral for the horizontal section of the wellbore to obtain the bottom hole pseudo pressure solution expression, as follows, that accounts for the effect of the threshold pressure gradient:
mD ¼
1 2s
Z1
qD ðsÞ ¼
1 s2 mwD ðsÞ
(69)
Following Eq. (67)e(69), we can determine the horizontal well dimensionless rate solution under the influence of the threshold pressure gradient for low permeability gas reservoirs. 4. Type curves and discussion 4.1. Rate decline performance type curve analysis Type curves of the production rate can indicate the characteristics of the flow behavior to some degree and can be used to apply a curve-fitting analysis to determine the property parameters of the well and the formation. The solutions of the dimensionless production rate in real space can be determined by applying the Stehfest numerical inversion algorithm (Stehfest, 1970). The whole production decline behavior of a horizontal well with low permeability is plotted in Fig. 3. The following five flow regimes can be recognized: Regime I: The pure wellbore storage effect flow regime. This stage occurs early, and the curves of the production rate and the derivative have the same trend, which is described by a downward sloping line, similar to a conventional reservoir. The gas stored in the wellbore controls this period.
Zþ∞ Z1 qffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj h2D L2D s da þ GðjxD aj; tÞdadt
1
0
1
1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∞ Z 1X þ K0 jxD aj n2 p2 L2D þ h2D L2D s cos npzwD cos npzD da s n¼1
(67)
1
Because the wellbore has a certain volume and natural gas has compressibility, the wellbore storage effect occurs when the well is open or shut and is represented by the wellbore coefficient. In addition, the skin effect around the bottom hole describes the situation of formation damage or stimulation treatment and is represented by the skin factor. When taking these two effects into account, the definite conditions will become nonhomogeneous. Duhamel’s theorem (or homogenized principle) can be used to divide the nonhomogeneous problem into homogeneous subproblems that are easily solved. The solutions to the subproblems are then combined to yield a solution to the original problem, which is named the superposition principle. Applying Duhamel’s theorem and the superposition principle to incorporate the wellbore storage coefficient and the skin factor into the well response, the pseudo pressure expression can be rewritten as:
mwD ¼
smD þ S s þ CD s2 ðsmD þ SÞ
(68)
The dimensionless wellbore flow rate equation for the constantpressure production case in low permeability gas reservoirs can be determined by the relationship between the dimensionless pressure and the rate in the Laplace space (Van Everdingen and Hurst, 1949):
Regime II: The skin effect transition flow regime. The drop in the degree of the production rate curve slowly changes, and the rate derivative curve continues to decline with a concave shape, which corresponds to the hump on the transient pressure curve. This stage reflects the combined action of the skin and wellbore storage. Regime III: The early radial flow regime in the vertical plane. The production rate curve continues to decrease, and the derivative curve is a level straight line, with a value related to LD. This period reflects the radial flow feature perpendicular to the horizontal axis before the pressure wave propagates to the top and bottom outer boundary (Fig. 4). Regime IV: The mid-term linear flow regime. The curves of the production rate and derivative are two parallel lines with the same slope. During this flow stage, the pressure wave has already propagated to the top and bottom boundary, and the natural gas flows linearly perpendicular to the horizontal section (Fig. 5). The influence of the threshold pressure gradient also appears in this period. Regime V: The late pseudo-radial flow regime affected by the TPG. The curves of the production rate and the derivative describe the pseudo-radial flow characteristics of the whole system. The effect of the threshold pressure gradient mainly exists in the fifth flow regime. It leads to a downwarping of the production curve and an upwarping of the rate derivative curve. This stage represents the radial flow characteristic of the entire system (Fig. 6).
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
711
Fig. 3. Type curves of the production rate decline for a horizontal well in a low permeability gas reservoir (solid lines represent production rate curves, dash lines represent the production rate derivative curves, red lines represent lmBD ¼ 0, blue lines represent lmBD ¼ 0.5). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. The early radial flow regime in the vertical plane of the horizontal well.
Fig. 5. The mid-term linear flow regime of the horizontal well.
4.2. Parameter influence on rate decline performance Fig. 7 shows the effect of the threshold pressure gradient on the production rate decline type curves. The figure suggests that the TPG mainly plays its role during the mid-term linear regime to the
Fig. 6. The late pseudo-radial flow regime of the horizontal well.
late pseudo-radial flow regime. The greater the TPG, the worse the reservoir properties, the more difficult it is for the fluid to flow, and the lower the production rate will be, resulting in a steeper downwarping feature on the production rate curve. Correspondingly, the faster the production curve decreases, the more obvious the upwarping of the rate derivative curve. Fig. 8 shows the effect of the wellbore storage coefficient on the production rate decline type curves. The wellbore storage coefficient has a major influence on the first two flow regimes: the pure wellbore storage effect regime and the skin effect transition regime. The greater the wellbore storage coefficient, the longer the wellbore storage period lasts and the less obvious the concavity of the derivative curves. When CD exceeds a certain value, the vertical radial flow regime might be concealed. Fig. 9 illustrates the effect of the skin factor on the production rate decline type curves. The skin factor reflects the amount of reservoir pollution near the wellbore area, and it mainly affects the skin transition flow regime. The incremental value of the skin factor results in increasing the additional filtration resistance. The bigger the skin factor, the longer the skin effect transition period lasts and the more obvious the concavity of the production rate derivative curves. After the skin effect transition regime, the impact of the skin effect is gradually removed.
712
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Fig. 7. Type curves of the production rate decline affected by the threshold pressure gradient (red line represents lmBD ¼ 0, blue line represents lmBD ¼ 0.1, green line represents lmBD ¼ 0.5).
Fig. 8. Type curves of the production rate decline affected by the wellbore storage coefficient (red line represents CD ¼ 104, blue line represents CD ¼ 103, green line represents CD ¼ 102).
Fig. 10 shows the effect of the horizontal section length of the well on the production rate decline type curves. Increasing the horizontal section length results in a shift upward of both the production rate curve and the rate derivative curve. This shift occurs because the longer the horizontal section, the larger the gas drainage area. The resistance of the gas flow becomes relatively small and the corresponding production rate will be faster. Moreover, the greater the value of LD, the longer the linear flow period lasts, and the next pseudo-radial flow regime will be postponed. Fig. 11 shows the effect of reservoir thickness on the production rate decline type curves. A change in the reservoir thickness has a significant influence on the type curves. When the reservoir thickness increases, the pressure wave needs more time to propagate to the top and bottom boundaries. Consequently, the skin transition flow period will have a more obvious V-shaped concavity, and the early vertical radial flow regime will longer. Moreover, the mid-term linear flow stage and the late pseudo-radial flow stage are delayed.
5. Model application with field data The JB4-10H well is a horizontal development well in the JingBian gas field in China. The total well depth is 3309 m, and the well is drilled in the Majiagou Group. After drilling is finished, the well is completed by perforation and acidification. The reservoir thickness is 8.6 m, and the formation porosity is 6.72%. The gas gravity is 0.6148, the gas viscosity is 0.02146 mPa s, and the gas deviation factor Z is 0.9821. The well radius is 0.0797 m. According to the deliverability test that occurred from 7 to 11 December 2013, the calculated absolute open flow was 41.93 104 m3/d. Before production, the test data suggested that the initial formation pressure was 30.34 MPa, the reservoir temperature was 103.9 C, and the tubing pressure and casing pressure were 24.20 MPa and 24.10 MPa, respectively. The commissioning date of the JB4-10H well was 15 December 2013. In the early stages of production, the proration flow was 12 104 m3/d, which allowed the pressure to drop quickly and the tubing pressure and casing pressure were
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
713
Fig. 9. Type curves of the production rate decline affected by the skin factor (blue line represents S ¼ 0.01, red line represents S ¼ 0.1, green line represents S ¼ 0.2).
Fig. 10. Type curves of the production rate decline affected by the horizontal section length (red line represents LD ¼ 5, blue line represents LD ¼ 10, green line represents LD ¼ 20).
reduced to 18.70 MPa and 19.60 MPa, respectively. Then, the JB410H well was shut in to allow the formation pressure to build up. Starting on 8 February 2014, the JB4-10H well produced under constant bottom hole pressure for 458 days. The production rate log-log plot of the JB4-10H well is shown in Fig. 12. Substituting Eqs. (67) and (68) into Eq. (69) and applying the Stehfest numerical inversion algorithm, the expression of the production type curve in the real domain is written as:
(" Z qD ¼
N ln 2 X V tD i¼1 i
1 2
1 þ CD s (" Z s 12
1
1
1 1
Z qffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj h2D L2D s da þ s
Z qffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj h2D L2D s da þ s
0
þ∞ Z 1 0
þ∞ Z 1 1
1
where N is an even number and s ¼ i$ln 2=tD . The Green function GðrD ; tÞ is shown as Eq. (64). The weight coefficient Vi is given by:
Vi ¼ ð1Þ 2 þi N
k¼
GðjxD aj; tÞdadt þ
GðjxD aj; tÞdadt þ
minði;N=2Þ X
∞ Z X
1
n¼1 1
(71)
# ) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj n2 p2 L2D þ h2D L2D s cos npzwD cos npzD da þ S
n¼1 1 1
∞ Z X
iþ1 2
k 2 þ1 ð2kÞ! N k !ðk!Þ2 ði kÞ!ð2k iÞ! 2 N
K0 jxD aj
# ) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 p2 L2D þ h2D L2D s cos npzwD cos npzD da þ S (70)
714
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Fig. 11. Type curves of the production rate decline affected by the reservoir thickness (red line represents hD ¼ 100, blue line represents hD ¼ 200, green line represents hD ¼ 400).
Fig. 12. Log-log plots of the production rate of the JB4-10H well.
Table 2 Rate transient analysis results of the JB4-10H well. Parameter
Value 3
Wellbore storage coefficient C (m /MPa) Skin factor S Horizontal permeability Kh (mm2) Vertical permeability Kv (mm2) Effective horizontal section length L (m) Threshold pressure gradient lB (MPa/m)
2.028 1.625 8.5 104 3.6 105 334.31 7.82 103
The parameters of the theoretic type curves are CD ¼ 104, S ¼ 0.1, LD ¼ 5, hD ¼ 100, and lmBD ¼ 0.5. By applying the type curves, the rate transient analysis results of the JB4-10H well are shown in Table 2. 6. Summary and conclusions In this paper, a mathematical model for a horizontal well in a low permeability gas reservoir that considers the threshold
pressure gradient is established. The production rate decline behavior under the effects of relative influence factors is analyzed. The following conclusions can be summarized: (1) The TPG mainly affects the mid-term linear flow regime and late pseudo-radial flow regime of horizontal well production decline type curves. A larger TPG can obviously cause a downward shift of the production rate curve and an upward shift of the rate derivative curve. (2) A large wellbore storage coefficient results in a longer pure wellbore storage period, even though the vertical radial flow regime might be hidden. Simultaneously, the concavity of the production rate derivative curve becomes more obvious when the skin factor is greater, which means the formation near the wellbore area is more polluted. (3) The increasing gas drainage area is caused by a larger horizontal section length, which can increase both the production rate curve and rate derivative curve. Type curves suggest that the linear flow period lasts longer and the pseudo-radial flow appears later.
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
(4) The reservoir thickness has little influence on the entire curves except in the pure wellbore storage regime. The greater the reservoir thickness, the deeper the V-shaped concavity of the derivative curve and the longer the early vertical radial flow. Acknowledgment The authors would like to acknowledge the support of the National Basic Research Program (973 Program) of China under Grant No. 2014CB239205. Nomenclature
Latin symbol C wellbore storage coefficient, m3/MPa CD dimensionless wellbore storage coefficient Ct total compressibility coefficient, MPa1 Cr fluid compressibility coefficient, MPa1 C4 rock compressibility coefficient, MPa1 h reservoir thickness, m hD dimensionless reservoir thickness Kh horizontal permeability, mm2 Kv vertical permeability, mm2 L horizontal section length, m LD dimensionless horizontal section length M molecular mass, kg/mol m pseudo pressure, MPa2/(mPa$s) mD dimensionless pseudo pressure mi initial pseudo pressure, MPa2/(mPa s) p pressure, MPa p0 reference pressure, MPa pi initial formation pressure, MPa q gas production rate, 104 m3/d qD dimensionless gas production rate qsc surface gas production rate, 104 m3/d R gas constant (0.008314 MPa m3/(kmol K) r radial distance, m rD dimensionless radial distance re outer boundary distance, m reD dimensionless outer boundary distance rw wellbore radius, m rwD dimensionless wellbore radius s Laplace transform variable S skin factor, dimensionless T absolute temperature, K t time, hours tD dimensionless time v velocity of gas flow, m/h Z gas deviation factor z vertical distance, m zD dimensionless vertical distance zw horizontal section position, m zwD dimensionless horizontal section position Greek symbols ε tiny variable l eigenvalue, dimensionless lB threshold pressure gradient, MPa/m lmB pseudo threshold pressure gradient, MPa2/(mPa s m) lmBD dimensionless pseudo threshold pressure gradient m gas viscosity, mPa s r gas density, kg/m3 n order number of modified Bessel equation
f
715
porosity of reservoir, fraction
Superscripts Laplace transform domain ¼ Orthogonal transform domain Subscripts D dimensionless h horizontal i initial m pseudo sc standard condition t total v vertical w wellbore Abbreviations TRDA transient rate decline analysis TPG threshold pressure gradient References Aminian, K., Ameri, S., 1989. Predicting horizontal well production performance using type curves. In: Paper SPE 19342 Presented at the SPE Eastern Regional Meeting, Morgantown, West Virginia, 24e27 October. Ayala, L., Ye, P., 2013. Unified decline type-curve analysis for natural gas wells in boundary-dominated flow. SPE J. 18 (1), 97e113. Bahadori, A., 2012. Analysing gas well production data using a simplified decline curve analysis method. Chem. Eng. Res. Des. 90 (4), 541e547. Bennion, D.B., Thomas, F.B., Ma, T., 2000. Recent advances in laboratory test protocols evaluate optimum drilling, completion and stimulation practices for low permeability gas reservoirs. In: Paper SPE 60324 Presented at SPE Rocky Mountain Regional/Low-permeability Reservoirs Symposium and Exhibition, Denver, Colorado, 12e15 March. Boukadi, F., Bemani, A., Rumhy, M., Kalbani, M., 1998. Threshold pressure as a measure of degree of rock wettability and diagenesis in consolidated Omani limestone cores. Mar. Petroleum Geol. 15 (1), 33e39. Cai, J., 2014. A fractal approach to low velocity non-Darcy flow in a low permeability porous medium. Chin. Phys. B 23 (4), 044701. Cai, J., Hu, X., Standnes, D.C., You, L., 2012. An analytical model for spontaneous imbibition in fractal porous media including gravity. Colloids Surfaces A 414, 228e233. Cai, J., Yu, B., 2011. A discussion of the effect of tortuosity on the capillary imbibition in porous media. Transp. Porous Media 89 (2), 251e263. Cai, J., Yu, B., Zou, M., Luo, Cai, Yu, L., 2010. Fractal characterization of spontaneous co-current imbibition in porous media. Energy Fuels 24 (3), 1860e1867. Das, A.K., 1997. Generalized Darcy’s law including source effect. J. Can. Petroleum Technol. 36 (6), 57e59. Dey, G.R., 2007. Nonlinear flow of gas at low velocity in porous media. J. Nat. Gas Chem. 16 (3), 217e226. Doublet, L.E., Pande, P.K., McCollum, T.J., Blasingame, T.A., 1994. Decline curve analysis using type curvesdanalysis of oil well production data using material balance time: application to field cases. In: Paper SPE 28688 Presented at the International Petroleum Conference and Exhibition of Mexico, Veracruz, Mexico, 10e13 October. Ehlig-Economides, C.A., Ramey Jr., H.J., 1981. Transient rate decline analysis for wells produced at constant pressure. Soc. Petroleum Eng. J. 21 (1), 98e104. Florin, B.A., 1965. Principles of Soil Mechanics (Tongji University). China’s industrial press, Beijing, pp. 38e39. Original work published 1959. Fraim, M.L., Wattenbarger, R.A., 1987. Gas reservoir decline-curve analysis using type curves with real gas pseudopressure and normalized time. SPE Form. Eval. 2 (4), 671e682. Guo, J., Zhang, S., Zhang, L., Qing, H., Liu, Q., 2012. Well testing analysis for horizontal well with consideration of threshold pressure gradient in tight gas reservoirs. J. Hydrodynamics 24 (4), 561e568. Kewen, L., Horne, R.N., 2003. A decline curve analysis model based on fluid flow mechanisms. In: Paper SPE 83470 Presented at the SPE Western Regional/AAPG Pacific Section Joint Meeting, Long Beach, California, 19e24 May. Li, X.-P., Yan, N.-P., Tan, X.-H., 2014. Characteristic value method of well test analysis for horizontal gas well. Math. Problems Eng. 2014. Article ID 472728. 10 pages. Liu, W., Yao, J., 2015. Numerical investigations of the effect of nonlinear quadratic pressure gradient term on a moving boundary problem of radial flow in lowpermeable reservoirs with threshold pressure gradient. Math. Problems Eng. 2015. Article ID 275057, 12 pages. Medeiros, F., Kurtoglu, B., Ozkan, E., Kazemi, H., 2007. Analysis of production data from hydraulically fractured horizontal wells in shale reservoirs. SPE Reserv. Eval. Eng. 13 (3), 559e568. Miller, R.J., Low, P.F., 1963. Threshold gradient for water flow in clay system. Soil Sci.
716
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Soc. Am. J. 27 (3), 605e609. Nie, R., Guo, J., Jia, Y., Zhu, Sh, Rao, Zh, Zhang, Ch, 2011. New modeling of transient well test and rate decline analysis for a horizontal well in a multiple-zone reservoir. J. Geophys. Eng. 8 (3), 464e476. Ozkan, E., Brown, M.L., Raghavan, R.S., 2009. Comparison of fractured-horizontalwell performance in conventional and unconventional reservoirs. In: Paper SPE 121290 Presented at the SPE Western Regional Meeting, San Jose, California ,USA, 24e26 March. Ozkan, E., Raghavan, R., 1991. New solutions for well-test-analysis problems: part 1analytical considerations. SPE Form. Eval. 6 (3), 359e368. Ozkan, E., Raghavan, R.S., Apaydin, O.G., 2010. M odeling of fl uid transfer from shale matrix to fracture network. In: SPE 134830, Presented at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19e22 September. Poon, D.C., 1991. Decline curves for predicting production performance from horizontal wells. J. Can. Petroleum Technol. 30 (1), 77e81. Song, H., Liu, Q., Yang, D., Yu, M., Lou, Y., Zhu, W., 2014. 2014. Productivity equation of fractured horizontal well in a water-bearing tight gas reservoir with lowvelocity non-Darcy flow. J. Nat. Gas Sci. Eng. 18, 467e473. Stehfest, H., 1970. Algorithm 368: numerical inversion of Laplace transforms. Commun. ACM 13 (1), 47e49. Tan, X.-H., Li, X.-P., Liu, J.-Y., Tang, C., Li, J., 2013. 2013. Pressure transient analysis of
dual fractal reservoir. J. Appl. Math. 2013. Article ID 137518, 9 pages. Tan, X.-H., Li, X.-P., Zhang, L.-H., Liu, J.-Y., Cai, J., 2015. Analysis of transient flow and starting pressure gradient of power-law fluid in fractal porous media. Int. J. Mod. Phys. C 26 (4), 1550045. Thomas, L.K., Katz, D.L., Tek, M.R., 1968. Threshold pressure phenomena in porous media. Soc. Petroleum Eng. J. 8 (02), 174e184. Van Everdingen, A.F., Hurst, W., 1949. The application of the Laplace transformation to flow problems in reservoirs. Pet. Trans. AIME 1 (12), 305e324. Xie, W., Li, X., Zhang, L., Wang, J., Cao, L., Yuan, L., 2014. Two-phase pressure transient analysis for multi-stage fractured horizontal well in shale gas reservoirs. J. Nat. Gas Sci. Eng. 21, 691e699. Xie, W.-Y., Li, X.-P., Zhang, L.-H., Tan, X.-H., Wang, J.-C., Wang, H.-T., 2015. Production decline analysis for two-phase flow in multifractured horizontal well in shale gas reservoirs. J. Chem. 2015. Article ID 212103, 10 pages. Xu, J., Jiang, R., Teng, W., 2015. Nonlinear flow characteristics and horizontal well pressure transient analysis for low-permeability offshore reservoirs. Math. Problems Eng. 2015. Article ID 387149, 13 pages. Zhao, Y., Zhang, L., He, Z., Zhang, B., 2014. Productivity for horizontal wells in lowpermeability reservoir with oil/water two-Phase flow. Math. Problems Eng. 2014. Article ID 364678, 9 pages.
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
Horizontal well transient rate decline analysis in low permeability gas reservoirs employing an orthogonal transformation method Li-Na Cao a, *, Xiao-Ping Li a, Cheng Luo b, Lin Yuan c, Ji-Qiang Zhang d, Xiao-Hua Tan a, ** a
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, No. 8 Xindu Road, Xindu Street, Chengdu, Sichuan 610500, China b Geological Exploration and Development Research Institute of CNPC Chuanqing Drilling Engineering Ltd, Chengdu, Sichuan 610056, China c Northeastern Sichuan Production Gas Plant of Sinopec Southwest Oil and Gas Field Co., Langzhong, Sichuan 637402, China d CNOOC Zhanjiang Co., Ltd., Zhanjiang, Guangdong 524057, China
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 December 2015 Received in revised form 26 May 2016 Accepted 1 June 2016 Available online 6 June 2016
Transient rate decline analysis (TRDA) is an important analysis in modern reservoir engineering and applications. Most studies have focused on the threshold pressure gradient (TPG) of low permeability reservoirs and have placed an emphasis on the transient pressure response or steady productivity, whereas transient production decline has not been a focus. This paper presents a transient production decline analysis model of a horizontal well in low permeability gas reservoirs considering the threshold pressure gradient. An accurate solution is determined by employing the Laplace transformation, SturmLiouville eigenvalue method, and orthogonal transformation. The bi-logarithmic type curves of dimensionless production and the derivative are plotted by the Stehfest numerical inversion method. Five different flow regimes were recognized, and the effects of the influencing factors, such as the threshold pressure gradient, wellbore storage coefficient, skin factor, horizontal section length, and formation thickness, are discussed. This research contributes to the understanding of transient flow behavior and provides the theoretical basis for efficiently exploiting low permeability reservoirs. © 2016 Elsevier B.V. All rights reserved.
Keywords: Transient rate decline analysis Horizontal well Orthogonal transformation Low permeability
1. Introduction As a hot research topic in modern reservoir engineering and applications, transient rate decline analysis (TRDA) has become the focus of study in recent years. TRDA can quantitatively calculate the physical parameters of gas wells and reservoirs and can also be used to monitor reservoir dynamics (Aminian and Ameri, 1989). Horizontal wells are suitable for the development of low permeability reservoirs. TRDA of horizontal wells in low permeability gas reservoirs has both theoretical and practical significance to understand fluid transport behavior characteristics and correctly guide reservoir production. Numerous studies on the production rate transient analysis of horizontal wells have been performed. Ehlig-Economides and Ramey Jr. (1981) studied the transient rate response of a well that
* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (L.-N. Cao), [email protected] (X.-H. Tan). http://dx.doi.org/10.1016/j.jngse.2016.06.001 1875-5100/© 2016 Elsevier B.V. All rights reserved.
produced at a constant pressure according to three outer boundary conditions and several values of the wellbore skin factor by using the numerical Laplace transformation inversion algorithm. Fraim and Wattenbarger (1987) introduced the concept of normalized time for a gas reservoir producing against a constant wellbore pressure during (external) boundary-dominated flow. Therefore, type curve matching can be applicable to a reservoir with any shape. Aminian and Ameri (1989) discussed the application of the production type curves and issues concerning the production performance of horizontal wells in low permeability gas reservoirs. They obtained consistent results between the pre-stimulation data and type curve. Poon (1991) developed production decline curves for horizontal wells under various reservoir boundary conditions using Green’s and source functions. He noted that the wellbore stimulation, well spacing, and well length are probably the most important factors determining the effectiveness of a horizontal well. Doublet et al. (1994) presented a type of method to perform a decline curve analysis for both the variable rate and variable bottom hole pressure cases, without regard to the structure of the reservoir (shape and size) or the reservoir drive mechanisms. The
704
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
parameters determined using this method include the skin factor for near-well damage or stimulation, formation permeability, reservoir drainage area, and original oil-in-place. Accounting for the relative permeability and capillary pressure, Kewen and Horne (2003) proposed a decline analysis model based on fluid flow mechanisms. The model reveals a linear relationship between the oil production rate and reciprocal of the oil recovery or accumulated oil production. Medeiros et al. (2007) employed a semi-analytical model that incorporates the key features of reservoir heterogeneity as well as the details of hydraulic fracture and the wellbore flow to compute the production decline. They proposed a transient productivity index to represent the production decline characteristics. Ozkan et al. (2009) presented a discussion of the fractured horizontal-well performance in conventional (milli-Darcy permeability) and unconventional (micro-to nanoDarcy permeability) reservoirs and provided an interpretation of the different production performances of fracturing horizontal wells in both types of formations. Cai et al. (2010, 2012) (Cai and Yu, 2011)analyzed the capillary-driven flow in gas-saturated porous media based on fractal theory and modified the classical LucasWashburn equation. The factors influencing the imbibition process upon approaching the equilibrium weight were also analyzed. Later, the hydraulic conductivity and effective porosity in porous media were investigated. Nie et al. (2011) proposed the negative skin problem in a transient well test and rate decline analysis; they established a new mathematical model for a horizontal well in a multiple-zone composite reservoir. The contribution of their work was the removal of the oscillation in type curves, thus solving the negative skin problem. Based on the widely familiar Arp’s equation, Bahadori (2012) proposed a simple-to-use method, which was formulated to arrive at an appropriate estimation of the nominal (initial) decline rate and Arp’s decline-curve exponent. This method is quite simple to apply and accurately generates the coefficients of the equations instead of using already generated coefficients with uncertainty. Without applying the empirical concepts from Arps decline models and the explicit calculations from pseudofunctions, Ayala and Ye (2013) suggested a new generation analytical decline equation for boundary-dominated flow in a gas well producing at constant pressure for both full and partial production potential. Tan et al. (2013) created a dual fractal reservoir transient flow model by embedding a fracture system simulated by a tree-shaped fractal network into a matrix system simulated by fractal porous media. They plotted the bilogarithmic type curves of the dual fractal reservoirs and discussed the influence of different fractal factors on the pressure transient responses. Xie et al. (2014, 2015) focused on the two-phase flow caused by flow back after hydrofracturing shale gas reservoirs and established a two-phase pressure transient analysis model of a multi-stage fractured horizontal well. They studied not only the transient pressure but also the transient production rate decline features. Li et al. (2014) established a characteristic value method of the well test analysis for horizontal gas wells. Following the characteristic lines that are manifested by the pseudo-pressure derivative curve of each flow period, they developed formulas to determine the horizontal permeability, vertical permeability, skin factor, reservoir pressure, and pore volume of the gas reservoir. However, flow in low permeability reservoirs does not obey Darcy’s law due to the threshold pressure gradient (TPG), which makes the flow mechanism of horizontal wells in low permeability reservoirs more complicated. The concept of a threshold pressure gradient was suggested by Florin in 1951 (Florin, 1965). In the middle of the 20th century, numerous scholars concluded that the non-Darcy percolation that exists in low permeability reservoirs is
caused by the threshold pressure gradient (TPG), and they investigated the flow problems related with TPG (Miller and Low, 1963; Thomas et al., 1968; Das, 1997; Boukadi et al., 1998; Bennion et al., 2000; Dey, 2007). Guo et al. (2012) presented a horizontal gas well model accounting for the threshold pressure gradient in tight gas reservoirs. They note that the effect of the threshold pressure gradient on the pressure behavior for horizontal gas wells mainly occurs during the mid-late time periods. Song et al. (2014) studied the steady productivity equations of fractured wells in waterbearing tight gas reservoirs under low-velocity non-Darcy flow. The influence of the threshold pressure gradient on the volumetric flow rate of a gas well is significant and should not be ignored. Zhao et al. (2014) established a nonlinear steady flow mathematical model of a horizontal well for oil-water two-phase flow, accounting for the permeability stress sensitivity and threshold pressure gradient. Xu et al. (2015) conducted TPG and stress sensitivity experiments and investigated the existing condition of the TPG and change rule of permeability modulus. Then, they presented a transient pressure analysis mathematical model for a horizontal well in low-permeability offshore reservoirs. Liu and Yao (2015) constructed a moving boundary model of the radial flow in lowpermeability reservoirs, considering the influence of the threshold pressure gradient and quadratic pressure gradient term. The model was solved using numerical methods because it is nonlinear. Cai (2014) and Tan et al. (2015) investigated the low velocity non-Darcy flow, which starts the pressure gradient based on fractal theory. Based on the above discussion, most studies on the TPG of low permeability reservoirs were mainly focused on the transient pressure or steady productivity, and solutions of an unsteady model were determined using point source functions theory following Ozkan and Raghavan (1991) and Ozkan et al. (2010). This paper presents a transient production decline analysis model of horizontal wells in low permeability gas reservoirs accounting for the threshold pressure gradient. An accurate solution is determined by applying the Laplace transformation, SturmLiouville eigenvalue method, and orthogonal transformation. The bilogarithmic type curves of the dimensionless production and derivative are plotted using the Stehfest numerical inversion method. The characteristics of the transient rate decline behavior are thoroughly analyzed, and different flow regimes are observed in each type of curve. In addition, the effects of the relevant parameters are discussed, which improves the understanding of transient flow behavior and efficient exploitation of low permeability reservoirs.
2. Mathematical model A physical model of a low permeability gas reservoir, as shown in Fig. 1, is based on fluid flow in a horizontal well. To make the mathematical model more tractable and easier to understand, the following hypothesis and descriptions are outlined: (1) The reservoir is homogeneous with impermeable top and bottom boundaries, and the lateral boundary is infinite. The initial formation pressure is pi and is equal everywhere; (2) The reservoir is anisotropic, and the horizontal and vertical permeability are Kh and Kv, respectively; (3) The horizontal section of the well is parallel to the closed top and bottom boundaries, and its location is unlimited and represented by zw; (4) The single-phase gas flow obeys a non-Darcy law that accounts for the threshold pressure gradient;
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
705
Fig. 1. Sketch map of a horizontal gas well in a formation (revised from Xiao-ping Li et al. (2014)).
(5) A slightly compressible gas has a constant viscosity and compressibility coefficient, and the gas flow is under isothermal conditions; (6) The influence of gravity and the capillary force are ignored, but both the wellbore storage effect and the skin effect are taken into consideration.
obeyed only if the pressure gradient reaches a certain limit value (threshold pressure gradient), lB (Fig. 2). There is no longer a linear relationship between the flow velocity and pressure gradient because of the threshold pressure gradient. The flow velocity can be expressed by (Zhao et al., 2014):
.
v ¼
2.1. Basic differential equation of flow 2.1.1. Mass conservation equation Based on the mass conservation principle, the generalized form of the continuity equation for single-phase fluid flow can be written as:
vðrfÞ . ¼0 V$ r v þ vt
(1)
In cylindrical coordinates, the partial differential form of Eq. (1) is:
! ! ! v ¼ vr r þ vq q þ vz z
(2)
vðrvr Þ rvr 1 vðrvq Þ vðrvz Þ vðrfÞ þ ¼ þ þ vr r vq vz vt r
(3)
.
8 > < K ðVp lB Þ; > :
m
0;
jVpj > lB
(6)
jVpj lB
where lB is the threshold pressure gradient, MPa/m. In cylindrical coordinates, the components of Eq. (6) is:
8 Kh vp > > ¼ 3:6 l v > r B < m vr > Kv vp > > : vz ¼ 3:6 lB m vz
(7)
Gas flow only exists along the radial r direction and vertical z direction around the horizontal well in the reservoir. There is no rotational motion along the q axis, which is represented by vq ¼ 0. Thus, the flow velocity vector (2) and the continuity equation (3) can be simplified to:
! ! v ¼ vr r þ vz z
(4)
vðrvr Þ rvr vðrvz Þ vðrfÞ þ ¼ þ vr vz vt r
(5)
.
2.1.2. Non-Darcy equation The pores of a low permeability reservoir are very small. The thickness of the fluid adsorbed layer on the wall of pores is as large as the size of the pore, and the flow resistance of the adsorbed layer is larger still, so fluid flow is difficult when the pressure gradient is not large enough. In a low permeability reservoir, Darcy’s law is
Fig. 2. Non-Darcy flow relationship with the threshold pressure gradient in a low permeability gas reservoir.
706
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Table 1 Definitions of the dimensionless variables. Dimensionless pseudo pressure
hh mD ¼ 78:489K ðmi mÞ Tqsc
Dimensionless pseudo threshold pressure gradient
h hL lmBD ¼ 78:489K lmB Tqsc
Dimensionless time
3:6Kh t tD ¼ 4m 2 Ct rw CD ¼ 2p4CC hL2 t qffiffiffiffi ffi LD ¼ hL KKhv qffiffiffiffiffi hD ¼ rhw KKhv
Dimensionless wellbore storage coefficient Dimensionless horizontal section length Dimensionless gas reservoir thickness Dimensionless radial distance Dimensionless outer boundary distance
rw L
Dimensionless horizontal section position
zwD ¼ zhw
f ¼ f0 1 þ Cf ðp p0 Þ
Mp RTZ
2p
mZ
lB
(18)
Eq. (17) can then be rewritten as:
(8)
(9)
All of the dimensionless definitions are shown in Table 1. The dimensionless flow differential equation, converted from Eq. (19), is:
The density of natural gas can be expressed in another form according to the equation of state for real gas:
r¼
lmB ¼
1 v vm l Kv v2 m fmCt vm r mB þ ¼ r vr vr r Kh vz2 3:6Kh vt
The gas state equation:
r ¼ r0 1 þ Cr ðp p0 Þ
(17)
The pseudo threshold pressure gradient is introduced as follows:
2.1.3. State equation Porous rock and gas are both elastic and compressible, so both the rock and gas state equations should be considered. The rock state equation:
The value of lB is commonly on the order of magnitude of 103. p v As a quadratic gradient term, the product term of lB and vr mZ or p v vz mZ can be ignored for the sake of simplicity. Therefore, Eq (16)
1 v vm 1 2p K v2 m fmCt vm r lB þ v 2 ¼ r vr vr r mZ Kh vz 3:6Kh vt
reD ¼ rLe rwD ¼ zD ¼ hz
(16)
can be rearranged into:
rD ¼ Lr
Dimensionless vertical distance
Dimensionless wellbore radius
1 v vm 2p Kv v vm 2p fmCt vm r lB lB ¼ þ r vr vr mZ Kh vz vz mZ 3:6Kh vt
1 v vm 1 v2 m vm rD D þ lmBD þ L2D 2D ¼ ðhD LD Þ2 D rD vrD rD vrD vtD vzD
(19)
(20)
(10) 2.2. Definite conditions
2.1.4. Dimensionless differential equation of flow Combining Eq. (5) and Eq. (7)e(10) and simplifying leads to the real gas flow differential equation that accounts for the TPG:
Kh v p vp v p vp fmCt p vp r lB lB þ Kv ¼ vz mZ vz r vr mZ vr 3:6 mZ vt (11) where the total compressibility coefficientCt ¼ Cr þ Cf , MPa1. The pseudo pressure is introduced in the following mathematical model:
Zp mðpÞ ¼ p0
2p dp mZ
Initial condition:
mD ðrD ; zD ; 0Þ ¼ 0
(21)
Inner boundary condition:
2
ε
lim 4 lim
εD /0
rD /0
εD 2
D zwD Zþ2
rD
ε zwD 2D
vmD þ lmBD vrD
3 1 dzwD 5 ¼ ; jzD zwD j 2
(22) Infinite lateral boundary condition:
(12)
mD ðrD /∞; zD ; tD Þ ¼ 0
(23)
Closed top and bottom boundaries: where p0 is the reference pressure and can be variable. Therefore:
vm 2p vp ¼ vr mZ vr
(13)
vm 2p vp ¼ vz mZ vz
(14)
vm 2p vp ¼ vt mZ vt
(15)
by
Combining Eq. (13)e(15) into Eq. (11) and multiplying both sides 2 Kh, we obtain:
vmD
vmD
¼ 0; ¼0 vzD zD ¼0 vzD zD ¼1
(24)
2.3. Problem for definite solutions By simultaneously solving Eq. (20)e(24), the mathematical flow model of a horizontal well in a low permeability gas reservoir can be written as follows:
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
707
8 1 v vm 1 v2 m vm > rD D þ lmBD þ L2D 2D ¼ ðhD LD Þ2 D > > rD vrD vtD > rD vrD vz > D > > > > > mD ðrD ; zD ; 0Þ ¼ 0 > > > > εD > 2 3 > zwD > Zþ2 < vmD 1 ε 4 þ lmBD rD dzwD 5 ¼ ; jzD zwD j D lim lim > 2 vr 2 ε /0 r /0 D D > D > ε > zwD 2D > > > > > > > mD ðrD /∞; zD ; tD Þ ¼ 0 > > >
> > vmD
> : vmD
¼ 0; ¼0 vzD zD ¼0 vzD zD ¼1
(25)
3.2. Sturm-Liouville eigenvalue problem
3. Solutions of the mathematical model Importantly, Eq. (20) is so nonlinear that it is hard to calculate an exact solution by conventional methods. For the initial boundary value problem (25), the first step is to apply the Laplace transform to the time variable. The second step is to apply the orthogonal transform to the spatial variable, thus reducing the nonlinear degree of the problem to help calculate the final solution.
The differential equation of model (28) can be rewritten as:
v2 mD 1 vmD 1 RmD þ þ f ðrD ; zD Þ ¼ 2 rD vrD gðzD Þ vrD
where R is the differential operator, and pðzD Þ, qðzD Þ, and gðzD Þ are all real functions on [0,1]. The expressions are:
3.1. Laplace transformation The Laplace transformation on mD ðrD ; zD ; tD Þ from tD to s defines mD as the Laplace transformation function of mD , that is:
Zþ∞ mD ðrD ; zD ; sÞ ¼ L½mD ðrD ; zD ; tD Þ ¼
mD ðrD ; zD ; tD ÞestD dtD
0
(26) Applying the initial conditions, the derivative of mD with respect to tD in the Laplace domain can be written as:
dmD ðrD ; zD ; tD Þ ¼ s$L½mD ðrD ; zD ; tD Þ mD ðrD ; zD ; tD ÞjtD ¼0 L dtD
(29)
R¼
v v pðzD Þ qðzD Þ vzD vzD
(30)
pðzD Þ ¼ L2D
(31)
qðzD Þ ¼ 0
(32)
gðzD Þ ¼ 1
(33)
f ðrD ; zD Þ ¼
1 l ðhD LD Þ2 smD srD mBD
(34)
The integral transformation on the space variable zD is as follows:
¼ smD (27) The dimensionless mathematical model (25) in the Laplace domain can be written as:
MðrD ; mÞ ¼ F½mD ðrD ; zD Þ ¼
Z1
gðzD ÞZðzD ; mÞmD ðrD ; zD ÞdzD
(35)
0
where gðzD ÞZðzD ; mÞ is known as the integral transformation kernel.
8 1 v vm 1 v2 m > rD D þ lmBD þ L2D 2D ¼ ðhD LD Þ2 smD > > r vr sr vr > D D vzD D D > > > > εD > 2 3 > zwD > Zþ2 > > vmD lmBD 1 ε > < lim 4 lim dzwD 5 ¼ ; jzD zwD j D rD þ 2s vrD s 2 εD /0 rD /0 ε > zwD 2D > > > > > > > mD ðrD /∞; zD ; sÞ ¼ 0 > > >
> > > vmD
: vmD
¼ 0; ¼0
vzD zD ¼0 vzD zD ¼1
(28)
708
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
The Sturm-Liouville theorem solves an eigenvalue problem that is composed of a second-order variable coefficient differential equation and corresponding boundary conditions. To determine the transformation kernel, we need to solve a Sturm-Liouville eigenvalue problem relating to ZðzD ; mÞ, which is:
lmBD
F
Z1 ¼
s
lmBD s
0
8 < lmBD ; n ¼ 0 s cos npzD dzD ¼ : 0 ; ns0
(41)
"
8 lgZ; 0 < zD < 1 > < RZ ¼
dZ
dZ
> ¼ ¼0 :
dzD zD ¼0 dzD zD ¼1
(36)
# Z1 v2 mD v2 mD F cos npzD dzD ¼ n2 p2 mD ¼ vz2D vz2D
(42)
1 v vmD 4 rD þ n n2 p2 L2D mD ¼ ðhD LD Þ2 smD rD vrD vrD rD
(43)
8 < lmBD ; n ¼ 0 s , the differential equation of Defining 4n ¼ : 0 ; ns0 model (28),can be orthogonally transformed to: 0
Merging the constant coefficient, Eq. (36) can be simplified to:
8 00 > < Z þ lZ ¼ 0; 0 < zD < 1
dZ
> dZ
¼ ¼0 : dzD zD ¼0 dzD zD ¼1
(37)
where l is the corresponding eigenvalue when Eq. (37) has a nonzero solution. Characteristic equation: r 2 þ l ¼ 0 ①When l ¼ 0, r1 ¼ r2 ¼ 0, then ZðzD Þ ¼ c1 þ c2 zD . Z’ð0Þ ¼ 00c2 ¼ 0, naturally, meet the requirement of Z’ð1Þ ¼ 0. So ZðzD Þ ¼ c1 ðc1 s0Þ. pffiffiffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffi lzD . ②When l < 0, r1;2 ¼ ± l, then ZðzD Þp¼ffiffiffiffiffic1 e lzDpþ ffiffiffiffiffic2 e Z’ð0Þ ¼ 00c1 c2 ¼ 0; Z’ð1Þ ¼ 00c1 e l c2 e l ¼ 0. Solve simultaneously and obtainc1 ¼ c2 ¼ 0, such that ZðzD Þ has a non-zero solution. pffiffiffi ③When pffiffiffi l > 0,pffiffiffi r1;2 ¼ ± li, then ZðzD Þ ¼ c1 cos lzD þ c2 sin lzD . Z’ð0Þ ¼ 00c2p¼ffiffiffi 0; pffiffiffi Z’ð1Þ ¼ 00c1 sin l ¼ 00 l ¼ np; n ¼ 1; 2; 3:::. So ln ¼ ðnpÞ2 ðn ¼ 1; 2; :::Þ, Zn ðzD Þ ¼ cn cos npzD . According to the above discussion on the S-L problem (37), the corresponding eigenvalues are ln ¼ ðnpÞ2 , and the corresponding eigenfunctions are Zn ðzD Þ ¼ cos npzD (where n ¼ 0; 1; 2:::), which are composed of a weighted function gðzD Þ into a complete orthogonal system in the [0,1] domain. Thus, the integral transformation kernel is:
gðzD ÞZðzD ; mÞ ¼ 1$Zn ðzD Þ ¼ cos npzD ; n ¼ 0; 1; 2:::; m ¼ n
(38)
The integral transformation whose kernel is a complete orthogonal function system can collectively be known as the orthogonal transformation.
The next step is to take the orthogonal transformation based on the inner boundary conditions. The right side of the inner boundary 1 , which is independent of ε . Therecondition of model (28) is 2s D fore, the following expression is established:
1 1 ¼ lim 2s εD /0 2s
(44)
Then, the inner boundary condition can be written as:
2
ε
D zwD Zþ2
lim 4 lim
εD /0
rD
rD /0
ε
¼ lim
εD /0
zwD 2D
vmD lmBD þ vrD s
3 dzwD 5
1 ε ; jzD zwD j D 2s 2
(45)
Simplifying Eq. (45), we obtain:
lim
lim εD rD
εD /0 rD /0
vmD 1 ε ; jzD zwD j D ¼ lim 2s vrD 2 εD /0
Because the integral calculation of the variable ZD is independent of variable εD, the orders of the orthogonal transformation and the limit calculation of εD can be exchanged, leading to:
ℱ
lim
εD /0
lim εD rD
rD /0
vmD vrD
¼ℱ
lim
εD /0
1 2s
i
vm 1 lim ℱ lim εD rD D ¼ lim ℱ 2s vrD εD /0 rD /0 εD /0
3.3. Orthogonal transformation
(46)
(47)
Because: The orthogonal transformation of variable zD , with eigenfunctions cos npzD as the kernel and mD as the orthogonal transformation function of mD , can be written as:
Z1 mD ðrD ; n; sÞ ¼ F½mD ðrD ; zD ; sÞ ¼
Based on the differentiability properties, we obtain:
Z1 vmD v vmD ¼ mD ðrD ; zD ; sÞcos npzD dzD ¼ vrD vrD vrD 0
rD /0
1 ℱ 2s (39)
lim εD rD
mD ðrD ; zD ; sÞcos npzD dzD
0
F
ℱ
(40)
vmD vrD
¼ lim εD rD rD /0
vmD vrD
(48)
ε
D zwD Zþ2
¼
ε
cos npzD dzD 2s
zwD 2D
ε ε sin np zwD D sin np zwD þ D 2 2 ¼ 2snp 3 2 cos npzwD sin np D 2 ¼ 2snp Thus:
(49)
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
lim εD rD
rD /0
cos npzwD sin np ε2D vmD ¼ vrD snp
(50)
Because variable εD is independent of rD, moving εD to the right side of Eq. (50), we obtain:
lim rD
rD /0
cos npzwD sin np ε2D vmD ¼ vrD snpεD
εD /0
lim rD
rD /0
vmD vrD
¼ lim
εD /0
cos npzwD sin np ε2D snpεD
(51)
rD /0
cos npzwD sin np ε2D vmD ¼ lim vrD snpεD εD /0
(52)
lim
εD /0
cos npzwD sin np snpεD
εD 2
ε D d cos npz wD sin np dðsnpεD Þ 2 ¼ lim dεD dεD εD /0 εD np cos npzwD cos np 2 ¼ cos npzwD ¼ lim 2s 2snp εD /0
(54)
lim rD
vmD cos npzwD ¼ vrD 2s
v2 mD 2 vrD
¼
v vrD
vmD vx $ vx vrD
(55)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v vmD $ n2 p2 L2D þ h2D L2D s vrD vx qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 v mD vx $ ¼ n2 p2 L2D þ h2D L2D s$ vx2 vrD v2 m D ¼ n2 p2 L2D þ h2D L2D s vx2
¼
(61) Combining Eq. (60) and Eq. (61) into Eq. (59), we have:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2 m vmD D 2 n2 p2 L2D þ h2D L2D s rD þ r n2 p2 L2D þ h2D L2D s D vx vx2 2 2 2 2 2 2 rD n p LD þ hD LD s mD ¼0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mD ¼ AI0 rD n2 p2 L2D þ h2D L2D s þ BK0 rD n2 p2 L2D þ h2D L2D s (62)
(56)
Simultaneously solving Eqs. (43), (55), and (56) for a mathematical model on mD can be written as:
8 > > > > > > > > <
4 vmD 1 v vmD þ n n2 p2 L2D mD ¼ ðhD LD Þ2 smD lim rD rD > rD vrD rD /0 r vr vr > D D D > > > > > > :
where I0 is the first zero order modified Bessel function, K0 is the second zero order modified Bessel function, and A and B are constants. If n ¼ 0, 4n ¼ lmBD s s0, the basic differential equation of model (57) is inhomogeneous because its free term is not equal to zero. The particular solution can be expressed by the Green functions. The general solution can be written as:
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi mD ¼ AI0 rD n2 p2 L2D þ h2D L2D s þ BK0 rD n2 p2 L2D þ h2D L2D s Zþ∞ þ
GðrD ; tÞdt
0
cos npzwD ¼ mD ðrD /∞; n; sÞ ¼ 0 2s
(63) (57)
If ns0 and 4n ¼ 0, the basic differential equation of model (57) can be converted into a zero order modified Bessel equation. The derivation process is shown in the following steps:
v2 mD 1 vmD 2 2 2 þ n p LD þ h2D L2D s mD ¼ 0 2 rD vrD vrD
(60)
2
Taking the orthogonal transformation on the outer boundary conditions of model (28):
mD ðrD /∞; n; sÞ ¼ 0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi vmD vmD vx vmD $ $ n2 p2 L2D þ h2D L2D s ¼ ¼ vrD vx vrD vx
2 2 D Namely: x2 vvxm2D þ x vm vx ðx þ y ÞmD ¼ 0, where y ¼ 0 and y represents the order number of the modified Bessel equation. The general solution of the equation is written as:
Therefore, the inner boundary condition is simplified to:
rD /0
n2 p2 LD þhD LD s
(53)
Applying the L’Hospital Rule to the right side of Eq. (53) is “00” uncertain limit type, that is:
(59)
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x ffi. Following the Let x ¼ rD n2 p2 L2D þ h2D L2D s, so rD ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2
In addition, because variable εD is independent of rD, Eq. (52) can be rewritten as:
lim rD
v2 mD vmD 2 2 2 2 2 2 þ r n p L þ h L s rD mD ¼ 0 D D D D 2 vrD vrD
compound function derivation rule, we obtain:
Taking the limit calculation of εD on both sides of Eq. (51):
lim
2 rD
709
(58)
2 , and Eq. (58) can be Both sides of Eq. (58) are multiplied byrD rewritten as:
where the Green function
R∞ 0
GðrD ; tÞ is:
qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi 8 lmBD 2 2 2 2 > > K ð0 < t < rD Þ > 0 rD hD LD s I0 t hD LD s < s GðrD ; tÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi > lmBD qffiffiffiffiffiffiffiffiffiffiffiffiffi > > : K0 t h2D L2D s I0 rD h2D L2D s ðrD < t < þ ∞Þ s (64) Combining Eq. (56) with an infinite outer boundary condition and I0 ð∞Þ ¼ ∞, we obtain A ¼ 0. So, Eq. (63) is rewritten as:
710
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
mD ¼ BK0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Zþ∞ rD n2 p2 L2D þ h2D L2D s þ GðrD ; tÞdt
(65)
0
Combining the inner boundary conditions of Eq. (55) and the properties of the Bessel functions, if z/0, then I1 ðzÞy0, so npzwD B ¼ cos 2s . Replacing B into Eq. (65) and taking the inverse orthogonal transformation, we determine the pressure solution in the Laplace domain:
qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 mD ¼ K0 rD h2D L2D s 2s qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∞ 1X K0 rD n2 p2 L2D þ h2D L2D s cos npzwD cos npzD þ s n¼1 Zþ∞ þ
GðrD ; tÞdt
0
(66) Following the source function idea method, let yD ¼ 0, and use Eq. (66) to mark the point-source solution. Then, compute the integral for the horizontal section of the wellbore to obtain the bottom hole pseudo pressure solution expression, as follows, that accounts for the effect of the threshold pressure gradient:
mD ¼
1 2s
Z1
qD ðsÞ ¼
1 s2 mwD ðsÞ
(69)
Following Eq. (67)e(69), we can determine the horizontal well dimensionless rate solution under the influence of the threshold pressure gradient for low permeability gas reservoirs. 4. Type curves and discussion 4.1. Rate decline performance type curve analysis Type curves of the production rate can indicate the characteristics of the flow behavior to some degree and can be used to apply a curve-fitting analysis to determine the property parameters of the well and the formation. The solutions of the dimensionless production rate in real space can be determined by applying the Stehfest numerical inversion algorithm (Stehfest, 1970). The whole production decline behavior of a horizontal well with low permeability is plotted in Fig. 3. The following five flow regimes can be recognized: Regime I: The pure wellbore storage effect flow regime. This stage occurs early, and the curves of the production rate and the derivative have the same trend, which is described by a downward sloping line, similar to a conventional reservoir. The gas stored in the wellbore controls this period.
Zþ∞ Z1 qffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj h2D L2D s da þ GðjxD aj; tÞdadt
1
0
1
1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∞ Z 1X þ K0 jxD aj n2 p2 L2D þ h2D L2D s cos npzwD cos npzD da s n¼1
(67)
1
Because the wellbore has a certain volume and natural gas has compressibility, the wellbore storage effect occurs when the well is open or shut and is represented by the wellbore coefficient. In addition, the skin effect around the bottom hole describes the situation of formation damage or stimulation treatment and is represented by the skin factor. When taking these two effects into account, the definite conditions will become nonhomogeneous. Duhamel’s theorem (or homogenized principle) can be used to divide the nonhomogeneous problem into homogeneous subproblems that are easily solved. The solutions to the subproblems are then combined to yield a solution to the original problem, which is named the superposition principle. Applying Duhamel’s theorem and the superposition principle to incorporate the wellbore storage coefficient and the skin factor into the well response, the pseudo pressure expression can be rewritten as:
mwD ¼
smD þ S s þ CD s2 ðsmD þ SÞ
(68)
The dimensionless wellbore flow rate equation for the constantpressure production case in low permeability gas reservoirs can be determined by the relationship between the dimensionless pressure and the rate in the Laplace space (Van Everdingen and Hurst, 1949):
Regime II: The skin effect transition flow regime. The drop in the degree of the production rate curve slowly changes, and the rate derivative curve continues to decline with a concave shape, which corresponds to the hump on the transient pressure curve. This stage reflects the combined action of the skin and wellbore storage. Regime III: The early radial flow regime in the vertical plane. The production rate curve continues to decrease, and the derivative curve is a level straight line, with a value related to LD. This period reflects the radial flow feature perpendicular to the horizontal axis before the pressure wave propagates to the top and bottom outer boundary (Fig. 4). Regime IV: The mid-term linear flow regime. The curves of the production rate and derivative are two parallel lines with the same slope. During this flow stage, the pressure wave has already propagated to the top and bottom boundary, and the natural gas flows linearly perpendicular to the horizontal section (Fig. 5). The influence of the threshold pressure gradient also appears in this period. Regime V: The late pseudo-radial flow regime affected by the TPG. The curves of the production rate and the derivative describe the pseudo-radial flow characteristics of the whole system. The effect of the threshold pressure gradient mainly exists in the fifth flow regime. It leads to a downwarping of the production curve and an upwarping of the rate derivative curve. This stage represents the radial flow characteristic of the entire system (Fig. 6).
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
711
Fig. 3. Type curves of the production rate decline for a horizontal well in a low permeability gas reservoir (solid lines represent production rate curves, dash lines represent the production rate derivative curves, red lines represent lmBD ¼ 0, blue lines represent lmBD ¼ 0.5). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. The early radial flow regime in the vertical plane of the horizontal well.
Fig. 5. The mid-term linear flow regime of the horizontal well.
4.2. Parameter influence on rate decline performance Fig. 7 shows the effect of the threshold pressure gradient on the production rate decline type curves. The figure suggests that the TPG mainly plays its role during the mid-term linear regime to the
Fig. 6. The late pseudo-radial flow regime of the horizontal well.
late pseudo-radial flow regime. The greater the TPG, the worse the reservoir properties, the more difficult it is for the fluid to flow, and the lower the production rate will be, resulting in a steeper downwarping feature on the production rate curve. Correspondingly, the faster the production curve decreases, the more obvious the upwarping of the rate derivative curve. Fig. 8 shows the effect of the wellbore storage coefficient on the production rate decline type curves. The wellbore storage coefficient has a major influence on the first two flow regimes: the pure wellbore storage effect regime and the skin effect transition regime. The greater the wellbore storage coefficient, the longer the wellbore storage period lasts and the less obvious the concavity of the derivative curves. When CD exceeds a certain value, the vertical radial flow regime might be concealed. Fig. 9 illustrates the effect of the skin factor on the production rate decline type curves. The skin factor reflects the amount of reservoir pollution near the wellbore area, and it mainly affects the skin transition flow regime. The incremental value of the skin factor results in increasing the additional filtration resistance. The bigger the skin factor, the longer the skin effect transition period lasts and the more obvious the concavity of the production rate derivative curves. After the skin effect transition regime, the impact of the skin effect is gradually removed.
712
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Fig. 7. Type curves of the production rate decline affected by the threshold pressure gradient (red line represents lmBD ¼ 0, blue line represents lmBD ¼ 0.1, green line represents lmBD ¼ 0.5).
Fig. 8. Type curves of the production rate decline affected by the wellbore storage coefficient (red line represents CD ¼ 104, blue line represents CD ¼ 103, green line represents CD ¼ 102).
Fig. 10 shows the effect of the horizontal section length of the well on the production rate decline type curves. Increasing the horizontal section length results in a shift upward of both the production rate curve and the rate derivative curve. This shift occurs because the longer the horizontal section, the larger the gas drainage area. The resistance of the gas flow becomes relatively small and the corresponding production rate will be faster. Moreover, the greater the value of LD, the longer the linear flow period lasts, and the next pseudo-radial flow regime will be postponed. Fig. 11 shows the effect of reservoir thickness on the production rate decline type curves. A change in the reservoir thickness has a significant influence on the type curves. When the reservoir thickness increases, the pressure wave needs more time to propagate to the top and bottom boundaries. Consequently, the skin transition flow period will have a more obvious V-shaped concavity, and the early vertical radial flow regime will longer. Moreover, the mid-term linear flow stage and the late pseudo-radial flow stage are delayed.
5. Model application with field data The JB4-10H well is a horizontal development well in the JingBian gas field in China. The total well depth is 3309 m, and the well is drilled in the Majiagou Group. After drilling is finished, the well is completed by perforation and acidification. The reservoir thickness is 8.6 m, and the formation porosity is 6.72%. The gas gravity is 0.6148, the gas viscosity is 0.02146 mPa s, and the gas deviation factor Z is 0.9821. The well radius is 0.0797 m. According to the deliverability test that occurred from 7 to 11 December 2013, the calculated absolute open flow was 41.93 104 m3/d. Before production, the test data suggested that the initial formation pressure was 30.34 MPa, the reservoir temperature was 103.9 C, and the tubing pressure and casing pressure were 24.20 MPa and 24.10 MPa, respectively. The commissioning date of the JB4-10H well was 15 December 2013. In the early stages of production, the proration flow was 12 104 m3/d, which allowed the pressure to drop quickly and the tubing pressure and casing pressure were
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
713
Fig. 9. Type curves of the production rate decline affected by the skin factor (blue line represents S ¼ 0.01, red line represents S ¼ 0.1, green line represents S ¼ 0.2).
Fig. 10. Type curves of the production rate decline affected by the horizontal section length (red line represents LD ¼ 5, blue line represents LD ¼ 10, green line represents LD ¼ 20).
reduced to 18.70 MPa and 19.60 MPa, respectively. Then, the JB410H well was shut in to allow the formation pressure to build up. Starting on 8 February 2014, the JB4-10H well produced under constant bottom hole pressure for 458 days. The production rate log-log plot of the JB4-10H well is shown in Fig. 12. Substituting Eqs. (67) and (68) into Eq. (69) and applying the Stehfest numerical inversion algorithm, the expression of the production type curve in the real domain is written as:
(" Z qD ¼
N ln 2 X V tD i¼1 i
1 2
1 þ CD s (" Z s 12
1
1
1 1
Z qffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj h2D L2D s da þ s
Z qffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj h2D L2D s da þ s
0
þ∞ Z 1 0
þ∞ Z 1 1
1
where N is an even number and s ¼ i$ln 2=tD . The Green function GðrD ; tÞ is shown as Eq. (64). The weight coefficient Vi is given by:
Vi ¼ ð1Þ 2 þi N
k¼
GðjxD aj; tÞdadt þ
GðjxD aj; tÞdadt þ
minði;N=2Þ X
∞ Z X
1
n¼1 1
(71)
# ) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K0 jxD aj n2 p2 L2D þ h2D L2D s cos npzwD cos npzD da þ S
n¼1 1 1
∞ Z X
iþ1 2
k 2 þ1 ð2kÞ! N k !ðk!Þ2 ði kÞ!ð2k iÞ! 2 N
K0 jxD aj
# ) qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 p2 L2D þ h2D L2D s cos npzwD cos npzD da þ S (70)
714
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Fig. 11. Type curves of the production rate decline affected by the reservoir thickness (red line represents hD ¼ 100, blue line represents hD ¼ 200, green line represents hD ¼ 400).
Fig. 12. Log-log plots of the production rate of the JB4-10H well.
Table 2 Rate transient analysis results of the JB4-10H well. Parameter
Value 3
Wellbore storage coefficient C (m /MPa) Skin factor S Horizontal permeability Kh (mm2) Vertical permeability Kv (mm2) Effective horizontal section length L (m) Threshold pressure gradient lB (MPa/m)
2.028 1.625 8.5 104 3.6 105 334.31 7.82 103
The parameters of the theoretic type curves are CD ¼ 104, S ¼ 0.1, LD ¼ 5, hD ¼ 100, and lmBD ¼ 0.5. By applying the type curves, the rate transient analysis results of the JB4-10H well are shown in Table 2. 6. Summary and conclusions In this paper, a mathematical model for a horizontal well in a low permeability gas reservoir that considers the threshold
pressure gradient is established. The production rate decline behavior under the effects of relative influence factors is analyzed. The following conclusions can be summarized: (1) The TPG mainly affects the mid-term linear flow regime and late pseudo-radial flow regime of horizontal well production decline type curves. A larger TPG can obviously cause a downward shift of the production rate curve and an upward shift of the rate derivative curve. (2) A large wellbore storage coefficient results in a longer pure wellbore storage period, even though the vertical radial flow regime might be hidden. Simultaneously, the concavity of the production rate derivative curve becomes more obvious when the skin factor is greater, which means the formation near the wellbore area is more polluted. (3) The increasing gas drainage area is caused by a larger horizontal section length, which can increase both the production rate curve and rate derivative curve. Type curves suggest that the linear flow period lasts longer and the pseudo-radial flow appears later.
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
(4) The reservoir thickness has little influence on the entire curves except in the pure wellbore storage regime. The greater the reservoir thickness, the deeper the V-shaped concavity of the derivative curve and the longer the early vertical radial flow. Acknowledgment The authors would like to acknowledge the support of the National Basic Research Program (973 Program) of China under Grant No. 2014CB239205. Nomenclature
Latin symbol C wellbore storage coefficient, m3/MPa CD dimensionless wellbore storage coefficient Ct total compressibility coefficient, MPa1 Cr fluid compressibility coefficient, MPa1 C4 rock compressibility coefficient, MPa1 h reservoir thickness, m hD dimensionless reservoir thickness Kh horizontal permeability, mm2 Kv vertical permeability, mm2 L horizontal section length, m LD dimensionless horizontal section length M molecular mass, kg/mol m pseudo pressure, MPa2/(mPa$s) mD dimensionless pseudo pressure mi initial pseudo pressure, MPa2/(mPa s) p pressure, MPa p0 reference pressure, MPa pi initial formation pressure, MPa q gas production rate, 104 m3/d qD dimensionless gas production rate qsc surface gas production rate, 104 m3/d R gas constant (0.008314 MPa m3/(kmol K) r radial distance, m rD dimensionless radial distance re outer boundary distance, m reD dimensionless outer boundary distance rw wellbore radius, m rwD dimensionless wellbore radius s Laplace transform variable S skin factor, dimensionless T absolute temperature, K t time, hours tD dimensionless time v velocity of gas flow, m/h Z gas deviation factor z vertical distance, m zD dimensionless vertical distance zw horizontal section position, m zwD dimensionless horizontal section position Greek symbols ε tiny variable l eigenvalue, dimensionless lB threshold pressure gradient, MPa/m lmB pseudo threshold pressure gradient, MPa2/(mPa s m) lmBD dimensionless pseudo threshold pressure gradient m gas viscosity, mPa s r gas density, kg/m3 n order number of modified Bessel equation
f
715
porosity of reservoir, fraction
Superscripts Laplace transform domain ¼ Orthogonal transform domain Subscripts D dimensionless h horizontal i initial m pseudo sc standard condition t total v vertical w wellbore Abbreviations TRDA transient rate decline analysis TPG threshold pressure gradient References Aminian, K., Ameri, S., 1989. Predicting horizontal well production performance using type curves. In: Paper SPE 19342 Presented at the SPE Eastern Regional Meeting, Morgantown, West Virginia, 24e27 October. Ayala, L., Ye, P., 2013. Unified decline type-curve analysis for natural gas wells in boundary-dominated flow. SPE J. 18 (1), 97e113. Bahadori, A., 2012. Analysing gas well production data using a simplified decline curve analysis method. Chem. Eng. Res. Des. 90 (4), 541e547. Bennion, D.B., Thomas, F.B., Ma, T., 2000. Recent advances in laboratory test protocols evaluate optimum drilling, completion and stimulation practices for low permeability gas reservoirs. In: Paper SPE 60324 Presented at SPE Rocky Mountain Regional/Low-permeability Reservoirs Symposium and Exhibition, Denver, Colorado, 12e15 March. Boukadi, F., Bemani, A., Rumhy, M., Kalbani, M., 1998. Threshold pressure as a measure of degree of rock wettability and diagenesis in consolidated Omani limestone cores. Mar. Petroleum Geol. 15 (1), 33e39. Cai, J., 2014. A fractal approach to low velocity non-Darcy flow in a low permeability porous medium. Chin. Phys. B 23 (4), 044701. Cai, J., Hu, X., Standnes, D.C., You, L., 2012. An analytical model for spontaneous imbibition in fractal porous media including gravity. Colloids Surfaces A 414, 228e233. Cai, J., Yu, B., 2011. A discussion of the effect of tortuosity on the capillary imbibition in porous media. Transp. Porous Media 89 (2), 251e263. Cai, J., Yu, B., Zou, M., Luo, Cai, Yu, L., 2010. Fractal characterization of spontaneous co-current imbibition in porous media. Energy Fuels 24 (3), 1860e1867. Das, A.K., 1997. Generalized Darcy’s law including source effect. J. Can. Petroleum Technol. 36 (6), 57e59. Dey, G.R., 2007. Nonlinear flow of gas at low velocity in porous media. J. Nat. Gas Chem. 16 (3), 217e226. Doublet, L.E., Pande, P.K., McCollum, T.J., Blasingame, T.A., 1994. Decline curve analysis using type curvesdanalysis of oil well production data using material balance time: application to field cases. In: Paper SPE 28688 Presented at the International Petroleum Conference and Exhibition of Mexico, Veracruz, Mexico, 10e13 October. Ehlig-Economides, C.A., Ramey Jr., H.J., 1981. Transient rate decline analysis for wells produced at constant pressure. Soc. Petroleum Eng. J. 21 (1), 98e104. Florin, B.A., 1965. Principles of Soil Mechanics (Tongji University). China’s industrial press, Beijing, pp. 38e39. Original work published 1959. Fraim, M.L., Wattenbarger, R.A., 1987. Gas reservoir decline-curve analysis using type curves with real gas pseudopressure and normalized time. SPE Form. Eval. 2 (4), 671e682. Guo, J., Zhang, S., Zhang, L., Qing, H., Liu, Q., 2012. Well testing analysis for horizontal well with consideration of threshold pressure gradient in tight gas reservoirs. J. Hydrodynamics 24 (4), 561e568. Kewen, L., Horne, R.N., 2003. A decline curve analysis model based on fluid flow mechanisms. In: Paper SPE 83470 Presented at the SPE Western Regional/AAPG Pacific Section Joint Meeting, Long Beach, California, 19e24 May. Li, X.-P., Yan, N.-P., Tan, X.-H., 2014. Characteristic value method of well test analysis for horizontal gas well. Math. Problems Eng. 2014. Article ID 472728. 10 pages. Liu, W., Yao, J., 2015. Numerical investigations of the effect of nonlinear quadratic pressure gradient term on a moving boundary problem of radial flow in lowpermeable reservoirs with threshold pressure gradient. Math. Problems Eng. 2015. Article ID 275057, 12 pages. Medeiros, F., Kurtoglu, B., Ozkan, E., Kazemi, H., 2007. Analysis of production data from hydraulically fractured horizontal wells in shale reservoirs. SPE Reserv. Eval. Eng. 13 (3), 559e568. Miller, R.J., Low, P.F., 1963. Threshold gradient for water flow in clay system. Soil Sci.
716
L.-N. Cao et al. / Journal of Natural Gas Science and Engineering 33 (2016) 703e716
Soc. Am. J. 27 (3), 605e609. Nie, R., Guo, J., Jia, Y., Zhu, Sh, Rao, Zh, Zhang, Ch, 2011. New modeling of transient well test and rate decline analysis for a horizontal well in a multiple-zone reservoir. J. Geophys. Eng. 8 (3), 464e476. Ozkan, E., Brown, M.L., Raghavan, R.S., 2009. Comparison of fractured-horizontalwell performance in conventional and unconventional reservoirs. In: Paper SPE 121290 Presented at the SPE Western Regional Meeting, San Jose, California ,USA, 24e26 March. Ozkan, E., Raghavan, R., 1991. New solutions for well-test-analysis problems: part 1analytical considerations. SPE Form. Eval. 6 (3), 359e368. Ozkan, E., Raghavan, R.S., Apaydin, O.G., 2010. M odeling of fl uid transfer from shale matrix to fracture network. In: SPE 134830, Presented at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19e22 September. Poon, D.C., 1991. Decline curves for predicting production performance from horizontal wells. J. Can. Petroleum Technol. 30 (1), 77e81. Song, H., Liu, Q., Yang, D., Yu, M., Lou, Y., Zhu, W., 2014. 2014. Productivity equation of fractured horizontal well in a water-bearing tight gas reservoir with lowvelocity non-Darcy flow. J. Nat. Gas Sci. Eng. 18, 467e473. Stehfest, H., 1970. Algorithm 368: numerical inversion of Laplace transforms. Commun. ACM 13 (1), 47e49. Tan, X.-H., Li, X.-P., Liu, J.-Y., Tang, C., Li, J., 2013. 2013. Pressure transient analysis of
dual fractal reservoir. J. Appl. Math. 2013. Article ID 137518, 9 pages. Tan, X.-H., Li, X.-P., Zhang, L.-H., Liu, J.-Y., Cai, J., 2015. Analysis of transient flow and starting pressure gradient of power-law fluid in fractal porous media. Int. J. Mod. Phys. C 26 (4), 1550045. Thomas, L.K., Katz, D.L., Tek, M.R., 1968. Threshold pressure phenomena in porous media. Soc. Petroleum Eng. J. 8 (02), 174e184. Van Everdingen, A.F., Hurst, W., 1949. The application of the Laplace transformation to flow problems in reservoirs. Pet. Trans. AIME 1 (12), 305e324. Xie, W., Li, X., Zhang, L., Wang, J., Cao, L., Yuan, L., 2014. Two-phase pressure transient analysis for multi-stage fractured horizontal well in shale gas reservoirs. J. Nat. Gas Sci. Eng. 21, 691e699. Xie, W.-Y., Li, X.-P., Zhang, L.-H., Tan, X.-H., Wang, J.-C., Wang, H.-T., 2015. Production decline analysis for two-phase flow in multifractured horizontal well in shale gas reservoirs. J. Chem. 2015. Article ID 212103, 10 pages. Xu, J., Jiang, R., Teng, W., 2015. Nonlinear flow characteristics and horizontal well pressure transient analysis for low-permeability offshore reservoirs. Math. Problems Eng. 2015. Article ID 387149, 13 pages. Zhao, Y., Zhang, L., He, Z., Zhang, B., 2014. Productivity for horizontal wells in lowpermeability reservoir with oil/water two-Phase flow. Math. Problems Eng. 2014. Article ID 364678, 9 pages.