# A Semianalytical Solution for Multifractured Horizontal Wells in Box-Shaped Reservoirs

Authors: Lei Wang; Xiaodong Wang; He Zhang; Yunpeng Hu; Chen Li
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper presented a 3D point sink model through using Dirac function. Then, 3D point sink solution in boxed reservoirs was obtained through using Laplace transform and Fourier transform methods. Based on the flux and pressure equivalent conditions in Laplace space, a semianalytical solution for multifractured horizontal wells was also proposed for the first time. The production rate distribution was discussed in detail for multifractured horizontal wells. The calculative results show the outermost fractures had higher production ratio due to larger drainage area and the inner fractures were lower due to the strong interface between fractures. Type curves were established to analyze the flow characteristics, which would be divided into six stages, for example, bilinear flow region, the first linear flow region, the first radial flow region, the second linear flow region, the second radial flow region, and the boundary dominated flow region, respectively. Finally, effects of some sensitive parameters on type curves were also analyzed in detail.


---

## Body

## 1. Introduction

The development of tight reservoirs has been paid more and more attention in China. Stimulation for a horizontal well in tight reservoirs may further enhance well productivity. Unlike a vertical well, a horizontal well may be more than one point along the well length. During the last two decades, horizontal wells have become a common applied completion in the petroleum industry. With a large reservoir contact area, horizontal wells can greatly improve well productivity; however, it is most advantageous to drill horizontal well in thin and tight reservoirs with vertical fractures [1, 2].The main difficulty of transient-pressure analysis for multifractured horizontal wells (MFHWs) is how to solve those problems of multiple-fracture interference, finite conductivity transverse fractures coupled with the formation flow. Some typical research papers are worth reviewing [2–9]. Since 1972, some attempts have been made to simulate the pressure transient behavior for either horizontal or vertical wells, with or without hydraulic fractures. Larsen and Hegre [3] rigorously presented a transient pressure solution of horizontal wells with circular finite-conductivity fractures in the three-dimensional unbounded formation. The results just showed the early-stage and middle-stage features of fracture system because the effects of outer boundary were not considered into their model. Guo et al. [2] presented the pressure transient behavior for a horizontal well with multiple randomly distributed vertical fractures in the infinite reservoir and bounded reservoir. Wan and Aziz [6] described a new semianalytical solution for horizontal wells with multiple hydraulic fractures. The fractures can be rotated at any horizontal angle to the well and need not fully penetrate the formation in the vertical direction. Al-Kobaisi et al. [7] presented a hybrid numerical analytical model with a finite-conductivity vertical fracture intercepted by a horizontal well, which dynamically couples a numerical fracture model with an analytical reservoir model. Their approach allows us to include finer details of the fractures characteristics while keeping the computational work manageable. Valkó and Amini [8] proposed a DVS (distributed volumetric sources) method to predict gas production for a horizontal well with multiple transverse fractures in a bounded reservoir. But it was only an approximate approach and the fracture conductivity was not considered in their model. AI Rbeawi and Tiab [9] introduced a new technique for interpreting the pressure transient behavior for a horizontal well with multiple infinite conductivity fractures which could be longitudinal or transverse, vertical or inclined, or symmetrical or asymmetrical. Recently, some meshless methods have been presented [10–13] and were possibly applied to solve the fracture problems.The main objective of this paper is to develop a computational model to investigate the transient-pressure behavior for multifractured horizontal wells (MFHWs) in box-shaped reservoirs. A semianalytical solution for MFHWs in boxed reservoirs is obtained through using Laplace transform and Fourier transform. The fluid flow in the fracture system and reservoir system is computed separately, and the flux and pressure equivalent conditions in Laplace space are applied in the fracture wall to couple the fluids flow in both systems. The result is validated accurately through comparing with previous results in the literature.


## 2. Physical Model

Figures1 and 2 show the physical model assumed and the assumptions of this model are as follows.(1)
A single-phase flow was assumed in the reservoir. There is a fractured horizontal well in the boxed, closed, and homogeneous formation whose length isx
e, width is y
e, and height is h (see Figure 1).
(2)
The well produces with constant viscosity and slightly compressible fluid at the total flow rate ofQ, but the flow rates from each fracture may change over time.
(3)
The oil production is assumed to be an isothermal process. Darcy flow is assumed in the fracture. The linear flow occurs in the fracture (see Figure2).
(4)
Two wings of each fracture may be of unequal length, and the tips of the fractures are assumed as no-flow boundaries (see Figure2).
(5)
The horizontal wellbore is assumed to be infinitely conductive and the fractures are assumed to be finitely conductive.
(6)
The effect of conductivity can be described by the common vertically fractured wells model.
(7)
At the starting time of production, the pressure is uniformly distributed and is equal to the initial pressure (p
i).Figure 1
Physical model for multiple fractured horizontal wells.Figure 2
Flow model within the fractures.


## 3. Mathematical Model

### 3.1. Reservoir Model

