Modeling the characteristics of Bingham porous-flow mechanics for a horizontal well in a heavy oil reservoir

Accepted Manuscript Modeling the characteristics of Bingham porous-flow mechanics for a horizontal well in a heavy oil reservoir Ren-Shi Nie, Yi-Min Wang, Yi-Li Kang, Yong-Lu Jia PII:

S0920-4105(18)30601-6

DOI:

10.1016/j.petrol.2018.07.026

Reference:

PETROL 5123

To appear in:

Journal of Petroleum Science and Engineering

Received Date: 7 March 2018 Revised Date:

29 June 2018

Accepted Date: 9 July 2018

Please cite this article as: Nie, R.-S., Wang, Y.-M., Kang, Y.-L., Jia, Y.-L., Modeling the characteristics of Bingham porous-flow mechanics for a horizontal well in a heavy oil reservoir, Journal of Petroleum Science and Engineering (2018), doi: 10.1016/j.petrol.2018.07.026. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Modelling the characteristics of Bingham porous-flow

2

mechanics for a horizontal well in a heavy oil reservoir

3

Ren-Shi Nie, Yi-Min Wang, Yi-Li Kang*, Yong-Lu Jia

4

State Key Laboratory of Oil & Gas Reservoir Geology and Exploitation, Southwest

5

Petroleum University, Chengdu 610500, Sichuan Province, People’s Republic of China.

6

* Yi-Li Kang is the corresponding author.

7

In order to investigate hydrodynamic characteristics of Bingham porous-flow mechanics in a

8

heavy oil reservoir, this paper originally proposes a physical model of Bingham fluid flow for

9

a horizontal well in a heavy oil reservoir. A corresponding mathematical model is established

10

and solved. A series of typical pressure and pressure derivative curves are drawn by the use of

11

the semi-analytical solution of the mathematical model. Hydrodynamic flow characteristics of

12

Bingham fluid are thoroughly analyzed. Through the recognition of flow regimes from the

13

typical curves, it is found that the flow process of Bingham fluid mainly includes the early

14

radial flow regime, linear flow regime and late pseudo-radial flow regime. The curves

15

controlled by different boundary conditions are especially studied, and it is found that the late

16

curve shape, especially the shape of pressure derivative curve, differs completely for the

17

various property of boundary. The curve sensitivity to different model parameters, such as

18

dimensionless threshold pressure gradient and horizontal wellbore length, is deeply discussed.

19

The findings are as follow: except the pure wellbore storage regime, the threshold pressure

20

gradient influences the hydrodynamic flow behavior in all the regimes; the curve of Bingham

21

fluid governed by the threshold pressure gradient upward deviates from that of Newtonian

22

fluid; and the deviation of the curves increases with the increase of the threshold pressure

23

gradient and the elapsed time. The research results of this paper sufficiently demonstrate the

AC C

EP

TE D

M AN U

SC

RI PT

1

1

ACCEPTED MANUSCRIPT 24 25

essential property of non-Newtonian porous-flow mechanics. Keywords: hydrodynamic characteristics; non-Newtonian fluid; Bingham fluid; pressure and pressure derivative curves; threshold pressure gradient.

27

1. Introduction

RI PT

26

Bingham fluid as a kind of non-Newtonian fluid has attracted more and more interests to

29

scientists and researchers (Frigaard and Nouar, 2003; Zhang, 2003; Osalusi et al., 2007;

30

Firouzi and Hashemabadi, 2008; Staron et al., 2013; Liu and Liu, 2014). The Bingham fluids

31

belong to viscous plastic fluids (Ancey et al., 2009; Chevalier et al., 2013; Chevalier et al.,

32

2015) and exhibit a yield stress (Abdali et al., 1992; Li and Yu, 2010). Only after the applied

33

stress exceeds the yield stress can Bingham fluids flow (Tichy, 1991). Bingham fluids

34

beginning to flow must overcome the threshold pressure gradient caused by the yield stress

35

(Oliveira et al., 2012). Besides, the flow of Bingham fluid considering the threshold pressure

36

gradient follows the modified Darcy’ law in porous media (Bear, 1972; Wu et al., 1992).

37

Based on the modified Darcy’ law, Wu (1990) developed a fully implicit, integral finite

38

difference model for simulation of Bingham fluid flow through porous media. Later, Martinez

39

et al. (2011) presented a technique for interpreting the behavior of pressure and pressure

40

derivative for a Bingham type fluid in a homogeneous reservoir drained by a vertical well

41

using the TDS technique.

AC C

EP

TE D

M AN U

SC

28

42

Despite the copious literature on the flow issues of Bingham fluid, little work is available in

43

petroleum literature on the hydrodynamic circumstance of Bingham fluid flow for horizontal

44

wells with different boundary properties. Most works concerned the hydrodynamic 2

ACCEPTED MANUSCRIPT circumstance of Newtonian fluid flow (Zerzar and Bettam, 2004; Guo et al., 2012; Nie et al.,

46

2012; Luo et al., 2014). The problem set of the flow to horizontal well in porous formations

47

have been studied during last decade. The flow problem set mainly includes two aspects: one

48

is the establishment of physical and mathematical model for various types of porous

49

formation; and the other is the types of fluids. Main types of porous formation for horizontal

50

well flow models were studied: the horizontal well flow model in a homogenous formation

51

(Kawecki, 2000); the horizontal well flow model in a composite formation that had two

52

regions with distinct properties traversed by a horizontal well (Raghavan, 2012); the

53

horizontal well flow model in a naturally fractured formation (Nie et al., 2012); and the

54

horizontal well flow model in a naturally fractured-vuggy carbonate formation (Guo et al.,

55

2012). Three typical types of fluids were considered in underground porous formations: the

56

first type was light oil considered as Newtonian fluid (Ozkan, 2001; Zerzar and Bettam, 2004;

57

Luo et al., 2014); the second type was water (Zhan and Zlotnik, 2002; Samani et al., 2004;

58

Sun and Zhan, 2006); and the third type was gas (Tabatabaei and Zhu, 2010; Zhu et al., 2015).

59

As stated above, the flow problem of horizontal well in porous formations focused on the

60

Newtonian fluids. Therefore, it is very interesting to study the flow problem of horizontal well

61

focused on the Bingham fluids. Heavy oil is a kind of Bingham fluid (Mirzajanzade et al.,

62

1971; Barenblatt et al., 1984; Oosthuizen et al., 2007). Studying characteristics of Bingham

63

porous-flow mechanics for a horizontal well in a heavy oil reservoir will be important for

64

heavy oil development. Therefore, a mathematical model of Bingham fluid flow for a

65

horizontal well in a heavy reservoir was proposed originally in this paper and then solved

66

under different boundary conditions. Besides, new log-log curves of p wD and p wD ′ were

AC C

EP

TE D

M AN U

SC

RI PT

45

3

ACCEPTED MANUSCRIPT drawn. At the end, the hydrodynamic behavior of Bingham fluid flow was analyzed

68

thoroughly. The research results of this paper will enable the public to comprehend the

69

characteristics of porous-flow mechanics for Bingham fluid flow in underground formations.

70

The presented model of Bingham fluid will enrich the theoretical model of non-Newtonian

71

fluid mechanics.

72

2. Modeling of bingham fluid flow

73

2.1. Physical Model

M AN U

SC

RI PT

67

A single horizontal well in an underground formation is established, whose production rate

75

is a constant (see Fig 1), when the wellbore and formation exists a pressure drop, the Bingham

76

fluid in the formation, such as the visco-plastic oil with slight elasticity under the high

77

pressure conditions of the formation, would flow into the wellbore.

TE D

74

Physical model assumptions:

79

(1) Bingham fluid flow with a threshold pressure gradient (Bear, 1972; Wu et al., 1992) is

