- Email: [email protected]
Nonsteady flow of non-Newtonian fluids through a porous medium
In/. 1 Engng Sci Vol. 21, No. 3, PP. W-210, Printed in Great Britain.
OOZO-7225/83/030199-l2W3.00/0 @ 1983 Pergarnon Press Ltd.
1983
NONSTEADY FLOW OF NON-NE~~NIAN THROUGH A POROUS TEDIUM
FLUIDS
H. PASCAL Department of Physics, Unkersity of Alberta, Edmonton, Alberta, Canada T6G 2Jl Abstract-In this paper, the nonsteady Row of power-law fluids with a yield stress through a porous medium is investigated. In order to illustrate the rheological behaviour effects on pressure and flow rate distributions, a flow system of practical interest was analysed. The approximate solutions in a closed form were obtained by using the integral method. For a short time, these have been formulated in terms of a moving boundary problem. The flow deviation from Newtonian behaviour was emphasized by means of several dimensionless groups. These have been found to be relevant in evaluating the rheological effects on steady and nonsteady flow behaviour through a porous medium.
INTRODUCTION
THETHEORIESused now in oil reservoir engineering to analyse the flow through a porous medium and to predict the reservoir performance are based on the assumption that crude oil is a Newtonian fluid and as a result Darcy’s law holds. While for light crudes, which are essentially Newtonian fluids, this assumption is justified, this is not so for a large class of complex fluids produced currently from certain oil fields. The heavy crudes, such as occur in Canada and Venezuela, are characterized by completely non-Newtonian behaviour as has been repeatedly demonstrated and reported in laboratory studies. Another illustrative example of non-Newtonian fluids are waxy crudes located in reservoirs of shallow depth. The reported results from the laboratory and field data have shown that heavy exudes at reservoir temperature display a rheological behaviour of power law with a yield stress. As a result of the presence of a gel structure, the shear stress-shear rate relationship was found to be a curve which does not pass through the origin. An example to support the experimental evidence of this behaviour was presented in papers [1, 21. Some heavy crudes show a reduction in apparent viscosity as the shear rate increases, others show an increase in apparent viscosity. From a rheological point of view these fluids are of pseudo-plastic or dilatant type. It should be mentioned that this behaviour has been observed above the yield stress, i.e. when this value was exceeded and presumably any structure that the fluid has is broken down. While for light crudes the viscosity is a parameter depending on temperature but not on conditions of measurement, the heavy crude viscosity depends strongly on shear rate. In order to fit the experimental data into a large range of shear rate variations, an analytical expression of the curves of rheological behaviour is required in solving the problems associated with flow trough porous medium. As a result, the knowledge of values of the apparent viscosity corresponding to conditions + + 0 and 9 + ~0is important. Most of the frequently used empirical relations to express analytically the apparent viscosity in terms of shear rate are various forms of the power law. For example, the most widely used for a fluid with a structure at zero shear rate is the empirical equation 7 =
H(j)” + 70;
7 >
70.
(0
By using this equation the apparent viscosity will be q,(i)
=
&y-’
+:
(2)
where n, H afid T,,are the rheological parameters to be determined from the rheological tests. These parameters depend strongly on temperature. Useful empiricism for temperature depen199
200
H.PASCAL
dence is presented in the literature as H(T)
= H(T,)
TV = ho
e-a(T-Td/Ta
(3)
+ NT - ToYTo
(4)
and n(T) = n(T,)+ c(T - T,)/T,
(5)
in which To is the reference temperature. The notations used are shown in the nomenclature. For n = 1, relation (1) reduces to the well known Bingham’s equation. This model was used in papers[3-51 to investigate the yield stress effect, i.e. threshold gradient, on the flow behaviour through a porous medium in several problems of fluid mechanics. Obviously the rheological model described by eqn (1) reflects the two most important aspects of the behaviour of heavy and waxy crudes observed from the experimental data; the existence of a yield stress and non-Newtonian flow curve above the yield stress. For example, at T > T,,the viscous flow becomes possible and its rheological behaviour may be described by eqn (1). It should be noted that at moderate and high shear rates the yield stress effect on apparent viscosity is insignificant. However, at low shear rate this effect becomes significant and cannot be neglected in eqn (2). In the production of heavy crudes, the shear rates which currently appear in flow through porous medium are expected to be in the range 0.1-1.0 set-‘. Consequently lower shear rate results are of practical interest, in which case a better rheological description should be the generalized Bingham’s model given by eqn (1). For example, it was found that, in general, the rheological models which included a yield stress better described the laboratory data for fluids with a gel structure than did the power law model. Another class of non-Newtonian fluids are polymer solutions, micro-macro emulsions. These fluids are often injected in oil reservoirs in order to improve the recovery efficiency. However this class of non-Newtonian fluids may be described by a power law equation without having to consider the yield stress. From the discussion presented above, it is evident that an understanding of non-Newtonian flow behaviour through porous medium is essential for practical purposes. The main objective of this investigation is to show the flow deviation from Newtonian behaviour expressed in terms of rheological parameters, using several examples of practical interest in engineering science. BASIC EQUATIONS
To consider the rheological effect in the flow description through a porous medium a modified Darcy’s law is required. For the rheological model given by eqn (1) the modtied Darcy’s law may be written as - & (grad p - a,,)]“”
v’= [
(6)
which is valid provided that lgrad pJ > cue;v # 0 [grad p[
(7)
Equation (6) was obtained by applying Kozeny’s approach to a capillary model of a power-law fluid in the presence of a yield stress. Using this approach, the following expression for pef was given in[6, 71 /.Lef= ; (9 + 3/n)“(lSOkf#-““2
(8)
Nonsteady flow of non-Newtonian fluids through a porous medium
201
Recently, the term k/per in (6) was found in[8] as (9) The threshold gradient in (6) may be expressed in terms of yield stress[3]
(10) where p is a dimensionless constant to be determined experimentally. As a result of the presence of a structure in heavy crudes which resist at a mechanical deformation, the modified Darcy’s law (6) implies a pressure gradient which increases or decreases monotonically to a given value. If there is successive increasing and decreasing of the pressure gradient it should be expected there will be a hysteresis loop effect, so that in this situation eqn (6) could not be valid. The material balance equation for a plan radial flow is f
& (Ru) =-(C&J + c@.
(11)
From equations (6) and (11) one obtains
$&[R($] =a2$
(12)
in which function x(R, t) is related to pressure by transformation
x(R t) = P@, t) - w,R.
(13)
For one dimensional flow, eqn (12) becomes
--a
ax ( > =aat
ax ax ax
I/”
2
(14)
where constant a2 in (12) and (14) is
a*= y
()
““(C&t
Cp).
(15)
Basic eqns (12) and (14) give an infinite velocity of wave front or pressure pulse propagation through a fluid flowing through a porous medium. To obtain a finite velocity the inertial term au/at should be considered in deriving the basic eqns (12) and (14). Although the infiniteness of the pressure wave velocity cannot be accepted on physical grounds, because it implies that the pressure perturbation generated in pore fluid by a wave front is felt instantaneously at any distance in an oil reservoir, the effect of finite velocity correction is extremely small and will be negligible in most problems of practical interest. Therefore, the unsteady flow through a porous medium may be described by using the parabolic eqns (12) and (14). The steady state solutions may be derived from (12) and (6) and expressed as P(R)=P,+~,(R-R,)+(~~~)~
(
$-
“(RI-” _ Rf;“) >
(16)
and me- pw -
m(Re- RJ
R;-“-R;-”
I
l/n
.
(17)
202
H.PASCAL
In practice we currently have RJR, = 10m3,so that in (17) and (16) one can use the following relation R w
l-n
(R >
-l+(l-n)ln+.
e
(18)
As a result, one obtains from (17)
(19)
FORMULATION
OFPROBLEMSANDSOLUTIONS
To illustrate the rheological behaviour effect on the unsteady flow through a porous medium, a classical problem of oil reservoir engineering will be considered. The flow system is a closed reservoir in which the porous medium is homogeneous. The solutions of practical interest are those corresponding to a centrally located well in the reservoir and producing at a constant pressure or constant flow rate. It is assumed that at initial moment t = 0, the pressure in the system is constant, say pe, and for t > 0, the system is depleted by a constant pressure of production, say pW As a result, the flow is initiated and maintained in the porous medium provided that pW< pe. A plan radial transient flow response will appear at t > 0 which may be analyzed by means of basic eqn (12) associated with certain boundary and initial conditions. Basic equations are nonlinear. Consequently, an approximate solution in closed form could be obtained by means of the integral method. Using this approach we have a moving boundary problem. This situation exists until the decompression front reaches the boundary of the reservoir located at distance R,. As a result, the appropriate boundary and initial conditions in this case will be POL, 0 = pw; 8P s R=R(t)=lYo; I PVW,
tl=
me;
R(O)=R,;
0 < t < 6~ O
(20)
O< t
t=O
where tM is the time required for the front to reach the closed frontier and R(t) represents the location of this front. The conditions (20) on function x(R, t), given by relation (13), become
xUL, 0 = pw - ~loRw ax 0 aR IR=R@)=
x[RW,
tl = me-
; O
(21)
4W)
R(0) = R,. The integral method has been used in our previous paper [3], in order to analyze an identical problem of one-dimensional flow for the particular case n = 1 and 7. # 0. This case corresponds to a Bingham fluid. Although no error analysis was made there, a comparison between the results obtained from numerical and exact solutions has shown a good agreement with those obtained from the integral method. The main advantage of the integral method appears to be the possibility of determining the pressure and flow rate-distributions without extensive computations and at a reasonable accuracy for engineering purposes. Using this approach, useful information on the flow behaviour has been obtained in[3]. However, it should be emphasized that the assumed pressure distributions must be selected with great care. For example, taking higher-order polynomial forms does not necessarily give better accuracy.
Nonsteadyflowof non-Newtonianfluidsthrougha porous medium
203
Taking into account the boundary conditions (21), the basic eqn (12) may be written as
a*
Rx(R, 0 dR - -_[P, - dW)l~ 2 A convenient consider
dR2
(22)
form of pressure distribution in the domain R, < R(t) < R,, for 0 < t < fM, is to
R x(R,t)=A+B~(t)+c~~(t)i
R2
O
(23)
Using the conditions (21) one obtains x(R, t) = me- LU,W) + C(t)( 1 - 2 &
(24)
+ &)
in which
C(f) =
W
(25)
The front location R(t) may be determined from (22) and (24) and expressed in terms of dimensionless time as (1 t $u)d’+“)‘” du (1 - RU)“n
T=
(26)
where
t(T)==;
T=
(27)
R,
From conditions (20) one obtains
(28) and one can observe that function (24) satisfies this equation. The particular case corresponding to a fluid of power law only, i.e. Cl= 0, may be determined from (26) with n # 1 and R = 0. In this situation, the front location becomes y
= [‘l t;n’T]“‘1+2”
(29)
and when n = 1 we have
( >
R(T) 3z3T R, ’
(30)
Once R(T) is known from (26), the pressure distribution in function of time and distance may be obtained from (24), (25) and (26). The time TM may be also determined from (26) with t(T) = 1. However, it should be mentioned that in deriving eqn (26), there was assumed 1 - RJR(t) = 1. This assumption is quite justified provided that RJR(t) CO.1; after a very short time of DES Vol. 21, No. 3-B
204
H.PASCAL
production the front location will reach such a distance in which condition RJR(t) < 0.1 is satisfied. Rather than to determine the pressure distribution in the system, the practical interest is to predict the flow rate behaviour in the well. For this purpose, eqns (6) and (24) should be used. For example, from these equations, the dimensionless flow rate in the well may be obtained and expressed as
Q*UL T) = [
1 - Q(7) 5(T)
1
I’”
(31)
with
Q*UL ‘J’)=
T) ($Y$ 2TR,h )I’”QoUL
(32)
and where t(T) one obtains from (26). To analyze the rheological effect on flow behaviour for t > t,, the following equation must be used _R
Rx(R, t) dR;
t > tM.
(33)
Obviously, the moving boundary problem has been eliminated, but the appropriate boundary and initial conditions in this case will be on function x(R, t)
xUL, t) = pw - aoh;
I
-dx aR
=o t> =
t&f
t>
tM
R=R,
xUL
t>
p,(t)
-
aoR,;
t >
(34)
ty
in which pe(tM) = pe = constant. As the numerical solution of the basic eqn (12) indicated, an appropriate form of the function x(R, t) for t > tM, associated with conditions (34), should be expressed as x(R, t) = A In F + RR + CR2; t > tM
(35)
w
which for conditions (34) becomes R(R - R,) t RJR, - 2R,) In $
x(R, t) = (pW- (Y,,R~)
1
W
(36)
where C(r)=n.(f)-
aoRe-
~(~~--a~R,)(1-ln~)
(37) RJR, - R,) t RJR, - 2R,) In +
W
Equation (37) contains the unknown function p,(t), i.e. pressure variation at closed frontier. This may be determined from (33) and (36). It should be noted that prior to the front reaching the closed frontier we have p(R,, t) = pe = constant, for t < tM, and p(R,, t) = p,(t) for t > tM, where pe = constant represents the pressure value in the system corresponding at initial moment t = 0. Using the condition pe(tM) = pe = constant, function p,(t) for t > tM may be obtained from
Nonsteady flow of non-Newtonian fluids through a porous medium
205
(33) and (36) and expressed as
2RtWf YOR,_ 2R:CM
2(n - l)R,R:(t - t,,,) na*~,(2RZ,C(t,) t y,,R,$-I)‘” I
+ Y~R, -
n/n-l (38)
where C(t) is given by relation (37) and yo =
(P, -
aoRd
3W
(39)
(40) In practice we currently have RJR, = 10m3,so that assumption 1 -RJR, = 1 was used in deriving the relations (38) and (40). C(t,) in (38) is obtained from (37) with pe(tM) = pe and fM from (26) with [(TM) = 1, i.e. R(T,) = R,. The flow rate at the sand face flow for t > tM is determined from the modified Darcy’s law (6) and previous relations. In dimensionless form we have Q*(R,, t) = [ 1t
y]““; t>
tM
(41)
in which Q*(Rw,
T) =
(42)
For one dimensional flow, the previous relations (38) and (41) become P,(T) - PW= R +
(1
_
Pe-Pw
n)
1
_
[
3(2)““(n - l)(T - TM) “‘(n-m” 4n(l - fl)(“_‘)‘” I
(43)
and (44) where T
M
‘(l+flu)u”“du.
Q*(o
[2(1-flu)]“”
’
; “=F:
t)=
’
ffoL pe-p,,,=Ap
(45)
(46)
Figures l-3 show the effects of the rheological parameters n and R on the flow rate behaviour in function of time for t > tM, considering one dimensional flow, while Figs. 4-6 show the rheological effect on dimensionless pressure drop obtained from relation (43). These plots indicate a significant deviation from Newtonian flow, i.e. n = 1 and Sz= 0. For example, for a given value of the dimensionless group 0 = cw,LlAp,i.e. the flow in the presence of a threshold gradient, the deviation from case n = 1 is more significant for a pseudo-plastic fluid compared with a dilatant fluid. On the other hand, for a given value of the rheological parameter n, Figs. l-3 indicate that the threshold gradient has a significant effect on the flow rate behaviour for a power law fluid. Similar effects may also be observed from Figs. 4-6 with regard to the pressure
H. PASCAL
206
0 0.1
0.5
0.9
1.3
Dimensionless
1.7
2.1
2.5
Time, T
Fig. 1. Effect of rheological parameter n on dimensionless flow rate for Q = 0.
behaviour. Obviously, Figs. 1-6 show that the depletion time for an oil reservoir which contains a non-Newtonian Auid of power law, having a gel structure at zero rate of shear, depends strongly of rheological parameters n, H and TV. These parameters are sensitive to the temperature changes. Consequently, a correct evaluation of the rheological effect on flow behaviour should consider the temperature effect on rheological parameters. For example, flow at the reservoir conditions is isothermic, since
1.0-
I
I
I
I
T > TM
0 0 A + v 0
1.3
1.7
0.8 2a, a
2 0.6 z
A
n n n n n n
= = = = = =
,
1
2.1
2.5
0.4 0.6 0.8 1.0 1.2 1.4
i-i
Q)
z
.o !L! 0.4 E a
0.2
0 0.1
0.5
09
Dimensionless
Time, T
Fig. 2. Eflect of rheological parameter n on dimensionless flow rate for II = 0.2.
207
Nonsteady flow of non-Newtonian fluids through a porous medium
‘.OI
T>TM
a, 0.8 E 3 ; 0.6
0 n = o n = A n = + n = v n = 0 n =
0.4 0.6 0.8 1.0 1.2 1.4
z a, .g 0.4 S .-E cl 0.2
0 0.2
0.6
1.o
1.4
Dimensionless
1.8
2.2
2.6
Time, T
Fig. 3. Effect of rheological parameter n on dimensionless flow rate for 0 = 0.4.
temperature is a constant along the flow path. As a result, the reservoir temperature must be known and useful empiricism expressed by relations (3)-U), could be used in evaluating the temperature effect on rheological parameters. We now ask, what is the rheological effect on flow behaviour when instead of a constant pressure at sand face flow a constant flow rate is imposed. In this situation, the following boundary condition is valid (47)
Dimensionless
Time, T
Fig. 4. Effect of rheological parameter n on dimensionless pressure drop for Cl= 0.
fI. PASCAL
0 l-l =
T > T,
0 n An + n vn 0 n
-0.1
0.5
0.9
Dimensionless Fig. 5. Effect of rheological parameter
I.7
1.3
= = = = =
0.4 0.6 0.8 1.0 1.2 1.4
2.1
;s
Time, T
II on dimensionless
prtswe
drvp fcx ti = 0.2.
The corresponding solutions in this case may be obtained by a similar approach to that used for
D 0 A + v 0
n n n n n n
= = = = = =
0.4 0.6 0.8 1.0 1.2 1.4 i
/
0
0.2
I
0.6
1
1,O
1A
Dimensionless Fig. 6. Effect rrf rhwiogical
parameter
I
I
1,8
2,2
Time, T
n on dimensionhs
pressure drop ior n=0_4.
2.6
Nonsteady flow of non-Newtonian fluids through a porous medium a
209
constant pressure. For example, the pressure distribution for t < fM may be written as R2 -J-+X’(t) I
AR, t) = me- aoRw-
*-2R(t)
(49)
and for t > tM
R,QQ(R, XV, t>=
Inf-/)
Re-Rw
+C(t)[R’-R,R(l+~)+2R,R,ln~].
(50)
where
c(t) =
(51) R:[2$$
I)-
11
*
The front tocation R(t) in (49) is determined from (EI&+)I+T
(52)
in which dimensionless time is expressed as
T=
(53)
The Term (R,/R,J3 in (52) may be neglected, so that tM is obtained from (52) with R(T,,,)/R, = I and expressed, taking into account (53), as
t&#=
~~~R~(co#C CJ ; + 5d ( ) 3kR 2rR,h (*-‘) ’ w ( 00 )
(54)
Once the front location is known, eqn (49) provides the pressure distribution in system and pressure variation at sand-face flow in function of time for t < tM This variation for t > tM requires the knowledge of pressure variation p,(t) at the closed frontier. This may be determined from relation p,(r) = Pe -
YoQo (JQ)““(t a2,y,2ah k
- f&f); t > fW
where yo= R:[2~(ln~-
1)- 1]
and x,=-&R:+2RRi
i!$lnRlfl w R,
4(
l_ Rf d
(55)
210
H. PASCAL
An examination of eqns (50), (51) and (55) indicates that the pressure response is essentially linear with respect to time. As a result, for t > tM a semisteady flow occurs in the system. During this period splat is constant at any distance, as it may be seen from (50), (51) and (55). CONCLUSIONS
The behaviour of the solutions in closed form developed in this investigation shows that the rheological parameters n, 7. and H have a significant effect on the flow behaviour through a porous medium. The rheological parameters are associated with physical properties of the porous medium and boundary conditions, imposed at sand-face flow, in several dimensionless groups. For example, for a reservoir producing at a constant pressure the dimensionless group cyoRJAp determines the influence of threshold gradient on flow behaviour. If the reservoir is produced at a constant flow rate, this group becomes
in which the threshold gradient is expressed in terms of yield stress by relation CY~ = @T,Jv%. When these dimensionless groups have small values, the threshold gradient may be neglected for a power-law fluid when a yield stress. Since the rheological parameters n, 7. and H depend strongly on temperature, then rheological effect on flow behaviour should be evaluated at a specified temperature corresponding at flow conditions. For an isothermic flow, temperature is a constant known along the flow path. If the flow is nonisothermic, the temperature distribution is an unknown function of time and distance. To close the mathematical model for this situation, the use of energy equation, along with basic eqns (6) and (12), is required. NOMENCLATURE constant oil compressibility coefficient porous medium compressibility coefficient consistency index reservoir thickness permeability rheological parameter pressure distribution in system pressure at internal boundary pressure at external boundary flow rate distribution in system flow rate at the outface flow, (internal boundary) radial distance external boundary radius well radius front location time of flowing dimensionless time time required for the depression front to reach the closed boundary velocity porosity threshold gradient effective viscosity (rheological parameter for a power-law fluid) shear rate shear stress shear stress at zero rate of shear REFERENCES [I] ROJAG et al. In The Oil Sands of Canada-Venezuela, pp. 284-302 (1977). [2] W. SEITZER and P. LOWELL, Paper S.P.E. 9500 presented at 55th Ann. Fall Conf. of the Society of Petroleum Engineers, Dallas, Sept. (1980j. [3] H. PASCAL, Acta Mechanica, 39, (3-4) (1981). [4] H. PASCAL, Rev. Instit. Franc. du Petrole, Tom 34, No. 3 (1979). [5] F. PASCAL, H. PASCAL and D. MURRAY, J. Numerical and Analytic Methods in Geomechanics, NO. 3 (1981). [6] R. H. CHRISTOPHER and S. MIDDLEMAN, Power-law flow through a packed tube. Ind. Eng. Chem. Fund (1965). [7] B. R. BIRD el al., Transport Phenomena. Wiley, New York (1960). [8] D. TEEUW and T. HESSELINK, Paper 8982presented at the 5th Int. Symp. on Oiljield and Geomethermal Chemistry, California (1980). (Received 4 March 1982)
OOZO-7225/83/030199-l2W3.00/0 @ 1983 Pergarnon Press Ltd.
1983
NONSTEADY FLOW OF NON-NE~~NIAN THROUGH A POROUS TEDIUM
FLUIDS
H. PASCAL Department of Physics, Unkersity of Alberta, Edmonton, Alberta, Canada T6G 2Jl Abstract-In this paper, the nonsteady Row of power-law fluids with a yield stress through a porous medium is investigated. In order to illustrate the rheological behaviour effects on pressure and flow rate distributions, a flow system of practical interest was analysed. The approximate solutions in a closed form were obtained by using the integral method. For a short time, these have been formulated in terms of a moving boundary problem. The flow deviation from Newtonian behaviour was emphasized by means of several dimensionless groups. These have been found to be relevant in evaluating the rheological effects on steady and nonsteady flow behaviour through a porous medium.
INTRODUCTION
THETHEORIESused now in oil reservoir engineering to analyse the flow through a porous medium and to predict the reservoir performance are based on the assumption that crude oil is a Newtonian fluid and as a result Darcy’s law holds. While for light crudes, which are essentially Newtonian fluids, this assumption is justified, this is not so for a large class of complex fluids produced currently from certain oil fields. The heavy crudes, such as occur in Canada and Venezuela, are characterized by completely non-Newtonian behaviour as has been repeatedly demonstrated and reported in laboratory studies. Another illustrative example of non-Newtonian fluids are waxy crudes located in reservoirs of shallow depth. The reported results from the laboratory and field data have shown that heavy exudes at reservoir temperature display a rheological behaviour of power law with a yield stress. As a result of the presence of a gel structure, the shear stress-shear rate relationship was found to be a curve which does not pass through the origin. An example to support the experimental evidence of this behaviour was presented in papers [1, 21. Some heavy crudes show a reduction in apparent viscosity as the shear rate increases, others show an increase in apparent viscosity. From a rheological point of view these fluids are of pseudo-plastic or dilatant type. It should be mentioned that this behaviour has been observed above the yield stress, i.e. when this value was exceeded and presumably any structure that the fluid has is broken down. While for light crudes the viscosity is a parameter depending on temperature but not on conditions of measurement, the heavy crude viscosity depends strongly on shear rate. In order to fit the experimental data into a large range of shear rate variations, an analytical expression of the curves of rheological behaviour is required in solving the problems associated with flow trough porous medium. As a result, the knowledge of values of the apparent viscosity corresponding to conditions + + 0 and 9 + ~0is important. Most of the frequently used empirical relations to express analytically the apparent viscosity in terms of shear rate are various forms of the power law. For example, the most widely used for a fluid with a structure at zero shear rate is the empirical equation 7 =
H(j)” + 70;
7 >
70.
(0
By using this equation the apparent viscosity will be q,(i)
=
&y-’
+:
(2)
where n, H afid T,,are the rheological parameters to be determined from the rheological tests. These parameters depend strongly on temperature. Useful empiricism for temperature depen199
200
H.PASCAL
dence is presented in the literature as H(T)
= H(T,)
TV = ho
e-a(T-Td/Ta
(3)
+ NT - ToYTo
(4)
and n(T) = n(T,)+ c(T - T,)/T,
(5)
in which To is the reference temperature. The notations used are shown in the nomenclature. For n = 1, relation (1) reduces to the well known Bingham’s equation. This model was used in papers[3-51 to investigate the yield stress effect, i.e. threshold gradient, on the flow behaviour through a porous medium in several problems of fluid mechanics. Obviously the rheological model described by eqn (1) reflects the two most important aspects of the behaviour of heavy and waxy crudes observed from the experimental data; the existence of a yield stress and non-Newtonian flow curve above the yield stress. For example, at T > T,,the viscous flow becomes possible and its rheological behaviour may be described by eqn (1). It should be noted that at moderate and high shear rates the yield stress effect on apparent viscosity is insignificant. However, at low shear rate this effect becomes significant and cannot be neglected in eqn (2). In the production of heavy crudes, the shear rates which currently appear in flow through porous medium are expected to be in the range 0.1-1.0 set-‘. Consequently lower shear rate results are of practical interest, in which case a better rheological description should be the generalized Bingham’s model given by eqn (1). For example, it was found that, in general, the rheological models which included a yield stress better described the laboratory data for fluids with a gel structure than did the power law model. Another class of non-Newtonian fluids are polymer solutions, micro-macro emulsions. These fluids are often injected in oil reservoirs in order to improve the recovery efficiency. However this class of non-Newtonian fluids may be described by a power law equation without having to consider the yield stress. From the discussion presented above, it is evident that an understanding of non-Newtonian flow behaviour through porous medium is essential for practical purposes. The main objective of this investigation is to show the flow deviation from Newtonian behaviour expressed in terms of rheological parameters, using several examples of practical interest in engineering science. BASIC EQUATIONS
To consider the rheological effect in the flow description through a porous medium a modified Darcy’s law is required. For the rheological model given by eqn (1) the modtied Darcy’s law may be written as - & (grad p - a,,)]“”
v’= [
(6)
which is valid provided that lgrad pJ > cue;v # 0 [grad p[
(7)
Equation (6) was obtained by applying Kozeny’s approach to a capillary model of a power-law fluid in the presence of a yield stress. Using this approach, the following expression for pef was given in[6, 71 /.Lef= ; (9 + 3/n)“(lSOkf#-““2
(8)
Nonsteady flow of non-Newtonian fluids through a porous medium
201
Recently, the term k/per in (6) was found in[8] as (9) The threshold gradient in (6) may be expressed in terms of yield stress[3]
(10) where p is a dimensionless constant to be determined experimentally. As a result of the presence of a structure in heavy crudes which resist at a mechanical deformation, the modified Darcy’s law (6) implies a pressure gradient which increases or decreases monotonically to a given value. If there is successive increasing and decreasing of the pressure gradient it should be expected there will be a hysteresis loop effect, so that in this situation eqn (6) could not be valid. The material balance equation for a plan radial flow is f
& (Ru) =-(C&J + c@.
(11)
From equations (6) and (11) one obtains
$&[R($] =a2$
(12)
in which function x(R, t) is related to pressure by transformation
x(R t) = P@, t) - w,R.
(13)
For one dimensional flow, eqn (12) becomes
--a
ax ( > =aat
ax ax ax
I/”
2
(14)
where constant a2 in (12) and (14) is
a*= y
()
““(C&t
Cp).
(15)
Basic eqns (12) and (14) give an infinite velocity of wave front or pressure pulse propagation through a fluid flowing through a porous medium. To obtain a finite velocity the inertial term au/at should be considered in deriving the basic eqns (12) and (14). Although the infiniteness of the pressure wave velocity cannot be accepted on physical grounds, because it implies that the pressure perturbation generated in pore fluid by a wave front is felt instantaneously at any distance in an oil reservoir, the effect of finite velocity correction is extremely small and will be negligible in most problems of practical interest. Therefore, the unsteady flow through a porous medium may be described by using the parabolic eqns (12) and (14). The steady state solutions may be derived from (12) and (6) and expressed as P(R)=P,+~,(R-R,)+(~~~)~
(
$-
“(RI-” _ Rf;“) >
(16)
and me- pw -
m(Re- RJ
R;-“-R;-”
I
l/n
.
(17)
202
H.PASCAL
In practice we currently have RJR, = 10m3,so that in (17) and (16) one can use the following relation R w
l-n
(R >
-l+(l-n)ln+.
e
(18)
As a result, one obtains from (17)
(19)
FORMULATION
OFPROBLEMSANDSOLUTIONS
To illustrate the rheological behaviour effect on the unsteady flow through a porous medium, a classical problem of oil reservoir engineering will be considered. The flow system is a closed reservoir in which the porous medium is homogeneous. The solutions of practical interest are those corresponding to a centrally located well in the reservoir and producing at a constant pressure or constant flow rate. It is assumed that at initial moment t = 0, the pressure in the system is constant, say pe, and for t > 0, the system is depleted by a constant pressure of production, say pW As a result, the flow is initiated and maintained in the porous medium provided that pW< pe. A plan radial transient flow response will appear at t > 0 which may be analyzed by means of basic eqn (12) associated with certain boundary and initial conditions. Basic equations are nonlinear. Consequently, an approximate solution in closed form could be obtained by means of the integral method. Using this approach we have a moving boundary problem. This situation exists until the decompression front reaches the boundary of the reservoir located at distance R,. As a result, the appropriate boundary and initial conditions in this case will be POL, 0 = pw; 8P s R=R(t)=lYo; I PVW,
tl=
me;
R(O)=R,;
0 < t < 6~ O
(20)
O< t
t=O
where tM is the time required for the front to reach the closed frontier and R(t) represents the location of this front. The conditions (20) on function x(R, t), given by relation (13), become
xUL, 0 = pw - ~loRw ax 0 aR IR=R@)=
x[RW,
tl = me-
; O
(21)
4W)
R(0) = R,. The integral method has been used in our previous paper [3], in order to analyze an identical problem of one-dimensional flow for the particular case n = 1 and 7. # 0. This case corresponds to a Bingham fluid. Although no error analysis was made there, a comparison between the results obtained from numerical and exact solutions has shown a good agreement with those obtained from the integral method. The main advantage of the integral method appears to be the possibility of determining the pressure and flow rate-distributions without extensive computations and at a reasonable accuracy for engineering purposes. Using this approach, useful information on the flow behaviour has been obtained in[3]. However, it should be emphasized that the assumed pressure distributions must be selected with great care. For example, taking higher-order polynomial forms does not necessarily give better accuracy.
Nonsteadyflowof non-Newtonianfluidsthrougha porous medium
203
Taking into account the boundary conditions (21), the basic eqn (12) may be written as
a*
Rx(R, 0 dR - -_[P, - dW)l~ 2 A convenient consider
dR2
(22)
form of pressure distribution in the domain R, < R(t) < R,, for 0 < t < fM, is to
R x(R,t)=A+B~(t)+c~~(t)i
R2
O
(23)
Using the conditions (21) one obtains x(R, t) = me- LU,W) + C(t)( 1 - 2 &
(24)
+ &)
in which
C(f) =
W
(25)
The front location R(t) may be determined from (22) and (24) and expressed in terms of dimensionless time as (1 t $u)d’+“)‘” du (1 - RU)“n
T=
(26)
where
t(T)==;
T=
(27)
R,
From conditions (20) one obtains
(28) and one can observe that function (24) satisfies this equation. The particular case corresponding to a fluid of power law only, i.e. Cl= 0, may be determined from (26) with n # 1 and R = 0. In this situation, the front location becomes y
= [‘l t;n’T]“‘1+2”
(29)
and when n = 1 we have
( >
R(T) 3z3T R, ’
(30)
Once R(T) is known from (26), the pressure distribution in function of time and distance may be obtained from (24), (25) and (26). The time TM may be also determined from (26) with t(T) = 1. However, it should be mentioned that in deriving eqn (26), there was assumed 1 - RJR(t) = 1. This assumption is quite justified provided that RJR(t) CO.1; after a very short time of DES Vol. 21, No. 3-B
204
H.PASCAL
production the front location will reach such a distance in which condition RJR(t) < 0.1 is satisfied. Rather than to determine the pressure distribution in the system, the practical interest is to predict the flow rate behaviour in the well. For this purpose, eqns (6) and (24) should be used. For example, from these equations, the dimensionless flow rate in the well may be obtained and expressed as
Q*UL T) = [
1 - Q(7) 5(T)
1
I’”
(31)
with
Q*UL ‘J’)=
T) ($Y$ 2TR,h )I’”QoUL
(32)
and where t(T) one obtains from (26). To analyze the rheological effect on flow behaviour for t > t,, the following equation must be used _R
Rx(R, t) dR;
t > tM.
(33)
Obviously, the moving boundary problem has been eliminated, but the appropriate boundary and initial conditions in this case will be on function x(R, t)
xUL, t) = pw - aoh;
I
-dx aR
=o t> =
t&f
t>
tM
R=R,
xUL
t>
p,(t)
-
aoR,;
t >
(34)
ty
in which pe(tM) = pe = constant. As the numerical solution of the basic eqn (12) indicated, an appropriate form of the function x(R, t) for t > tM, associated with conditions (34), should be expressed as x(R, t) = A In F + RR + CR2; t > tM
(35)
w
which for conditions (34) becomes R(R - R,) t RJR, - 2R,) In $
x(R, t) = (pW- (Y,,R~)
1
W
(36)
where C(r)=n.(f)-
aoRe-
~(~~--a~R,)(1-ln~)
(37) RJR, - R,) t RJR, - 2R,) In +
W
Equation (37) contains the unknown function p,(t), i.e. pressure variation at closed frontier. This may be determined from (33) and (36). It should be noted that prior to the front reaching the closed frontier we have p(R,, t) = pe = constant, for t < tM, and p(R,, t) = p,(t) for t > tM, where pe = constant represents the pressure value in the system corresponding at initial moment t = 0. Using the condition pe(tM) = pe = constant, function p,(t) for t > tM may be obtained from
Nonsteady flow of non-Newtonian fluids through a porous medium
205
(33) and (36) and expressed as
2RtWf YOR,_ 2R:CM
2(n - l)R,R:(t - t,,,) na*~,(2RZ,C(t,) t y,,R,$-I)‘” I
+ Y~R, -
n/n-l (38)
where C(t) is given by relation (37) and yo =
(P, -
aoRd
3W
(39)
(40) In practice we currently have RJR, = 10m3,so that assumption 1 -RJR, = 1 was used in deriving the relations (38) and (40). C(t,) in (38) is obtained from (37) with pe(tM) = pe and fM from (26) with [(TM) = 1, i.e. R(T,) = R,. The flow rate at the sand face flow for t > tM is determined from the modified Darcy’s law (6) and previous relations. In dimensionless form we have Q*(R,, t) = [ 1t
y]““; t>
tM
(41)
in which Q*(Rw,
T) =
(42)
For one dimensional flow, the previous relations (38) and (41) become P,(T) - PW= R +
(1
_
Pe-Pw
n)
1
_
[
3(2)““(n - l)(T - TM) “‘(n-m” 4n(l - fl)(“_‘)‘” I
(43)
and (44) where T
M
‘(l+flu)u”“du.
Q*(o
[2(1-flu)]“”
’
; “=F:
t)=
’
ffoL pe-p,,,=Ap
(45)
(46)
Figures l-3 show the effects of the rheological parameters n and R on the flow rate behaviour in function of time for t > tM, considering one dimensional flow, while Figs. 4-6 show the rheological effect on dimensionless pressure drop obtained from relation (43). These plots indicate a significant deviation from Newtonian flow, i.e. n = 1 and Sz= 0. For example, for a given value of the dimensionless group 0 = cw,LlAp,i.e. the flow in the presence of a threshold gradient, the deviation from case n = 1 is more significant for a pseudo-plastic fluid compared with a dilatant fluid. On the other hand, for a given value of the rheological parameter n, Figs. l-3 indicate that the threshold gradient has a significant effect on the flow rate behaviour for a power law fluid. Similar effects may also be observed from Figs. 4-6 with regard to the pressure
H. PASCAL
206
0 0.1
0.5
0.9
1.3
Dimensionless
1.7
2.1
2.5
Time, T
Fig. 1. Effect of rheological parameter n on dimensionless flow rate for Q = 0.
behaviour. Obviously, Figs. 1-6 show that the depletion time for an oil reservoir which contains a non-Newtonian Auid of power law, having a gel structure at zero rate of shear, depends strongly of rheological parameters n, H and TV. These parameters are sensitive to the temperature changes. Consequently, a correct evaluation of the rheological effect on flow behaviour should consider the temperature effect on rheological parameters. For example, flow at the reservoir conditions is isothermic, since
1.0-
I
I
I
I
T > TM
0 0 A + v 0
1.3
1.7
0.8 2a, a
2 0.6 z
A
n n n n n n
= = = = = =
,
1
2.1
2.5
0.4 0.6 0.8 1.0 1.2 1.4
i-i
Q)
z
.o !L! 0.4 E a
0.2
0 0.1
0.5
09
Dimensionless
Time, T
Fig. 2. Eflect of rheological parameter n on dimensionless flow rate for II = 0.2.
207
Nonsteady flow of non-Newtonian fluids through a porous medium
‘.OI
T>TM
a, 0.8 E 3 ; 0.6
0 n = o n = A n = + n = v n = 0 n =
0.4 0.6 0.8 1.0 1.2 1.4
z a, .g 0.4 S .-E cl 0.2
0 0.2
0.6
1.o
1.4
Dimensionless
1.8
2.2
2.6
Time, T
Fig. 3. Effect of rheological parameter n on dimensionless flow rate for 0 = 0.4.
temperature is a constant along the flow path. As a result, the reservoir temperature must be known and useful empiricism expressed by relations (3)-U), could be used in evaluating the temperature effect on rheological parameters. We now ask, what is the rheological effect on flow behaviour when instead of a constant pressure at sand face flow a constant flow rate is imposed. In this situation, the following boundary condition is valid (47)
Dimensionless
Time, T
Fig. 4. Effect of rheological parameter n on dimensionless pressure drop for Cl= 0.
fI. PASCAL
0 l-l =
T > T,
0 n An + n vn 0 n
-0.1
0.5
0.9
Dimensionless Fig. 5. Effect of rheological parameter
I.7
1.3
= = = = =
0.4 0.6 0.8 1.0 1.2 1.4
2.1
;s
Time, T
II on dimensionless
prtswe
drvp fcx ti = 0.2.
The corresponding solutions in this case may be obtained by a similar approach to that used for
D 0 A + v 0
n n n n n n
= = = = = =
0.4 0.6 0.8 1.0 1.2 1.4 i
/
0
0.2
I
0.6
1
1,O
1A
Dimensionless Fig. 6. Effect rrf rhwiogical
parameter
I
I
1,8
2,2
Time, T
n on dimensionhs
pressure drop ior n=0_4.
2.6
Nonsteady flow of non-Newtonian fluids through a porous medium a
209
constant pressure. For example, the pressure distribution for t < fM may be written as R2 -J-+X’(t) I
AR, t) = me- aoRw-
*-2R(t)
(49)
and for t > tM
R,QQ(R, XV, t>=
Inf-/)
Re-Rw
+C(t)[R’-R,R(l+~)+2R,R,ln~].
(50)
where
c(t) =
(51) R:[2$$
I)-
11
*
The front tocation R(t) in (49) is determined from (EI&+)I+T
(52)
in which dimensionless time is expressed as
T=
(53)
The Term (R,/R,J3 in (52) may be neglected, so that tM is obtained from (52) with R(T,,,)/R, = I and expressed, taking into account (53), as
t&#=
~~~R~(co#C CJ ; + 5d ( ) 3kR 2rR,h (*-‘) ’ w ( 00 )
(54)
Once the front location is known, eqn (49) provides the pressure distribution in system and pressure variation at sand-face flow in function of time for t < tM This variation for t > tM requires the knowledge of pressure variation p,(t) at the closed frontier. This may be determined from relation p,(r) = Pe -
YoQo (JQ)““(t a2,y,2ah k
- f&f); t > fW
where yo= R:[2~(ln~-
1)- 1]
and x,=-&R:+2RRi
i!$lnRlfl w R,
4(
l_ Rf d
(55)
210
H. PASCAL
An examination of eqns (50), (51) and (55) indicates that the pressure response is essentially linear with respect to time. As a result, for t > tM a semisteady flow occurs in the system. During this period splat is constant at any distance, as it may be seen from (50), (51) and (55). CONCLUSIONS
The behaviour of the solutions in closed form developed in this investigation shows that the rheological parameters n, 7. and H have a significant effect on the flow behaviour through a porous medium. The rheological parameters are associated with physical properties of the porous medium and boundary conditions, imposed at sand-face flow, in several dimensionless groups. For example, for a reservoir producing at a constant pressure the dimensionless group cyoRJAp determines the influence of threshold gradient on flow behaviour. If the reservoir is produced at a constant flow rate, this group becomes
in which the threshold gradient is expressed in terms of yield stress by relation CY~ = @T,Jv%. When these dimensionless groups have small values, the threshold gradient may be neglected for a power-law fluid when a yield stress. Since the rheological parameters n, 7. and H depend strongly on temperature, then rheological effect on flow behaviour should be evaluated at a specified temperature corresponding at flow conditions. For an isothermic flow, temperature is a constant known along the flow path. If the flow is nonisothermic, the temperature distribution is an unknown function of time and distance. To close the mathematical model for this situation, the use of energy equation, along with basic eqns (6) and (12), is required. NOMENCLATURE constant oil compressibility coefficient porous medium compressibility coefficient consistency index reservoir thickness permeability rheological parameter pressure distribution in system pressure at internal boundary pressure at external boundary flow rate distribution in system flow rate at the outface flow, (internal boundary) radial distance external boundary radius well radius front location time of flowing dimensionless time time required for the depression front to reach the closed boundary velocity porosity threshold gradient effective viscosity (rheological parameter for a power-law fluid) shear rate shear stress shear stress at zero rate of shear REFERENCES [I] ROJAG et al. In The Oil Sands of Canada-Venezuela, pp. 284-302 (1977). [2] W. SEITZER and P. LOWELL, Paper S.P.E. 9500 presented at 55th Ann. Fall Conf. of the Society of Petroleum Engineers, Dallas, Sept. (1980j. [3] H. PASCAL, Acta Mechanica, 39, (3-4) (1981). [4] H. PASCAL, Rev. Instit. Franc. du Petrole, Tom 34, No. 3 (1979). [5] F. PASCAL, H. PASCAL and D. MURRAY, J. Numerical and Analytic Methods in Geomechanics, NO. 3 (1981). [6] R. H. CHRISTOPHER and S. MIDDLEMAN, Power-law flow through a packed tube. Ind. Eng. Chem. Fund (1965). [7] B. R. BIRD el al., Transport Phenomena. Wiley, New York (1960). [8] D. TEEUW and T. HESSELINK, Paper 8982presented at the 5th Int. Symp. on Oiljield and Geomethermal Chemistry, California (1980). (Received 4 March 1982)