To obtain the fracture solution in the boxed reservoirs, we must firstly obtain the point sink solution. The partial differential equation which is closed in all directions in the Cartesian coordinate system is given by(1)
ϕ
c
t
∂
p
∂
t
-
(
∂
(
(
k
/
μ
)
(
∂
p
/
∂
x
)
)
∂
x
+
∂
(
(
k
/
μ
)
(
∂
p
/
∂
y
)
)
∂
y
+
∂
(
(
k
/
μ
)
(
∂
p
/
∂
z
)
)
∂
z
)
=
q
f
(
γ
,
t
)
δ
(
x
-
x
w
)
δ
(
y
-
y
w
)
δ
(
z
-
z
w
)
,
where ϕ is porosity, c
t is total compressibility, t is the time, p is reservoir pressure, k is reservoir permeability, q
f is point flux, δ is Dirac function, (x
w, y
w, z
w) is the point position, and x
e and y
e are the boundary length and width. All the boundaries of the box-shaped drainage volume are closed; therefore, boundary conditions can be given by
(2)
∂
p
∂
x
|
x
=
0
,
x
e
=
0
∂
p
∂
y
|
y
=
0
,
y
e
=
0
∂
p
∂
z
|
z
=
0
,
z
e
=
0
.
And the initial condition can be given by
(3)
p
(
x
,
y
,
z
,
0
)
=
p
i
.
In order to simplify (1)–(3), we take the following dimensionless transforms:
(4)
p
D
=
2
π
k
h
(
p
i
-
p
)
Q
μ
B
;
t
D
=
k
t
ϕ
μ
c
t
L
2
;
y
D
=
y
L
;
y
w
D
=
y
w
L
;
y
e
D
=
y
e
L
;
x
D
=
x
L
;
x
w
D
=
x
w
L
;
x
e
D
=
x
e
L
;
z
D
=
z
L
;
z
w
D
=
z
w
L
;
z
e
D
=
z
e
L
,
where L is the characteristic length, q
f is the point convergence intensity at the point sink (x
w, y
w, z
w), and the total productivity of the well is Q; there holds
(5)
q
f
D
=
q
f
B
Q
.
Then, through using dimensionless transforms, (1) becomes
(6)
∂
p
D
∂
t
D
-
(
∂
2
p
D
∂
x
D
2
+
∂
2
p
D
∂
y
D
2
+
∂
2
p
D
∂
z
D
2
)
=
2
π
q
f
D
(
γ
D
,
t
D
)
δ
(
x
D
-
x
w
D
)
×
δ
(
y
D
-
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Equation (2) becomes
(7)
∂
p
D
∂
x
D
|
x
D
=
0
,
x
e
D
=
0
∂
p
D
∂
y
D
|
y
D
=
0
,
y
e
D
=
0
∂
p
D
∂
z
D
|
z
D
=
0
,
z
e
D
=
0
.
And initial condition becomes
(8)
p
D
(
x
D
,
y
D
,
z
D
,
0
)
=
0
.
The Laplace transform is based on t
D and functions as follows:
(9)
p
~
D
(
r
D
,
s
)
=
∫
0
∞
p
D
(
r
D
,
s
)
e
-
s
t
D
d
t
D
.
Applying the Laplace transforms to (6)–(8), we have
(10)
s
p
~
D
-
(
∂
2
p
~
D
∂
x
D
2
+
∂
2
p
~
D
∂
y
D
2
+
∂
2
p
~
D
∂
z
D
2
)
=
2
π
q
~
f
D
(
s
)
δ
(
x
D
-
x
w
D
)
×
δ
(
y
D
-
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Outer boundary conditions in Laplace space are
(11)
∂
p
~
D
∂
x
D
|
x
D
=
0
,
x
e
D
=
0
∂
p
~
D
∂
y
D
|
y
D
=
0
,
y
e
D
=
0
∂
p
~
D
∂
z
D
|
z
D
=
0
,
z
e
D
=
0
.
Fourier cosine transform based on x
D can be defined as
(12)
p
-
D
(
u
n
)
=
∫
0
x
e
D
p
D
cos
⁡
(
u
n
x
D
)
d
x
D
.
The characteristic equation is
(13)
sin
(
u
m
x
e
D
)
=
0
.
Through solving (13), characteristic number is obtained by
(14)
u
m
=
m
π
x
e
D
.
Imposing the Fourier cosine transform on the variable x
D in (10), we have
(15)
s
p
~
-
D
-
(
-
u
m
2
p
~
-
D
+
∂
2
p
~
-
D
∂
y
D
2
+
∂
2
p
~
-
D
∂
z
D
2
)
=
2
π
q
~
f
D
(
s
)
cos
⁡
(
u
m
x
w
D
)
×
δ
(
y
D
-
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Similarly, imposing the Fourier cosine transform on the variable y
D in (15), we have
(16)
s
p
~
-
-
D
-
(
-
u
m
2
p
~
-
¯
D
-
v
n
2
p
~
-
-
D
+
∂
2
p
~
-
-
D
∂
z
D
2
)
=
2
π
q
~
f
D
(
s
)
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Similarly, imposing the Fourier cosine transform on the variable z
D in (16), we have(17a)
p
~
-
-
-
D
=
2
π
q
~
f
D
(
s
)
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
,
where
(17b)
u
m
=
m
π
x
e
D
v
n
=
n
π
y
e
D
w
p
=
p
π
z
e
D
.Imposing the Fourier inverse transform on the variablez
D in  (17a) and (17b), we obtain
(18)
p
~
-
-
D
q
~
f
D
(
s
)
=
1
z
e
D
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
+
2
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
.
Imposing the Fourier inverse transform on the variable y
D in (18), we have
(19)
p
~
-
D
q
~
f
D
(
s
)
=
1
y
e
D
1
z
e
D
2
π
cos
⁡
(
u
m
x
w
D
)
s
+
u
m
2
+
2
y
e
D
∑
n
=
1
∞
1
z
e
D
cos
⁡
(
v
n
y
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
+
1
y
e
D
2
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
w
p
2
+
2
y
e
D
2
z
e
D
∑
n
=
1
∞
∑
p
=
1
∞
cos
⁡
(
v
n
y
D
)
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
.
Imposing the Fourier inverse transform on the variable x
D in (19), we have
(20)
p
~
D
q
~
f
D
(
s
)
=
2
π
x
e
D
y
e
D
z
e
D
s
+
2
x
e
D
y
e
D
z
e
D
×
∑
m
=
1
∞
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
s
+
u
m
2
+
2
x
e
D
y
e
D
1
z
e
D
∑
n
=
1
∞
cos
⁡
⁡
(
v
n
y
D
)
×
2
π
cos
⁡
(
v
n
y
w
D
)
s
+
v
n
2
+
4
x
e
D
y
e
D
z
e
D
×
∑
m
=
1
∞
∑
n
=
1
∞
cos
⁡
(
v
n
y
D
)
cos
⁡
(
u
m
x
D
)
×
2
π
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
u
m
x
w
D
)
s
+
u
m
2
+
v
n
2
+
1
x
e
D
y
e
D
2
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
2
π
cos
⁡
(
w
p
z
w
D
)
s
+
w
p
2
+
2
x
e
D
y
e
D
2
z
e
D
×
∑
m
=
1
∞
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
u
m
x
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
w
p
2
+
4
x
e
D
y
e
D
z
e
D
∑
n
=
1
∞
∑
p
=
1
∞
cos
⁡
(
v
n
y
D
)
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
v
n
2
+
w
p
2
+
8
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
∑
n
=
1
∞
∑
p
=
1
∞
cos
⁡
(
u
m
x
D
)
×
cos
⁡
(
v
n
y
D
)
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
,
where we use the following formulas:
(21)
∑
k
=
1
∞
cos
⁡
k
x
k
2
+
α
2
=
π
2
α
cosh
⁡
α
(
π
-
x
)
sinh
⁡
α
π
-
1
2
α
2
.
Therefore, we have
(22)
∑
n
=
1
∞
cos
⁡
(
(
n
π
/
y
e
D
)
(
y
D
-
y
w
D
)
)
s
+
(
n
π
/
y
e
D
)
2
=
y
e
D
2
s
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
sinh
⁡
s
y
e
D
-
1
2
s
and the cosine product formula can be given by
(23)
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
=
cos
⁡
(
v
n
(
y
D
-
y
w
D
)
)
+
cos
⁡
(
v
n
(
y
D
+
y
w
D
)
)
=
cos
⁡
(
n
π
y
e
D
(
y
D
-
y
w
D
)
)
+
cos
⁡
(
n
π
y
e
D
(
y
D
+
y
w
D
)
)
.
Through using (22)-(23), the following formulas can be obtained:
(24)
∑
n
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
v
n
2
=
∑
n
=
1
∞
(
cos
⁡
(
(
n
π
/
y
e
D
)
(
y
D
-
y
w
D
)
)
+
cos
⁡
(
(
n
π
/
y
e
D
)
(
y
D
+
y
w
D
)
)
)
×
(
s
+
v
n
2
)
-
1
=
y
e
D
s
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
sinh
⁡
s
y
e
D
-
1
2
s
+
y
e
D
s
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
sinh
⁡
s
y
e
D
-
1
2
s
=
y
e
D
s
×
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
y
e
D
)
-
1
-
1
s
,
∑
n
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
=
y
e
D
2
s
+
u
m
2
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
-
1
s
+
u
m
2
,
∑
p
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
v
n
2
+
w
p
2
=
y
e
D
2
s
+
w
p
2
×
(
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
w
p
2
y
e
D
)
-
1
-
1
s
+
w
p
2
,
∑
p
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
=
y
e
D
2
s
+
u
m
2
+
w
p
2
×
(
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
+
w
p
2
y
e
D
)
-
1
-
1
s
+
u
m
2
+
w
p
2
.
Substituting (21)–(24) into (20), we can have
(25)
p
~
D
q
~
f
D
(
s
)
=
2
π
x
e
D
y
e
D
z
e
D
s
+
2
π
x
e
D
y
e
D
z
e
D
×
∑
m
=
1
∞
2
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
s
+
u
m
2
+
2
π
x
e
D
y
e
D
z
e
D
×
[
y
e
D
2
s
×
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
y
e
D
)
-
1
-
1
s
]
+
4
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
[
y
e
D
2
s
+
u
m
2
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
-
1
s
+
u
m
2
y
e
D
2
s
+
u
m
2
]
+
2
x
e
D
y
e
D
z
e
D
×
∑
p
=
1
∞
2
π
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
w
p
2
+
4
x
e
D
y
e
D
z
e
D
×
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
u
m
x
D
)
×
∑
m
=
1
∞
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
w
p
2
+
4
π
x
e
D
y
e
D
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
×
[
y
e
D
2
s
+
w
p
2
h
×
(
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
hh
+
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
h
×
(
sinh
⁡
s
+
w
p
2
y
e
D
)
-
1
-
1
s
+
w
p
2
y
e
D
2
s
+
u
m
2
+
w
p
2
]
+
8
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
×
[
y
e
D
2
s
+
u
m
2
+
w
p
2
×
(
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
+
w
p
2
y
e
D
)
-
1
-
1
s
+
u
m
2
+
w
p
2
y
e
D
2
s
+
u
m
2
+
w
p
2
]
.
Further, (25) could be simplified as the following equation:
(26)
p
~
D
q
~
f
D
(
s
)
=
2
π
x
e
D
y
e
D
z
e
D
×
[
y
e
D
2
s
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
h
hh
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
y
e
D
)
-
1
y
e
D
2
s
]
+
4
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
y
e
D
2
s
+
u
m
2
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
+
4
π
x
e
D
y
e
D
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
y
e
D
2
s
+
w
p
2
×
(
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
w
p
2
y
e
D
)
-
1
+
8
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
y
e
D
2
s
+
u
m
2
+
w
p
2
×
(
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
+
w
p
2
y
e
D
)
-
1
.
Equation (26) is the point sink solution in box-shaped reservoirs.


