- Email: [email protected]
The Solution of the Riemann Problem for a Hyperbolic System Modeling Polymer Flooding with Hysteresis
JOURNAL OF MATHEMATICAL ANALYSIS AND APPLICATIONS ARTICLE NO.
206, 205]233 Ž1997.
AY975222
The Solution of the Riemann Problem for a Hyperbolic System Modeling Polymer Flooding with Hysteresis Khaled M. Furati Department of Mathematical Sciences, King Fahd Uni¨ ersity of Petroleum and Minerals, Dhahran 31261, Saudi Arabia Submitted by Barbara Lee Keyfitz Received January 23, 1995
It is well known that multiphase flow in porous media exhibits hysteresis. This is typically modeled by modifying the saturation dependence of the relative permeabilities. In this paper, a model for hysteretic relative permeabilities is built into the polymer flooding model and the analytical solution to the corresponding Riemann problem is constructed. This produces a nonstrictly hyperbolic system of conservation laws with a history-dependent flux function. Because the polymer model without hysteresis possesses Riemann problem solutions that are not monotonic, the introduction of hysteresis necessarily produces structurally different solutions. We show that hysteresis produces more complicated solutions with more fronts and expansions; and removes some nonuniqueness of solutions. Q 1997 Academic Press
1. INTRODUCTION In this paper we solve the Riemann problem for the system of conservation laws t s q x F s 0, Ž 1. t Ž cs . q x Ž cF . s 0, with the flux function F of the form Fs
½
f Ž s, c .
if s s p , t s ) 0
h Ž s, p , c .
otherwise
Ž primary flow . Ž secondary flow .
Ž 2.
and initial conditions
Ž s, p , c . Ž x, 0 . s
½
Ž sL , p L , cL . , Ž sR , p R , cR . ,
x F 0, x ) 0.
Ž 3.
205 0022-247Xr97 $25.00 Copyright Q 1997 by Academic Press All rights of reproduction in any form reserved.
206
KHALED M. FURATI
where, t ) 0, x g R. The variables s, p , c are in I and s F p , where I is the interval w0, 1x. The system Ž1. is used in enhanced oil recovery to model polymer flooding in which oil is displaced by water containing dissolved polymer. The variable s is the saturation of the aqueous phase and c is the concentration of polymer. The flux function F is the fractional flow function for the aqueous phase. By analogy to the elastic-plastic bar problem in w13x, where the history of the material is accounted for by the plastic displacement gradient p , we let p be the hysteresis parameter. For the derivation of Ž1. and physical assumptions see w12x. Experiments have shown that the flow functions are not only functions of current fluid saturations but also functions of the saturation history. Moreover, the solution of the Riemann problem for the polymer flooding model without hysteresis w5x shows cases with nonmonotonic behavior in saturation. The form of F in Ž2. is based on remembering the saturation history through parameterizing the states at which the imbibition process is reversed. The function F s f Ž s, c . describes the flow when the saturation increases monotonically from the irreducible aqueous phase saturation sra to the maximum saturation 1}srl , where srl is the residual liquid phase saturation. We consider this imbibition of the aqueous phase as the primary flow. When this imbibition process is reversed at s s p , a new flow function F s hŽ s, p , c . appears and the flow is called a secondary flow. The main contribution of this paper is incorporating the flow history into the flux function and constructing the unique global solution of the Riemann problem. The Oleinik entropy condition w11x for scalar equations and the generalized Lax entropy condition w7x are used to distinguish the physically meaningful discontinuities. The solution is constructed by piecing together a compatible sequence of elementary waves Žshocks, rarefactions, and contact discontinuities . to connect constant states. Including hysteresis produces more complicated solutions. Another effect of including hysteresis is that it removes some of the nonuniqueness cases arises in w5x, where the solution is unique in the x, t-space but not in the state space. However, for some p R large enough the same nonuniqueness problem persists. Although for certain critical values of the left and right states the solution is not pointwise continuous with respect to the left and right states, it is continuous in L1 norm. Buckley and Leverett w1x constructed the analytical solution for immiscible waterflooding. System Ž1. reduces to their single equation model when c is constant. Pope w12x generalized the fractional flow theory to more complicated floodings and used the method of characteristics to construct the corresponding solutions. The polymer flooding model without hysteresis was analyzed and solved by Isaacson w5x and Johansen and Winther w6x.
POLYMER FLOODING WITH HYSTERESIS
207
The problem is similar to the Riemann problem for vibration of an elastic string considered by Keyfitz and Kranzer w7x. Hysteresis phenomena in multiphase flows in porous media has been recognized by many researchers. For example, Colonna et al. w2x, Killough w8x, Gladfelter and Gupta w4x, and Lenhard w9x. Marchesin et al. w10x used a double-valued flux function to model two-phase flow with hysteresis and constructed, graphically, the complete solution of the associated Riemann problem. However, in this paper we consider a three-component two-phase flow with a more complicated model of hysteresis. The treatment of the history dependence in this paper is similar to that of motion in an elastic-plastic bar studied by Trangenstein and Pember in w13x. The outline of the paper is as follows. In Section 2 we present a brief derivation of the model. In Section 3 we state the mathematical assumptions and properties of the flux functions and analyze the hyperbolicity of the model. In Section 4 we discuss the different wave families. In Sections 5]8 we construct the global solution for the Riemann problem.
2. DERIVATION OF THE MODEL In polymer flooding, a solution of water and polymer is injected into an oil reservoir to push the oil through the rock. The reservoir fluid is considered to consist of three components: oil, water, and polymer. The oil forms its own liquid phase, and water with polymer form the aqueous phase. The polymer transports only in the aqueous phase and does not partition to oil. The purpose of the polymer is to increase the viscosity of water, thereby enhancing its ability to push the oil. System Ž1. is derived under many assumptions. We assume that the reservoir rock is homogeneous, so that the porosity and total permeability of the rock are constant. For the phases, we assume incompressibility and that the liquid phase viscosity is constant, but the aqueous phase viscosity increases as the concentration of polymer increases. Moreover, we assume that there are no diffusive forces, such as capillary pressure, and no gravity forces. The mass conservation equation for water, oil, and polymer, respectively, can be formulated as
f t sa q x¨ a s 0, f t s l q x¨ l s 0, f t Ž csa . q x Ž c¨ a . s 0,
Ž 4.
208
KHALED M. FURATI
where f is the rock porosity, c is the polymer concentration, and sa and s l are the aqueous and liquid phases saturations, respectively. The fluxes ¨ a and ¨ l are the volumetric flow rates of the aqueous and liquid phases, respectively, and are given by Darcy’s law: ¨a s y
Kk ra
ma Ž c .
x P ,
¨l s y
Kk rl
ml
x P.
Ž 5.
Here, P is the fluid pressure, K is the absolute permeability of the rock, k ra and k rl are the relative permeabilities of the aqueous and liquid phases, respectively, and ma , m l are the corresponding viscosities. Next, we use the assumption that the total velocity ¨ t s ¨ a q ¨ l is constant to eliminate the pressure gradient. Adding the first two equations in Ž4., using the assumption that s l q sa s 1 and changing x to f xr¨ t we obtain Ž1. with s s sa and Fs
k ra k ra q k rl m Ž c .
,
Ž 6.
where m s marm l is the viscosity ratio. To complete the description of the conservation laws we need to describe the relative permeability curves. In this paper, we consider the typical relative permeability curves shown in Fig. 1, with the aqueous phase being the wetting phase. The solid curves give the relative permeability values when the saturation s is monotonically increasing from sra . We call such curves the primary curves. When the aqueous phase saturation decreases, the relative permeability is smaller and follows the dashed
FIG. 1. Hysteretic relative permeability curves.
209
POLYMER FLOODING WITH HYSTERESIS
curves. We call such curves the secondary curves. When the saturation increases again, we assume that the relative permeability traces back the secondary curve until it meets the primary curve, then follows the primary curve. We parameterize the secondary curves by the symbol p , which we refer to as the hysteresis parameter. The primary and secondary curves are denoted by k r j Ž s . and k rhj Ž s, p ., respectively, where j s a, l. Thus we have
krj s
½
krj Ž s. k rhj
Ž s, p .
if s s p , t s ) 0, otherwise,
5
j s a, l.
Ž 7.
With the above forms of the relative permeability curves, the fractional flow function F takes the form Ž2. with f Ž s, c . s h Ž s, p , c . s
k ra Ž s . k ra Ž s . q k rl Ž s . m Ž c .
,
h k ra Ž s, p . h k ra Ž s, p . q k rlh Ž s, p . m Ž c .
Ž 8. .
Figure 2 shows typical fractional flow curves. Note that since the polymer increases water viscosity, the ratio m s m Ž c . is an increasing function of c. This implies that the fractional flow function F is decreasing with respect to c. Physically, this means that adding polymer will increase the oil flow rate which is equal to 1 y F.
FIG. 2. Fractional flow curves, c 2 ) c1.
210
KHALED M. FURATI
3. THE MATHEMATICAL DESCRIPTION OF THE MODEL In this section we present the mathematical assumptions and properties of the relative permeabilities and, hence, the flux functions. Accordingly, we analyze the hyperbolicity of the system of conservation laws Ž1.. Let I s w0, 1x, Is s w sra , 1 y s sl x, and Ip s w sp , p x, where sp s max s : hŽ s, p , c . s 04 . We assume that the relative permeability curves satisfy the following conditions, as shown in Fig. 1: A1. j s a, l;
k r j Ž s . g C 2 Ž Is ., k rhj Ž?, p . g C 2 Ž Ip ., and k rhj Ž s, ? . g C 1 Ž Is ., where
A2.
k rhj Ž s, s . s k r j Ž s . ;s g Is , j s a, l;
A3.
kXra ) 0, kYra ) 0, kXrl - 0, and kYrl ) 0 ;s g Is ;
h h h A4. s k ra ) 0, p k ra - 0, s s k ra ) 0, s k rlh - 0, p k rlh ) 0, and s s k rlh ) 0 in Ip = Is ; A5. kXra Ž sra . q m Ž c . kXrl Ž sra . F 0 and kXraŽ1 y srl . q m Ž c . kXrl Ž1 y srl . G 0 ;c g I; hŽ A6. s k ra sp , p . q m Ž c . s k rlh Ž sp , p . F 0 ;c g I.
Assumption A5 is a sufficient Žnot necessary . condition for the primary fractional flow function, f Ž., c ., to be convex in a neighborhood of s s sra and concave in a neighborhood of s s 1 y srl . This implies the existence of an inflection point of f. h We assumed p k ra - 0 and p k rlh ) 0 so that the secondary relative permeability curves are disjoint for distinct values of p . Assumption A6 is a sufficient Žnot necessary . condition for the secondary fractional flow function hŽ s, p , c . to be convex in some neighborhood of the irreducible saturation sp . With direct examinations of the derivatives, the above assumptions imply the following properties for the flux functions f and h as shown in Fig. 2: P1. f Ž s, c . is a C 2 function in Is = I, and f Ž sra , c . s 0 and f Ž1 y srl , c . s 1 ;c g I; P2. s f Ž s, c . ) 0 and c f Ž s, c . - 0 in Is = I; P3. For each c g I, there is an inflection point, s I s s I Ž c . g Ž sra , 1 y srl ., such that s s f Ž s, c . ) 0 for s - s I ; P4. hŽ s, p , c . is a C 2 function in s g Ip , C 1 in p g Is and C 2 in c g I; P5. hŽ sp , p , c . s 0 and hŽ s, s, c . s f Ž s, c . for each p , s g Is and c g I;
POLYMER FLOODING WITH HYSTERESIS
P6. = I;
211
s hŽ s, p , c . ) 0, p hŽ s, p , c . - 0 and c hŽ s, p , c . F 0 in Ip = Is
P7. s s hŽ sp , p , c . ) 0 for each c g I and p g Is ; P8. s f Žp , c . F sy hŽp , p , c ., for each c g I, where sy is the left derivative. Note that the assumptions on relative permeability curves imply existence Žproperty P3. but not necessarily uniqueness of the inflection points of f. However, for most applications, experimental data indicates that the fraction flow function does not have more than one inflection point. Therefore, for each c g I we assume that s I Ž c . is the only point of inflection, and consequently, f Ž., c . is convex in Ž sra , s I Ž c .. and concave in Ž s I Ž c ., 1 y srl .. Property P6 implies that the secondary fractional flow curves for a fixed c are disjoint for distinct values of p . Also, it implies that for a fixed p the curves are disjoint for distinct values of c. Unlike the primary fractional flow function, f Ž s, c ., assumptions on hŽ k ra s, p . and k rlh Ž s, p . do not guarantee the existence of inflection points for the secondary fractional flow function, h. Moreover, the assumptions do not guarantee the uniqueness of such points for the curves of h, as for f. The assumptions only guarantee property P7 which implies that all the curves of h are convex in some neighborhood of sp . We assume that for each c g I and p g Is , hŽ s, p , c . has at most one inflection point, s I Žp , c ., for s g Ip . For the primary flow, system Ž1. can be written in the quasilinear form
t cs q A x cs s 0,
As
s f 0
c f . frs
Ž 9.
Since the eigenvalues of A are real, the system of conservation laws is hyperbolic and has the two characteristic speeds and the corresponding eigenvectors:
l b s s f ,
rb s
1 , 0
lp s
f s
,
rp s
c f . l p y lb
Ž 10 .
The characteristic speed l b is the Buckley]Leverett rarefaction wave speed, which is equal to the slope of f. The second characteristic speed, l p , is the aqueous particle velocity, which is equal to the slope of the chord from the origin to the fractional flow function.
212
KHALED M. FURATI
If the flow is secondary, F s hŽ s, p , c ., then the system of conservation laws Ž1. can be closed by adding the constraint t p s 0 and written in the quasilinear form s s t c q A˜x c s 0, p p
A˜ s
s h 0 0
c h hrs 0
p h 0 . 0
Ž 11 .
The characteristic speeds and their corresponding eigenvectors are
l b s s s h,
1 rb s s 0 , 0
lp s 0,
lps s
h s
,
p h rp s 0 yl b s
c h rp s s l p s y lb s , 0
Ž 12 .
Note that all the characteristic speeds are nonnegative. This implies that all the wave families are either stationary or traveling to the right. For the border states u s Žp , p , c ., we define the Buckley]Leverett characteristic speeds by l b Žp , p , c . s s f Žp , c . and l b s Žp , p , c . s sy hŽp , p , c .. The type of the flow Žsecondary or primary. determines which velocity is to be used in determining the Riemann problem solutions. Note that for a given c the Buckley]Leverett and particle velocities might coalesce as shown in Fig. 3. However, this is not always the case, as shown in Fig. 4. We let utan Ž c . and utans Žp , c . be the states such that l b Ž utan . s l p Ž utan . and l b s Ž utans . s l p s Ž utans .. For the fractional flow curves described above we have the following results. LEMMA 3.1. For each c g I, f Ž., c . has at most one tangential state utan Ž c .. If such a state exists then l b Ž s, c . ) l p Ž s, c . for sr a - s - s tan Ž c . and l b Ž s, c . - l p Ž s, c . for s tan Ž c . - s. Proof. Note that by definition, the tangential states utan are the solutions of the equation l b y l p s s f y frs s 0. Given c g I, let g Ž s . s s s f Ž s, c . y f Ž s, c . for s g Ž sra , 1 y srl .. Since g X Ž s . s s s s f Ž s, c . and f Ž s, c . is assumed to have a unique inflection point u I s Ž s I Ž c ., c ., g Ž s . has only one local maximum at s I in Ž sra , 1 y srl .. Therefore, g s 0 can have at most one solution since g Ž sra . G 0. This proves the uniqueness because any s tan is by definition a solution of g s 0. The second part follows from the observation that g changes from positive to negative at s s s tan .
POLYMER FLOODING WITH HYSTERESIS
213
FIG. 3. Tangential states.
The next lemma shows for which values of p the tangential state utans might exist. LEMMA 3.2. For a gi¨ en c g I and p - s tan Ž c ., l b s Ž s, p , c . ) l p s Ž s, p , c . for s g Ž sp , p x. Proof. For a fixed c g I and p g Ž sra , s tan Ž c .., let g Ž s . s s s hŽ s, p , c . y hŽ s, p , c .. From Property P7 of the secondary fractional flow curves, we
FIG. 4. No tangential states.
214
KHALED M. FURATI
have g X Ž sp . s sp s s hŽ sp , p , c . ) 0. In addition, g satisfies the following inequalities: g Ž sp . s sp s h Ž sp , p , c . G 0, g Ž p . s p syh Ž p , p , c . y h Ž p , p , c .
Ž 13 .
) p s f Ž p , c . y f Ž p , c . ) 0. From those inequalities and the assumption that h has at most one inflection state, it follows that g is positive for s g Ž sp , p x. This proves the result. The above lemma implies that the tangential state, utans s Ž s tans Žp , c ., p , c ., may only exists if p G s tan Ž c .. The following lemma is about the uniqueness of utans . LEMMA 3.3. Gi¨ en c g I and p G s tan Ž c ., there exists at most one tangential state utans Žp , c . s Ž s tans Žp , c ., p , c ., sra - s tans - p . Proof. As in the proof of Lemma 3.1, the result follows from the uniqueness of the inflection point of h. In the above proofs, note that the uniqueness of the tangential states follows from the uniqueness of the inflection points. However, this does not guarantee their existence. In this paper we assume the existence of utan but not necessarily utans . However, if for p G s tan the state utans does not exist we have the following lemma. LEMMA 3.4. Gi¨ en c g I and p G s tan Ž c ., if there is no tangential state u , then the state up s Žp , p , c . satisfies tans
l b s Ž up . y l p s Ž up .
l b Ž up . y l p Ž up . F 0.
Ž 14 .
Proof. Again, for fixed c and p consider the function g Ž s . s s s h y h. By assumption, g s 0 has no solution for s g Ž sp , p .. But g Ž sp . G 0 and g X Ž sp . ) 0. Therefore, g Ž s . ) 0 for sp - s - p . By continuity, this implies that g Žp . G 0, and thus the first bracket in the inequality Ž14. is nonnegative. The second bracket is nonpositive by Lemma 3.1. This proves the result. To summarize, for each c g I we assume that there exists a state utan s Ž s tan Ž c ., c . such that l p Ž utan . s l b Ž utan .. Moreover, for each c g I and p g Is , if p - s tan Ž c . then l p s Ž u. - l b s Ž u. for all u s Ž s, p , c . with sp - s F p . If p G s tan Ž c ., then either there exists a state utans s Ž s tans Žp , c ., p , c . such that l p s Ž utans . s l b s Ž utans . or the state Žp , p , c . is
POLYMER FLOODING WITH HYSTERESIS
215
a transitional state across which the relation between the two speeds is reversed. Both states utan and utans are unique. We partition the set of all states into the following sets: P s u : l p Ž u . ) lb Ž u . 4 ,
B s u : lb Ž u . ) l p Ž u. 4 ,
T s u : l p Ž u . s lb Ž u . 4 , Ps s u : l p s Ž u . ) l b s Ž u . 4 ,
Bs s u : l b s Ž u . ) l p s Ž u . 4 ,
Ts s u : l p s Ž u . s l b s Ž u . 4 . We will use these sets to characterize the different cases for the solutions.
4. WAVE FAMILIES In this section we determine all wave families associated with the Riemann problem Ž1. ] Ž3.. Those wave families consist of the elementary waves: rarefactions, shocks, and contact discontinuities. To find all the possible waves we need to examine the nonlinearity of the characteristic fields, the conditions for discontinuities and the integral curves of the characteristic eigenvectors. For the nonlinearity we need to check the projection of the characteristic speed gradient into the corresponding eigenvector. For our problem we have =l p ? r p s =p l p s ? r p s s =p lp ? rp s 0,
=l b ? r b s s s f ,
=p l b s ? r b s s s s h, where = s Ž s , c . and =p s Ž s , c , p .. Consequently, the p-field, ps-field, and p-field are linearly degenerate and therefore the corresponding families consist only of contact discontinuities. On the other hand, due to the inflection points of f and h, the b-family and bs-family are not typically genuinely nonlinear. We expect both families to consist of rarefactions and shocks. As we know, any discontinuity solution of speed s must satisfy the Rankine]Hugoniot conditions. For the primary flow those conditions take the form fRyfL c f yc f R
R
L
L
s
sR y sL s. c s y cL sL R R
Ž 15 .
216
KHALED M. FURATI
One solution for Ž15. is s s l p s frs. Another solution is determined by taking c to be constant and s s w f xrw s x, where w?x denotes the jump of the argument across the discontinuity. For the secondary flow, the Rankine]Hugoniot conditions take the form sR y sL hR y h L R L L s cR sR y cL sL s . c h yc h 0 pRypL R
Ž 16 .
The conditions in Ž16. have the three solutions: 1. 2. 3.
s s w h xrw s x and w c x s wp x s 0; s s l p s s hrs and wp x s 0; s s 0 and w h x s w c x s 0.
To search for rarefactions we need to consider the integral curves of the eigenvectors r b and r b s . The integral curves of r b are the curves where c is constant. Similarly, the integral curves of r b s are the curves where c and p are constant. Accordingly, the wave families for the Riemann problem for system Ž1. ] Ž3. can be classified as follows. b-wa¨ es and bs-wa¨ es. The two families correspond to the characteristic speeds l b and l b s , respectively, and consist of rarefactions and shocks. Note that c is a Riemann invariant of both families and p is a Riemann invariant of the bs-family. Thus, from the solutions of the Rankine]Hugoniot conditions, the wave curves for the shocks and rarefactions of both families are identical. The b-family curves are the ones along which c is constant. The bs-family curves are the ones along which both c and p are constants. As a result, both families can be considered as solutions of the scalar equation t s q x F s 0, F s f, h. To distinguish the physically meaningful shocks, we need to impose an entropy condition. For our purpose, a shock is considered admissible if it satisfies the Oleinik condition w11x for scalar equations. This condition requires that F Ž uR . y F Ž u . sR y s
Fss
F Ž uR . y F Ž u L . sR y sL
Ž 17 .
for all s between s L and s R . In other words, the chord from Ž s R , F Ž u R .. to Ž s L , F Ž u L .. must lie above the chord from Ž s R , F Ž u R .. to Ž s, F Ž u.. for any s in between. Since both families are not genuinely nonlinear, the solution may consist of a combination of a shock and rarefaction. For the rarefactions, the
POLYMER FLOODING WITH HYSTERESIS
217
characteristic speeds should be monotonically increasing from left to right. This monotonicity and Oleinik condition are equivalent to constructing the concave hull for s L ) s R and the convex hull for s L - s R as shown in Fig. 5. A straight line segment represents a shock and a curve segment represents a rarefaction. The shock speed is equal to the slope of the straight line. In petroleum engineering literature, the above solutions are usually explained in terms of Welge tangents w3, 14x. The Welge tangents are equivalent to the Oleinik chords for such problems. We use k-wave, k s b, bs, to denote a k-shock, k-rarefaction, or any consecutive combination of them. Note that to have a b-wave we must have s L ) s R so it corresponds to a primary flow. p-wa¨ es and ps-wa¨ es. The two families correspond to the characteristic speed l p and l p s , respectively. The two families are linearly degenerate and consist of contact discontinuities since l p and l p s are Riemann invariants. The wave curves for the p-family are the curves c s cŽ s . along which l p is constant. The slope of such curves depends on the relation between l p and l b since dcrds s Ž l p y l b .rf c . The same argument is true for the ps-family curves c s cŽ s, p .. The states u L s Ž s L , c L . and u R s Ž s R , c R ., s L ) s R , can be connected by a contact discontinuity if l p Ž u L . s l p Ž u R . and they both belong to either P j T or B j T. The last condition is to fulfill the generalized Lax’s entropy condition w7x which requires that exactly one of the b-characteristics leave the discontinuity. With p L s p R , similar results are true for the ps-family.
FIG. 5. Concave and convex hulls.
218
KHALED M. FURATI
st-wa¨ es. This family consists of stationary contact discontinuities satisfying Rankine]Hugoniot conditions with speed s s 0. Across this discontinuity, w h x s w c x s 0.
5. GLOBAL SOLUTION OF THE RIEMANN PROBLEM The global solution to the Riemann problem can be constructed graphically using the wave families mentioned above. The uniqueness of the solutions can easily be verified by imposing the compatibility condition, i.e., the initial speed of each wave is greater than or equal to the final speed of the preceding wave. We start with the following remarks. Remark 1 ŽNotation.. Although the problem during a primary flow is independent of the history parameter p , we write any state u s Ž s, c . as u s Ž s, s, c .. Thus, for primary and secondary flows we can let u L s Ž s L , p L , c L . and u R s Ž s R , p R , c R . be the left and right states of the Riemann problem. For the rest of this paper we use hŽ s, s, c . to denote f Ž s, c .. For simplicity, we use h a to denote hŽ u a . and l a to denote lŽ u a .. k
Moreover, we use ª to denote a k-wave. Remark 2 ŽReduction.. For any given two states u L s Ž s L , p L , c L . and u s Ž s R , p R , c R . let u Q s Žp R , p R , c L .. Accordingly, there are two cases as shown in Fig. 6. If hQ G hL , then there is a unique state u A s Ž s A , p R , c L . such that h A s hL , and u L can be connected to u A by a R
FIG. 6. Left: hQ G hL . Right: hQ - h L .
POLYMER FLOODING WITH HYSTERESIS
219
st
stationary wave u L ª u A . Otherwise, there is a unique state u P s Ž s P , s P , c L ., s P F s L , such that h P s hL and u L can be connected to u P by st
a stationary wave u L ª u P . Thus, the problem can be reduced to solving either for u L s Ž s L , p R , c L . or u L s Ž s L , s L , c L ., p R F s L , for any arbitrary u R . To account for all possible cases, we treat the three different cases separately: c L s c R , c L ) c R , and c L - c R . Solutions of the Riemann problem will be composed of compatible sequences of the st-waves, bswaves, b-waves, p-waves, and ps-waves.
6. THE CASE c L s c R ŽBUCKLEY]LEVERETT EQUATION WITH HYSTERESIS. When c L s c R , the system Ž1. degenerates to the Buckley]Leverett equation and the solution consists of b-waves and bs-waves. For this special case we omit the concentration argument and write any state u as u s Ž s, p .. From Remark 2, we only need to consider the following two cases: 1. u L s Ž s L , p R .. In this case, u L and u R lie on the same secondary curve h s hŽ s, p R .. The solution is a bs-wave determined by the concave hull if s R - s L or by the convex hull if s R ) s L . Thus, we have the wave bs
u L ª u R across which p is constant and is equal to p R . See Fig. 7 for the case s R - s L .
FIG. 7. Solution for c L s c R , p L s p R and s L ) s R .
220
KHALED M. FURATI
2. u L s uŽ s L , s L . and p R - s L . Consider first the intermediate states u s Ž s, s ., p R F s - s L . The left state u L can be connected to any of those states by a unique b-wave determined by the concave hull of the primary curve joining the two states. In particular, for the state Žp R , p R . we have b the connection u L ª Žp R , p R .. On the other hand, the state Žp R , p R . can be connected to the right state u R by a unique bs-wave determined by the bs concave hull of the secondary curve. Thus, we have Žp R , p R . ª u R . If the final speed of the b-wave is less than the initial speed of the bs-wave, then the two waves are compatible and can be composed to give the unique b bs solution u L ª Žp R , p R . ª u R as shown in Fig. 8a. Note that, in this case, the union of the two concave hulls is also concave. Otherwise, we need to construct the concave hull of the union of the two hulls as in Fig. 8b. The saturation s M satisfies f X Ž sM . s
f Ž sM . y hŽ sR , p R . sM y sR
.
Ž 18 .
Thus, the solution starts with a rarefaction from u L to u M along the primary curve, followed by a shock to u R . This shock is actually a composition of a b-shock and a bs-shock traveling with the same speed. 7. THE CASE c L ) c R Considering the left state u L s Ž s L , p L , c L ., it is enough to solve for the following two cases as pointed out in Remark 2.
FIG. 8. Solution for c L s c R and s L s p L ) p R .
POLYMER FLOODING WITH HYSTERESIS
221
7.1. The Case p L s p R The different possibilities for this Riemann problem are listed in Table I. Note that when u L g Bs j Ts , u R g Ps , and lLp s s lRp s , the solution is not unique in the state space but it is unique in the physical space. Similar nonuniqueness occurs when u L g Ps , u R g Ps , and lRp s s lTp s , where
TABLE I Solutions for c L ) c R and p L s p R
222
KHALED M. FURATI
uT s Ž sT, p R , c L . g Ts is the tangential state as in Fig. 3. Moreover, both special cases are examples of solutions that are not pointwise continuous with respect to the initial data but continuous in L1 norm. Note also that all the flow is secondary for all different cases. 7.2. The Case s L s p L ) p R With such initial conditions, part of the flow is primary and the other is secondary. This is because the saturations of some intermediate states exceed the history parameter p s p R . Consider the following two cases listed in Table II. Case 1. u L g B j T. As shown in Fig. 9, there are two possibilities for u L with respect to u F s Žp R , p R , c R .. For both, any compatible sequence consists of a contact discontinuity leaving the left state followed by a b-wave or a bs-wave. Case 2. u L g P. The slowest wave that leaves u L is a b-wave. How far this wave can be extended depends on the position of uT s Ž sT, sT, c L . g T with respect to u R . If sT F p R , then both u L and u A s Žp R , p R , c L . lie on the concave part of the primary curve and can be connected by the rarefaction b
u L ª u A with final speed equal to lAb . Afterwards, u A can be connected to u R as in Section 7.1. There, the wave that leaves u A is either a bs-rarefaction or a ps-wave. The initial speed of both waves is greater than or equal to lAb . Thus, the whole connection is compatible. When sT ) p R , uT considered as the left state, can be connected to u R as in Case 1 in which the initial speed of the p-wave leaving uT is equal to lTp . Moreover, u L and uT can be connected through a b-rarefaction. From the definition of uT, the two waves are compatible.
8. THE CASE c L - c R Considering the characteristic speeds at the right state, we treat the following two situations separately. 8.1. The Case u R g Ps j Ts As shown in Fig. 10, for a given c L , the state u R determines uniquely different states of the same particle velocity. The first one is the state u A s Ž s A , p R , c L . such that lAp s s lRp s and u A g Bs . If lRp s G lFp s , where u F s Žp R , p R , c L ., then there exists a state u B s Ž s B , p R , c L . g Ps such that lBp s s lRp s . Otherwise, there exist the two states u M s Žp R , p R , c M . g
POLYMER FLOODING WITH HYSTERESIS
TABLE II Solutions for c L ) c R and s L s p L ) p R
223
224
KHALED M. FURATI
FIG. 9. The state u F and the intermediate states.
B P, c L - c M - c R , and u B s Ž s B , s B , c L . g P, s B ) p R such that lM p s lp R s l p s as in Fig. 11. If lRp s - lFp s , consider the case shown in Fig. 11. In that case, we consider a left state u L s Ž s L , p L , c L . such that lLp s ) lRp s and either p L s p R or s L s p R . Let u L s Ž s L , s B , c L . be the state such that hL s hL . Thus, u L can be connected to u L by a stationary contact discontinuity and to u B by a bs-wave along the convex hull with constant p s s B. Since
FIG. 10. The states u A and u B.
225
POLYMER FLOODING WITH HYSTERESIS
FIG. 11. The cd-wave
u B lies on a concave segment, the bs-wave ends with a shock. To have a compatible sequence, we must have the shock speed less than or equal to lBp . By inspecting the convex hull joining u L and u B , this is equivalent to the condition that lLp s G lRp s and the bs-wave is necessarily a shock. bs
However, when lLp s - lRp s , the final speed of the wave u L ª u B is greater than lBp and the above sequence is incompatible. Figure 11 suggests the following lemma for an alternative connection. LEMMA 8.1. Let a be a gi¨ en particle ¨ elocity. Let u L s Ž s L , p L , c L . and u s Ž s B , s B , c L . g P be such that lLp s ) lBp s a . Let u L s Ž s L , s B , c L . be such that hL s h L . Then lLp s - a if and only if there exist the unique states ˆ ˆ ˆ ˆ u J s Ž s J, s J, c J ., c L - c J, and u L s Ž s L , s J, c L ., hL s hL , such that lLp s s J lp s a. B
Proof. The necessity follows from the fact that we have three condiˆ tions for the three unknowns s J, c J, and s L . For sufficiency, notice that for ˆL ˆL J L such u and u , a s l p s ) l p s . Note that if u I s Ž s I , s J, c J . g Bs is the unique state for which lIp s s lJp s , then u I and u J can be connected by a bs-shock of speed lJp . Afterwards, u J can be connected to u M by a p-wave. Hence, the sequence ˆ ps
bs
p
uL ª uI ª uJ ª uM
Ž 19 .
is a compatible sequence that consists of a single discontinuity. For ˆ cd
simplicity we denote such a discontinuity by u L ª u M . Finally, u M and u R
226
KHALED M. FURATI
can be connected by a ps-wave. As a result, we have proven the following lemma. LEMMA 8.2. Let u R g Ps j Ts and u L be such that lFp s , lLp s ) lRp s . Let u L , u B , and u M be as in Fig. 11. If lLp s G lRp s then the solution is gi¨ en by the unique compatible sequence st
p
bs
ps
u L ª u L ª u B ª u M ª uR .
Ž 20 .
Otherwise, the solution is gi¨ en by the unique sequence st
ˆ cd
ps
u L ª u L ª u M ª uR ,
Ž 21 .
cd
ˆ where u L is as in Lemma 8.1 and ª is as in Ž19..
Note that all the waves in Ž21. are discontinuities that travel with speed lRp s . Thus the connection uˆL ª u R in the physical space is just one contact discontinuity. Next, to find the complete set of solutions we consider the following two cases for the left state: Case 1. p L s p R . For such left states, there are different possibilities described in Table III. Note that in the first two solutions, all the flow is secondary with p s p R . However, when s L ) s A and lRp s - lFp s , any compatible sequence passes through some intermediate states with larger p . Solutions for such Riemann problems are given by Lemma 8.2. Note that when s L s s A and lRp s G lFp s , the solution is unique in the physical space but not in the state space and is L1 continuous with respect to the initial data but not pointwise. Case 2. s L s p L ) p R . The different possibilities for this case are shown in Table IV. Note that when lRp s G lFp s both states u B and u L lie on a concave segment of the flow curve and their Buckley]Leverett speeds are less than the particle velocities. Hence, both states can be connected by a primary rarefaction followed by a secondary one with final speed equals to lBb s . Thereafter, a faster ps-wave takes over to u R . When lRp s - lFp s we need to distinguish two possibilities for s L with respect to s B. First, if s L G s B , then u L can be connected to u B with a b-rarefaction of final speed equals to lBb . Since u B g P, any p-wave passes through u B is faster and thus the two waves are compatible. Second, if s L - s B , then the solutions are given by Lemma 8.2. 8.2. The Case u R g Bs Unlike in Section 8.1, where p R G s tan Ž c R ., in this section p R can also be less than s tan Ž c R .. For the first case the same connections in Section 8.1 can be used to connect u L to uT followed by a bs-wave to u R . For the
POLYMER FLOODING WITH HYSTERESIS
227
TABLE III Solutions for c L - c R , uR g Ps j Ts and p L s p R
second case, we need to construct other different compatible connections. To account for all different possibilities, we characterize the left states according to the relation between p R and s tan Ž c R . and some relating states as follows: Case 1. p R G s tan Ž c R .. Let uT s Ž s tans Žp R , c R ., p R , c R . if it exists, otherwise let uT s Žp R , p R , c R .. Note that for u L s Ž s L , p R , c L . g Bs
228
KHALED M. FURATI
TABLE IV Solutions for c L - c R , u R g Ps j Ts and s L s p L ) p R
such that lLp s F lTp s , the solution is the same as in Case 1 of Section 8.1 for s L F s A and shown in Fig. 12a. Otherwise, uT can be connected to u L as in Section 8.1 with uT in place of u R and to u R by a bs-wave. See Figs. 12 and 13 for the different solutions. Case 2. p R - s tan Ž c R .. The solution varies according to the position L of u with respect to some states associated with uT s Ž sT, sT, c R ., sT s
POLYMER FLOODING WITH HYSTERESIS
229
FIG. 12. Solutions for c L - c R , u R g Bs , p R G s tan Ž c R . and p L s p R .
s tan Ž c R ., and u F s Žp R , p R , c R .. Let u A s Ž s A , p R , c L ., u E s Ž s E , sT, c L ., and u D s Ž s D , s D , c L . be the states that satisfy lAp s s lFp and lEp s s lDp s lTp as shown in Fig. 14. Considering the left states with p L s p R , we can characterize three different possibilities listed in Table V. If hL F h A , the solution is as in Case 1 of Section 8.1 with s L F s A and the flow is secondary. When h A - hL F h E , a geometric inspection shows that there exists a unique pair of states u I s Ž s I , s I , c R . and u L s Ž s L , s I , c L . that satisfy the two conditions: lLp s s lIp and h L s h L . These intermediate states are used to form the compatible sequence. Finally, when h E - h L , let u L s Ž s L , s D , c L . be the state such that hL s h L and the solution can be constructed as in Lemma 8.2. Similarly, for left states with s L s p L ) p R , solutions vary according to the values of hL with respect to hE as listed in Table VI.
SUMMARY In this paper we studied the effect of hysteresis on the multiphase flow in porous media. We presented a model for the hysteresis in relative
230
KHALED M. FURATI
FIG. 13. Solutions for c L - c R , u R g Bs , p R G s tan Ž c R . and s L s p L ) p R .
FIG. 14. The states uT, u D , u F and u A .
POLYMER FLOODING WITH HYSTERESIS
231
TABLE V Solutions for c L - c R , uR g Bs , p R - s tan Ž c R . and p L s p R
permeabilities and discussed the corresponding fractional flows. The model is based on considering the initial imbibition of the aqueous phase as the primary flow and any other flow as secondary. We assumed that secondary curves for imbibition and drainage are identical. We analyzed a nonstrictly hyperbolic system modeling polymer flooding with hysteretic fractional
232
KHALED M. FURATI
TABLE VI Solutions for c L - c R , u R g Bs , p R - s tan Ž c R . and s L s p L ) p R
flows and constructed the global solution for the corresponding Riemann problem. Due to the hysteresis and generality of the model, it was not practical to deal with a fixed partition of the state space. Instead, for each case we looked at a suitable criterion to account for all possibilities. Although many possibilities end up having similar solutions, we did not group them together because of their different criteria and to make the construction process easy to follow. Our purpose in presenting the different cases is to provide useful test problems for numerical methods. Including hysteresis produced more complicated flows. This indicates that hysteresis may effect the oil recovery efficiency. In a later paper we will verify numerically the solutions above using second-order Godunov method and present the effect of hysteresis on oil recovery efficiency.
ACKNOWLEDGMENTS I thank Professor John Trangenstein for his help and guidance in this research. His invaluable advice, direction, and encouragement are greatly appreciated. I also thank Dr. Richard Hornung for his useful comments. The author is grateful for the financial support provided by King Fahd University of Petroleum and Minerals.
POLYMER FLOODING WITH HYSTERESIS
233
REFERENCES 1. S. E. Buckley and M. C. Leverett, Mechanisms of fluid displacement in sands, Trans. Amer. Inst. Min. Meta. Engrs. 146 Ž1942., 107]116. 2. J. Colonna, F. Brissaud, and J. L. Millet, Evolution of capillarity and relative permeability hysteresis, SPEJ Trans. 253 Ž1972., 28]38. 3. F. F. Craig, ‘‘The Reservoir Engineering Aspects of Waterflooding,’’ Monograph Series, Society of Petroleum Engineers, Dallas, 1971. 4. R. E. Gladfelter and S. P. Gupta, Effect of fractional flow hysteresis on recovery of tertiary oil, Soc. Pet. Engrs. J. 20 Ž1980., 508]520. 5. E. L. Isaacson, Global Solution of a Riemann problem for a non-strictly hyperbolic system of conservation laws arising in enhanced oil recovery, unpublished. 6. T. Johansen and R. Winther, The solution of the Riemann problem for a hyperbolic system of conservation laws modeling polymer flooding, SIAM J. Math. Anal. 19, No. 3 Ž1988., 541]566. 7. B. L. Keyfitz and H. C. Kranzer, A system of non-strictly hyperbolic conservation laws arising in elasticity, Arch. Rational Mech. Anal. 72 Ž1980., 220]241. 8. J. E. Killough, Reservoir simulation with history-dependent saturation functions, Soc. Pet. Engrs. J. 16 Ž1976., 37]48. 9. R. J. Lenhard and J. C. Parker, A model for hysteretic constitutive relations governing multiphase flow. 2. permeability-saturation relations, Water Resources Res. 23, No. 12 Ž1987., 2197]2206. 10. D. Marchesin, H. B. Medeiros, and P. J. Paes-Leme, A model for two phase flow with hysteresis, Contemp. Math. 60 Ž1987., 89]107. 11. O. A. Oleinik, Uniqueness and a stability of the generalized solution of the Cauchy problem for a quasilinear equation, Uspekhi Mat. Nauk 14 Ž1959., 165]170; English transl., Amer. Math. Soc. Trans. Ser. 2, 33 Ž1964., 285]290. 12. G. A. Pope, The application of fractional flow theory to enhanced oil recovery, Soc. Pet. Engrs. J. 20 Ž1980., 191]205. 13. J. A. Trangenstein and R. B. Pember, The Riemann problem for longitudinal motion in an elastic-plastic bar, SIAM J. Sci. Statist. Comput. 12, No. 1 Ž1991., 180]207. 14. H. J. Welge, A simplified method for computing oil recovery by gas or water drive, Trans. Amer. Inst. Min. Meta. Engrs. 195 Ž1952., 91]98.
206, 205]233 Ž1997.
AY975222
The Solution of the Riemann Problem for a Hyperbolic System Modeling Polymer Flooding with Hysteresis Khaled M. Furati Department of Mathematical Sciences, King Fahd Uni¨ ersity of Petroleum and Minerals, Dhahran 31261, Saudi Arabia Submitted by Barbara Lee Keyfitz Received January 23, 1995
It is well known that multiphase flow in porous media exhibits hysteresis. This is typically modeled by modifying the saturation dependence of the relative permeabilities. In this paper, a model for hysteretic relative permeabilities is built into the polymer flooding model and the analytical solution to the corresponding Riemann problem is constructed. This produces a nonstrictly hyperbolic system of conservation laws with a history-dependent flux function. Because the polymer model without hysteresis possesses Riemann problem solutions that are not monotonic, the introduction of hysteresis necessarily produces structurally different solutions. We show that hysteresis produces more complicated solutions with more fronts and expansions; and removes some nonuniqueness of solutions. Q 1997 Academic Press
1. INTRODUCTION In this paper we solve the Riemann problem for the system of conservation laws t s q x F s 0, Ž 1. t Ž cs . q x Ž cF . s 0, with the flux function F of the form Fs
½
f Ž s, c .
if s s p , t s ) 0
h Ž s, p , c .
otherwise
Ž primary flow . Ž secondary flow .
Ž 2.
and initial conditions
Ž s, p , c . Ž x, 0 . s
½
Ž sL , p L , cL . , Ž sR , p R , cR . ,
x F 0, x ) 0.
Ž 3.
205 0022-247Xr97 $25.00 Copyright Q 1997 by Academic Press All rights of reproduction in any form reserved.
206
KHALED M. FURATI
where, t ) 0, x g R. The variables s, p , c are in I and s F p , where I is the interval w0, 1x. The system Ž1. is used in enhanced oil recovery to model polymer flooding in which oil is displaced by water containing dissolved polymer. The variable s is the saturation of the aqueous phase and c is the concentration of polymer. The flux function F is the fractional flow function for the aqueous phase. By analogy to the elastic-plastic bar problem in w13x, where the history of the material is accounted for by the plastic displacement gradient p , we let p be the hysteresis parameter. For the derivation of Ž1. and physical assumptions see w12x. Experiments have shown that the flow functions are not only functions of current fluid saturations but also functions of the saturation history. Moreover, the solution of the Riemann problem for the polymer flooding model without hysteresis w5x shows cases with nonmonotonic behavior in saturation. The form of F in Ž2. is based on remembering the saturation history through parameterizing the states at which the imbibition process is reversed. The function F s f Ž s, c . describes the flow when the saturation increases monotonically from the irreducible aqueous phase saturation sra to the maximum saturation 1}srl , where srl is the residual liquid phase saturation. We consider this imbibition of the aqueous phase as the primary flow. When this imbibition process is reversed at s s p , a new flow function F s hŽ s, p , c . appears and the flow is called a secondary flow. The main contribution of this paper is incorporating the flow history into the flux function and constructing the unique global solution of the Riemann problem. The Oleinik entropy condition w11x for scalar equations and the generalized Lax entropy condition w7x are used to distinguish the physically meaningful discontinuities. The solution is constructed by piecing together a compatible sequence of elementary waves Žshocks, rarefactions, and contact discontinuities . to connect constant states. Including hysteresis produces more complicated solutions. Another effect of including hysteresis is that it removes some of the nonuniqueness cases arises in w5x, where the solution is unique in the x, t-space but not in the state space. However, for some p R large enough the same nonuniqueness problem persists. Although for certain critical values of the left and right states the solution is not pointwise continuous with respect to the left and right states, it is continuous in L1 norm. Buckley and Leverett w1x constructed the analytical solution for immiscible waterflooding. System Ž1. reduces to their single equation model when c is constant. Pope w12x generalized the fractional flow theory to more complicated floodings and used the method of characteristics to construct the corresponding solutions. The polymer flooding model without hysteresis was analyzed and solved by Isaacson w5x and Johansen and Winther w6x.
POLYMER FLOODING WITH HYSTERESIS
207
The problem is similar to the Riemann problem for vibration of an elastic string considered by Keyfitz and Kranzer w7x. Hysteresis phenomena in multiphase flows in porous media has been recognized by many researchers. For example, Colonna et al. w2x, Killough w8x, Gladfelter and Gupta w4x, and Lenhard w9x. Marchesin et al. w10x used a double-valued flux function to model two-phase flow with hysteresis and constructed, graphically, the complete solution of the associated Riemann problem. However, in this paper we consider a three-component two-phase flow with a more complicated model of hysteresis. The treatment of the history dependence in this paper is similar to that of motion in an elastic-plastic bar studied by Trangenstein and Pember in w13x. The outline of the paper is as follows. In Section 2 we present a brief derivation of the model. In Section 3 we state the mathematical assumptions and properties of the flux functions and analyze the hyperbolicity of the model. In Section 4 we discuss the different wave families. In Sections 5]8 we construct the global solution for the Riemann problem.
2. DERIVATION OF THE MODEL In polymer flooding, a solution of water and polymer is injected into an oil reservoir to push the oil through the rock. The reservoir fluid is considered to consist of three components: oil, water, and polymer. The oil forms its own liquid phase, and water with polymer form the aqueous phase. The polymer transports only in the aqueous phase and does not partition to oil. The purpose of the polymer is to increase the viscosity of water, thereby enhancing its ability to push the oil. System Ž1. is derived under many assumptions. We assume that the reservoir rock is homogeneous, so that the porosity and total permeability of the rock are constant. For the phases, we assume incompressibility and that the liquid phase viscosity is constant, but the aqueous phase viscosity increases as the concentration of polymer increases. Moreover, we assume that there are no diffusive forces, such as capillary pressure, and no gravity forces. The mass conservation equation for water, oil, and polymer, respectively, can be formulated as
f t sa q x¨ a s 0, f t s l q x¨ l s 0, f t Ž csa . q x Ž c¨ a . s 0,
Ž 4.
208
KHALED M. FURATI
where f is the rock porosity, c is the polymer concentration, and sa and s l are the aqueous and liquid phases saturations, respectively. The fluxes ¨ a and ¨ l are the volumetric flow rates of the aqueous and liquid phases, respectively, and are given by Darcy’s law: ¨a s y
Kk ra
ma Ž c .
x P ,
¨l s y
Kk rl
ml
x P.
Ž 5.
Here, P is the fluid pressure, K is the absolute permeability of the rock, k ra and k rl are the relative permeabilities of the aqueous and liquid phases, respectively, and ma , m l are the corresponding viscosities. Next, we use the assumption that the total velocity ¨ t s ¨ a q ¨ l is constant to eliminate the pressure gradient. Adding the first two equations in Ž4., using the assumption that s l q sa s 1 and changing x to f xr¨ t we obtain Ž1. with s s sa and Fs
k ra k ra q k rl m Ž c .
,
Ž 6.
where m s marm l is the viscosity ratio. To complete the description of the conservation laws we need to describe the relative permeability curves. In this paper, we consider the typical relative permeability curves shown in Fig. 1, with the aqueous phase being the wetting phase. The solid curves give the relative permeability values when the saturation s is monotonically increasing from sra . We call such curves the primary curves. When the aqueous phase saturation decreases, the relative permeability is smaller and follows the dashed
FIG. 1. Hysteretic relative permeability curves.
209
POLYMER FLOODING WITH HYSTERESIS
curves. We call such curves the secondary curves. When the saturation increases again, we assume that the relative permeability traces back the secondary curve until it meets the primary curve, then follows the primary curve. We parameterize the secondary curves by the symbol p , which we refer to as the hysteresis parameter. The primary and secondary curves are denoted by k r j Ž s . and k rhj Ž s, p ., respectively, where j s a, l. Thus we have
krj s
½
krj Ž s. k rhj
Ž s, p .
if s s p , t s ) 0, otherwise,
5
j s a, l.
Ž 7.
With the above forms of the relative permeability curves, the fractional flow function F takes the form Ž2. with f Ž s, c . s h Ž s, p , c . s
k ra Ž s . k ra Ž s . q k rl Ž s . m Ž c .
,
h k ra Ž s, p . h k ra Ž s, p . q k rlh Ž s, p . m Ž c .
Ž 8. .
Figure 2 shows typical fractional flow curves. Note that since the polymer increases water viscosity, the ratio m s m Ž c . is an increasing function of c. This implies that the fractional flow function F is decreasing with respect to c. Physically, this means that adding polymer will increase the oil flow rate which is equal to 1 y F.
FIG. 2. Fractional flow curves, c 2 ) c1.
210
KHALED M. FURATI
3. THE MATHEMATICAL DESCRIPTION OF THE MODEL In this section we present the mathematical assumptions and properties of the relative permeabilities and, hence, the flux functions. Accordingly, we analyze the hyperbolicity of the system of conservation laws Ž1.. Let I s w0, 1x, Is s w sra , 1 y s sl x, and Ip s w sp , p x, where sp s max s : hŽ s, p , c . s 04 . We assume that the relative permeability curves satisfy the following conditions, as shown in Fig. 1: A1. j s a, l;
k r j Ž s . g C 2 Ž Is ., k rhj Ž?, p . g C 2 Ž Ip ., and k rhj Ž s, ? . g C 1 Ž Is ., where
A2.
k rhj Ž s, s . s k r j Ž s . ;s g Is , j s a, l;
A3.
kXra ) 0, kYra ) 0, kXrl - 0, and kYrl ) 0 ;s g Is ;
h h h A4. s k ra ) 0, p k ra - 0, s s k ra ) 0, s k rlh - 0, p k rlh ) 0, and s s k rlh ) 0 in Ip = Is ; A5. kXra Ž sra . q m Ž c . kXrl Ž sra . F 0 and kXraŽ1 y srl . q m Ž c . kXrl Ž1 y srl . G 0 ;c g I; hŽ A6. s k ra sp , p . q m Ž c . s k rlh Ž sp , p . F 0 ;c g I.
Assumption A5 is a sufficient Žnot necessary . condition for the primary fractional flow function, f Ž., c ., to be convex in a neighborhood of s s sra and concave in a neighborhood of s s 1 y srl . This implies the existence of an inflection point of f. h We assumed p k ra - 0 and p k rlh ) 0 so that the secondary relative permeability curves are disjoint for distinct values of p . Assumption A6 is a sufficient Žnot necessary . condition for the secondary fractional flow function hŽ s, p , c . to be convex in some neighborhood of the irreducible saturation sp . With direct examinations of the derivatives, the above assumptions imply the following properties for the flux functions f and h as shown in Fig. 2: P1. f Ž s, c . is a C 2 function in Is = I, and f Ž sra , c . s 0 and f Ž1 y srl , c . s 1 ;c g I; P2. s f Ž s, c . ) 0 and c f Ž s, c . - 0 in Is = I; P3. For each c g I, there is an inflection point, s I s s I Ž c . g Ž sra , 1 y srl ., such that s s f Ž s, c . ) 0 for s - s I ; P4. hŽ s, p , c . is a C 2 function in s g Ip , C 1 in p g Is and C 2 in c g I; P5. hŽ sp , p , c . s 0 and hŽ s, s, c . s f Ž s, c . for each p , s g Is and c g I;
POLYMER FLOODING WITH HYSTERESIS
P6. = I;
211
s hŽ s, p , c . ) 0, p hŽ s, p , c . - 0 and c hŽ s, p , c . F 0 in Ip = Is
P7. s s hŽ sp , p , c . ) 0 for each c g I and p g Is ; P8. s f Žp , c . F sy hŽp , p , c ., for each c g I, where sy is the left derivative. Note that the assumptions on relative permeability curves imply existence Žproperty P3. but not necessarily uniqueness of the inflection points of f. However, for most applications, experimental data indicates that the fraction flow function does not have more than one inflection point. Therefore, for each c g I we assume that s I Ž c . is the only point of inflection, and consequently, f Ž., c . is convex in Ž sra , s I Ž c .. and concave in Ž s I Ž c ., 1 y srl .. Property P6 implies that the secondary fractional flow curves for a fixed c are disjoint for distinct values of p . Also, it implies that for a fixed p the curves are disjoint for distinct values of c. Unlike the primary fractional flow function, f Ž s, c ., assumptions on hŽ k ra s, p . and k rlh Ž s, p . do not guarantee the existence of inflection points for the secondary fractional flow function, h. Moreover, the assumptions do not guarantee the uniqueness of such points for the curves of h, as for f. The assumptions only guarantee property P7 which implies that all the curves of h are convex in some neighborhood of sp . We assume that for each c g I and p g Is , hŽ s, p , c . has at most one inflection point, s I Žp , c ., for s g Ip . For the primary flow, system Ž1. can be written in the quasilinear form
t cs q A x cs s 0,
As
s f 0
c f . frs
Ž 9.
Since the eigenvalues of A are real, the system of conservation laws is hyperbolic and has the two characteristic speeds and the corresponding eigenvectors:
l b s s f ,
rb s
1 , 0
lp s
f s
,
rp s
c f . l p y lb
Ž 10 .
The characteristic speed l b is the Buckley]Leverett rarefaction wave speed, which is equal to the slope of f. The second characteristic speed, l p , is the aqueous particle velocity, which is equal to the slope of the chord from the origin to the fractional flow function.
212
KHALED M. FURATI
If the flow is secondary, F s hŽ s, p , c ., then the system of conservation laws Ž1. can be closed by adding the constraint t p s 0 and written in the quasilinear form s s t c q A˜x c s 0, p p
A˜ s
s h 0 0
c h hrs 0
p h 0 . 0
Ž 11 .
The characteristic speeds and their corresponding eigenvectors are
l b s s s h,
1 rb s s 0 , 0
lp s 0,
lps s
h s
,
p h rp s 0 yl b s
c h rp s s l p s y lb s , 0
Ž 12 .
Note that all the characteristic speeds are nonnegative. This implies that all the wave families are either stationary or traveling to the right. For the border states u s Žp , p , c ., we define the Buckley]Leverett characteristic speeds by l b Žp , p , c . s s f Žp , c . and l b s Žp , p , c . s sy hŽp , p , c .. The type of the flow Žsecondary or primary. determines which velocity is to be used in determining the Riemann problem solutions. Note that for a given c the Buckley]Leverett and particle velocities might coalesce as shown in Fig. 3. However, this is not always the case, as shown in Fig. 4. We let utan Ž c . and utans Žp , c . be the states such that l b Ž utan . s l p Ž utan . and l b s Ž utans . s l p s Ž utans .. For the fractional flow curves described above we have the following results. LEMMA 3.1. For each c g I, f Ž., c . has at most one tangential state utan Ž c .. If such a state exists then l b Ž s, c . ) l p Ž s, c . for sr a - s - s tan Ž c . and l b Ž s, c . - l p Ž s, c . for s tan Ž c . - s. Proof. Note that by definition, the tangential states utan are the solutions of the equation l b y l p s s f y frs s 0. Given c g I, let g Ž s . s s s f Ž s, c . y f Ž s, c . for s g Ž sra , 1 y srl .. Since g X Ž s . s s s s f Ž s, c . and f Ž s, c . is assumed to have a unique inflection point u I s Ž s I Ž c ., c ., g Ž s . has only one local maximum at s I in Ž sra , 1 y srl .. Therefore, g s 0 can have at most one solution since g Ž sra . G 0. This proves the uniqueness because any s tan is by definition a solution of g s 0. The second part follows from the observation that g changes from positive to negative at s s s tan .
POLYMER FLOODING WITH HYSTERESIS
213
FIG. 3. Tangential states.
The next lemma shows for which values of p the tangential state utans might exist. LEMMA 3.2. For a gi¨ en c g I and p - s tan Ž c ., l b s Ž s, p , c . ) l p s Ž s, p , c . for s g Ž sp , p x. Proof. For a fixed c g I and p g Ž sra , s tan Ž c .., let g Ž s . s s s hŽ s, p , c . y hŽ s, p , c .. From Property P7 of the secondary fractional flow curves, we
FIG. 4. No tangential states.
214
KHALED M. FURATI
have g X Ž sp . s sp s s hŽ sp , p , c . ) 0. In addition, g satisfies the following inequalities: g Ž sp . s sp s h Ž sp , p , c . G 0, g Ž p . s p syh Ž p , p , c . y h Ž p , p , c .
Ž 13 .
) p s f Ž p , c . y f Ž p , c . ) 0. From those inequalities and the assumption that h has at most one inflection state, it follows that g is positive for s g Ž sp , p x. This proves the result. The above lemma implies that the tangential state, utans s Ž s tans Žp , c ., p , c ., may only exists if p G s tan Ž c .. The following lemma is about the uniqueness of utans . LEMMA 3.3. Gi¨ en c g I and p G s tan Ž c ., there exists at most one tangential state utans Žp , c . s Ž s tans Žp , c ., p , c ., sra - s tans - p . Proof. As in the proof of Lemma 3.1, the result follows from the uniqueness of the inflection point of h. In the above proofs, note that the uniqueness of the tangential states follows from the uniqueness of the inflection points. However, this does not guarantee their existence. In this paper we assume the existence of utan but not necessarily utans . However, if for p G s tan the state utans does not exist we have the following lemma. LEMMA 3.4. Gi¨ en c g I and p G s tan Ž c ., if there is no tangential state u , then the state up s Žp , p , c . satisfies tans
l b s Ž up . y l p s Ž up .
l b Ž up . y l p Ž up . F 0.
Ž 14 .
Proof. Again, for fixed c and p consider the function g Ž s . s s s h y h. By assumption, g s 0 has no solution for s g Ž sp , p .. But g Ž sp . G 0 and g X Ž sp . ) 0. Therefore, g Ž s . ) 0 for sp - s - p . By continuity, this implies that g Žp . G 0, and thus the first bracket in the inequality Ž14. is nonnegative. The second bracket is nonpositive by Lemma 3.1. This proves the result. To summarize, for each c g I we assume that there exists a state utan s Ž s tan Ž c ., c . such that l p Ž utan . s l b Ž utan .. Moreover, for each c g I and p g Is , if p - s tan Ž c . then l p s Ž u. - l b s Ž u. for all u s Ž s, p , c . with sp - s F p . If p G s tan Ž c ., then either there exists a state utans s Ž s tans Žp , c ., p , c . such that l p s Ž utans . s l b s Ž utans . or the state Žp , p , c . is
POLYMER FLOODING WITH HYSTERESIS
215
a transitional state across which the relation between the two speeds is reversed. Both states utan and utans are unique. We partition the set of all states into the following sets: P s u : l p Ž u . ) lb Ž u . 4 ,
B s u : lb Ž u . ) l p Ž u. 4 ,
T s u : l p Ž u . s lb Ž u . 4 , Ps s u : l p s Ž u . ) l b s Ž u . 4 ,
Bs s u : l b s Ž u . ) l p s Ž u . 4 ,
Ts s u : l p s Ž u . s l b s Ž u . 4 . We will use these sets to characterize the different cases for the solutions.
4. WAVE FAMILIES In this section we determine all wave families associated with the Riemann problem Ž1. ] Ž3.. Those wave families consist of the elementary waves: rarefactions, shocks, and contact discontinuities. To find all the possible waves we need to examine the nonlinearity of the characteristic fields, the conditions for discontinuities and the integral curves of the characteristic eigenvectors. For the nonlinearity we need to check the projection of the characteristic speed gradient into the corresponding eigenvector. For our problem we have =l p ? r p s =p l p s ? r p s s =p lp ? rp s 0,
=l b ? r b s s s f ,
=p l b s ? r b s s s s h, where = s Ž s , c . and =p s Ž s , c , p .. Consequently, the p-field, ps-field, and p-field are linearly degenerate and therefore the corresponding families consist only of contact discontinuities. On the other hand, due to the inflection points of f and h, the b-family and bs-family are not typically genuinely nonlinear. We expect both families to consist of rarefactions and shocks. As we know, any discontinuity solution of speed s must satisfy the Rankine]Hugoniot conditions. For the primary flow those conditions take the form fRyfL c f yc f R
R
L
L
s
sR y sL s. c s y cL sL R R
Ž 15 .
216
KHALED M. FURATI
One solution for Ž15. is s s l p s frs. Another solution is determined by taking c to be constant and s s w f xrw s x, where w?x denotes the jump of the argument across the discontinuity. For the secondary flow, the Rankine]Hugoniot conditions take the form sR y sL hR y h L R L L s cR sR y cL sL s . c h yc h 0 pRypL R
Ž 16 .
The conditions in Ž16. have the three solutions: 1. 2. 3.
s s w h xrw s x and w c x s wp x s 0; s s l p s s hrs and wp x s 0; s s 0 and w h x s w c x s 0.
To search for rarefactions we need to consider the integral curves of the eigenvectors r b and r b s . The integral curves of r b are the curves where c is constant. Similarly, the integral curves of r b s are the curves where c and p are constant. Accordingly, the wave families for the Riemann problem for system Ž1. ] Ž3. can be classified as follows. b-wa¨ es and bs-wa¨ es. The two families correspond to the characteristic speeds l b and l b s , respectively, and consist of rarefactions and shocks. Note that c is a Riemann invariant of both families and p is a Riemann invariant of the bs-family. Thus, from the solutions of the Rankine]Hugoniot conditions, the wave curves for the shocks and rarefactions of both families are identical. The b-family curves are the ones along which c is constant. The bs-family curves are the ones along which both c and p are constants. As a result, both families can be considered as solutions of the scalar equation t s q x F s 0, F s f, h. To distinguish the physically meaningful shocks, we need to impose an entropy condition. For our purpose, a shock is considered admissible if it satisfies the Oleinik condition w11x for scalar equations. This condition requires that F Ž uR . y F Ž u . sR y s
Fss
F Ž uR . y F Ž u L . sR y sL
Ž 17 .
for all s between s L and s R . In other words, the chord from Ž s R , F Ž u R .. to Ž s L , F Ž u L .. must lie above the chord from Ž s R , F Ž u R .. to Ž s, F Ž u.. for any s in between. Since both families are not genuinely nonlinear, the solution may consist of a combination of a shock and rarefaction. For the rarefactions, the
POLYMER FLOODING WITH HYSTERESIS
217
characteristic speeds should be monotonically increasing from left to right. This monotonicity and Oleinik condition are equivalent to constructing the concave hull for s L ) s R and the convex hull for s L - s R as shown in Fig. 5. A straight line segment represents a shock and a curve segment represents a rarefaction. The shock speed is equal to the slope of the straight line. In petroleum engineering literature, the above solutions are usually explained in terms of Welge tangents w3, 14x. The Welge tangents are equivalent to the Oleinik chords for such problems. We use k-wave, k s b, bs, to denote a k-shock, k-rarefaction, or any consecutive combination of them. Note that to have a b-wave we must have s L ) s R so it corresponds to a primary flow. p-wa¨ es and ps-wa¨ es. The two families correspond to the characteristic speed l p and l p s , respectively. The two families are linearly degenerate and consist of contact discontinuities since l p and l p s are Riemann invariants. The wave curves for the p-family are the curves c s cŽ s . along which l p is constant. The slope of such curves depends on the relation between l p and l b since dcrds s Ž l p y l b .rf c . The same argument is true for the ps-family curves c s cŽ s, p .. The states u L s Ž s L , c L . and u R s Ž s R , c R ., s L ) s R , can be connected by a contact discontinuity if l p Ž u L . s l p Ž u R . and they both belong to either P j T or B j T. The last condition is to fulfill the generalized Lax’s entropy condition w7x which requires that exactly one of the b-characteristics leave the discontinuity. With p L s p R , similar results are true for the ps-family.
FIG. 5. Concave and convex hulls.
218
KHALED M. FURATI
st-wa¨ es. This family consists of stationary contact discontinuities satisfying Rankine]Hugoniot conditions with speed s s 0. Across this discontinuity, w h x s w c x s 0.
5. GLOBAL SOLUTION OF THE RIEMANN PROBLEM The global solution to the Riemann problem can be constructed graphically using the wave families mentioned above. The uniqueness of the solutions can easily be verified by imposing the compatibility condition, i.e., the initial speed of each wave is greater than or equal to the final speed of the preceding wave. We start with the following remarks. Remark 1 ŽNotation.. Although the problem during a primary flow is independent of the history parameter p , we write any state u s Ž s, c . as u s Ž s, s, c .. Thus, for primary and secondary flows we can let u L s Ž s L , p L , c L . and u R s Ž s R , p R , c R . be the left and right states of the Riemann problem. For the rest of this paper we use hŽ s, s, c . to denote f Ž s, c .. For simplicity, we use h a to denote hŽ u a . and l a to denote lŽ u a .. k
Moreover, we use ª to denote a k-wave. Remark 2 ŽReduction.. For any given two states u L s Ž s L , p L , c L . and u s Ž s R , p R , c R . let u Q s Žp R , p R , c L .. Accordingly, there are two cases as shown in Fig. 6. If hQ G hL , then there is a unique state u A s Ž s A , p R , c L . such that h A s hL , and u L can be connected to u A by a R
FIG. 6. Left: hQ G hL . Right: hQ - h L .
POLYMER FLOODING WITH HYSTERESIS
219
st
stationary wave u L ª u A . Otherwise, there is a unique state u P s Ž s P , s P , c L ., s P F s L , such that h P s hL and u L can be connected to u P by st
a stationary wave u L ª u P . Thus, the problem can be reduced to solving either for u L s Ž s L , p R , c L . or u L s Ž s L , s L , c L ., p R F s L , for any arbitrary u R . To account for all possible cases, we treat the three different cases separately: c L s c R , c L ) c R , and c L - c R . Solutions of the Riemann problem will be composed of compatible sequences of the st-waves, bswaves, b-waves, p-waves, and ps-waves.
6. THE CASE c L s c R ŽBUCKLEY]LEVERETT EQUATION WITH HYSTERESIS. When c L s c R , the system Ž1. degenerates to the Buckley]Leverett equation and the solution consists of b-waves and bs-waves. For this special case we omit the concentration argument and write any state u as u s Ž s, p .. From Remark 2, we only need to consider the following two cases: 1. u L s Ž s L , p R .. In this case, u L and u R lie on the same secondary curve h s hŽ s, p R .. The solution is a bs-wave determined by the concave hull if s R - s L or by the convex hull if s R ) s L . Thus, we have the wave bs
u L ª u R across which p is constant and is equal to p R . See Fig. 7 for the case s R - s L .
FIG. 7. Solution for c L s c R , p L s p R and s L ) s R .
220
KHALED M. FURATI
2. u L s uŽ s L , s L . and p R - s L . Consider first the intermediate states u s Ž s, s ., p R F s - s L . The left state u L can be connected to any of those states by a unique b-wave determined by the concave hull of the primary curve joining the two states. In particular, for the state Žp R , p R . we have b the connection u L ª Žp R , p R .. On the other hand, the state Žp R , p R . can be connected to the right state u R by a unique bs-wave determined by the bs concave hull of the secondary curve. Thus, we have Žp R , p R . ª u R . If the final speed of the b-wave is less than the initial speed of the bs-wave, then the two waves are compatible and can be composed to give the unique b bs solution u L ª Žp R , p R . ª u R as shown in Fig. 8a. Note that, in this case, the union of the two concave hulls is also concave. Otherwise, we need to construct the concave hull of the union of the two hulls as in Fig. 8b. The saturation s M satisfies f X Ž sM . s
f Ž sM . y hŽ sR , p R . sM y sR
.
Ž 18 .
Thus, the solution starts with a rarefaction from u L to u M along the primary curve, followed by a shock to u R . This shock is actually a composition of a b-shock and a bs-shock traveling with the same speed. 7. THE CASE c L ) c R Considering the left state u L s Ž s L , p L , c L ., it is enough to solve for the following two cases as pointed out in Remark 2.
FIG. 8. Solution for c L s c R and s L s p L ) p R .
POLYMER FLOODING WITH HYSTERESIS
221
7.1. The Case p L s p R The different possibilities for this Riemann problem are listed in Table I. Note that when u L g Bs j Ts , u R g Ps , and lLp s s lRp s , the solution is not unique in the state space but it is unique in the physical space. Similar nonuniqueness occurs when u L g Ps , u R g Ps , and lRp s s lTp s , where
TABLE I Solutions for c L ) c R and p L s p R
222
KHALED M. FURATI
uT s Ž sT, p R , c L . g Ts is the tangential state as in Fig. 3. Moreover, both special cases are examples of solutions that are not pointwise continuous with respect to the initial data but continuous in L1 norm. Note also that all the flow is secondary for all different cases. 7.2. The Case s L s p L ) p R With such initial conditions, part of the flow is primary and the other is secondary. This is because the saturations of some intermediate states exceed the history parameter p s p R . Consider the following two cases listed in Table II. Case 1. u L g B j T. As shown in Fig. 9, there are two possibilities for u L with respect to u F s Žp R , p R , c R .. For both, any compatible sequence consists of a contact discontinuity leaving the left state followed by a b-wave or a bs-wave. Case 2. u L g P. The slowest wave that leaves u L is a b-wave. How far this wave can be extended depends on the position of uT s Ž sT, sT, c L . g T with respect to u R . If sT F p R , then both u L and u A s Žp R , p R , c L . lie on the concave part of the primary curve and can be connected by the rarefaction b
u L ª u A with final speed equal to lAb . Afterwards, u A can be connected to u R as in Section 7.1. There, the wave that leaves u A is either a bs-rarefaction or a ps-wave. The initial speed of both waves is greater than or equal to lAb . Thus, the whole connection is compatible. When sT ) p R , uT considered as the left state, can be connected to u R as in Case 1 in which the initial speed of the p-wave leaving uT is equal to lTp . Moreover, u L and uT can be connected through a b-rarefaction. From the definition of uT, the two waves are compatible.
8. THE CASE c L - c R Considering the characteristic speeds at the right state, we treat the following two situations separately. 8.1. The Case u R g Ps j Ts As shown in Fig. 10, for a given c L , the state u R determines uniquely different states of the same particle velocity. The first one is the state u A s Ž s A , p R , c L . such that lAp s s lRp s and u A g Bs . If lRp s G lFp s , where u F s Žp R , p R , c L ., then there exists a state u B s Ž s B , p R , c L . g Ps such that lBp s s lRp s . Otherwise, there exist the two states u M s Žp R , p R , c M . g
POLYMER FLOODING WITH HYSTERESIS
TABLE II Solutions for c L ) c R and s L s p L ) p R
223
224
KHALED M. FURATI
FIG. 9. The state u F and the intermediate states.
B P, c L - c M - c R , and u B s Ž s B , s B , c L . g P, s B ) p R such that lM p s lp R s l p s as in Fig. 11. If lRp s - lFp s , consider the case shown in Fig. 11. In that case, we consider a left state u L s Ž s L , p L , c L . such that lLp s ) lRp s and either p L s p R or s L s p R . Let u L s Ž s L , s B , c L . be the state such that hL s hL . Thus, u L can be connected to u L by a stationary contact discontinuity and to u B by a bs-wave along the convex hull with constant p s s B. Since
FIG. 10. The states u A and u B.
225
POLYMER FLOODING WITH HYSTERESIS
FIG. 11. The cd-wave
u B lies on a concave segment, the bs-wave ends with a shock. To have a compatible sequence, we must have the shock speed less than or equal to lBp . By inspecting the convex hull joining u L and u B , this is equivalent to the condition that lLp s G lRp s and the bs-wave is necessarily a shock. bs
However, when lLp s - lRp s , the final speed of the wave u L ª u B is greater than lBp and the above sequence is incompatible. Figure 11 suggests the following lemma for an alternative connection. LEMMA 8.1. Let a be a gi¨ en particle ¨ elocity. Let u L s Ž s L , p L , c L . and u s Ž s B , s B , c L . g P be such that lLp s ) lBp s a . Let u L s Ž s L , s B , c L . be such that hL s h L . Then lLp s - a if and only if there exist the unique states ˆ ˆ ˆ ˆ u J s Ž s J, s J, c J ., c L - c J, and u L s Ž s L , s J, c L ., hL s hL , such that lLp s s J lp s a. B
Proof. The necessity follows from the fact that we have three condiˆ tions for the three unknowns s J, c J, and s L . For sufficiency, notice that for ˆL ˆL J L such u and u , a s l p s ) l p s . Note that if u I s Ž s I , s J, c J . g Bs is the unique state for which lIp s s lJp s , then u I and u J can be connected by a bs-shock of speed lJp . Afterwards, u J can be connected to u M by a p-wave. Hence, the sequence ˆ ps
bs
p
uL ª uI ª uJ ª uM
Ž 19 .
is a compatible sequence that consists of a single discontinuity. For ˆ cd
simplicity we denote such a discontinuity by u L ª u M . Finally, u M and u R
226
KHALED M. FURATI
can be connected by a ps-wave. As a result, we have proven the following lemma. LEMMA 8.2. Let u R g Ps j Ts and u L be such that lFp s , lLp s ) lRp s . Let u L , u B , and u M be as in Fig. 11. If lLp s G lRp s then the solution is gi¨ en by the unique compatible sequence st
p
bs
ps
u L ª u L ª u B ª u M ª uR .
Ž 20 .
Otherwise, the solution is gi¨ en by the unique sequence st
ˆ cd
ps
u L ª u L ª u M ª uR ,
Ž 21 .
cd
ˆ where u L is as in Lemma 8.1 and ª is as in Ž19..
Note that all the waves in Ž21. are discontinuities that travel with speed lRp s . Thus the connection uˆL ª u R in the physical space is just one contact discontinuity. Next, to find the complete set of solutions we consider the following two cases for the left state: Case 1. p L s p R . For such left states, there are different possibilities described in Table III. Note that in the first two solutions, all the flow is secondary with p s p R . However, when s L ) s A and lRp s - lFp s , any compatible sequence passes through some intermediate states with larger p . Solutions for such Riemann problems are given by Lemma 8.2. Note that when s L s s A and lRp s G lFp s , the solution is unique in the physical space but not in the state space and is L1 continuous with respect to the initial data but not pointwise. Case 2. s L s p L ) p R . The different possibilities for this case are shown in Table IV. Note that when lRp s G lFp s both states u B and u L lie on a concave segment of the flow curve and their Buckley]Leverett speeds are less than the particle velocities. Hence, both states can be connected by a primary rarefaction followed by a secondary one with final speed equals to lBb s . Thereafter, a faster ps-wave takes over to u R . When lRp s - lFp s we need to distinguish two possibilities for s L with respect to s B. First, if s L G s B , then u L can be connected to u B with a b-rarefaction of final speed equals to lBb . Since u B g P, any p-wave passes through u B is faster and thus the two waves are compatible. Second, if s L - s B , then the solutions are given by Lemma 8.2. 8.2. The Case u R g Bs Unlike in Section 8.1, where p R G s tan Ž c R ., in this section p R can also be less than s tan Ž c R .. For the first case the same connections in Section 8.1 can be used to connect u L to uT followed by a bs-wave to u R . For the
POLYMER FLOODING WITH HYSTERESIS
227
TABLE III Solutions for c L - c R , uR g Ps j Ts and p L s p R
second case, we need to construct other different compatible connections. To account for all different possibilities, we characterize the left states according to the relation between p R and s tan Ž c R . and some relating states as follows: Case 1. p R G s tan Ž c R .. Let uT s Ž s tans Žp R , c R ., p R , c R . if it exists, otherwise let uT s Žp R , p R , c R .. Note that for u L s Ž s L , p R , c L . g Bs
228
KHALED M. FURATI
TABLE IV Solutions for c L - c R , u R g Ps j Ts and s L s p L ) p R
such that lLp s F lTp s , the solution is the same as in Case 1 of Section 8.1 for s L F s A and shown in Fig. 12a. Otherwise, uT can be connected to u L as in Section 8.1 with uT in place of u R and to u R by a bs-wave. See Figs. 12 and 13 for the different solutions. Case 2. p R - s tan Ž c R .. The solution varies according to the position L of u with respect to some states associated with uT s Ž sT, sT, c R ., sT s
POLYMER FLOODING WITH HYSTERESIS
229
FIG. 12. Solutions for c L - c R , u R g Bs , p R G s tan Ž c R . and p L s p R .
s tan Ž c R ., and u F s Žp R , p R , c R .. Let u A s Ž s A , p R , c L ., u E s Ž s E , sT, c L ., and u D s Ž s D , s D , c L . be the states that satisfy lAp s s lFp and lEp s s lDp s lTp as shown in Fig. 14. Considering the left states with p L s p R , we can characterize three different possibilities listed in Table V. If hL F h A , the solution is as in Case 1 of Section 8.1 with s L F s A and the flow is secondary. When h A - hL F h E , a geometric inspection shows that there exists a unique pair of states u I s Ž s I , s I , c R . and u L s Ž s L , s I , c L . that satisfy the two conditions: lLp s s lIp and h L s h L . These intermediate states are used to form the compatible sequence. Finally, when h E - h L , let u L s Ž s L , s D , c L . be the state such that hL s h L and the solution can be constructed as in Lemma 8.2. Similarly, for left states with s L s p L ) p R , solutions vary according to the values of hL with respect to hE as listed in Table VI.
SUMMARY In this paper we studied the effect of hysteresis on the multiphase flow in porous media. We presented a model for the hysteresis in relative
230
KHALED M. FURATI
FIG. 13. Solutions for c L - c R , u R g Bs , p R G s tan Ž c R . and s L s p L ) p R .
FIG. 14. The states uT, u D , u F and u A .
POLYMER FLOODING WITH HYSTERESIS
231
TABLE V Solutions for c L - c R , uR g Bs , p R - s tan Ž c R . and p L s p R
permeabilities and discussed the corresponding fractional flows. The model is based on considering the initial imbibition of the aqueous phase as the primary flow and any other flow as secondary. We assumed that secondary curves for imbibition and drainage are identical. We analyzed a nonstrictly hyperbolic system modeling polymer flooding with hysteretic fractional
232
KHALED M. FURATI
TABLE VI Solutions for c L - c R , u R g Bs , p R - s tan Ž c R . and s L s p L ) p R
flows and constructed the global solution for the corresponding Riemann problem. Due to the hysteresis and generality of the model, it was not practical to deal with a fixed partition of the state space. Instead, for each case we looked at a suitable criterion to account for all possibilities. Although many possibilities end up having similar solutions, we did not group them together because of their different criteria and to make the construction process easy to follow. Our purpose in presenting the different cases is to provide useful test problems for numerical methods. Including hysteresis produced more complicated flows. This indicates that hysteresis may effect the oil recovery efficiency. In a later paper we will verify numerically the solutions above using second-order Godunov method and present the effect of hysteresis on oil recovery efficiency.
ACKNOWLEDGMENTS I thank Professor John Trangenstein for his help and guidance in this research. His invaluable advice, direction, and encouragement are greatly appreciated. I also thank Dr. Richard Hornung for his useful comments. The author is grateful for the financial support provided by King Fahd University of Petroleum and Minerals.
POLYMER FLOODING WITH HYSTERESIS
233
REFERENCES 1. S. E. Buckley and M. C. Leverett, Mechanisms of fluid displacement in sands, Trans. Amer. Inst. Min. Meta. Engrs. 146 Ž1942., 107]116. 2. J. Colonna, F. Brissaud, and J. L. Millet, Evolution of capillarity and relative permeability hysteresis, SPEJ Trans. 253 Ž1972., 28]38. 3. F. F. Craig, ‘‘The Reservoir Engineering Aspects of Waterflooding,’’ Monograph Series, Society of Petroleum Engineers, Dallas, 1971. 4. R. E. Gladfelter and S. P. Gupta, Effect of fractional flow hysteresis on recovery of tertiary oil, Soc. Pet. Engrs. J. 20 Ž1980., 508]520. 5. E. L. Isaacson, Global Solution of a Riemann problem for a non-strictly hyperbolic system of conservation laws arising in enhanced oil recovery, unpublished. 6. T. Johansen and R. Winther, The solution of the Riemann problem for a hyperbolic system of conservation laws modeling polymer flooding, SIAM J. Math. Anal. 19, No. 3 Ž1988., 541]566. 7. B. L. Keyfitz and H. C. Kranzer, A system of non-strictly hyperbolic conservation laws arising in elasticity, Arch. Rational Mech. Anal. 72 Ž1980., 220]241. 8. J. E. Killough, Reservoir simulation with history-dependent saturation functions, Soc. Pet. Engrs. J. 16 Ž1976., 37]48. 9. R. J. Lenhard and J. C. Parker, A model for hysteretic constitutive relations governing multiphase flow. 2. permeability-saturation relations, Water Resources Res. 23, No. 12 Ž1987., 2197]2206. 10. D. Marchesin, H. B. Medeiros, and P. J. Paes-Leme, A model for two phase flow with hysteresis, Contemp. Math. 60 Ž1987., 89]107. 11. O. A. Oleinik, Uniqueness and a stability of the generalized solution of the Cauchy problem for a quasilinear equation, Uspekhi Mat. Nauk 14 Ž1959., 165]170; English transl., Amer. Math. Soc. Trans. Ser. 2, 33 Ž1964., 285]290. 12. G. A. Pope, The application of fractional flow theory to enhanced oil recovery, Soc. Pet. Engrs. J. 20 Ž1980., 191]205. 13. J. A. Trangenstein and R. B. Pember, The Riemann problem for longitudinal motion in an elastic-plastic bar, SIAM J. Sci. Statist. Comput. 12, No. 1 Ž1991., 180]207. 14. H. J. Welge, A simplified method for computing oil recovery by gas or water drive, Trans. Amer. Inst. Min. Meta. Engrs. 195 Ž1952., 91]98.