81 82

considered;

(2) The formation is assumed as homogeneity, so the formation permeability has no

AC C

80

EP

78

dependence on the formation location;

83

(3) The formation is assumed as the isotropy only in the horizontal plane, so the horizontal

84

permeability of the formation has no dependence on the direction angle and the horizontal

85

permeability of the formation is not equal to the perpendicular permeability;

86 87

(4) The external top boundary and bottom boundary may be constant pressure or closed, the external side boundary may be constant pressure or closed or infinite (see Fig 1); 4

ACCEPTED MANUSCRIPT 88

(5) It is assumed that the Bingham fluid in the formation flows into the horizontal wellbore

89

is uniform and along the wellbore each point source is same considering open-hole

90

completion; (6) Wellbore storage effect is considered (Mathias and Butler, 2007; Park and Zhan, 2002);

92

(7) Skin effect is considered (Park and Zhan, 2002).

94

2.2. Mathematical Model

SC

93

RI PT

91

To solve the mathematical model (see in Fig 1), we establish a set of coordinate system (see Fig 2).

96

The conversion relationships of all coordinate are given below:

r=

x2 + y2

98

r′ =

y 2 + z ′2

99

z ′ = z − zw

TE D

97

M AN U

95

(1) (2) (3)

In a radial cylindrical system, the governing differential equation for a horizontal well (Guo

101

et al., 2012; Nie et al., 2012; Wu et al., 1992) is given by Eq. 4. The derivation process of the

102

Eq. 4 is given in Appendix. It is notable that the governing equation in this cylindrical system

103

does not have the polar dependence because the formation was assumed as isotropy in the

104

horizontal plane.

AC C

105

EP

100

∂2 p 1 ∂p λB kp ∂2 p µφCt ∂p + − + = ∂r 2 r ∂r r kh ∂z 2 3.6kh ∂t

(4)

106

where Cρ is fluid compressibility, 1/MPa; λB is threshold pressure gradient, MPa/m; r is radial

107

coordinate, m; Ct is total compressibility, 1/MPa; kp is perpendicular permeability, µm2; kh is

108

horizontal permeability, µm2; ϕ is porosity, fraction; µ is viscosity, mPa s; p is pressure, MPa; 5

ACCEPTED MANUSCRIPT

112 113 114

standard international units. (1) Inner boundary condition:

RI PT

111

Eq. 4 has a coefficient equal to 3.6 because it uses a set of engineering units instead of

The rate of well production is a constant (Guo et al., 2012; Nie et al., 2012; Wu et al., 1992): kh r → 0 1.842 × 10 −3 µε

115

q% = lim[lim

116

qBρ = ∫

117

pw

ε →0

zw −ε / 2

r(

∂p − λB )dz ] ∂r

L /2

q% ( r , t )dx

− L /2

r →0

∫

zw +ε / 2

=∫

L/2

−L/ 2

SC

110

z is perpendicular coordinate, m; t is time, h.

M AN U

109

p(r , t )dx

(5)

(6) (7)

where q% is production rate of a point source, m3/d; ε is a variable in z direction, m; q is

119

production rate, m3/d; L is horizontal wellbore length, m; Bρ is fluid volume factor,

120

dimensionless; x is the x direction that is along the horizontal wellbore, m; pw is wellbore

121

pressure, MPa.

124

p

t =0

= pi

where pi is initial formation pressure, MPa;

125

(3) External boundary conditions:

126

①For bottom

127

128 129

(8)

EP

123

(2) Initial conditions:

AC C

122

TE D

118

closed boundary:

∂p ∂z

z =0

=0

constant pressure boundary: p

(9)

z =0

= pi

② For top

6

(10)

ACCEPTED MANUSCRIPT ∂p ∂z

=0

130

closed boundary:

131

constant pressure boundary: p

constant pressure boundary: p

141 142

closed boundary:

∂p ∂r

r = re

=0

2.3. Dimensionless mathematical model Dimensionless z coordinate

zD = z / h

Dimensionless perpendicular distance

zwD = zw / h

144

Total skin factor, dimensionless

AC C

where zw is perpendicular distance, m.

S=

148

kh h ∆ps 1.842 × 10 −3 q µ Bρ

St =

146

147

= pi

where re is radial radius, m.

143

145

r = re

TE D

140

r →∞

EP

139

RI PT

135

138

(12)

③ For side infinite boundary: lim p = pi

137

= pi

where h is formation thickness, m.

134

136

z =h

SC

133

(11)

M AN U

132

z =h

St =

L kp ⋅S 2h k h

( L / 2) k p k h 1.842 × 10 −3 qµ Bρ

where ∆ps is additional pressure drop near wellbore, MPa. 7

∆ps

(13) (14)

(15)

ACCEPTED MANUSCRIPT 149

Dimensionless radial radius of side boundary

150

reD =

Dimensionless radial distance in horizontal plane

154

Dimensionless time tD =

hD =

pD =

EP AC C

163

166

khh ( pi − p ) 1.842 × 10 −3 qBρ µ

CD =

Cs 6.2832φ C t hrw 2

where Cs is wellbore storage coefficient, m3/MPa. Dimensionless threshold pressure gradient

164

165

kh kp

Dimensionless wellbore storage coefficient

161

162

h rw e− S

TE D

Dimensionless pressure

159

160

3.6kh t φµ Ct rw 2

M AN U

Dimensionless formation thickness

157

158

−S

where rw is real wellbore radius, m.

155

156

r e rw

RI PT

rD =

152

153

e rw

SC

151

r

e −S

λBD =

khh λB rw 1.842 × 10 −3 Bρ q µ

In a radial cylindrical system, the dimensionless governing equation for a horizontal well is given by:

8

ACCEPTED MANUSCRIPT 167

∂2 pD 1 ∂pD 1 ∂2 pD λBDe− S ∂p + + 2 + = e−2S D 2 2 ∂rD rD ∂rD hD ∂zD rD ∂tD

168

(1) Inner boundary condition:

169

The rate of well production is a constant

RI PT

172

rD → 0

∂pD ) = −(1 + λBD e− S ) ∂rD

(2) Initial conditions: pD

tD = 0

=0

173

(3) External boundary conditions:

174

① For bottom

176 177

closed boundary:

∂p D ∂z D

zD = 0

=0

(20)

(21)

∂p D ∂z D

closed boundary:

179

constant pressure boundary: p D

EP

AC C

③ For side

z D =1

181

infinite boundary: lim pD = 0

182

constant pressure boundary: p D

183

zD =0

=0

178

180

z D =1

=0

∂p D ∂rD

rD = reD

(22)

(23)

rD →∞

closed boundary:

(19)

=0

constant pressure boundary: p D

② For top

(18)

TE D

175

(17)

SC

171