### 3.2. Fracture Model

Fluid flow insidenth fracture in the system may be defined by the following set of equations:
(27)
∂
(
(
k
f
n
/
μ
)
(
∂
p
f
n
/
∂
z
n
)
)
∂
z
n
+
q
w
h
f
n
w
f
n
h
δ
(
z
n
-
z
w
h
f
n
)
-
Ω
f
n
(
z
n
,
t
)
=
0
,
where
(28)
Ω
f
n
(
z
n
,
t
)
=
∫
0
L
f
n
q
z
f
n
(
z
n
,
t
)
w
f
n
h
L
f
n
δ
(
z
n
-
γ
n
)
d
γ
n
,
where z
w
h
f
n is the well position in the fracture, k
f
n is the fracture permeability, q
w
h
f
n is well rate in the fracture, L
f
n is the fracture length, γ
n is the source position, h is the reservoir thick, and q
z
f
n is the fracture flux. In order to simplify (27), we take the following dimensionless transforms:
(29)
P
f
D
=
2
π
k
h
(
p
i
-
p
f
n
)
B
Q
μ
,
C
f
D
n
=
k
f
n
w
f
n
k
L
f
n
,
q
w
h
f
D
n
=
q
w
h
f
n
Q
,
L
f
D
n
=
L
f
n
L
,
γ
D
n
=
γ
n
L
,
z
D
n
=
z
n
L
,
z
w
h
f
D
n
=
z
w
h
f
n
L
,
q
z
f
D
n
=
q
z
f
n
Q
.
Therefore, the following equations can be obtained by
(30)
∂
2
p
f
D
n
∂
z
D
n
2
+
2
π
q
w
h
f
D
n
L
f
D
n
C
f
D
n
δ
(
z
D
n
-
z
w
h
f
D
n
)
-
∫
0
L
f
D
n
2
π
q
z
f
D
n
(
z
D
n
,
t
D
)
L
f
D
n
2
C
f
D
n
δ
(
z
D
n
-
γ
D
n
)
d
γ
D
n
=
0
.
The initial condition can be given by
(31)
p
f
D
n
(
z
D
n
,
0
)
=
0
.
The boundary conditions are given by
(32)
∂
p
f
D
n
(
0
,
t
D
)
∂
z
D
n
=
0
,
∂
p
f
D
n
(
L
f
D
n
,
t
D
)
∂
z
D
n
=
0
.
Imposing Laplace transforms on t
D in (30)–(32), we have
(33)
∂
2
p
~
f
D
n
∂
z
D
n
2
+
2
π
q
~
w
h
f
D
n
L
f
D
n
C
f
D
n
δ
(
z
D
n
-
z
w
h
f
D
n
)
-
∫
0
L
f
D
n
2
π
q
~
z
f
D
n
(
z
D
n
,
s
)
L
f
D
n
2
C
f
D
n
δ
(
z
D
n
-
γ
D
n
)
d
γ
D
n
=
0
.
The boundary conditions become
(34)
∂
p
~
f
D
n
(
0
,
s
)
∂
z
D
n
=
0
,
∂
p
~
f
D
n
(
L
f
D
n
,
s
)
∂
z
D
n
=
0
.


### 3.3. Semianalytical Solution

The reservoir- and fracture-flow equations are described in Sections3.1 and 3.2. To build the reservoir-flow model, Fourier transformation on space and Laplace transformation on time have been applied and point sink solution has been also obtained in Section 3.1. Thus, line sink solution can be derived by integrating the solution given in (26) with respect to z
w
D over 0 to z
e
D. (Note that only the integral of cos (w
p
z
w
D) needs to be determined.) The result can be given by
(35)
p
~
D
=
q
~
f
D
(
γ
D
,
s
)
F
(
x
D
,
y
D
,
x
w
D
,
y
w
D
,
x
e
D
,
y
e
D
,
s
)
,
where
(36)
F
(
x
D
,
y
D
,
x
w
D
,
y
w
D
,
x
e
D
,
y
e
D
,
s
)
=
π
x
e
D
×
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
h
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
s
sinh
⁡
s
y
e
D
)
-
1
+
2
π
x
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
s
+
u
m
2
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
.Fracture number for MFHW isN
f; the pressure drop formulas at each fracture segment for MFHW can be obtained through using superposition principle and point sink integral method,
(37)
p
~
D
=
∑
n
=
1
N
f
∫
0
L
f
D
n
q
~
f
D
n
(
γ
D
,
s
)
×
F
(
x
D
,
y
D
,
x
f
D
n
+
γ
D
n
cos
⁡
θ
f
n
,
y
f
D
n
+
γ
D
n
sin
θ
f
n
,
x
e
D
,
y
e
D
,
s
)
d
γ
D
,
whereθ
f
n is the angle between nth fracture and x-axis, F is defined in (36), and γ
D
n is source position in the fracture. Now, we divide per fracture into N
s
n segments; thus, pressure drop formulas for jth segment of nth fracture can be given by
(38)
p
~
D
j
=
∑
n
=
1
N
f
∑
i
=
1
N
s
n
q
~
f
D
n
i
(
γ
D
,
s
)
×
F
(
x
D
j
,
y
D
j
,
x
f
D
n
+
γ
D
n
cos
⁡
θ
f
n
,
y
f
D
n
+
γ
D
n
sin
θ
f
n
,
x
e
D
,
y
e
D
,
s
)
d
γ
D
.
To obtain the solution for the fracture-flow models, Laplace transformations on both space and time have been carried out. The pressure-drop evaluation for a fracture is presented in the following equation. The pressure drop at kth segment on the nth fracture is given by
(39)
p
~
f
D
n
k
(
z
D
n
k
,
s
)
=
p
~
f
D
n
(
0
,
s
)
+
1
L
f
D
n
2
C
f
D
n
R
1
-
2
π
L
f
D
n
C
f
D
n
R
2
,
where
(40)
R
1
=
∑
i
=
1
k
-
1
q
~
f
D
n
i
(
1
2
+
k
-
i
)
Δ
z
D
n
2
+
q
~
f
D
n
k
Δ
z
D
n
2
8
,
R
2
=
q
~
w
h
f
D
n
(
z
w
h
f
D
n
,
s
)
[
z
f
D
n
k
-
z
w
h
f
D
n
]
×
H
(
z
f
D
n
k
-
z
w
h
f
D
n
)
,
where H is the Heaviside unit step function.For each fracture, the total flux equation can be given by(41)
∑
i
=
1
N
s
n
q
~
f
D
n
i
Δ
x
D
n
=
q
~
w
h
f
D
n
and the pressure for each fracture is equal to the pressure of horizontal well; it must have
(42)
p
~
f
D
1
(
z
w
h
f
D
1
,
s
)
=
p
~
f
D
2
(
z
w
h
f
D
2
,
s
)
⋯
=
p
~
f
D
N
f
-
1
(
z
w
h
f
D
N
f
-
1
,
s
)
=
p
~
f
D
N
f
(
z
w
h
f
D
N
f
,
s
)
.
In addition, total rate of horizontal wells should be equal to sum of flux for each fracture; thus,
(43)
∑
i
=
1
N
f
q
~
w
h
f
D
n
=
1
s
.
Coupling the fracture- and reservoir-flow models at the fracture faces, we obtain a system of N
t
s
+
N
f
+
N
f equations with N
t
s
+
N
f
+
N
f unknowns. N
t
s is the total number of fracture segments, and N
f is the total number of fractures. This system of equations could be rearranged and written in the matrix format. Once the matrix is solved for unknowns, the rate at each fracture segment and the pressure drop at each fracture will be obtained. These rates and pressure drops are used to calculate the pressure drop at the wellbore for MFHW.


