Fluid Ovna,,us, t61 . 30. No. l, 1995
SOLUTION OF SOME PROBLEMS OF FLOW THROUGH FRACTURED POROUS MEDIA R. Sh . Mardanov, F M . Mukhametzyanov, A. G . Fatykhov
UDC 532 .546
Problems of fluid flow through a fractured porous medium consisting of fractures and blocks with different filter characteristics are solved . The mass exchange between fractures and blocks is assumed to be proportional to the pressure difference between them . The porosity in the fractures is assumed to be negligibly small . Under these assumptions the determination of the pressure fields reduces to the integration of a system of linear differential equations . The solution is found by the operational method using the Efros theorem . The cases of oil reservoir operation by means of both galleries and wells are considered . The solutions are obtained in an analytical form convenient for calculations .
1 . In [1] it was shown that the distribution of the pressuresp l and p 2 in the fractures and blocks of fractured media, respectively, can be determined from the following system of linear differential equations describing fluid flow through porous fractures and the mass exchange between the latter and the blocks : 2 -
¢
t + A(p2 -
ak2 A= N 12m1A ( F'22 + ¢ .)~
p1)
xApl
=0,
-
A(p2
kt x= µmzo(¢22 +
p .) '
-
p1) =0
¢= ¢t2
¢ 22
+ ¢.
where A is the Laplace operator, and ¢, is the compressibility coefficient of the fluid . With allowance for the dependence of the block porosity m 2 on the pressures p 1 andp 2, the constant quantities m m , and are related by the formula ¢21, ate = m2ol ¢Zt
dt +
¢22
atl )
The mass exchange between the fractures and the blocks is characterized by the amount of fluid leaking per unit time t per unit volume q q=a
pk2 p 2
µ
-p
1
IZ
where µ is the viscosity of the fluid, p is its density, k 1 is the permeability of the fractures, k2 is the permeability of the blocks, 1 is the average dimension of the latter, and a is a dimensionless constant associated with the geometry of the medium . The fracture voidage m 1 (ratio of the fracture volume to the total volume of the medium) is assumed to be negligibly small . A feature of the solution of the system (1.1) is the fact that the initial values of pl and p2 cannot be arbitrary assigned [1] . Given the initial value of p 1 it is necessary to find the initial value of the reduced pressure p=p 2 - ¢p 1 from the condition xApt - (1 - ¢)Ap 1 =Ap(x, y, 0) and then to find the initial value of p 2 as follows p2(x, y, 0) =p(x, y, 0) + ¢pt(x, y, 0) From the system (11) we can readily eliminate p 2 , expressing it in terms of the second equation [1] :
Kazan' . Translated from Izvestiya Rossiiskoi Akademii Nauk, Mekhanika Zhidkosti i Gaza, No . 1, pp . 94102, January-February, 1995 . Original article submitted October 14, 1993 . 0015-4628/95/3001-0077$12 .50 ® 1995 Plenum Publishing Corporation
77
~,,
c7t
x
c7
t'1 -
3t
71 12 x _ rl = A(1 - 13) ak2 (1 -
~P
P1_
A
1' 1 - 13
p)
When a << 1 the quantity a can be neglected as compared with unity . As i -. 0 equation (L2) goes over into the usual equation of the elastic regime with the piezoconductivity coefficient x/(1 - a) . On the basis of this model vee will consider the one-dimensional flow between two galleriesy=0 and y=L where the pressure p 1 (y, t) is equal to p 11 and p 12 , respectively. The initial condition is p 1 (y, 0)=0 . In this case Eq .(L2) takes the form: a
"p1
32P1
2p1 a
Applying to Eq .(L2) the Laplace transform P 1 (a, o) = fe -otp1(y,
t)
dt=L(p 1 )
0
in place of (L3) for p1 we obtain : d2p1 _ wp1, dy
_
a x + qa
The boundary conditions for pl can be written in the form : P11
PjIy=0-
U
P12
- P11 ,
P1Iy=L=
Q
-P12
The general solution of (1S) is given by the expression p1 =A 1e~y + B1 e
- may
where A 1 and B 1 are constants to be found from (L6) . Finally, for p l we obtain -1 P1= [ff, sink /jy + P11 sink V Y' (L - Y))( h f L) We now represent (L7) in the form : p1 =F(r1 Vl)~(a) -1 F(a) = P12 sink
y + p 11 sink
° (L - y) a sink
(1 .8)
L
`
x + ii a
In order to fmd the inverse transform of (L7), initially we consider the function -1 F(a) = P12 sink
a Y + P11 sink
° (L - y) a sink
aL
(1 .9)
The inverse transform f(r) of this transform can readily be found . In fact, in accordance with the inversion theorem
78
Y * ic.)
fit) = 1 lim f
2ni
e ° `F(a) da
where y > y o > 0, and determines the domain of analyticity of F(a) in the half-plane to the right of the straight line
Rea Z yo , so that all singular points of F(a) must lie to the left of this straight line . At the same time, F(a) must tend uniformly to zero as ~a~ - ~ . The function (1 .9) satisfies these conditions and has simple poles at the points a=0, i a/rl L=nn (n E N) or a=0, a"=n'n 2 /L 2 (n=1, 2) . Introducing into consideration the closed contour consisting of the interval of integration (y - ic., y - icc) and a semicircle t of radius R 1 =nZ(n + 1/2) 211/L 2 , so that the pole does not lie on I', and taking into account the fact that the integral f e°`F(a) do along this closed contour is equal to the sum of the residues of F(a) with respect to the above simple poles multiplied by e ° ", , we obtain : Y .im
f e°LF(a)da= fe°TF(a)da + ~ResF(a)e ° ' t Then, letting
R1
(or n) tend to infinity and taking into consideration the fact that in this case
[2]
fe°TF(a)da - 0 r we find (generalized Vashchenko-Zakharchenko theorem) : fir)
=E Res
F(a)
e ° "`
In relation to (1.9) this formula gives : fir) =p 12
p
y
L
2(-1) "'1 sin nny/L Q -(xn/L)2T1 !< n_1 n
_ 1
11~(1 - y) - 1E2(_l)n~1sinan(1 L n n=1 n
-
+
(1 .10)
y/L) e -cxnm2 Tn
The series entering into (110) converge uniformly and rapidly . Then, in accordance with the Efros theorem [2], we find that the inverse transform p 1 of the transform (L8) is given by the expression p 1 = ffit)Y(t, t)
dt
0
in which
¶(t,
t) is the solution of the integral equation e-zq(°>~(a) =fe-°T(t, t) dt, 0
x + rl a
or, in this case, e-`exTm(°(a)= fe - °`~( t, t) dt
(1 .12)
0
From (L12) it is obvious that the function e(0)/'lc(a) is the transform of the original function It is well known [ 2] that the original function of the transform 1/aex`/^° is
e-°` Y(t,
t) .
79
where to is a Bessel function of zeroth order depending on an imaginary argument . In accordance with the delay theorem [2], the expression e-"4'Io2 xtt
is the inverse transform of Thus
e°> 1" .
~(t, t)=e -(nI71 Io 2,
(1 .13)
t
Substituting this expression in (111) and taking into account (110), we fmd the pressure distribution in the fractures : P1(Y, t)=e - "~
if
[P1Z{2
-
(_1)n+ 1 sin
? n i=1
anYIL e - (nn/L)Zns}
+
nn
(1 .14) p 11 {(1
-
y)
-
L
?~(-1)n+1smnn(1 - YIL) e_(sn/L)=rq n n .1 n
l0
2
xtt =
dt
~
Expanding to in a series and taking into account the fact that
E (k!)2(11) xit
dt i10(2 xtt rl _,
k
(1 .15)
t ke - T dt
=E
we can write formula (114) in the form convenient for calculations : P 1 (y,
t) =pu
+ p
L
11 (1 - ~) + e'11 x (1 .16)
E
(-l)n
P12 sm anyIL + P11 sm itn(1 -y/L) xtt
k
n[(x/L) 2n 2ti + 1]
n=1 k0
Denoting the unknown function for nonzero initial values by means of pi and settingp 1 =p ; -p 10 , we arrive at a zero initial value, and in the final formula (1 IS) the numbersp 21 -p u, andp12 -p10 will appear 1n place of p01 andp12. Knowing p 1(y, t), from the second equation of system (U) we also find the pressure distribution in the blocks p z (y, t) . 2. In the case of circular galleries the problem is similarly solved . Here, the pressure distribution in the fractures p 1 must satisfy the equation: a?1 a l a aP1 l a apt - rl- --r- =x--r-at at r ar ar r ar ar
Given constant pressures on the circles
r=r and r=R, the boundary conditions are P1I r=r~ = Pc '
P1I r=R
- Pk
where r is the radius of the well . Using the Laplace transform (L7), for zero initial conditions in place of (21) and (2 .2) we obtain the relations :
80
ld(dPll
r dr P1 I r=r~
= PaC
r
dr
= Pc,
Vfp1 =0
Pl
Ir=R
Pk
a
= Pk
The general solution of (2.3) and (2 .4) can be written in the form : p1
=A 1I0(y"jr) + B1Ko(/r)
where to and Ko are zeroth-order Bessel functions of the first and second kind . Satisfying boundary conditions (2.4), we find : A 1 =tP~Ko(~R) B1
={pklo(Jr~)
-
A =Io(~r~Ko(~R) - lo(~R)Ko(~r~ Thus, for p, consistent with (2.4) we obtain the expression P1= jp~[lo(~r)Ko(~R) - I0(j R)K0 (Ir)+
Pk[IO(VrC)KO(Jr) - I0(f r)K0(/r)]}A-i Again, as in linear galleries, in order to find the inverse transform of (2.5) we consider the transform F(a) =
1
a
fp~[lo(ar)K0(vK) - lo(aR)K0 (ar)] +
p k[I0 (ar~)K0 (ar)
- I0(ar)K0(ar)]}A -1
A =1o(ar)Ko (aR) - I0 (aR)K0(ar~),
a =
(2 .6) v ti 0
Proceeding as before in (19), we find that the inverse transform of the transform (2 .6) is the function [3] : fl) =p~
+ pk
v
law
+ 2~ [p J0 (a~u) - p .J0(a„)] x
[J0(a nv)Y0(a,) - J0 (a,)Y0(a v)]Jo(anw)[Jo(a,,w) - J(a)1 l exp(-a~r
where Jo and to are zeroth-order Bessel functions of the first and second kind, respectively, and a, are roots of the transcendental equation [3] J0(a)Y0(a 1w) - J0 (a,~w)Y0 (a) =0 In accordance with the same Efros theorem and (113), for the required inverse transform of the transform (2 .5) describing the pressure distribution in the fractures in the ring in question for a zero initial value we obtain the expression: + + 2E {p~J0(a~u) - peJo (a,,)} x pkw ~w J lnw e •1 o [J0 (aw)Yo(a„) - Jo(a,)Yo(a,,v)]Jo(a 1w)exp(-( r4jr~rt+1)t)[Jo(a~w) - JI (2%/xrt/rl)Jdt
P1(r>r) = e
xdn (~
81
or, simplifying this expression as in (114), using (115) and (116), we finally arrive at the following formula convenient for calculations : P1(r,
t)
=P
`lnw
+
lnv Pk
lnw
- 2e"
a^E Jo(a„u) „_ 1
[P J0(a„w) -
pkJ0(a„)]
x k
[Jo(aw)Yo(an) - J0(a„)Y0(a„v)][Jo(a~w) -
Jo(a„)] -'
E [1 + aRr~(r~ kk!] r1 1
( x t)
k-I
As before, for an initial value equal to a constant quantity p 10 , in formula (2 .7) this quantity is added additively and the quantitiesp~ -P10 andpk -P10 appear in place of p~ andpk , because the replacement of the unknown function p,' by P, =P1 -P10 leads to a zero initial value . 3 . The problem is similarly solved in the case of systems of wells . For example, if a well operates at the center of a circular reservoir of radius R bounded by a supply contour on which a constant pressure Pk is maintained, i.e ., if at the center there is a point sink or source with a constant flow rate Q, then for zero initial pressure the boundary conditions for Eq .(21) will be, respectively, P 1 (R,
aPl = µQ =Q . r ar r=O 2nk1 h
t) =p Ik ,
or in Laplace transforms Pl(R, t) = P 'k a
- ~'=~O =Q . (r)
Here, the solution is the Laplace transform of Eq.(21) and conditions (2.2), written in the cylindrical coordinate system (r, ~). We seek this solution in the form : pl =A I I0 (t/jr) + B 1Ko(Ir)
where Io and Ka are wroth-order modified Bessel functions of the first and second kind . It is known that for r=0 K°(4rr)" 2 has a logarithmic singularity and, consequently, takes into account the presence of a source (sink) at the pole . Determining the constants A 1 and B from conditions (31) and (3 .2), we find
a I0(jjR)
( I R)I°( Q~ I°(~R)K°(/ r) a I0(R)
)
(3.3)
To calculate the inverse transform (3 .3) initially we determine the inverse transform of the expression - p lo(ar) + Q, lo(aR)Ko(ar) - lo(ar)K0(aR) F(a)
a 4(aR)
(3 .4)
J0(aR)
o
Using the inversion theorem [2], all of whose conditions of applicability are satisfied in the case in question, and choosing the same contour of integration as before with allowance for the fact that the integrand (3.4) has the simple poles a=0, a„= iv 2/R4z (n=1, 2, » .), we find J0(a u) J°(an ~e E E Ik 2Plk J° ,(a ~ e + Q In u + 2Q~ =1 1(T) =P E E a i (a„) R2
where a„ are roots of the equation 82
r _R
J0(an)=0, n=1, 2 . . . In fording the residues with respect to a„ we used both the identities 10(z)K0 (z) - K° (z)1° ' (z) =1 / z and the fact that 1o(iz) =J°(z)1°'(tz) _ -1J0 1 (z) Representing (3 .3), as for galleries, in the form : p1 =F(o(a)) I (a)
where w(a) and ~(a) are the notation adopted previously and again using the Efros theorem as well as making the same simplifications as before, with allowance for the fact that
fo e -r ° r k dt = k!Y~ we finally obtain:
E
~ J°(ar/R) 1 ( xtl k + n - - 2p1ke-"~^ P1(r, t)=p er + Q 1n! a,,J1 (a n) k-° ( 1 + anti/R2ykk1 r) r n-1 2Q~e _, 4
~~
1
JO(anr/R)
xt k
a i(a n) k-° (1 + a~T)/R2)kk!
n-1
rl
If in a circular reservoir the well is arranged off center at a point with coordinates (p„ c ;), then, as follows from (L2), p1 must satisfy the equation T(r
p1 + 1 a 1J r 2 a p2
I
-
V1P1
=0
We seek the solution of this equation in the form : p 1 = Q'Ko(%/iR') + EA,,1~(Vr) cos (c
-
cp,)
R' = Jr2 - 2ps cos(cp - (p,) + Pi where I„ (x) is a modified Bessel function of n-th order and A„ are constants to be determined from condition (31) . Satisfying this condition, we use the well-known relation for r > K0(fjR') =I0(~P,)Ko(*,,/ijr) + 2~ I / p,)K~(~r) cos (p n=1
where I„ is a modified Bessel function of n-th order. We find Ao
=Lp°Y A =-2Q •1 (%P .)K (~R)In1(.,/jR) a
Thus, for p1 we obtain the expression
83
. I°(/r) P1= p
I°(jp i ) + Q• I°(JR)K°(Jr) - I°(vr)K(JR)
2Q ' E [If (~f KA(%r) - I,,(~r)K"(sR)] A( cos( cp
°
-
I
i=1
Again with the intention of using the Efros theorem, we first calculate the inverse transform f(r) of the function F(a) I°(ar)
F(a) = P a I° (aR)
+ Q• [I°(aR)K°(ar) - I°(ar)K°(aR)] I°(a a
p')
+
I°(aR)
2Q ' __[In(vI~K"(ar) - I(r)K(R)]
j"~o
cos (gyp - cp s)
n
In so doing, proceeding as in the case of a central well, we find : J°(a"r~ E .1(t) =P1t - 2p ~ et + Q,In u + 2Q ~~ ,E e + n=1 anl1(an) n=1 a I2(a n) " J° (akp~R}I° (aku) t 4Q,~ a cos (tp n-I t-I (ak)2J,2„1(at)
( -
p,)
where a n and ak are the known roots of the equations : J°(a n)=0,
J° (ak) _O
Now, knowing f(r) and having ~P(r, t) in (L13), and taking into account (3 .5) and the expansion of to in its Efros theorem power series, we can readily obtain the following expression for the unknown pressure field p t (r,tp,t) due to the operation of an off-center well in a circular reservoir :
J0(a nr/R) 1 (xt\ t 1 2pue xanE ~- + P1(r, ~, t)=pu "=1 anJ1(an) t~ k! 1~ (1 + tt T)/R2)t
r 4Q, e-
J°(a",p/R) '
+ 2Q a-=an
Q In.
ni
a12( a")
E E J°(at prr)I°(atr/R) r=1
t~
1
Ixt~t 1
+
k! r~ (1 + a/R 2 )t
cos n(p -
(at')2J2 (at•)
tp,)E- xt~t 1 - 2)'° m~
I
(1 - an r~/R
REFERENCES 1.
2. 3.
84
G. I . Barenblatt, V M. Ento and V M . Ryzhyk, Motion of Fluids and Gases in Porous Strata [in Russian], Nedra, Moscow (1984) . A . V Lyko The Theory of Heat Conduction [in Russian], Gostoptekhizdat, Moscow (1952) . A. Gray and G . B . Mathews, A Treatise on Bessel Functions and Their Applications to Physics, Macmillan, London (1922).