lim (rD

M AN U

170

(16)

rD = reD

=0

=0

(24) (25)

9

ACCEPTED MANUSCRIPT 184

3. Solution to mathematical modeling

185

3.1. Laplace transform Based on tD, introducing the Laplace transform, Eq. 26 can be given as follow:

187

L[ pD ( rD , tD )] = p D ( rD , u ) = ∫ pD ( rD , tD )e − utD dtD

∞

0

190

∂ 2 pD 1 ∂ pD λBDe− S 1 ∂2 pD −2 S + + − u e p = − D ∂rD 2 rD ∂rD hD 2 ∂zD 2 urD

191

(1) Inner boundary conditions

192

The rate of well production is a constant

∂ pD (1 + λBD e− S ) lim(rD )=− rD →0 ∂rD u

194

(2) External boundary conditions:

195

① For bottom

∂ pD ∂z D

=0

closed boundary:

197

constant pressure boundary: p D

AC C

EP

196

198

zD =0

∂ pD ∂z D

(29)

zD = 0

=0

=0

closed boundary:

200

constant pressure boundary: p D

202

(28)

(30)

② For top

199

201

(27)

TE D

193

(26)

Thus, the governing differential equation after Laplace transform is given by:

SC

189

where u is time, dimensionless; L[] is Laplace operator.

M AN U

188

RI PT

186

zD =1

(31)

z D =1

=0

(32)

③ For side infinite boundary: lim p D = 0

(33)

rD →∞

10

ACCEPTED MANUSCRIPT 203

constant pressure boundary: p D

204

closed boundary:

208

209 210

=0

(35)

Eq. 27 is a nonhomogeneous differential equation because of the right term in Eq. 27, and the corresponding homogeneous equation is written by

∂2 pD 1 ∂ pD 1 ∂2 pD + + − ue−2 S p D = 0 ∂rD 2 rD ∂rD hD 2 ∂zD 2

RI PT

207

(34)

(36)

SC

206

rD = reD

=0

3.2. Separation of variables

In Laplace space, the dimensionless pressure can be separated by

p D = Z ( zD ) R(rD )

M AN U

205

∂ pD ∂rD

rD = reD

(37)

where Z is separated variable in vertical direction, dimensionless; R is separated variable

212

in horizontal direction, dimensionless.

213

Substituting Eq. 37 into Eq. 36 :

1 ′ R + R ′′ − ue −2 S R Z ′′ r hD2 D =− =λ R Z

EP

215

1 1 Z R ′ + Z R ′′ + 2 Z ′′ R − f ( u ) Z R = 0 rD hD

216

So

217

R′′ +

218

ξ = ue −2 S +

219

(38)

(39)

AC C

214

TE D

211

1 ′ R −ξ R = 0 rD

(40)

λ

(41)

hD2

Z ′′ + λ Z = 0

(42)

220

Solutions in vertical direction

221

The general solution of Eq. 42 is given as follow: 11

ACCEPTED MANUSCRIPT 222 223 224

Z = C cos( λ zD ) + D sin( λ zD )

(43)

where C and D are undetermined coefficients. Substituting Eq. 43 into the equations of external boundary conditions of top and bottom, we can solve the coefficients C and D by using the eigenfunction and eigenvalue method

226

(Guo et al., 2012; Nie et al., 2012):

229

1 C = cos( λn zwD ), λn = (n π)2 , n = 0,1, 2 L 2

D=0,λ = λn = (n π)2 , n = 0,1,2 L

SC

228

(1) Top and bottom boundaries are both closed boundary:

M AN U

227

RI PT

225

230

where λ is equal to cos( λ zD ) which is a eigenvalue of eigenfunction.

231

(2) Top and bottom boundaries are both constant pressure boundary:

232

C =0,λ = λn = (n π)2 , n = 0,1,2 L

235 236

TE D

(45)

(46) (47)

where λ is equal to sin( λ zD ) which is a eigenvalue of eigenfunction. (3) The top boundary is closed boundary and bottom boundary is constant pressure boundary

EP

234

1 D= sin( λn zwD ), λn = (n π)2 , n = 0,1,2 L 2

AC C

233

(44)

237

1 C =0,λn = [(n- )π]2 , n = 0,1, 2 L 2

(48)

238

1 1 D= sin( λn zwD ), λn = [(n- )π]2 , n = 0,1, 2 L 2 2

(49)

239 240 241

(4) The top boundary is constant pressure boundary and bottom boundary is closed boundary

1 1 C = cos( λn zwD ), λn = [(n- )π]2 , n = 0,1, 2 L 2 2 12

(50)

ACCEPTED MANUSCRIPT (51)

Establish a new perpendicular coordinate z′ and radial coordinate r′ (see figure 2d):

244

r′ = y 2 + z′2 = y 2 + ( z − zw )2

245

rD′ = zD′ 2 + y D 2 = (

246

zD =

(52)

zD h − z wD h 2 h ) + yD 2 = ( zD − zwD ) 2 ( )2 + yD 2 rw rw

RI PT

243

1 D=0,λn = [(n- )π]2 , n = 0,1, 2 L 2

rw rD′ 2 − yD 2 + zwD h

At the wall of wellbore (Guo et al., 2012):

248

r = rw , rD′ = 1, yD = 0

250

251

zD =

rw + zwD h

(53)

(54)

(55)

(56)

The solution in vertical direction at the wall of wellbore is Eq. 57

Z = C cos[ λ (

rw r + zwD )] + D sin[ λ ( w + zwD )] h h

TE D

249

M AN U

247

SC

242

(57)

Solutions in horizontal direction

253

Eq. 40 is a standard Bessel function, thus the general solution of Eq. 40 can be expressed

255

by (Amos, 1974; Asmar, 2005)

R = AI0 (rD ξ ) + BK0 (rD ξ )

(58)

AC C

254

EP

252

256

where I0( ) is the first kind zero-order modified Bessel function; K0( ) is the second kind

257

zero-order modified Bessel function; A and B are Unknown coefficients.

258

(1)

259

The nonhomogeneous Eq. 27 has a particular solution for infinite boundary (Agarwal and

260

Infinite boundary

Regan, 2009; Asmar, 2005):

13

ACCEPTED MANUSCRIPT ∗

R =

261

1

G ( rD , σ )d σ

(59) 1 < σ < rD (60) rD < σ < ∞

where G(rD,σ) is Green function.

RI PT

264

∞

 λBDe − S  u K 0 ( ξ rD )I 0 ( ξσ ) G ( rD , σ ) =  −S  λBDe K ( ξσ )I ( ξ r ) 0 0 D  u

262

263

∫

At the bottom hole of the horizontal well, Eq. 59 can be expressed by

πλBDe− S R = I0 ( ξ ) 2u ξ

265

R = AI 0 ( ξ rD ) + BK 0 ( ξ rD ) + R

268 269

*

(62)

Substituting Eq. 62 into Eq. 33, the coefficient A can be obtained by

A= 0

270 271

M AN U

by

TE D

267

The general solution of the nonhomogeneous Eq. 27 in horizontal plane can be expressed

(63)

Substituting Eq. 62 and Eq. 63 into Eq. 28, the coefficient B can be obtained by

(1+ λBDe- S ) B= u

272

(64)

EP

266

(61)

SC

∗

(2)

274

The nonhomogeneous Eq. 27 has a particular solution for a finite boundary (Agarwal and

275

Regan, 2009; Asmar, 2005): ∗

R =

276

277 278 279

Constant pressure or closed boundary

AC C

273

λBD e − S u

∫

rD

1

[K 0 ( ξ rD )I 0 ( ξ ε ) −K 0 ( ξ ε )I0 ( ξ rD )]dε

(65)

The general solution of the nonhomogeneous Eq. 27 in horizontal plane can be expressed by R = AI 0 ( ξ rD ) + BK 0 ( ξ rD ) + R

*

(66) 14

ACCEPTED MANUSCRIPT 280 281

Substitute Eq. 66 into Eq. 28: lim[ rD ξ I1 ( ξ rD ) A − rD ξ K1 ( ξ rD ) B ] +

rD → 0

lim

rD → 0

λBD e − S rD u

reD

∫

1

[K 0 ( ξ reD )I 0 ( ξσ ) −K 0 ( ξσ )I 0 ( ξ reD )]dσ = −

283

For constant pressure boundary, substitute Eq. 66 into Eq. 34:

284

AI0 ( ξ reD ) + BK0 ( ξ reD ) + reD

1

[K 0 ( ξ reD )I 0 ( ξσ ) −K 0 ( ξσ )I 0 ( ξ reD )]dσ = 0

286

Substitute Eq. 66 into Eq. 35:

287

A ξ I1(reD ξ ) − ξ BK1(reD ξ ) −

288

λBD e − S u

I1 ( ξ reD ) ∫

ξ reD

ξ

K 0 (σ )dσ + K1 ( ξ reD ) ∫

ξ reD ξ

I 0 (σ )dσ = 0

where I1( ) is a first kind first-order modified Bessel function; K1( ) is a second kind

290

first-order modified Bessel function.

293

3.3. Solution

AC C

In Laplace space, the solution at the wall of the wellbore is given as follows: ∞

p sD = ∑ ∫ n=0

296

EP

Eq. 67 and Eq. 69.

295

L / 2 / rw

− L / 2 / rw

R ( xD , λn )dxD ⋅Z w (λn )

(70)

where p sD is dimensionless wellbore pressure in Laplace space. ∗

297

298

(69)

Coefficients A and B can be calculated by combining Eq. 67 and Eq. 68 and by combining

292

294

(68)

TE D

289

291

(67)

SC

u

∫

M AN U

λBD e − S

285

(1 + λBD e − S ) u

RI PT

282

R( xD , ξn ) = An I0 ( xD ξn ) + Bn K0 ( xD ξn ) + R ξ n = ue −2 S +

λn hD2

, n = 0,1, 2 L

(71) (72)

15

ACCEPTED MANUSCRIPT

300

301

Considering the effect of wellbore storage, the solution can be obtained by Eq. 73 based on Duhamel’s principle: p wD =

p sD CDu p sD + 1

(73)

2

where pwD is dimensionless wellbore pressure.

303

4. Hydrodynamic characteristics

SC

302

RI PT

299

The dimensionless bottomhole pressure ( pwD ) in real space can be obtained by using

305

Stehfest numerical inversion (Stehfest, 1970) to convert pwD back to pwD . Therefore the

306

standard bi-logarithmic curves of pressure and pressure derivative can be drawn using some

307

′ for Bingham computer program. Fig 3 and Figs 5-10 are log-log curves of pwD and pwD

308

fluid flow in an underground porous formation.

Through flow regime recognition, hydrodynamic flow behavior can be perspicuously

TE D

309

M AN U

304

′ . To conveniently observe the pwD and pwD

reflected from the log-log curves of

311

hydrodynamic flow behavior, a couple of pressure and pressure derivative curves were

312

especially drawn using the mathematical model under a group of model parameters – “St=0,

313

CD=500, h=30m, L=300m, kh/kp=10, zwD=0.5, λBD=10-4”, as shown in Fig 3. The following

314

five main flow regimes for a horizontal well in an underground formation with an infinite

315

boundary, a closed top boundary and a closed bottom boundary are recognized from Fig 3:

AC C

EP

310

316

(I) pure wellbore storage regime: in the formation the fluid does not flow, however the fluid

317

′ are overlapped as a straight line stored in wellbore flows; and the curves of pwD and pwD

318

with a unit slope.

319

(II) skin effect regime: the formation fluid near the horizontal wellbore flows toward to the

320

wellbore, however the formation fluid far away from the horizontal wellbore does not flow; 16

ACCEPTED MANUSCRIPT 321

′ diverges gradually from the curve of pwD and the curve of pwD

with the elapsed time.

(III) early radial flow regime (see Fig 4a): in vertical plane, the fluid radially flows around

323

the axis of the horizontal wellbore, which is visually exhibited with the streamlines that are

324

′ is a horizontal line; it orthogonal to the horizontal wellbore axis in Fig 4a; the curve of pwD

325

is not long for the propagation of the pressure wave to the top and bottom closed boundaries

326

because the formation thickness (only 30m for Fig 3) is usually very small compared with the

327

horizontal distance of infinite side boundary; Thus this regime lasts very short time, which is

328

directly manifested in the regime III of Fig 3. A larger formation thickness would lead to

329

relative longer time duration of this regime.

M AN U

SC

RI PT

322

(IV) linear flow regime (see Fig 4b): the fluid linearly flows into the horizontal wellbore in

331

horizontal plane, which is visually shown with a group of parallel streamlines that are

332

′ take on orthogonal to the horizontal wellbore axis in Fig 4b; the curves of pwD and pwD

333

upswept lines; the time duration of this regime is longer than that of the regime (III), which is

334

also directly manifested in Fig 3; and the duration of this regime mainly depends on the length

335

of horizontal wellbore.

TE D

330

(V) late pseudo-radial flow regime (see Fig 4c): the fluid far away from the horizontal

337

wellbore pseudo-radially flows in the formation area near the wellbore, and then flows into

338

the horizontal wellbore, which is visually shown with a group of radial streamlines in

339

′ also take on upswept lines; the horizontal plane in Fig 4c; the curves of pwD and pwD

340

pressure derivative curve gradually tilts up and ultimately runs parallel to the pressure curve;

341

′ demonstrate the hydrodynamic characteristics of and the parallel curves of pwD and pwD

342

infinite acting radial flow.

343 344

AC C

EP

336

The above flow regimes, especially regime (III) and (IV), reflect the typical hydrodynamic flow behavior of horizontal wells in an underground porous formation. 17

ACCEPTED MANUSCRIPT

′ controlled by different side boundary conditions The log-log curves of pwD and pwD

346

were simulated using mathematical model, as shown in Fig 5. It can be seen form Fig 5 that

347

different boundaries lead to different curve characteristics. For the constant pressure boundary,

348

′ the curve of pwD ultimately becomes a horizontal straight line, besides the curve of pwD

349

swiftly decrease (curves “②” in Fig 5); The horizontal straight line of dimensionless

350

pressure represents the pressure does not vary with time any more, which means that the fluid

351

flow has reached a steady state in which the derivative of pressure to time is 0 for the

352

sufficient energy supply from the constant pressure boundary. For the closed boundary, the

353

′ gradually converge to a upswept overlapped line (see curves “③” curves of pwD and pwD

354

in Fig 5). In this regime, the fluid flow has reached a pseudo-steady state in which the second

355

derivative of pressure to time is a constant for no energy supply from the closed pressure

356

boundary.

M AN U

SC

RI PT

345

′ were simulated using the mathematical model under The log-log curves of pwD and pwD

358

constant pressure boundary of bottom, as shown in Fig 6. Fig 6 reflects the curve

359

characteristics influenced by different top boundaries. After regime (III), the curve of pwD

360

′ swiftly decreases for both tends to become a horizontal straight line and the curve of pwD

361

closed and constant pressure boundaries of top. The horizontal straight line of dimensionless

362

pressure reflects that the fluid flow has reached a steady state because of the sufficient energy

363

′ begin to change for the constant supply. However, the time that the curves of pwD and pwD

364

pressure boundary of top is earlier than that for the closed top boundary because top and

365

bottom boundaries both provide energy supply.

AC C

EP

TE D

357

366

The hydrodynamic characteristics of Bingham fluid flow in an underground porous

367

′ affected formation can be analyzed through simulating a group of curves of pwD and pwD

368

by the dimensionless threshold pressure gradient. Under a group of fixed model parameters 18

ACCEPTED MANUSCRIPT (St=0, CD=500, h=30m, L=300m, kh/kp=10, zwD=0.5), the curves were simulated by setting

370

dimensionless threshold pressure gradient (λBD) as 10-3, 10-4, 10-5, and 0, as shown in Fig 7.

371

When the dimensionless threshold pressure gradient is 0, the Bingham fluid model is reduced

372

to the conventional model of Newtonian fluid. The curves of p wD and p wD ′ are influenced by

373

the dimensionless threshold pressure gradient significantly except regime (I), for the Bingham

374

fluid does not flow in the porous formation. The curves of Bingham fluid (see curves ①, ②

375

and ③ in Fig 7) upwardly deviate from the curves of Newtonian fluid (see curves ④ in Fig

376

7) and the deviation of the curves increases with the elapsed time. The gradual deviation of

377

the pressure and pressure derivative curves of Bingham fluid from those of Newtonian fluid is

378

the typical hydrodynamic characteristic of Bingham fluid flow. The larger the dimensionless

379

threshold pressure gradient is, the more conspicuous the hydrodynamic characteristic is. For

380

example, the curve location of “λBD=10-3” is much higher than that of “λBD=10-5” (see curves

381

① and ③ in Fig 7). Another conspicuous hydrodynamic characteristic is that the pressure

382

derivative curve of Bingham fluid is no longer a horizontal straight line in regime (V), which

383

is completely different from the hydrodynamic characteristic of Newtonian fluid where the

384

pressure derivative curve is a horizontal straight line (see Fig 7).

EP

TE D

M AN U

SC

RI PT

369

Well production causes the depletion of formation pressure (pressure drop) and the

386

propagation of pressure wave. The formation can be divided into two areas: one is the swept

387

area that has been swept by the pressure wave; the other is the unswept area that has not been

388

swept by the pressure wave. The pressure drop only exists in the swept area. Fluid does not

389

flow in the unswept area. For Newtonian fluid, fluid flows in the whole swept area. However,

390

for Bingham fluid, the swept area can be divided into two subareas: one is the mobile area

391

where the fluid can flow because threshold pressure gradient is smaller than formation

392

pressure gradient; the other is the stagnant area where the fluid cannot flow because formation

393

pressure gradient is smaller than threshold pressure gradient. Compared with Newtonian fluid,

AC C

385

19

ACCEPTED MANUSCRIPT the flow area is relatively smaller for Bingham fluid under the same conditions of model

395

parameters, such as fluid rate, permeability, and horizontal wellbore length. Therefore, under

396

the same conditions of model parameters, the pressure depletion of Bingham fluid is more

397

quickly than that of Newtonian fluid. According to the definition of dimensionless pressure,

398

the dimensionless pressure is directly proportional to the pressure drop. A larger pressure drop

399

results in a higher curves location of p wD and p wD ′ . The discussion above is the explanation

400

that why the Bingham fluid curves of p wD and p wD ′ (curves ①, ② and ③) are above

401

those of Newtonian fluid (curves ④) respectively, as shown in Fig 7.

SC

RI PT

394

It is the threshold pressure gradient that differentiates Bingham fluid with Newtonian fluid

403

in the porous flow model. Threshold pressure gradient gives rise to an additional pressure

404

drop for Bingham fluid flow in the formation. The larger the threshold pressure gradient is,

405

the greater the value of dimensionless pressure is. Therefore, when compared with the

406

location of curves ② and ③ in Fig 7, the location of curves ① is the highest.

M AN U

402

Under a group of fixed model parameters (St=0, CD=500, h=30m, λBD =10-5, kh/kp=10,

408

zwD=0.5), the pressure and pressure derivative curves were simulated by setting horizontal

409

wellbore length as 300m, 500m and 700m, as shown in Fig 8. A longer horizontal wellbore

410

length means a larger contacted area of the horizontal wellbore with the formation and,

411

accordingly, a strong supply ability of fluids in the formation. Under the same rate production

412

condition, the longer horizontal wellbore needs less pressure depletion. Therefore, the longer

413

the horizontal wellbore length is, the lower the location of curves is, and the longer the

414

duration of linear flow regime is. Fig 9 and Fig 10 show the curves sensitivity of

415

′ to skin factor and the ratio of horizontal and vertical permeability respectively. p wD and p wD

416

A relative larger skin factor means the formation has a more serious pressure drop near the

417

wellbore. Therefore, the greater skin factor causes higher curve location of p wD in the

AC C

EP

TE D

407

20

ACCEPTED MANUSCRIPT whole flow regimes. The dimensionless pressure derivative curves are also affected similarly

419

by skin factor in the regimes (I–IV), and the derivative curves of different skin factors

420

converge to a curved line in regime (V). The ratio of horizontal and vertical permeability

421

represents the anisotropy degree of the horizontal direction with the vertical direction. The

422

larger ratio means a more serious anisotropy for the underground porous formation and a

423

higher curve location mainly in regime (II) and (III).

RI PT

418

In sum, the hydrodynamic characteristics of Bingham fluid for a horizontal well in an

425

underground porous formation are obviously different from those of Newtonian fluid in the

426

whole flow process. The above figures of Bingham fluid flow sufficiently demonstrate the

427

essential property of non-Newtonian porous-flow mechanics.

428

5. Conclusions

M AN U

SC

424

In this paper, we established the Bingham flow model of horizontal well and revealed the

430

hydrodynamic characteristics thoroughly. It was found that the model parameters, such as

431

horizontal length and skin factor, influenced the flow characteristics differently. It was also

432

found that the pressure curve of Bingham fluid deviated from the curve of Newtonian fluid

433

mainly for the influence of threshold pressure gradient corresponding to the yield stress. A

434

larger threshold pressure gradient led to a heavier curve deviation between Bingham and

435

Newtonian fluids. When the threshold pressure gradient is equal to zero, the model of

436

Bingham fluid is reduced to the model of Newtonian fluid and the curve deviation is

437

disappeared accordingly. Under the same conditions of model parameters, the curves of

438

Bingham fluid are above those of Newtonian fluid, which implies that the pressure depletion

439

of Bingham fluid is faster than that of Newtonian fluid. In conclusion, the hydrodynamic

AC C

EP

TE D

429

21

ACCEPTED MANUSCRIPT characteristics of Bingham fluid are apparently different from Newtonian fluid. The results

441

obtained with Bingham fluid promote the model development of non-Newtonian fluid and

442

give some references to the related researchers.

443

Acknowledgments

RI PT

440

The authors would like to thank the NSFC (National Natural Science Foundation of China)

445

for supporting this article through two projects: one is the National Science Fund for Young

446

Scholars of China (Grant No. 51304164)—“Research on the pressure dynamics of

447

multiple-acidized-fractured horizontal wells in fractured-vuggy carbonate formations”; and

448

the other is the National Science Fund for Distinguished Young Scholars of China (Grant No.

449

51525404) — “Fracturing and acidizing in low permeability and tight reservoirs”. The paper

450

was financially supported by the Fok Ying Tung Education Foundation for Young Teachers in

451

the Higher Education Institutions of China (Grant No. 151050). The paper also was also

452

financially supported by a basic research project under Grant No. 2015JY0132 from the

453

Science and Technology Department of Sichuan Province. The authors would also like to

454

thank the reviewers and editors, whose critical comments were very helpful in preparing this

455

article.

456

References

457

Abdali, S.S., Mitsoulis, E., Markatos, N.C., 1992. Entry and Exit Flows of Bingham Fluids.

458 459

AC C

EP

TE D

M AN U

SC

444

Journal of Rheology. 36, 389-407. Agarwal, R.P., Regan, D.O., 2009. Ordinary and partial differential equations with special 22

ACCEPTED MANUSCRIPT functions, Fourier series, and boundary value problems. Publisher: Springer New York.

461

Amadei, B., 2000. A mathematical model for flow of Bingham materials in fractures. 4th

462

North American Rock Mechanics Symposium. American Rock Mechanics Association.

463

Amadei, B., Savage, W.Z., 2001. An analytical solution for transient flow of Bingham

464

viscoplastic materials in rock fractures. International Journal of Rock Mechanics &

465

Mining Sciences. 38, 285–296.

470 471 472 473 474 475

SC

M AN U

469

Ancey, C., Balmforth, N.J., Frigaard, I., 2009. Visco-plastic fluids: from theory to application. Journal of Non-Newtonian Fluid Mechanics. 158, 1-3.

Asmar, N.H., 2005. Partial differential equations with Fourier series and boundary value problems. Publisher: Upper Saddle River, N. J. by Pearson Prentice Hall.

TE D

468

Computation. 28, 239-51.

Barenblatt, G.E., Entov, B.M., Rizhik, B.M., 1984. Flow of liquids and gases in natural formations. Nedra, Moscow.

Bear, J., 1972. Dynamics of fluids in porous media. Publisher: Elsevier Science. New York City.

EP

467

Amos, D.E., 1974. Computation of modified Bessel functions and their ratios. Mathematics of

AC C

466

RI PT

460

476

Chevalier, T., Chevalier, C., Clain, X., Dupla, J.C., Canou, J., Rodts, S., Coussot, P., 2013.

477

Darcy’s law for yield stress fluid flowing through a porous medium. Journal of

478

Non-Newtonian Fluid Mechanics. 195, 57–66.

479

Chevalier, T., Rodts, S., Chateau, X., Chevalier, C., Coussot, P., 2014. Breaking of

480

non-Newtonian character in flows through a porous medium. Physical Review E. 89,

481

023002. 23

ACCEPTED MANUSCRIPT Chevalier, T., Talon, L., 2015. Generalization of Darcy's law for Bingham fluids in porous

483

media: From flow-field statistics to the flow-rate regimes. Physical Review E. 91, 023011.

484

Chien, S.F. 1970. Laminar flow pressure loss and flow pattern transition of Bingham plastics

485

in pipes and annuli. International Journal of Rock Mechanics and Mining Sciences &

486

Geomechanics Abstracts. Pergamon. 7, 339-356.

RI PT

482

Firouzi, M., Hashemabadi, S.H., 2008. Analytical solution for Newtonian–Bingham plastic

488

two-phase pressure driven stratified flow through the circular ducts. International

489

Communications in Heat & Mass Transfer. 35, 666-673.

491

M AN U

490

SC

487

Frigaard, I., Nouar, C., 2003. On three-dimensional linear stability of Poiseuille flow of Bingham fluids. Physics of Fluids. 15, 2843-2851.

Gjerstad, K., 2012. Simplified modelling of downhole pressure surges for Bingham fluids in

493

tripping operation. SPE Annual Technical Conference and Exhibition. Society of

494

Petroleum Engineers.

TE D

492

Guo, J.C., Nie, R.S., Jia, Y.L., 2012. Dual permeability flow behavior for modeling horizontal

496

well production in fractured-vuggy carbonate reservoirs, Journal of Hydrology. s464–465,

497

281–293.

499

AC C

498

EP

495

Guo, J.L., 2005. The modern mechanics of fluids flow in oil reservoir, Publisher: Petroleum Industry Press. Beijing City.

500

Kawecki, M.W., 2000. Transient flow to a horizontal water well, Ground water. 38, 842-850.

501

Kong, X.Y., 2010. Advanced mechanics of fluid flow in porous media. Publisher: University

502 503

of Science and Technology of China. Hefei City. Li, Y., Yu, B.M., 2010. Study of the starting pressure gradient in branching network. Science 24

ACCEPTED MANUSCRIPT

505 506 507 508

China Technological Sciences. 9, 2397-2403. Liu, K.F., Mei, C.C., 1990. Approximate equations for the slow spreading of a thin sheet of Bingham plastic fluid. Physics of Fluids A: Fluid Dynamics (1989-1993). 2, 30-36. Liu, R., Liu, Q.S., 2014. Non-modal stability in Hagen-Poiseuille flow of a Bingham fluid.

RI PT

504

Physics of Fluids. 26, 502-543.

Luo, W., Tang, C., Wang, X., 2014. Pressure transient analysis of a horizontal well intercepted

510

by multiple non-planar vertical fractures. Journal of Petroleum Science & Engineering.

511

124, 232–242.

513

M AN U

512

SC

509

Mathias, S.A., Butler, A.P., 2007. Flow to a finite diameter well in a horizontally anisotropic aquifer with wellbore storage. Water Resources Research. 43, 116-130. Martinez, J.A., Escobar, F.H., Montealegre-M, M., 2011. Vertical Well Pressure and Pressure

515

Derivative Analysis for Bingham Fluids in a Homogeneous Reservoirs. Dyna, Year 78,

516

Num.166, p.21-28.

TE D

514

Mirzajanzade, A.K.H., et al., 1971. On the special features of oil and gas field development

518

due to effects of initial pressure gradient. Proc., 8th world Pet.Cong., Special Papers.

519

Elsevier Science Publishers, London.

AC C

EP

517

520

Nie, R.S., Meng, Y.F., Jia, Y.L., Zhang, F.X., Yang, X.T., Niu, X. N., 2012. Dual porosity and

521

dual permeability modeling of horizontal well in naturally fractured reservoir. Transport in

522

Porous Media. 92, 213-235.

523 524 525

Nie, R.S., Meng, Y.F., Guo, J.C., Jia, Y.L., 2012. Modeling transient flow behavior of a horizontal well in a coal seam. International Journal of Coal Geology. 92, 54–68. Oliveira, G.M., Negrdão, C.O.R., Franco, A.T., 2012. Pressure transmission in Bingham fluids 25

ACCEPTED MANUSCRIPT 526

compressed within a closed pipe. Journal of Non-Newtonian Fluid Mechanics. 169–170,

527

121–125.

528

Oosthuizen, R., Al Naqi, A., Al-Anzi, K., Gok, I., Zeybek, M., Cig, K., Al Hashim, H., 2007. Horizontal-well-production logging experience in heavy-oil environment with sand screen:

530

a case study from Kuwait. SPE 105327. 15th SPE Middle East Oil & Gas Show and Conf.,

531

Bahrain, Mar. 11-14.

RI PT

529

Osalusi, E., Side, J., Harris, R., Johnston, B., 2007. On the effectiveness of viscous dissipation

533

and Joule heating on steady MHD flow and heat transfer of a Bingham fluid over a porous

534

rotating disk in the presence of Hall and ion-slip currents ☆ . International

535

Communications in Heat & Mass Transfer. 34, 1030–1040.

539 540 541

M AN U

TE D

538

Reserv. Eval. Eng. 4, 260-269.

Park, E., Zhan, H., 2002. Hydraulics of a finite-diameter horizontal well with wellbore storage and skin effect. Advances in Water Resources. 25, 389–400. Raghavan, R., 2012. A horizontal well in a composite system with planar interfaces. Advances

EP

537

Ozkan, E., 2001. Analysis of horizontal-well responses: contemporary vs. conventional. SPE

in Water Resources. 38, 38–47.

AC C

536

SC

532

542

Samani, N., Kompani-Zare, M., Seyyedian, H., Barry, D.A., 2004. Flow to horizontal and

543

slanted drains in anisotropic unconfined aquifers. Developments in Water Science. 55,

544

427-440.

545 546 547

Staron, L., Lagrée, P.Y., Ray, P., Ray, P., 2013. Scaling laws for the slumping of a Bingham plastic fluid. Journal of Rheology. 57, 1265-1280. Stehfest, H., 1970. Numerical inversion of Laplace transform - algorithm 368. Commun. 26

ACCEPTED MANUSCRIPT

550 551 552 553 554

Sun, D., Zhan, H., 2006. Flow to a horizontal well in an aquitard-aquifer system. Journal of Hydrology. 321, 364-376. Tabatabaei, M., Zhu, D., 2010. Generalized inflow performance relationships for horizontal

RI PT

549

ACM. 13, 47-49.

gas wells. Journal of Natural Gas Science & Engineering. 2, 132-142.

Tichy, J.A., 1991. Hydrodynamic lubrication theory for the Bingham plastic flow model. Journal of Rheology. 35, 477-496.

SC

548

Walicki, E., Walicka, A., Makhaniok, A., 2001. Pressure distribution in a curvilinear thrust

556

bearing with one porous wall lubricated by a Bingham fluid. Meccanica. 36, 709-716.

557

Wang, S., Yu, B., Zheng, Q., Duan, Y., Fang, Q., 2011. A fractal model for the starting

558

pressure gradient for Bingham fluids in porous media embedded with randomly

559

distributed fractal-like tree networks. Advances in Water Resources. 34, 1574–1580.

562 563

TE D

Porous Media. Ph.D. dissertation, U. of California, Berkeley. Wu, Y.S., Pruess, K., Witherspoon, P.A., 1992. Flow and displacement of bingham

EP

561

Wu, Y.S., 1990. Theoretical Studies of Non-Newtonian and Newtonian Fluid Flow Through

non-Newtonian fluids in porous media. SPE Reservoir Engineering. 7, 369-376.

AC C

560

M AN U

555

564

Yun, M., Yu, B., Cai, J., 2008. A fractal model for the starting pressure gradient for Bingham

565

fluids in porous media. International Journal of Heat & Mass Transfer. 51, 1402–1408.

566

Zerzar, A., Bettam, Y., 2004. Interpretation of multiple hydraulically fractured horizontal

567

wells in closed systems. Canadian International Petroleum Conference. Petroleum Society

568

of Canada.

569

Zhan, H., Zlotnik, V., 2002. A. Ground water flow to horizontal and slanted wells in 27

ACCEPTED MANUSCRIPT 570

unconfined aquifers. Water Resource Research. 38, 1108. Zhang, Y., 2003. Error estimates for the numerical approximation of time-dependent flow of

572

Bingham fluid in cylindrical pipes by the regularization method. Numerische Mathematik.

573

96, 153-184.

RI PT

571

Zhu, G., Yao, J., Fan, D., Zeng, H., 2015. Pressure transient analysis of fractured horizontal

575

well in shale gas reservoir. Chinese Journal of Theoretical and Applied Mechanics. 47,

576

945-954.

577

Appendix

578

1. Material balance equation

The material balance equation for fluid flow in porous formations is as follow (Guo 2005;

580

Kong 2010):

581

∇ ⋅ ( ρv ) +

∂ ( ρφ ) =0 ∂t

TE D

579

M AN U

SC

574

(A-1)

where ϕ is porosity, fraction; ρ is fluid density, kg/m3; t is time, s; v is velocity, m/s; q is rate,

583

m3/s.

584

2. State equations

AC C

EP

582

585

Under the high pressure conditions of underground formation, the fluid and rock have a

586

slightly compressible feature, so the porosity of rock and the density of viscoelastic Bingham

587

fluid depend on the formation pressure. The state equations of fluid and rock considering

588

slight compressibility when fluid density and porosity depend exponentially on pressure (Guo

589

2005; Kong 2010) can be given by

28

ACCEPTED MANUSCRIPT Cρ ( p − p0 )

590

ρ = ρ0e

591

φ = φ0eC ( p − p ) f

(A-2)

(A-3)

0

where Cρ is liquid compressibility, Pa-1; Cf is formation rock compressibility, Pa-1; p is

593

formation pressure, Pa; ρ0, φ0 , p0 are arbitrary reference values.

594

3. Motion equation

RI PT

592

For visco-plastic Bingham fluids, the rheological equation of shear stress with shear rate

596

can be represented by (Chevalier et al. 2013; Chevalier et al. 2014; Chevalier et al. 2015; Wu

597

et al. 1992)

M AN U

598

SC

595

τ =τ c + µγ&

(A-4)

599

where τc is yield stress, Pa; γ& is shear rate, s-1; τ is shear stress, Pa; µ is consistency factor

600

(Chevalier et al. 2013), namely, Bingham plastic factor (Wu et al. 1992), Pa s. Based on experimental observations (Chevalier et al. 2013), Bingham fluid can

602

significantly flow in an applied pressure that must exceed a constant pressure loss (a pressure

603

threshold) coming from the fluid yield stress (Chevalier et al. 2015). The non-Newtonian fluid

604

has a more complex relationship between pressure gradient and velocity, the formulation of

605

Darcy's law has been modified for non-Newtonian Bingham fluid flow in porous media (Bear

606

1972; Wu et al. 1992):