## 3.1. Reservoir Model

To obtain the fracture solution in the boxed reservoirs, we must firstly obtain the point sink solution. The partial differential equation which is closed in all directions in the Cartesian coordinate system is given by(1)
ϕ
c
t
∂
p
∂
t
-
(
∂
(
(
k
/
μ
)
(
∂
p
/
∂
x
)
)
∂
x
+
∂
(
(
k
/
μ
)
(
∂
p
/
∂
y
)
)
∂
y
+
∂
(
(
k
/
μ
)
(
∂
p
/
∂
z
)
)
∂
z
)
=
q
f
(
γ
,
t
)
δ
(
x
-
x
w
)
δ
(
y
-
y
w
)
δ
(
z
-
z
w
)
,
where ϕ is porosity, c
t is total compressibility, t is the time, p is reservoir pressure, k is reservoir permeability, q
f is point flux, δ is Dirac function, (x
w, y
w, z
w) is the point position, and x
e and y
e are the boundary length and width. All the boundaries of the box-shaped drainage volume are closed; therefore, boundary conditions can be given by
(2)
∂
p
∂
x
|
x
=
0
,
x
e
=
0
∂
p
∂
y
|
y
=
0
,
y
e
=
0
∂
p
∂
z
|
z
=
0
,
z
e
=
0
.
And the initial condition can be given by
(3)
p
(
x
,
y
,
z
,
0
)
=
p
i
.
In order to simplify (1)–(3), we take the following dimensionless transforms:
(4)
p
D
=
2
π
k
h
(
p
i
-
p
)
Q
μ
B
;
t
D
=
k
t
ϕ
μ
c
t
L
2
;
y
D
=
y
L
;
y
w
D
=
y
w
L
;
y
e
D
=
y
e
L
;
x
D
=
x
L
;
x
w
D
=
x
w
L
;
x
e
D
=
x
e
L
;
z
D
=
z
L
;
z
w
D
=
z
w
L
;
z
e
D
=
z
e
L
,
where L is the characteristic length, q
f is the point convergence intensity at the point sink (x
w, y
w, z
w), and the total productivity of the well is Q; there holds
(5)
q
f
D
=
q
f
B
Q
.
Then, through using dimensionless transforms, (1) becomes
(6)
∂
p
D
∂
t
D
-
(
∂
2
p
D
∂
x
D
2
+
∂
2
p
D
∂
y
D
2
+
∂
2
p
D
∂
z
D
2
)
=
2
π
q
f
D
(
γ
D
,
t
D
)
δ
(
x
D
-
x
w
D
)
×
δ
(
y
D
-
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Equation (2) becomes
(7)
∂
p
D
∂
x
D
|
x
D
=
0
,
x
e
D
=
0
∂
p
D
∂
y
D
|
y
D
=
0
,
y
e
D
=
0
∂
p
D
∂
z
D
|
z
D
=
0
,
z
e
D
=
0
.
And initial condition becomes
(8)
p
D
(
x
D
,
y
D
,
z
D
,
0
)
=
0
.
The Laplace transform is based on t
D and functions as follows:
(9)
p
~
D
(
r
D
,
s
)
=
∫
0
∞
p
D
(
r
D
,
s
)
e
-
s
t
D
d
t
D
.
Applying the Laplace transforms to (6)–(8), we have
(10)
s
p
~
D
-
(
∂
2
p
~
D
∂
x
D
2
+
∂
2
p
~
D
∂
y
D
2
+
∂
2
p
~
D
∂
z
D
2
)
=
2
π
q
~
f
D
(
s
)
δ
(
x
D
-
x
w
D
)
×
δ
(
y
D
-
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Outer boundary conditions in Laplace space are
(11)
∂
p
~
D
∂
x
D
|
x
D
=
0
,
x
e
D
=
0
∂
p
~
D
∂
y
D
|
y
D
=
0
,
y
e
D
=
0
∂
p
~
D
∂
z
D
|
z
D
=
0
,
z
e
D
=
0
.
Fourier cosine transform based on x
D can be defined as
(12)
p
-
D
(
u
n
)
=
∫
0
x
e
D
p
D
cos
⁡
(
u
n
x
D
)
d
x
D
.
The characteristic equation is
(13)
sin
(
u
m
x
e
D
)
=
0
.
Through solving (13), characteristic number is obtained by
(14)
u
m
=
m
π
x
e
D
.
Imposing the Fourier cosine transform on the variable x
D in (10), we have
(15)
s
p
~
-
D
-
(
-
u
m
2
p
~
-
D
+
∂
2
p
~
-
D
∂
y
D
2
+
∂
2
p
~
-
D
∂
z
D
2
)
=
2
π
q
~
f
D
(
s
)
cos
⁡
(
u
m
x
w
D
)
×
δ
(
y
D
-
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Similarly, imposing the Fourier cosine transform on the variable y
D in (15), we have
(16)
s
p
~
-
-
D
-
(
-
u
m
2
p
~
-
¯
D
-
v
n
2
p
~
-
-
D
+
∂
2
p
~
-
-
D
∂
z
D
2
)
=
2
π
q
~
f
D
(
s
)
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
δ
(
z
D
-
z
w
D
)
.
Similarly, imposing the Fourier cosine transform on the variable z
D in (16), we have(17a)
p
~
-
-
-
D
=
2
π
q
~
f
D
(
s
)
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
,
where
(17b)
u
m
=
m
π
x
e
D
v
n
=
n
π
y
e
D
w
p
=
p
π
z
e
D
.Imposing the Fourier inverse transform on the variablez
D in  (17a) and (17b), we obtain
(18)
p
~
-
-
D
q
~
f
D
(
s
)
=
1
z
e
D
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
+
2
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
.
Imposing the Fourier inverse transform on the variable y
D in (18), we have
(19)
p
~
-
D
q
~
f
D
(
s
)
=
1
y
e
D
1
z
e
D
2
π
cos
⁡
(
u
m
x
w
D
)
s
+
u
m
2
+
2
y
e
D
∑
n
=
1
∞
1
z
e
D
cos
⁡
(
v
n
y
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
+
1
y
e
D
2
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
w
p
2
+
2
y
e
D
2
z
e
D
∑
n
=
1
∞
∑
p
=
1
∞
cos
⁡
(
v
n
y
D
)
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
.
Imposing the Fourier inverse transform on the variable x
D in (19), we have
(20)
p
~
D
q
~
f
D
(
s
)
=
2
π
x
e
D
y
e
D
z
e
D
s
+
2
x
e
D
y
e
D
z
e
D
×
∑
m
=
1
∞
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
s
+
u
m
2
+
2
x
e
D
y
e
D
1
z
e
D
∑
n
=
1
∞
cos
⁡
⁡
(
v
n
y
D
)
×
2
π
cos
⁡
(
v
n
y
w
D
)
s
+
v
n
2
+
4
x
e
D
y
e
D
z
e
D
×
∑
m
=
1
∞
∑
n
=
1
∞
cos
⁡
(
v
n
y
D
)
cos
⁡
(
u
m
x
D
)
×
2
π
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
u
m
x
w
D
)
s
+
u
m
2
+
v
n
2
+
1
x
e
D
y
e
D
2
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
2
π
cos
⁡
(
w
p
z
w
D
)
s
+
w
p
2
+
2
x
e
D
y
e
D
2
z
e
D
×
∑
m
=
1
∞
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
u
m
x
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
w
p
2
+
4
x
e
D
y
e
D
z
e
D
∑
n
=
1
∞
∑
p
=
1
∞
cos
⁡
(
v
n
y
D
)
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
v
n
2
+
w
p
2
+
8
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
∑
n
=
1
∞
∑
p
=
1
∞
cos
⁡
(
u
m
x
D
)
×
cos
⁡
(
v
n
y
D
)
cos
⁡
(
w
p
z
D
)
×
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
v
n
y
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
,
where we use the following formulas:
(21)
∑
k
=
1
∞
cos
⁡
k
x
k
2
+
α
2
=
π
2
α
cosh
⁡
α
(
π
-
x
)
sinh
⁡
α
π
-
1
2
α
2
.
Therefore, we have
(22)
∑
n
=
1
∞
cos
⁡
(
(
n
π
/
y
e
D
)
(
y
D
-
y
w
D
)
)
s
+
(
n
π
/
y
e
D
)
2
=
y
e
D
2
s
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
sinh
⁡
s
y
e
D
-
1
2
s
and the cosine product formula can be given by
(23)
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
=
cos
⁡
(
v
n
(
y
D
-
y
w
D
)
)
+
cos
⁡
(
v
n
(
y
D
+
y
w
D
)
)
=
cos
⁡
(
n
π
y
e
D
(
y
D
-
y
w
D
)
)
+
cos
⁡
(
n
π
y
e
D
(
y
D
+
y
w
D
)
)
.
Through using (22)-(23), the following formulas can be obtained:
(24)
∑
n
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
v
n
2
=
∑
n
=
1
∞
(
cos
⁡
(
(
n
π
/
y
e
D
)
(
y
D
-
y
w
D
)
)
+
cos
⁡
(
(
n
π
/
y
e
D
)
(
y
D
+
y
w
D
)
)
)
×
(
s
+
v
n
2
)
-
1
=
y
e
D
s
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
sinh
⁡
s
y
e
D
-
1
2
s
+
y
e
D
s
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
sinh
⁡
s
y
e
D
-
1
2
s
=
y
e
D
s
×
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
y
e
D
)
-
1
-
1
s
,
∑
n
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
=
y
e
D
2
s
+
u
m
2
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
-
1
s
+
u
m
2
,
∑
p
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
v
n
2
+
w
p
2
=
y
e
D
2
s
+
w
p
2
×
(
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
w
p
2
y
e
D
)
-
1
-
1
s
+
w
p
2
,
∑
p
=
1
∞
2
cos
⁡
(
v
n
y
D
)
cos
⁡
(
v
n
y
w
D
)
s
+
u
m
2
+
v
n
2
+
w
p
2
=
y
e
D
2
s
+
u
m
2
+
w
p
2
×
(
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
+
w
p
2
y
e
D
)
-
1
-
1
s
+
u
m
2
+
w
p
2
.
Substituting (21)–(24) into (20), we can have
(25)
p
~
D
q
~
f
D
(
s
)
=
2
π
x
e
D
y
e
D
z
e
D
s
+
2
π
x
e
D
y
e
D
z
e
D
×
∑
m
=
1
∞
2
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
s
+
u
m
2
+
2
π
x
e
D
y
e
D
z
e
D
×
[
y
e
D
2
s
×
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
y
e
D
)
-
1
-
1
s
]
+
4
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
[
y
e
D
2
s
+
u
m
2
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
-
1
s
+
u
m
2
y
e
D
2
s
+
u
m
2
]
+
2
x
e
D
y
e
D
z
e
D
×
∑
p
=
1
∞
2
π
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
w
p
2
+
4
x
e
D
y
e
D
z
e
D
×
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
u
m
x
D
)
×
∑
m
=
1
∞
2
π
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
w
p
z
w
D
)
s
+
u
m
2
+
w
p
2
+
4
π
x
e
D
y
e
D
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
×
[
y
e
D
2
s
+
w
p
2
h
×
(
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
hh
+
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
h
×
(
sinh
⁡
s
+
w
p
2
y
e
D
)
-
1
-
1
s
+
w
p
2
y
e
D
2
s
+
u
m
2
+
w
p
2
]
+
8
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
×
[
y
e
D
2
s
+
u
m
2
+
w
p
2
×
(
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
+
w
p
2
y
e
D
)
-
1
-
1
s
+
u
m
2
+
w
p
2
y
e
D
2
s
+
u
m
2
+
w
p
2
]
.
Further, (25) could be simplified as the following equation:
(26)
p
~
D
q
~
f
D
(
s
)
=
2
π
x
e
D
y
e
D
z
e
D
×
[
y
e
D
2
s
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
h
hh
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
y
e
D
)
-
1
y
e
D
2
s
]
+
4
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
y
e
D
2
s
+
u
m
2
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
+
4
π
x
e
D
y
e
D
z
e
D
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
y
e
D
2
s
+
w
p
2
×
(
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
w
p
2
y
e
D
)
-
1
+
8
π
x
e
D
y
e
D
z
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
∑
p
=
1
∞
cos
⁡
(
w
p
z
D
)
cos
⁡
(
w
p
z
w
D
)
y
e
D
2
s
+
u
m
2
+
w
p
2
×
(
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
+
w
p
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
sinh
⁡
s
+
u
m
2
+
w
p
2
y
e
D
)
-
1
.
Equation (26) is the point sink solution in box-shaped reservoirs.


## 3.2. Fracture Model

Fluid flow insidenth fracture in the system may be defined by the following set of equations:
(27)
∂
(
(
k
f
n
/
μ
)
(
∂
p
f
n
/
∂
z
n
)
)
∂
z
n
+
q
w
h
f
n
w
f
n
h
δ
(
z
n
-
z
w
h
f
n
)
-
Ω
f
n
(
z
n
,
t
)
=
0
,
where
(28)
Ω
f
n
(
z
n
,
t
)
=
∫
0
L
f
n
q
z
f
n
(
z
n
,
t
)
w
f
n
h
L
f
n
δ
(
z
n
-
γ
n
)
d
γ
n
,
where z
w
h
f
n is the well position in the fracture, k
f
n is the fracture permeability, q
w
h
f
n is well rate in the fracture, L
f
n is the fracture length, γ
n is the source position, h is the reservoir thick, and q
z
f
n is the fracture flux. In order to simplify (27), we take the following dimensionless transforms:
(29)
P
f
D
=
2
π
k
h
(
p
i
-
p
f
n
)
B
Q
μ
,
C
f
D
n
=
k
f
n
w
f
n
k
L
f
n
,
q
w
h
f
D
n
=
q
w
h
f
n
Q
,
L
f
D
n
=
L
f
n
L
,
γ
D
n
=
γ
n
L
,
z
D
n
=
z
n
L
,
z
w
h
f
D
n
=
z
w
h
f
n
L
,
q
z
f
D
n
=
q
z
f
n
Q
.
Therefore, the following equations can be obtained by
(30)
∂
2
p
f
D
n
∂
z
D
n
2
+
2
π
q
w
h
f
D
n
L
f
D
n
C
f
D
n
δ
(
z
D
n
-
z
w
h
f
D
n
)
-
∫
0
L
f
D
n
2
π
q
z
f
D
n
(
z
D
n
,
t
D
)
L
f
D
n
2
C
f
D
n
δ
(
z
D
n
-
γ
D
n
)
d
γ
D
n
=
0
.
The initial condition can be given by
(31)
p
f
D
n
(
z
D
n
,
0
)
=
0
.
The boundary conditions are given by
(32)
∂
p
f
D
n
(
0
,
t
D
)
∂
z
D
n
=
0
,
∂
p
f
D
n
(
L
f
D
n
,
t
D
)
∂
z
D
n
=
0
.
Imposing Laplace transforms on t
D in (30)–(32), we have
(33)
∂
2
p
~
f
D
n
∂
z
D
n
2
+
2
π
q
~
w
h
f
D
n
L
f
D
n
C
f
D
n
δ
(
z
D
n
-
z
w
h
f
D
n
)
-
∫
0
L
f
D
n
2
π
q
~
z
f
D
n
(
z
D
n
,
s
)
L
f
D
n
2
C
f
D
n
δ
(
z
D
n
-
γ
D
n
)
d
γ
D
n
=
0
.
The boundary conditions become
(34)
∂
p
~
f
D
n
(
0
,
s
)
∂
z
D
n
=
0
,
∂
p
~
f
D
n
(
L
f
D
n
,
s
)
∂
z
D
n
=
0
.