EP

AC C

607

TE D

601

k v  v = − µ (∇p − λ ), r v = 0, 

(∇p > λ )

(A-5)

(∇p ≤ λ )

608

v where ∇p is pressure gradient, Pa/m; v is velocity, m/s; µ is viscosity (consistency factor) of

609

Bingham fluid, Pa s; λ is threshold pressure gradient corresponding to the yield stress in a

610

porous medium, Pa/m.

611

Therefore, the motion equation of fluid for r and z coordinates can be expressed by 29

ACCEPTED MANUSCRIPT

612

kh ∂p  vr = − µ ( ∂r − λB ),   v = 0,  r

∂p > λB ) ∂r ∂p ( ≤ λB ) ∂r

(A-6)

613

kp ∂p  ∂p vz = − µ ( ∂z − λz ), ( ∂z > λz )  ∂p  v = 0, ( ≤ λz ) z ∂z 

(A-7)

RI PT

(

where kh is horizontal permeability; kp is perpendicular permeability; z is vertical coordinate,

615

m; r is radial coordinate, m; vz is flow velocity in the vertical plane, m/s; vr is flow velocity in

616

the horizontal plane, m/s; λB is threshold pressure gradient for r coordinate, Pa/m; λz is

617

threshold pressure gradient for z coordinate, Pa/m.

618

4. Governing equation

M AN U

SC

614

619

In radial cylindrical system, the first item in Eq. (A-1) can be written by

620

∇ ⋅ ( ρv) =

(A-8)

TE D

621

∂ ( ρ vz ) 1 ∂ ( ρ rvr ) + r ∂r ∂z

where r is radial coordinate, m; vr is radial velocity, m/s; vz is velocity in z direction, m/s. Substitute Eq. (A-2) and Eq. (A-6) into the first item in the right side of Eq. (A-8):

623

1 ∂ ρ k 1 ∂ Cρ ( p − p0 ) ∂p [re ( − λB )] ( r ρ vr ) = − 0 h r ∂r µ r ∂r ∂r

625

626

AC C

624

EP

622

Cρ ( p − p0 )

ρ k 1 ∂ 1 ∂e =− 0 h [r µ r ∂r Cρ ∂r

− λBre

Cρ ( p − p0 )

]

(A-9)

Using Maclaurin series expansion, the function e x can be converted into ey = 1 + y +

y2 yn +L + +L 2 n!

(A-10)

627

Using Maclaurin series expansion for Eq.(A-9), Eq.(A-9) can be rewritten by

628

2 2 1 ∂ 1 ρ 0k h ∂ 1 ∂[1 + Cρ ( p − p0 ) + Cρ ( p − p0 ) / 2 + L] { } ρ r v = − r ( r) r ∂r r µ ∂r Cρ ∂r

30

ACCEPTED MANUSCRIPT +

629

630

ρ 0k h 1 ∂ {λBr[1 + Cρ ( p − p0 ) + Cρ2 ( p − p0 )2 / 2 + L]} µ r ∂r

(A-11)

Reasonable values of Cρ and (p-p0) are in the orders of 10-10Pa-1 and 107Pa respectively, so the square of Cρ(p-p0) is in the order of 10-6. Accordingly, the value of Cρ(p-p0) is far larger

632

than that of [Cρ(p-p0)]2/2. Through neglecting the second and higher order items in Eq.(A-11),

633

Eq.(A-11) can be simplified to

RI PT

631

ρk 1 ∂ 1 ∂ 1 ∂[1 + Cρ ( p − p0 )] − λBr[1 + Cρ ( p − p0 )]} {r ( r ρ vr ) = − 0 h µ r ∂r Cρ ∂r r ∂r

635

Eq.(A-12) can be changed to

638

M AN U

637

1 ∂ ρ k 1 ∂ ∂p {r − λ r[1 + Cρ ( p − p0 )]} ( r ρ vr ) = − 0 h µ r ∂r ∂r B r ∂r ρ k 1 ∂ ∂p ρ0kh 1 ∂ (r ) + {λ r[1 + Cρ ( p − p0 )]} =− 0 h µ r ∂r ∂r µ r ∂r B

(A-13)

When compared with 1 the term of Cρ(p-p0) can also be neglected (Kong 2010), therefore,

TE D

636

(A-12)

SC

634

Eq.(A-13) can be further simplified to

1 ∂ ρ 0k h ∂ 2 p 1 ∂p λB r ρ v = − [ + − ] ( r) r ∂r µ ∂r 2 r ∂r r

640

Substitute Eq. (A-2) and Eq. (A-7) into the second item in the right side of Eq. (A-8):

641

ρ 0kp ∂ Cρ ( p − p0 ) ∂p ρ 0k p ∂ 1 ∂eCρ ( p − p0 ) ∂ ( ρ vz ) C ( p− p ) =− [e ( − λz )]= − [ − λz e ρ 0 ] ∂z ∂z ∂z µ ∂z µ ∂z Cρ

642

Using Maclaurin series expansion for Eq.(A-15) and ignoring the higher order items,

644

645

646

AC C

643

EP

639

(A-14)

(A-15)

Eq.(A-15) can be rewritten by

ρ k ∂ 1 ∂[1 + Cρ ( p − p0 )] ∂ ( ρ vz ) =− 0 p { − λz [1 + Cρ ( p − p0 )]} ∂z ∂z µ ∂z Cρ =−

ρ0kp ∂ ∂p { − λ [1 + Cρ ( p − p0 )]} µ ∂z ∂z z

(A-16)

When compared with 1 the term of Cρ(p-p0) can also be neglected (Kong 2010), therefore, 31

ACCEPTED MANUSCRIPT Eq.(A-16) can be further simplified to

648

ρ0kp ∂ ∂p ρ0kp ∂ 2 p ∂ ( ρ vz ) =− ( −λ ) = − ∂z µ ∂z ∂z z µ ∂z 2

649

Substitute Eq. (A-14) and Eq. (A-17) into Eq. (A-8):

653 654

655 656

657

(A-18)

Substitute Eq. (A-2) and Eq. (A-3) into the second item in Eq. (A-1):

∂ ∂ ∂ C ( p− p ) ( C +C )( p − p0 ) ( ρφ ) = ( ρ0φ0e ρ 0 eCf ( p − p0 ) ) = [ ρ 0φ0e ρ f ] ∂t ∂t ∂t

SC

652

ρ 0k h ∂ 2 p 1 ∂p λB ρ 0k p ∂ 2 p [ + − ]− µ ∂r 2 r ∂r r µ ∂z 2

(A-19)

Using Maclaurin series expansion for Eq.(A-19) and ignoring the higher order items, Eq.(A-19) can be rewritten by

M AN U

651

∇ ⋅ ( ρv) = −

∂ ∂ ∂p ( ρφ ) = ρ 0φ0 [1+(Cρ + Cf )( p − p0 )] = ρ 0φ0Ct ∂t ∂t ∂t

(A-20)

where Ct is total compressibility, Pa-1.

Ct = Cf + Cρ

TE D

650

(A-17)

RI PT

647

(A-21)

Combine Eq. (A-18) , Eq. (A-20) and Eq. (A-1):

659

∂p ρ 0kh ∂ 2 p 1 ∂p λB ρ 0kp ∂ 2 p − [ 2 + − ]− + ρ 0φ0Ct =0 2 r ∂r r ∂t µ ∂r µ ∂z

660

φ0 is usually replaced by φ for it is an arbitrary reference value. Thus Eq. (A-22) can be

AC C

661

EP

658

(A-22)

changed to

662

∂ 2 p 1 ∂p λB kp ∂ 2 p φµCt ∂p + − + = ∂r 2 r ∂r r kh ∂z 2 k h ∂t

663

Eq. (A-23) is the governing differential equation of Non-Newtonian Bingham fluid. The

664

three parameters, µ, λB and Ct, represent viscous, plastic and slightly elastic property of

665

Bingham fluid in underground porous formations respectively.

(A-23)

666 32

ACCEPTED MANUSCRIPT List of figures:

667

Fig. 1. Horizontal well scheme in an underground formation.

669

Fig. 2. The views of coordinate system (Nie et al., 2012).

670

Fig. 3. Log-log curves of pwD and p wD ′ for flow regime recognition. The curves were simulated under a

671

group of model parameters as the legends in the figure. An infinite boundary, a closed top boundary and a

672

closed bottom boundary were considered in the model.

673

Fig. 4. Schemes of main flow regimes

674

Fig. 5. Log-log curves of p wD and p wD ′ controlled by different side boundary conditions. Closed top

675

and bottom boundaries were considered in the model. The curves were simulated under a group of model

676

parameters as the legends in the figure. Curves ①, ② and ③ represented the infinite boundary, the

677

constant pressure boundary and the closed boundary respectively.

678

′ under constant pressure boundary of bottom. Constant pressure Fig. 6. Log-log curves of p wD and p wD

679

boundary and closed boundary of top were considered in the model respectively. Under a group of model

680

parameters as the legends in the figure, the curves were simulated. The closed boundary and the constant

681

pressure boundary were respectively represented by Curves ① and ②.

682

Fig. 7. Log-log curves of p wD and p wD ′ affected by dimensionless threshold pressure gradient. Under a

683

group of fixed model parameters as the legends in the figure, the curves were simulated by setting

684

dimensionless threshold pressure gradient as 10-3, 10-4, 10-5 and 0. Curves ①, ② and ③ represented the

685

curves of Bingham fluid and curves ④ represented the curves of Newtonian fluid. The figure

686

demonstrated the curves sensitivity of p w D a n d p w′ D to the dimensionless threshold pressure gradient.

687

Fig. 8. Log-log curves of p w D a n d p w′ D affected by horizontal wellbore length. Under a group of fixed

688

model parameters as the legends in the figure, the curves were simulated by setting the horizontal wellbore

689

length as 300m, 500m and 700m. The figure demonstrated the curves sensitivity of p w D an d p w′ D to

690

the horizontal wellbore length.

691

Fig. 9. Log-log curves of p w D a n d p w′ D affected by skin factor. Under a group of fixed model

692

parameters as the legends in the figure, the curves were simulated by setting the skin factor as 0, -0.5 and -1.

693

The figure demonstrated the curves sensitivity of p w D a n d p w′ D to the skin factor.

AC C

EP

TE D

M AN U

SC

RI PT

668

33

ACCEPTED MANUSCRIPT Fig. 10. Log-log curves of p w D a n d p w′ D affected by the ratio of horizontal and vertical permeability.

695

Under a group of fixed model parameters as the legends in the figure, the curves were simulated by setting

696

the ratio of horizontal and vertical permeability as 50, 20 and 10. The figure demonstrated the curves

697

sensitivity of p w D a n d p w′ D to the ratio of horizontal and vertical permeability.

AC C

EP

TE D

M AN U

SC

RI PT

694

34

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT Research highlights A new porous flow model of Bingham fluid for horizontal well production was originally investigated. The viscous, plastic and slightly elastic properties of Bingham fluid in a heavy oil reservoir were considered in the model. The hydrodynamic behavior of Bingham fluid was modeled and the typical flow characteristics

RI PT

were thoroughly analyzed.

A series of pressure dynamic curves controlled by model parameters were created and studied.

AC C

EP

TE D

M AN U

SC

The influence of the threshold pressure gradient on the Bingham flow traits was deeply discussed.