## 3.3. Semianalytical Solution

The reservoir- and fracture-flow equations are described in Sections3.1 and 3.2. To build the reservoir-flow model, Fourier transformation on space and Laplace transformation on time have been applied and point sink solution has been also obtained in Section 3.1. Thus, line sink solution can be derived by integrating the solution given in (26) with respect to z
w
D over 0 to z
e
D. (Note that only the integral of cos (w
p
z
w
D) needs to be determined.) The result can be given by
(35)
p
~
D
=
q
~
f
D
(
γ
D
,
s
)
F
(
x
D
,
y
D
,
x
w
D
,
y
w
D
,
x
e
D
,
y
e
D
,
s
)
,
where
(36)
F
(
x
D
,
y
D
,
x
w
D
,
y
w
D
,
x
e
D
,
y
e
D
,
s
)
=
π
x
e
D
×
(
cosh
⁡
s
(
y
e
D
-
(
y
D
-
y
w
D
)
)
h
+
cosh
⁡
s
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
s
sinh
⁡
s
y
e
D
)
-
1
+
2
π
x
e
D
∑
m
=
1
∞
cos
⁡
(
u
m
x
w
D
)
cos
⁡
(
u
m
x
D
)
×
(
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
-
y
w
D
)
)
+
cosh
⁡
s
+
u
m
2
(
y
e
D
-
(
y
D
+
y
w
D
)
)
)
×
(
s
+
u
m
2
sinh
⁡
s
+
u
m
2
y
e
D
)
-
1
.Fracture number for MFHW isN
f; the pressure drop formulas at each fracture segment for MFHW can be obtained through using superposition principle and point sink integral method,
(37)
p
~
D
=
∑
n
=
1
N
f
∫
0
L
f
D
n
q
~
f
D
n
(
γ
D
,
s
)
×
F
(
x
D
,
y
D
,
x
f
D
n
+
γ
D
n
cos
⁡
θ
f
n
,
y
f
D
n
+
γ
D
n
sin
θ
f
n
,
x
e
D
,
y
e
D
,
s
)
d
γ
D
,
whereθ
f
n is the angle between nth fracture and x-axis, F is defined in (36), and γ
D
n is source position in the fracture. Now, we divide per fracture into N
s
n segments; thus, pressure drop formulas for jth segment of nth fracture can be given by
(38)
p
~
D
j
=
∑
n
=
1
N
f
∑
i
=
1
N
s
n
q
~
f
D
n
i
(
γ
D
,
s
)
×
F
(
x
D
j
,
y
D
j
,
x
f
D
n
+
γ
D
n
cos
⁡
θ
f
n
,
y
f
D
n
+
γ
D
n
sin
θ
f
n
,
x
e
D
,
y
e
D
,
s
)
d
γ
D
.
To obtain the solution for the fracture-flow models, Laplace transformations on both space and time have been carried out. The pressure-drop evaluation for a fracture is presented in the following equation. The pressure drop at kth segment on the nth fracture is given by
(39)
p
~
f
D
n
k
(
z
D
n
k
,
s
)
=
p
~
f
D
n
(
0
,
s
)
+
1
L
f
D
n
2
C
f
D
n
R
1
-
2
π
L
f
D
n
C
f
D
n
R
2
,
where
(40)
R
1
=
∑
i
=
1
k
-
1
q
~
f
D
n
i
(
1
2
+
k
-
i
)
Δ
z
D
n
2
+
q
~
f
D
n
k
Δ
z
D
n
2
8
,
R
2
=
q
~
w
h
f
D
n
(
z
w
h
f
D
n
,
s
)
[
z
f
D
n
k
-
z
w
h
f
D
n
]
×
H
(
z
f
D
n
k
-
z
w
h
f
D
n
)
,
where H is the Heaviside unit step function.For each fracture, the total flux equation can be given by(41)
∑
i
=
1
N
s
n
q
~
f
D
n
i
Δ
x
D
n
=
q
~
w
h
f
D
n
and the pressure for each fracture is equal to the pressure of horizontal well; it must have
(42)
p
~
f
D
1
(
z
w
h
f
D
1
,
s
)
=
p
~
f
D
2
(
z
w
h
f
D
2
,
s
)
⋯
=
p
~
f
D
N
f
-
1
(
z
w
h
f
D
N
f
-
1
,
s
)
=
p
~
f
D
N
f
(
z
w
h
f
D
N
f
,
s
)
.
In addition, total rate of horizontal wells should be equal to sum of flux for each fracture; thus,
(43)
∑
i
=
1
N
f
q
~
w
h
f
D
n
=
1
s
.
Coupling the fracture- and reservoir-flow models at the fracture faces, we obtain a system of N
t
s
+
N
f
+
N
f equations with N
t
s
+
N
f
+
N
f unknowns. N
t
s is the total number of fracture segments, and N
f is the total number of fractures. This system of equations could be rearranged and written in the matrix format. Once the matrix is solved for unknowns, the rate at each fracture segment and the pressure drop at each fracture will be obtained. These rates and pressure drops are used to calculate the pressure drop at the wellbore for MFHW.


## 4. Results and Discussion

### 4.1. Validation of Solutions in This Paper

Riley et al. [14] presented an analytical solution for a single fracture case. To validate the solutions proposed in this paper, some data were utilized (see Table 1) for comparing the solutions with those from Riley et al. [14] (see Figure 3). For the convenience of the comparison, only a fracture is considered into the MFHWs. As shown in Figure 3, there is a good agreement between the solutions obtained in this work and the results presented by Riley. This might indicate that previous models are special cases of our MFHWs model and also implies that our model is validated accurately.Table 1
Basic data.
Parameters
Dimensionless value
Fracture numberN
f
1
Fracture segment numberN
sn
16
Fracture positionx
f
D
n
2000
Fracture positiony
f
D
n
2000
Fracture lengthL
f
D
n
1
Well positionZ
w
h
f
D
n
0.5
Region lengthx
e
D
5000
Region widthy
e
D
5000
Fracture angleθ
f
n
π/2Figure 3
Validation of solutions presented in this paper.


### 4.2. The Production Rate Distribution

The production rate distribution is presented in Figure4 and the data is shown in Table 2. In this case, five transversal fractures have identical length and uniform distribution along a horizontal well located in a box-shaped reservoir. As shown in Figure 4, the outermost fractures have higher production ratio due to larger drainage area and the inner fractures are lower due to strong interference between fractures. For a specific fracture, the flux distribution along the discrete elements has the same property as the whole fracture system.Table 2
Basic data.
Parameters
Frac1
Frac2
Frac3
Frac4
Frac5
Fracture numberN
f
5
5
5
5
5
Fracture segment numberN
sn
20
20
20
20
20
Fracture positionx
f
D
n
2000
2000
2000
2000
2000
Fracture positiony
f
D
n
1900
1950
2000
2050
2100
Fracture lengthL
f
D
n
1
1
1
1
1
Well positionZ
w
h
f
D
n
0.5
0.5
0.5
0.5
0.5
Region lengthx
e
D
4000
4000
4000
4000
4000
Region widthy
e
D
4000
4000
4000
4000
4000
Fracture angleθ
f
n
π/2
π/2
π/2
π/2
π/2
Fracture conductivityC
f
D
n
6π
6π
6π
6π
6πFigure 4
The production rate distribution of MFHWs with five transversal fractures.


### 4.3. Flow Characteristics Analysis

The transient transport characteristics are graphically showed by type curves, which can be used to analyzing transient pressure and rate decline so as to recognize the flow characteristics of fluids in the reservoir. In addition, by type curves matching some reservoir property parameters, such as permeability, skin factor, gas in place, fracture half-length, and gas reservoir drainage area, can be obtained [15, 16].The pressure and its derivative curves are presented in Figure5, which shows flow characteristics for a fractured horizontal well with 5 identical lengths and uniform-distribution fractures, and the used data is shown in Table 2. As shown in Figure 5, the flow can be divided into six stages, which could be corresponding with Figure 6 and described as follows.Figure 5
Type curves used to analyze the flow characteristics.Six kinds of different flow patterns for MFHWs.
(a)
(b)
(c)
(d)
(e)
(f)Stage 1. In this stage, the segment has a straight line with 1/4 slope, reflecting the bilinear flow region (see Figure 6(a)). In this region, fluids flow from fracture to wellbore and from reservoirs to fracture at the same time. This region only occurs when the fracture conductivity is relatively small.Stage 2. The segment has a 1/2 slope straight line on both pressure and pressure derivative curves, namely, the first linear flow region. In this region, flow occurs linearly and directly from formation to individual fractures and each fracture behaves independently of the other fractures (see Figure 6(b)). The linear flow is the optimal flow mode because it can reduce the seepage resistance.Stage 3. If the pressure derivative curve shows the 1/(2
N) constant, the first radial flow region will be observed. In this region, radial flow occurs directly (see Figure 6(c)) from formation to individual fractures and the continuous time of radial flow for each fracture depends on the fracture length and fracture spacing. A big fracture length or big spacing between the fractures will lead to a short time for radial flow.Stage 4. The segment has a 0.36 slope straight line on both pressure and pressure derivative curves, namely, the second linear flow region. The reservoir linear flow occurs in this region (Figure 6(d)). Fluids flow will start to occur when the pressure wave spreads farther. The linear flow starts from reservoir to the vertical plane which contains the horizontal wellbore.Stage 5. The segment has a straight line with 0.5 constant, namely, the second radial flow region (Figure 6(e)). During this flow period, flow across the outer most producing elements becomes dominant.Stage 6. The segment has a unit slope straight line on both pressure and pressure derivative curves, namely, boundary dominated flow region. This flow period occurs when the reservoir boundary is reached (Figure 6(f)). Due to the fractured-horizontal wells centered in the reservoir, the fluids flow reached the six boundaries nearly at the same time. If the fractured-horizontal wells were not centered in the reservoir, the characteristics of boundary dominated flow may be not a unit slope straight line on both pressure and pressure derivative curves.The proposed flow regions above may not exist in a single test. Depending on the specific properties of fracture and reservoir, some of the flow regions may be absent. For the moderate ratios of fracture length and spacing, the boundary dominated flow region may be nonexistent or replaced by a transitional flow region. It is also possible that radial flow region may not exist because of well interference or boundary effects. Admittedly, affected by fracture conductivity and relative position, some characteristics of flow pattern may not clearly emerge at times. In the real case, engineers should analyze the results of flow regions according to the concrete case.


### 4.4. Effects of Sensitive Parameters on Type Curves

Based on the presented model, the main factors affecting transient-pressure characteristics of fractured horizontal wells in coal reservoirs are presented in Figures7–9. Those parameters are including fracture number, fracture conductivity C
f
D
n, and reservoir width y
e
D. The basic used parameters are listed in Table 2.Figure 7
The effect of fracture number on type curve.Figure 8
The effect of fracture conductivity on type curve.Figure 9
The effect of reservoir width on type curve.The effect of fracture number is shown in Figure7 and fracture number can be set to 3 and 7, respectively. It appears that the fracture number mainly affects the flow characteristics of bilinear flow region, the first linear flow region, and the first radial flow region, which could be seen on the pressure derivative curve. It can be seen from pressure derivative curve that the time of the first radial flow region with big fracture number is longer than that with small fracture number and the time of the second linear flow region with small fracture number is longer than that with big fracture number. It is known that the production mainly benefits from the second linear flow region, and the longer the second linear flow region is, the better the production is. Therefore, through the analysis above, increasing the fracture number is helpful to production, but it is not the more the fracture, the better the development. It can be analyzed from the pressure curve that pressure curve with big fracture number is always higher than the one with small fracture number, which means the increasing fractures number will cause small pressure depletion.The effect of fracture conductivity is shown in Figure8 and fracture conductivity can be set to 6 and 7, respectively. The derivative curve shows that the effect of the fracture conductivity on flow characteristics occurs obviously in the bilinear flow region. As shown in the figure, the bilinear flow characteristic will gradually disappear as the conductivity increases; finally, linear flow will appear when C
f
D
>
300. It can be also found from pressure curve that the pressure with small conductivity is higher than the one with big conductivity, which means fracture with high conductivity will cause big pressure depletion at early time.The effect of reservoir width is shown in Figure9 and reservoir width can be set to 1000 and 4000, respectively. Reservoir width mainly affects the boundary dominated flow region. It can be found from the pressure curve that pressure with big reservoir width is higher than the one with small reservoir width, which means that a small reservoir will cause big pressure depletion in the boundary dominated flow region. We also see from the pressure derivative curve that a small reservoir width value will shorten the time of the second radial flow and reach the boundary dominated flow region early.


## 4.1. Validation of Solutions in This Paper

Riley et al. [14] presented an analytical solution for a single fracture case. To validate the solutions proposed in this paper, some data were utilized (see Table 1) for comparing the solutions with those from Riley et al. [14] (see Figure 3). For the convenience of the comparison, only a fracture is considered into the MFHWs. As shown in Figure 3, there is a good agreement between the solutions obtained in this work and the results presented by Riley. This might indicate that previous models are special cases of our MFHWs model and also implies that our model is validated accurately.Table 1
Basic data.
Parameters
Dimensionless value
Fracture numberN
f
1
Fracture segment numberN
sn
16
Fracture positionx
f
D
n
2000
Fracture positiony
f
D
n
2000
Fracture lengthL
f
D
n
1
Well positionZ
w
h
f
D
n
0.5
Region lengthx
e
D
5000
Region widthy
e
D
5000
Fracture angleθ
f
n
π/2Figure 3
Validation of solutions presented in this paper.


## 4.2. The Production Rate Distribution

The production rate distribution is presented in Figure4 and the data is shown in Table 2. In this case, five transversal fractures have identical length and uniform distribution along a horizontal well located in a box-shaped reservoir. As shown in Figure 4, the outermost fractures have higher production ratio due to larger drainage area and the inner fractures are lower due to strong interference between fractures. For a specific fracture, the flux distribution along the discrete elements has the same property as the whole fracture system.Table 2
Basic data.
Parameters
Frac1
Frac2
Frac3
Frac4
Frac5
Fracture numberN
f
5
5
5
5
5
Fracture segment numberN
sn
20
20
20
20
20
Fracture positionx
f
D
n
2000
2000
2000
2000
2000
Fracture positiony
f
D
n
1900
1950
2000
2050
2100
Fracture lengthL
f
D
n
1
1
1
1
1
Well positionZ
w
h
f
D
n
0.5
0.5
0.5
0.5
0.5
Region lengthx
e
D
4000
4000
4000
4000
4000
Region widthy
e
D
4000
4000
4000
4000
4000
Fracture angleθ
f
n
π/2
π/2
π/2
π/2
π/2
Fracture conductivityC
f
D
n
6π
6π
6π
6π
6πFigure 4
The production rate distribution of MFHWs with five transversal fractures.


## 4.3. Flow Characteristics Analysis

The transient transport characteristics are graphically showed by type curves, which can be used to analyzing transient pressure and rate decline so as to recognize the flow characteristics of fluids in the reservoir. In addition, by type curves matching some reservoir property parameters, such as permeability, skin factor, gas in place, fracture half-length, and gas reservoir drainage area, can be obtained [15, 16].The pressure and its derivative curves are presented in Figure5, which shows flow characteristics for a fractured horizontal well with 5 identical lengths and uniform-distribution fractures, and the used data is shown in Table 2. As shown in Figure 5, the flow can be divided into six stages, which could be corresponding with Figure 6 and described as follows.Figure 5
Type curves used to analyze the flow characteristics.Six kinds of different flow patterns for MFHWs.
(a)
(b)
(c)
(d)
(e)
(f)Stage 1. In this stage, the segment has a straight line with 1/4 slope, reflecting the bilinear flow region (see Figure 6(a)). In this region, fluids flow from fracture to wellbore and from reservoirs to fracture at the same time. This region only occurs when the fracture conductivity is relatively small.Stage 2. The segment has a 1/2 slope straight line on both pressure and pressure derivative curves, namely, the first linear flow region. In this region, flow occurs linearly and directly from formation to individual fractures and each fracture behaves independently of the other fractures (see Figure 6(b)). The linear flow is the optimal flow mode because it can reduce the seepage resistance.Stage 3. If the pressure derivative curve shows the 1/(2
N) constant, the first radial flow region will be observed. In this region, radial flow occurs directly (see Figure 6(c)) from formation to individual fractures and the continuous time of radial flow for each fracture depends on the fracture length and fracture spacing. A big fracture length or big spacing between the fractures will lead to a short time for radial flow.Stage 4. The segment has a 0.36 slope straight line on both pressure and pressure derivative curves, namely, the second linear flow region. The reservoir linear flow occurs in this region (Figure 6(d)). Fluids flow will start to occur when the pressure wave spreads farther. The linear flow starts from reservoir to the vertical plane which contains the horizontal wellbore.Stage 5. The segment has a straight line with 0.5 constant, namely, the second radial flow region (Figure 6(e)). During this flow period, flow across the outer most producing elements becomes dominant.Stage 6. The segment has a unit slope straight line on both pressure and pressure derivative curves, namely, boundary dominated flow region. This flow period occurs when the reservoir boundary is reached (Figure 6(f)). Due to the fractured-horizontal wells centered in the reservoir, the fluids flow reached the six boundaries nearly at the same time. If the fractured-horizontal wells were not centered in the reservoir, the characteristics of boundary dominated flow may be not a unit slope straight line on both pressure and pressure derivative curves.The proposed flow regions above may not exist in a single test. Depending on the specific properties of fracture and reservoir, some of the flow regions may be absent. For the moderate ratios of fracture length and spacing, the boundary dominated flow region may be nonexistent or replaced by a transitional flow region. It is also possible that radial flow region may not exist because of well interference or boundary effects. Admittedly, affected by fracture conductivity and relative position, some characteristics of flow pattern may not clearly emerge at times. In the real case, engineers should analyze the results of flow regions according to the concrete case.


## 4.4. Effects of Sensitive Parameters on Type Curves

Based on the presented model, the main factors affecting transient-pressure characteristics of fractured horizontal wells in coal reservoirs are presented in Figures7–9. Those parameters are including fracture number, fracture conductivity C
f
D
n, and reservoir width y
e
D. The basic used parameters are listed in Table 2.Figure 7
The effect of fracture number on type curve.Figure 8
The effect of fracture conductivity on type curve.Figure 9
The effect of reservoir width on type curve.The effect of fracture number is shown in Figure7 and fracture number can be set to 3 and 7, respectively. It appears that the fracture number mainly affects the flow characteristics of bilinear flow region, the first linear flow region, and the first radial flow region, which could be seen on the pressure derivative curve. It can be seen from pressure derivative curve that the time of the first radial flow region with big fracture number is longer than that with small fracture number and the time of the second linear flow region with small fracture number is longer than that with big fracture number. It is known that the production mainly benefits from the second linear flow region, and the longer the second linear flow region is, the better the production is. Therefore, through the analysis above, increasing the fracture number is helpful to production, but it is not the more the fracture, the better the development. It can be analyzed from the pressure curve that pressure curve with big fracture number is always higher than the one with small fracture number, which means the increasing fractures number will cause small pressure depletion.The effect of fracture conductivity is shown in Figure8 and fracture conductivity can be set to 6 and 7, respectively. The derivative curve shows that the effect of the fracture conductivity on flow characteristics occurs obviously in the bilinear flow region. As shown in the figure, the bilinear flow characteristic will gradually disappear as the conductivity increases; finally, linear flow will appear when C
f
D
>
300. It can be also found from pressure curve that the pressure with small conductivity is higher than the one with big conductivity, which means fracture with high conductivity will cause big pressure depletion at early time.The effect of reservoir width is shown in Figure9 and reservoir width can be set to 1000 and 4000, respectively. Reservoir width mainly affects the boundary dominated flow region. It can be found from the pressure curve that pressure with big reservoir width is higher than the one with small reservoir width, which means that a small reservoir will cause big pressure depletion in the boundary dominated flow region. We also see from the pressure derivative curve that a small reservoir width value will shorten the time of the second radial flow and reach the boundary dominated flow region early.


## 5. Conclusions

Through using Fourier cosine transform and Laplace transform method, a semianalytical solution for multifractured horizontal wells in box-shaped reservoirs is presented to investigate the pressure transient behavior. The accuracy of the semianalytical solutions is verified by graphically comparing them with the published analytical solutions which are special cases of this paper. According to the presented solutions, the production rate distribution is also proposed for MFHWs, which show the outermost fractures have higher production ratio due to larger drainage area and the inner fractures are lower due to the strong interface between fractures. Type curves are established to analyze the flow characteristics, which can be divided into six stages. Effects of main sensitive parameters on type curves are analyzed in detail including fracture number, fracture conductivity, and reservoir width. The presented model could be applicable to develop well test soft for MFHWs, which will be our future research topic.



---

*Source: 716390-2014-04-03.xml*