A study of spatial approximations for simulating fluid displacements in petroleum reservoirs

c

N

number

of grid

intemaIs

or eb

men&,

I.3 phase dispersion tensor defined by (6) or coefficients cuf (A.21 approximate which the artificial dispersion term, *raetiwal AQWof phase 611* CQM* f ponent, k%Ega&ies sf ~~rnp~n~~~ i in olis and gas phases, coeRkient due to capiltary press g sure wr artificial dispsrsion, parameters in artificial dispersisrr terms, ~~~it~t~~~~~ @onsfant, depth, absofute permeability, relative permeability of phase, characteristic length, number of nonaqueow ci”)mponents,

P PC capiIJary pressure difference, J?e Peclet number, Llcu, PC” effective value crf PecJet number in numericai ~~~roximati~n~ s owzatf vd-ume fraczion of camponeat or phk%e satnratirsn, saturation at kading edge of dispJacement front, bubble and dew point coned position, B time, T cw%cients ol [A_2) wkich apm proximate ~~~~~c~ion term, coefkients in approximation of pressure equation (A.31, u phase or totaJ velocity, xi,Yi mass fraction Of component i in cd and gas phzw, respectiveIy, JXeSSWtS,

4

L.C. Young, Spatial approximations for simulating fluid displacements

x, y

Cartesian coordinates, grid interval or element size, dispersivity, Y function describing the effect of pressure on density, Kronecker delta, &j & maximum numerical error, coordinate which 71 dimensionless moves at the frontal advance rate, total fluid mobility, A, + A, + A,, phase mobilities, k,/F, viscosity of phase, dimensionless coordinate, dimensionless coordinate location of front,

AX,AY a

p PO T 4

phase or component density, reference value of density, dimensionless time, fractional porosity,

Subscripts and superscripts c denotes connate water, g value associated with gas phase, e value for longitudinal direction, o value associated with oil phase, r denotes value for residual and nonflowing conditions, t value for transverse direction, w value associated with water phase, 4 Y value for X or Y direction.

1. Introduction

All secondary and tertiary oil recovery processes involve the injection of fluids to displace oil towards production wells. To describe these processes the distinction is made between immiscible and miscible flow. Strictly immiscible flow is the flow of two or more phases with constant phase compositions, while strictly miscible flow is the flow of two or more components within a single phase. Commercial recovery methods involve the flow of multiple phases at reservoir conditions. Secondary recovery methods, e.g. waterflooding, involve primarily immiscible flow, while tertiary recovery methods, e.g. CO, flooding, involve aspects of both miscible and immiscible flow. This difference causes tertiary recovery processes to be more difficult to simulate numerically than secondary processes. When a solvent is used to displace reservoir fluid, the reservoir fluid and solvent are said to be first contact miscible if the fluids remain single phase when mixed in all proportions. Solvents of commercial interest, e.g., CO* or hydrocarbon gases, are not generally contact miscible with reservoir oils. Instead, they develop a solvent-like character as a result of component transfer between phases [l]. This displacement is termed ‘multiple-contact-miscible’ (MCM). These processes involve not only miscible and immiscible flow, but also the transition from a single nonaqueous phase (oil or gas) ahead of the displacement front to two nonaqueous phases (oil and gas) behind the displacement front. This phase transition adds an additional mathematical complexity to the problem. Phase transitions also occur in other tertiary recovery processes, e.g., surfactant flooding and steam flooding. This paper discusses numerical methods for simulating fluid displacements in petroleum reservoirs. The objective is to develop methods for tertiary processes, i.e., multiphase multicomponent flows. Several specific difficulties associated with the simulation of tertiary recovery processes have previously been discovered. The most widely recognized problem is the grid orientation sensitivity which is observed for adverse mobility ratio displacements [2].

L.C. Young, Spatial approximationsfor simulating fluid displacements

5

Although several solutions to this problem have been claimed, the only method in widespread use is the upstream weighted nine-point difference method [2]. Most studies of this problem have considered the orientation sensitivity but not convergence of the numerical ap proximation. The orientation sensitivity has been attributed to the functional form of the dispersion introduced by the use of upstream weighting in a five-point difference scheme [3]. When the physical dispersion term is included in a central difference approximation, the grid orientation sensitivity is not observed for a convergent solution. However, the magnitude of the physical dispersion term has a strong influence on the results [3]. We may anticipate that formulations which neglect physical dispersion, and instead introduce an artificial grid size dependent dispersion through the use of upstream weighting, may not converge as the grid is refined. This effect is demonstrated in two of the example problems presented in the body of this paper. In addition to its influence on the problem of grid orientation sensitivity, the following effects of dispersion have been identified: (1) Dispersion controls the length of mixing or transition zones between various fluid banks

PI*

(2) For adverse mobility ratio displacements, dispersion enhances sweep efficiency [3]. (3) For multiple-contact-miscible displacements, dispersion reduces displacement efficiency

Due to these observations, it is an objective of this study to develop and test numerical techniques for simulating multiphase multicomponent flows, including the effects of physical dispersion. Several example problems of increasing complexity are considered. To develop the necessary background, the simulation of miscible and immiscible flows are first considered separately. The multiphase multicomponent flow problems considered are for the displacement of oil by solvent. Numerical methods are tested to determine if they converge as the spatial grid is refined and whether convergent solutions have a sensitivity to grid orientation. Since spatial approximations are the primary subject of this paper, all examples presented used time steps small enough to make the temporal truncation errors small relative to the spatial truncation errors. A variety of numerical approaches for this class of problems is considered, Both upstream weighting and artificial dispersion techniques are studied. Finite difference and finite element methods are used. In two dimensions both five-point and nine-point difference schemes are applied. The finite element methods studied are based on the Galerkin method with Lagrange trial functions and Lobatto quadrature [3], These methods give a lumped capacity matrix and computational molecules with star structures. In fact, when linear polynomial trial functions are used, the method is a five-point central difference approximation. The lumped capacity matrix is essential for some solution techniques [6]. Previous work with these methods was promising [3]. Although this family of Lobatto-Galerkin methods is specific, we expect its numerical behavior to be similar to that of other higher-order finite difference and finite element methods. This paper is divided into several sections. In Section 2, the governing equations are described for three-phase flow with mass transfer between two of the phases. In Section 3, the miscible displacement problem is discussed. Several different discretizations are considered and their relative efficiency is compared. Next, in Section 4, the immiscible displacement

problem is considered. Both conventional upstream weighting methods and an artificial dispersion approach are applied to the immiscible problem. After analyzing methods for purely miscible flow and purely immiscible flows, problems are considered which contain both types of ftow. In Section 5, numericat methods for this type of problem are developed by generalizing the methods used for the miscible and immiscible problems. The contact miscibie displacement of oil by solvent is considered in Section 6, under conditions when mobile water exists. Then, in Section 7, the displacement of a volatile liquid component by a condensible gas component is considered. This problem not only contains aspects of miscible and immiscible flow, but also demonstrates the additional complexity introduced when phase transitions occur. Finally, the multiple-contact-miscible displacement of a gas condensate by nitrogen is considered in Section 8. This problem also involves phase transitions and illustrates the importance of physical dispersion in this type of displacement process.

2. Governing equations Three fluid phases (oil, gas and water) are considered, Only water is assumed to exist in the aqueous phase, while the rt nonaqueous components may transfer between the oil and gas phases under the assumption of local equilibria. The governing equations are of three types: (1) differential material balances, (2) constraints, and (3) phase equilibria relations. The material balance on the aqueous component is

while the material bafances on the nonaqueous

components

The velocities are given by Darcy’s Law extended

uo= -k

k (VP,- p&w),

are

to multiphase

flow,

(2)

The porosity 4 and absolute permeability k are primarily functions of spatial position; however, in certain cases their dependence on pressure also may be important. The constraint equations require the phase saturations (volume fractions) and the mass fractions in each phase to sum to unity,

L.C. Young, Spatial approximationsfor simuluting fruid dispiaceme~ts

s,+s,+s,=1,

5 i=l

Xi

=

1,

f:

yi =

1%

7

(3)

i=l

Phase equilibria is established by the requirement

of equal component

fugacities,

foi=fgi, i = 1, . . . , n.

(4)

The component fugacities, phase densities, and phase viscosities are provided by equations of state. Each of these properties are functions of pressure and the respective phase compositions. The phase pressures are related by the capillary pressure functions, (5) The capillary pressures and relative permeabilities, k,, krg, k,, are expressed as functions of the phase saturations. These functions are history dependent and may also depend on spatial position. In addition, for some recovery processes the dependence of these variables on capillary number must be considered. The dispersion tensors are given by [7]

where i = X, y, z and i = X, y, z. In most cases of practical interest, the molecular diffusivities, D, and D,, are negligible relative to the hydrodynamic dispersion terms. Tbe dispersivities, CY,~, (Y,~,agt, Q, are dependent on phase saturations [S, 91. For single phase contact miscible displacement of one fluid by another, the degree of mixing is a function of the viscosity ratio or mobility ratio of the two fluids [4, lo]. It may be possible to account for this phenomenon by allowing the dispersivities to be dependent on viscosity gradient. This approach has apparently not been considered. In some cases the functional dependencies, just discussed, have not been described quantitatively. This is not a drawback here, since the objective of this work is to investigate numerical methods. For this purpose, only the essential features of these relationships are maintained. For example, the dispersion relationship, (6), is simplified by taking (Y,(, o!013CQ-,apt to be equal in value and constant. For simplicity only horizontal flows in one and two dimensions are considered, i.e., VH = 0 in (2). The capillary pressure differences, (5), are neglected in some of the example problems. In order to gain insight into the physical and numerical aspects of displacement problems, it is useful to examine the special case for which ideal mixing of fluids, i.e. additive volumes, is assumed,

and for which all components

are assumed to have the same compressibility,

3

Pi = PPY(P)

Pw = P%P)

where p’j and p$ are constants. fractions, c OI . = ~~X~~i,

9

(8)

This assumption

permits

the definition

of phase volume

cgi = ~~yi/~i f

(9

and overall volume fraction, Si = SoCof+

SgCgi

(10)

e

The total fluid mobility is defined by A = h,+A,+A,,

(11)

and the phase and component fW= AWlA,

fo

=

fractional flow functions are h,ih

,

f,

=

h,th

,

W fi

=

focoi

+

f&i

-

The fractional flows and total mobility are ultimately functions of the fluid overall composition, S; however, the equations of state are not explicit in overall composition. Under the assumptions of ideal mixing, equal component compressibilities, and negligible capillary pressure, the mass balance equations are expressed in terms of these newly defined variables. The material balances are, for water,

and for the nonaqueous

components,

#~(ySi)=-V~yf,u+v‘a~~*~~Vx,+V.*~/u,/Oyj I

(13b)

i

for i = 1,. . . , n, where CY= aOe = sot = age = agt. In (13) the overall velocity is given by u = u,+ u,+ u,= -kAVp, The summation of the water and all of the nonaqueous following overall balance: dY 4 -=-Vyu=VykAVp. at

(14) component

balances, (13), gives the

(1%

9

L.C. Young, Spatial approximationsfor simulating fluid displacements

Although the assumption of ideal mixing is not always valid, it does not appear to destroy the essential features of the problem, either physically or numerically. The ideal mixing assumption allows the basic structure of the equations to be more clearly visualized. Equations (13) are a set of nonlinear convective-diffusion equations, where the velocity field is determined by the solution of (15). If physical dispersion is neglected, (13) are hyperbolic and shock fronts appear in many solutions. For most cases of practical interest, the dispersion terms are small and the solutions exhibit sharp fronts. The behavior of solutions to these equations are more clearly illustrated by the examples that follow.

3. Single phase miscible displacements This problem is one for which only a single phase is flowing. The displacement consists of one component (solvent) miscibly displaceing another (oil). Under the assumption of constant density, the material balance, (13b), for the solvent is

(16)

+p=-v.u=o,

(17)

where c is the solvent fraction and the subscripts on c and p have been omitted for simplicity. The boundary and initial conditions to be used with these equations are as follows: (1) initially c = 0; (2) at injection boundaries the injected component is flowing at a rate equal to the total flow rate; (3) at production boundaries, the dispersive flux is zero; (4) all other boundaries, no flow occurs. This problem is considered first in one dimension and then in two dimensions. 3.1. Miscible displacement in one dimension In one dimension balance is

25, a7

the velocity, u, is constant. The dimensionless

form of the component

ac- 1 a2c a(

Pe a["'

where 7 = tu/(4L), 5 = X/L and Pe = L/a. T is equivalent to the hydrocarbon fluid injected (HCPVI). The boundary and initial conditions are c(& 0) = y

= 0)

(18) pore volumes of

(1%

-1 ac(0, T) + c(0, T) = 1, Pe a#$

(19b)

For large values of Fe, the solution to (18) and (19) can be approximated

by [4]

C = $ erfc[(c - T) ZiPej471 . This solution possesses the property that the dispersed zone grows at a rate proportional to the square root of time. At any given time the length of the dispersed zone is inversely proportional to the square root of Pe. These features of the analytical solution have an effect on how numerical methods behave. The Lobatto-Gale&in method with ‘Linearpolynomial trial functions is equivalent to central differences. The approximation of (18) and (19) using this technique with uniform grid spacing is 4AX 2 AX

2-

The approximation differences is

- :@I- cz) -t &

$(c~-~-

(CI - cz) + (CI - 1)-O,

Cir~)-~(cj_~-2cjfcii~)=0 fori=2,...,N,

of (18) and (19) using conventional

upstream

weighting

or one-sided

~AX~+-&$-C,)+(Q-l)=O, AX g-

(cjml- ci) -&(Ci-~-2Ci+Ci+3=0

fori=2,,..,N,

(22)

By comparison of (21) and (22) it is apparent that if an artificially small Peclet number, Pe*, (artificially large dispersion coefficient) is used in (21) such that 1 -=&AX=&-+&, Pe*

then (21) becomes equivalent to (22) The numerical error with the upstream weighting method is due almost exclusively to the difference between the actual Peclet number and the effective value, Pe*, given by (23). For small differences the error can be approximated by a Taylor series expansion using (20). This type of analysis suggests that the maximum error at

L.C. Young, Spatial approximations for simulating fluid displacements

any value of 5 can be approximated E r= -&

11

by

Pe - Pe* e-l’* Fe J F:: = 0.242 Pe + Pe* ,

(24)

Fig. 1 shows the ratio of the actual numerical error to that predicted by (24). The actual numerical error was computed by using linear interpolation of the approximate solution within each grid interval and comparing the approximate and analytical solutions at points throughout each interval. Since an exact analytical solution does not exist for the boundary condition (19b), this condition was modified to ~(0, T) = 1. Equation (20) is also a good approximation of the solution with this modified boundary condition when Pe is large. It is apparent from the results in Fig. 1 that (24) gives a good approximation of the actual error as N becomes large and thus Pe* approaches Pe. Although the calculated errors are not constant with time as predicted by (24) they do approach an asymptotic value for large time. The agreeement with (24) is best for large Pe, since (20) becomes a better approximation to the true solution. By substituting (23) into (24) we have 0.242 ‘=1+4N/Pe

(25)

~0.06%.

Equation (25) indicates that the upstream weighting method is first-order correct, and for a given error tolerance the number of grid intervals must be proportional to Pe. In recent years a great deal of attention has been devoted to the so called ‘optimal upwind’ method (see [ll] and references therein). This method is also equivalent to replacing the Peclet number in (21) by an artificially small value, Pe* = 2N tanh (PeI2N).

_ kh

(26)

1.4-

O

-

x

l&p

A

o 0

0

I!!

0

0

;

0

X

0

u,

ON= w,Fe=250 XN- 200,Fe=250 AN= 500,F~=250 0 N=1000,F'e=1000

A

ya

::

;

2

.

1.0_____:_‘“_-~__a--:__:___

0.60

'

I 0.2

I

I 0.4

I

I 0.6

I 0.6

T

Fig. 1. Ratio of maximum

error to the value given by (24) for upstream weighting.

12

L.

c.

Young, Spatial upproxi~ati~ns for ~~~~lati~g j&id di~p~a&e~ent~

n

0

0.2

0.4

0.6

0.8

Fig. 2. Ratio of maximum error to the value given by (24) for the optimal upwind method.

Equation (26) is derived by forcing the approximate solution to fit the steady state convectivediffusion equation. Since this method also alters the Peclet number, (24) can be used to estimate the maximum numerical error. Fig. 2 shows the ratio of the maximum error to that predicted by (24). From Fig. 2, it is apparent that (24) adequately represents the error of the optimal upwind method. The deviations from (24) are similar to those observed for standard upstream weighting. The asymptotic error for this method is found by substituting (26) into (24) and expanding for large N, E = 0.01(Pe/N)2 .

(27)

Equation (27) indicates that the optimal upwind method is ultimately second-order correct and for an error of two percent this method requires one-fourth the number of grid points needed for conventional upstream weighting. However, this method is similar to standard upstream weighting, since the error is essentially constant with time and for a given accuracy the number of grid points must be proportional to the Peclet number. With conventional central differences, (21), and higher-order methods the behavior of the approximate solutions is different from that observed for upstream weighting techniques. To analyze these methods, first consider two systems of different lengths, but with the same value of LYand AX. Prior to the front reaching the outlet boundary, the maximum error for these two systems should be the same at a given value of dimensional time. This reasoning suggests that in terms of dimensionless variables the error should be the same for given values of Pe T and N/Pe. Since the convergence rate of these methods is known, the asymptotic error should be E = (Pe/N)“F(Pe

7) ,

where pn is the order of the method and F is some function.

(28)

13

L.C. Young, Spatial approximationsfor simulating fluid displacements 10

0-2

I

a LINEARS (CENTRALDIFFERENCES): x N=120 Pe= 250, AN=lSDF'e= 500, ~N=2DDh=lOOO,

x 8 R

x

.

%x

. QUADRATICS: AN=%I,Pe= mN=40,Ps=

OA

.

500 250

0

0

.N=80,F'e=lDOO

0-3

00

.A . -.

u

;d -

. 3

kl

.A

l

A % . l

0-4 l

0

I 100

I

WT

Fig. 3. Effect of grid size, Peclet number, functions) and quadratic trial functions.

and time on the maximum error with central differences (linear trial

Fig. 3 contains a plot of e(N/Pe)2 versus Pe T for central differences or linear trial functions. These results can be approximated by E =

0.023(Pe/N)2 (Pe 7)-O.’ .

(29)

Unlike the upstream weighting methods, with this method the error decreases with time. As a result, for an error specified at a given time, the number of grid intervals do not have to be directly proportional to the Peclet number. The approximation to (18) and (19) using the Lobatto-Galerkin method with quadratic polynomial trial functions and uniform grid spacing is fAXdC’+‘(-3c

dr

$AX

1

6

%-$(Ci_,-

$AXh+‘(c-_ dr ’

+4c

2

-c)+l3

Ci+l)-g

I2

,~(7C,-8C2+C~)+(C1-1)=0,

for i=2,4,...,2N,

&(ci-~-2ci+ci+~)=O

-4c._1 1 +4c. r+1-

ILki -

8Ci-l+

Ci+2)+3&(Ci-2-

(30)

8Ci+l+

for i = 3,5, . . . ,2N - 1, ZiAX

dC;+l +

&+I

- ‘&v

+ ~CZN+I) + ; &

(CAN-1

-

8C2r4

+

7C2~+1)

=

0

.

Ci+2)

=

0

14

L.C. Young, S~tiu~ upproxi~a~o~

Table 1 Number of grid intervals equation (18)

for maximum

for simulating j&.&d~~sp~ace~ents

error of 0.02 at T = 0.5 for

Pe Method

100

250

500

1000

Upstream weighting Optimal upstream Central differences Quadratic finite elementa

300 71 40 9

750 177 80 16

1500 354 135 24

3000 707 227 35

“The number of grid points are 2N + 1 for quadratics and N + 1 for the other methods.

For this method bY E =

e(N/Pe)3 is plotted versus Pe T in Fig. 3. These results can be approximated

0.0021~e/~)3

(Pe 7)-1,25.

(31)

For this method the error also decreases with time and at a faster rate than with the linear trial functions (central differences). Like the linear case, for an error specified at a given time, the number of grid intervals do not have to increase in direct proportion to the Peclet number. For the quadratic trial functions when the Peclet number is doubled, the number of grid intervals must only increase by fifty percent. For still higher-order methods the error decreases with time at an even faster rate. By differentiating (20), the maximum value of &/a[ is found to be proportional to T-‘,~. The higher derivatives of c decrease with time at a faster rate. The behavior of the higher-order methods is consistent with the notion that with additional orders of accuracy the error is proportional to higher derivatives of c. The results just described have important consequences for flows that are strongly dominated by convection, i.e., with very large values of Pe or small values of ty. For the various methods, Table 1 lists the number of grid intervals required for a maximum error of 0.02 at 7 = 0.5. These results suggest that, for this problem, if one is truly interested in resolving the dispersive effects in the solution, then conventional methods are superior to upstream weighting methods, especially for large values of Pe. Also, for large values of Pe, higher-order methods are favored. 3.2. vincible displacement in two dimensions The miscible displacement problem is considered for a two-dimensional five-spot pattern. Fig. 4 illustrates the pattern geometry and shows that two different square elements of symmetry can be used. For the diagonal grid, the line connecting injector and producer goes across the grid lines. For the parallel grid, the injector and producer are connected by a grid line. This problem has been studied using upstream weighting techniques [2] when the

L.C. Young, Spatial approximationsfor simulating fluid displacements

A

A

15

A

A

,DIAGONAL GRID PARALLEL A 0

A

l-4A

Fig. 4. Geometry with symmetry elements for a repeated five-spot pattern, 0 Production well; A Injection well.

Fig. 5. Contours of solvent fraction for the upstream weighted five-point method. M = 10; 20 X 20 Diagonal grid; ---28 x 28 Parallel grid.

dispersion term is neglected and with central differences and finite element methods [3,12] when the dispersion term is included in the numerical approximation. In one dimension, upstream weighting is equivalent to central differences with an artificially small Peclet number (large dispersion coefficient). A similar relationship does not exist in two dimensions, Equation (16) indicates that the dispersion coefficient is proportional to the magnitude of the velocity vector. A truncation error analysis of upstream weighted five-point differences gives a second derivative in the X-direction with a coefficient of $AXlu,l and a second derivative in the Y-direction with a coefficient $AYlu,\. These leading terms in the truncation error are not invariant to a rotation of the coordinate system; whereas the physical dispersion term in (16) is invariant to a coordinate rotation. Upstream weighted five-point differences give results that have a dependence on grid orientation when the viscosity of the solvent is less than the oil viscosity (an adverse mobility ratio). Fig. 5 shows concentration contours for this method for a viscosity ratio (oil : solvent) of ten. In this simulation mixture viscosities were calculated by the quarter power mixing rule [3]. From Fig. 5 it is clear that this method gives unsuitable results. This sensitivity to grid orientation is not diminished as the grid is refined. Fig. 6 shows concentration contours for the central difference method which corresponds to the use of linear trial functions with the Lobatto-Galerkin method [3]. For these results, the method was applied to the nondivergence form of (16) [3]. In this case, physical dispersion is included. The Peclet number, L/a, is defined using the characteristic length indicated in Fig. 4. These results have very little sensitivity to grid orientation. The results using Lobatto-Galerkin methods were tested for their sensitivity to grid size. Fig. 7 is a plot of recovery at 1.2 hydrocarbon pore volumes injected versus grid size. These results demonstrate that as the grids are refined, both the diagonal and parallel cases converge to the same answer. Fig. 8 contains a similar plot for the quadratic Lobatto-Galerkin method. With this method the solutions for both cases also converge to the same answer as in Fig. 7.

Fig. 6. Contours of solvent fraction for the menteal difference mathod or linear trial functions~ M = 10, Pe = 40; ----20x 20 Dia&mal grid; ---28 x 28 ParalleI grid.

0.82

0.76 2

With regard to the efficiency of the linear method versus the quadratic method, Figs. 7 and 8 indicate that for a given number of grid points the accuracy of the two methods is similar. lf other criteria are used, the quadratic method would be favored. For AX/L. = 0.10 the finear case gave osc~~~at~o~s of conce~t~at~o~ in excess of 0.20; whereas with quadratics and AX& = 0.2QS the rn~~um osc~~~at~o~~were on the order of 0234 The sensitivity of the results to dispersion level was atso investigattid, Fig. 9 shows for various Pe the cumulative oil recovered as a function of hydrocarbon pore volumes of fluid injected (HCPVI), Prom Fig. 9 it is clear that recovery is strongly influenced by Pe. The recovery is as sensitive to Pe as it is to the viscosity ratio [?I.

17

L.C. Young, Spatial approximationsfor simulating fluid displacements

i > a ?

0.82 -

0.76

0

I

, 0.04

I

I 0.06

I

I 0.12

I

I 0.16

I

I 0.20

I

_ 0.24

AX/L Fig. 8. Convergence grids.

of quadratic

0

trial functions for M = 10 and Pe = 40. -0-

0.2

0.4

0.6

HYDROCARBON

0.6

1.0

PORE VOLUMES

Diagonal grids; --a--

1.2

1.4

Parallel

1.6

INJECTED

Fig. 9. Effect of Pe on cumulative recovery for M = 10.

The results in Fig. 9 were computed for various grid sizes. The results in Figs. 7 and 8 indicate that the recovery curve for Pe = 40 in Fig. 9 is accurate to within 0.001. Similar calculations were performed for the other values of Pe on diagonal grids only with the finest grids being 25 X 25 with quadratics and 12 X 12 with cubits. These calculations suggest that all of the recoveries in Fig. 9 are accurate to within 0.01. Figs. 5 and 6 indicate that the form of the truncation error for the upstream weighted five-point method causes the grid orientation sensitivity. Although the truncation error is diffusive in nature, it is not invariant to grid rotation and it cannot be related to the physical dispersion term. Fig. 9 indicates that the solution is sensitive to the magnitude of the physical

18

L.C. Young, Spatial ~~proxi~atio~

for si~~~a~~g fruid ~~pla~e~e~ts

dispersion term. The results in Figs. 5 and 6 indicate that the solution is also very sensitive to the functional form of the dispersion. When the upstream weighted five-point method exhibits grid orientation sensitivity, the most widely used alternative is the upstream weighted nine-point method [2]. Fig. 10 shows concentration contours for this method. From Fig, IO it is apparent that very little orientation sensitivity exists. Apparently, with this method the dispersion introduced by upstream weighting is essentially invariant to grid rotation. Fig. 9 indicates that the results are sensitive to the magnitude of the dispersion term. Since upstream weighting introduces an artificial dispersion which is proportional to grid size, the sensitivity of the solution to grid size was investigated. Fig. 11 shows a plot of recovery versus grid size for the nine-point method. In Fig. 11 the results from Fig. 9 are also plotted versus 2/Pe, since in one dimension (23) indicates this to be the relationship between grid size and the effective dispersion level of upstream weighting. The effect of grid size on the nine-point results is qualitatively similar to the effect of varying the Peclet number. As the grid size (or 2/Pe) is decreased an unchanging answer is not ap proached. This result is not unreasonable, because in the absence of dispersion the physical problem is not well posed. From Fig. 11 one might attempt to conclude that the effective dispersion level with the nine-point method is less than (23) suggests. However, if Figs. 6 and 10 are compared it is apparent that lengths of the dispersed zones are similar with the two methods, while the recovery at 1.2 hydrocarbon pore volumes injected are 0.80 and 0.76, respectively. Thus, the numerical dispersion introduced by upstream weighting the nine-point difference method bears only a qualitative similarity to physical dispersion. These results show that for adverse mobility ratio miscible displacements both the functional form and magnitude of dispersion affect the solution. Upstream weighting techniques in

0.65

0

8 0.05

I 0.10

I 0.15

I 0.20

0.25

AWLor 2/l’@ Fig. 10. Contours of solvent fraction for the upstream weighted nine-point 20 x 20 method. M = 0; Diagonal grid; ---28 x 28 Parallel grid.

Fig. 11. The effect of grid size (AX/L) on recovery predicted by the upstream weighted nine-point method. M = 10, -0Diagonal grids; --a-Parallel grids. Effect of physical dispersion (2/Pe) on predicted recovery, -- X --.

effect introduce a grid size dependent artificial dispersion which cannot be quantitatively related to physical dispersion. For the upstream weighted five-point method, the artificial dispersion is not invariant to rotation of the coordinate system and the results are sensitive to grid orientation. With the upstream weighted nine-point method, no grid orientation sensitivity occurs; however, since the magnitude of dispersion decreases as the grid is refined, the results fail to converge to an unchanging answer. The most desirable solution to the numerical difficulties associated with upstream weighting is to instead include physical dispersion in the approximating equations. As in one dimension, when physical dispersion is included in the approximation, central difference methods or conventional finite element methods are superior to upstream weighting. For this particular problem the higher-order Lobatto-Galerkin methods are more efficient than central differences (linear trial functions).

4. Immiscible displacements This problem is one for which two completely immiscible components are flowing. The displacement consists of the component occupying one phase (water) displacing the other (oil). Assuming constant densities for this problem, the material balance equations reduce to

4

aS,_. - -v a7

9 ufw -

V - kfwk,Vpcw

,

v*u=o,

(33

where u = u, + u, = - khVp, + kh,Vp,,,

.

The initial and boundary conditions used with (32) are similar to those used with the miscible displacement problem. (1) initially S, = SwC; (2) at injection boundaries, the flow rate of water is equal to the total flow rate; (3) at production boundaries, the capillary pressure gradient is zero; (4) at all other boundaries, no flow occurs. This problem is considered first in one dimension and then in two dimensions. 4.1.

~~~isci~~e

dis~lace~en~

In one dimension, equation is

as af(q

a7c ag

in one dimension

the velocity, u, is constant.

=$[m$],

The dimensionless

form of the saturation

(33)

20

L.C. Young, Spatial approximation

for simulating @id di~piaceme~~s

where S = (S, - &x)/Q - Six - &Ix), 7 = tu/[&q1-

so,-

SW)] ,

5=mf, and

g=-(1-s..-&.)%d@ In (33) the subscript on f has been dropped for convenience. T is equivalent to the displaceable hydrocarbon pore volumes of injection (DI-ICPVI). The initial and boundary conditions are

f-g%=1CXS

att=O,

(34

J!z$.&O.

Equation (33) is the well-known Buckley-Leverett equation. This problem is equivalent to the miscible displacement problem when f = S and g = l/Pe. For an immiscible displacement the specific shape of the fractional flow curve can affect the performance of numerical methods. Two different cases are considered for this reason. The properties for Case 1 are listed in Table 2. Case 2 uses the properties from the immiscible displacement problem considered by Yanosik and McCracken [2]. The fractional flow curves for these two problems are plotted in Figs. 12 and 13, respectively. Equation (33) possesses analytical solutions for certain limiting cases (for a review of this problem, see [13]). Buckley and Leverett [14] and Welge [15] considered this problem for the case of negligible capillary pressure, i.e. g = 0. Applying the method of characteristics to this simplified equation gives

+$. Table 2 Physical properties for immiscible displacement, k,, = 0.5 S’ k,=(1-S)3, ~oIl.hv = 3 When capillary pressure is considered:

k(l-Sor-Swc)=003 pwuL .

(35)

Case 1

21

L.C. Young, Spatial approximationsfor simulating fluid displacements

0.8 -

0.8 -

A @P-l

f

)

f

i

0

1 6.2

; Sb = 0.287 , ’ 11 ’ 0.4

1 0.8

1

’ 0.8

’

0.001 1.0

S

Fig. 12. Fractional

Equation

flow curve for Case 1 (Table 2).

(35) applies for saturations

Fig. 13. Fractional Case 2 from [2].

flow and total mobility curves for

greater than S,,, which is determined

by

s,qp-f&)=0. For saturations

(36)

less than this value the Rankine-Hugoniot

jump condition [16] applies,

sb=T-yL b

(37)

Buckley and Leverett determined the relationship for & (36), by requiring the solution to conserve mass. Welge interpreted (36) graphically as a tangent to the fractional flow curve as illustrated in Figs. 12 and 13. For Case 1 the zero capillary pressure solution is plotted on Fig. 14. The concept of a tangent to the fractional flow curve was first introduced by Terwilliger et al. [17]. In addition, Terwilliger et al. considered the problem with nonzero capillary pressure. Their results can best be explained by first transforming to a coordinate system which moves with the velocity of the front given by (37) 17 = t--fl(Sb)T.

With this transformation,

(33) becomes

$+-&,,-f’(sb)s]=$(@

For small values of g, (38) is characterized

by two limiting conditions:

0.8

Fig. 14. Comparison of analytical solutions to exact solution for Case 1 at T = 0.4 with capillary pressure. Buckley-Leverett, g = 0; 000 Stabilized zone solution, g = 0.03(1- fi; ---Exact (fine grid) solution, g = 0.03(1 - f).

(1) For S>& the first and second terms are dominant and the zero capillary pressure solution is applicable. Terwilliger et al. called this the unstabilized zone, since the saturation profile continuously spreads with time by (35). (2) For S < Sb the time derivative is small relative to the last two terms, such that

[f-f”(Sl$q = g g

’

Terwilliger et al. called this the stabilized zone since the shape of the saturation profile is unchanging with time. Near S = S, all three terms in (38) are important and a transition zone results. Terwilliger et al. [17] verified these concepts experimentally and observed that the transition zone about Sb diminished with time. The data in Table 2 results in the following function: g(S) = 0.03(1- j-)

I

(4~~

The stabilized zone solution, i.e., the solution to (391, is plotted on Fig. 14. Also plotted an Fig. 14 is a numerical solution to (33) and (34) that was calculated using forty quadratic elements (N = 40). For practical purposes, this numerical solution may be regarded as exact. It should be noted that (39) provides only the shape of the stabilized zone and not its lateral location. In Fig. I4 the lateral location was determined by matching the solution to the numerical solution. From Fig. 14 it is apparent that the saturation profile is accurately described by the Buckley-Leverett solution for S > Sb and by the stabilized zone solution for S < Sb, A transition region exists near S = &,, and the size of this region decreases with time. These obse~ations are in agreement with those of Terwilliger et al. 1171, Most field scale problems of practical interest give parameter values such that g(S) is much smaller than that given by (40). Therefore, the capillary pressure effects are so small that the

of these effects are neither practical nor ae~~%sary.An artificial capillary pressure (artificial dispersion) term much larger than that given by (333 must be us@dto obtain suitable results with grids uf practical size, For these reasons, the physical value of this term is not considered in the remainder of this ~~~~~~~~~~ a farger ~~~~~c~a~ term is used instead, The idea af a stabifized zone and (3) is stiIf valid eye5 thorns the dispersion ternI is ~tj~~~a~ ra#her than F~~s~ca~= The stabilized zone concept is use&X for de~e~o~~ngan ~~~~~stand~ngof the performance of vmrious numerical methods, Upstream weighting is the norm&l means of incorporating artificial dispersion into the numerical approximation. Applying upstream weighting to (33) with g = 0 results in the foliowing ditierence aFpro~imat~un~ resolution

By analogy with the miscible displacement problem, (45) is equivalent to solving fhe following problem using central differences:

where go = SAX Rg. 15 illustrates s~~~t~~~sobtained using upstream w~~~~~~~~ for Case 1. From Fig. 15 it is apparent that with upstream weighting the numerical solutian converges to the Buckley-Lever&t solution.

If central differences are applied to (42) and go is fixed in value, then the saturation prafile

24

Fig. 16. Convergence of linear finite elements to stabilized zone solution of 1421, at T = O,d with N=20, go=;AX;xxxN=50, I. -Buckley-Leverett; ---Stabilized ZQIW solution; 000 OOON= 100, go=2<5AX.

for Case g,=L25AX;

go = 0.025

should converge ttrr the stabilized zone solution. Fig. I6 shows that indeed this is the case. It should be noted that for most relative permeability data

so that at S = 0, (43) reduces to

For the relative permeability data in TabIe 2, r(O) = 0, so for this problem the saturation gradient at S = 0 has an abrupt change from negative infinity to zero. For most problems of practical interest, the solution to (42) is not smooth at S = 0. Despite this lack of smoothness, upstream weighting performs reasonably well because the approximation possesses a maximum principal, viz., the vaiue of saturation at a particular grid point cannot increase (decrease) unless the value at a neighboring grid point is greater (less) than the value at the grid point of interest. This property fohows from the farm of (41) and the fact that f is a monotonically increasing function of S. The maximum principal insures that the upstream weighting method will not oscillate despite the lack of smoothness in the solution, Higher-order methods such as quadratic finite elements do not possess a maximum principal when artificial dispersion like that in (42) is used. The lack of smoothness about S = 0 causes these methods to oscillate about this point, even for large values of g,. For higher-urder methods, artificial dispersion of a more appropriate form can be determined by examination of (38). For example, the following ~nctional form cuuld be used:

and

25

L.C. Young, Spatial approximationsfor simulating fluid displacements

This artificial dispersion term was constructed stabilized zone:

to give the following saturation

gradient in the

where O(gI) denotes terms that are on the order of gl. For gl = 0, the maximum value of dS/dq = l/go, and the stabilized zone solution is not dependent on the fractional flow curve. Values of go = 2AX and gI = 0.04 are approximately optimal for the quadratic finite element method. These values were determined by numerical experimentation with several different fractional flow curves. Results calculated using this approach are compared to the analytical solution and to the upstream weighting results on Figs. 17 and 18 for Case 1 and Case 2. The

0.6s

0.4-

0.2

00

0.4

0.6

0.6

10

X/L

Fig. 17. Saturation profile from (44). -Buckley-Leverett;

at r = 0.4 predicted using quadratic finite elements O--O Upstream weighting, N = 20; x . .

l.O-

I

for Case 1 with artificial dispersion Quadratic finite element, N = 10.

. .x

1

0.8 -

0.2+.yq 0 0.0

0.2

\ .I.. I b__, 0.6

0.4

0.6

1.0

X/L

Fig. 18. Saturation profile from (44).-Buckley-Leverett;

at T = 0.2 predicted using quadratic finite elements for Case 2 with artificial dispersion O--O Upstream weighting, N = 20; x . . . . . . x Quadratic finite element, N = 10.

L.C. Young, Spatial upproxjrn~~o~ for simulating f&d ~i~p~~c~rn~n~

26

average absolute error is defined as

where SBL and SNUMare the Buckley-Leverett and numerical solutions, respectively. The average absolute error has been calculated and plotted against grid size for Case 1 and Case 2 in Figs. 19 and 20. For the upstream weighting method the numerical solution was represented by linear interpolation in order to calculate the error. The quadratic trial functions were used in the error calculations for the finite element method. The results in Figs. 17-20 show that when an appropriate artificial dispersion term, such as given by (44), is used with quadratic finite elements, the method is superior to upstream weighting. The only difficulty with this approach is that it relies heavily on knowing the characteristics of the solution when a given artificial dispersion term is used. For more general problems, such as those in the following sections, this information may not be available. The quadratic method can be extended to these cases only if one uses parameters like go which can be adjusted for any given problem based on numerical experiments. This approach is not considered worthwhile. Higher-order methods need a more universal artificial dispersion term. The artificial dispersion term in (42) with go = $AX, can be used with linear finite elements or finite differences for a wide range of problems. This form may not be optimal for a specific problem, but it generally gives satisfactory results without adjustable parameters that must be determined on a case-by-case basis.

0.005

0.01

0.02

0.05

I

0.1

l/N OR 1/(2N)

Fig. 19. Effect of grid size on average absolute for Case 1 at T = 0.4. M Upstream weighting versus l/N; X-X Quadratic finite elements versus 1/2N.

0.001

*

L

0.005

*

.

I 0.01

0.02

0.05

I

0.1

l/N OR 1/(2N)

Fig. 20. Effect of grid size on average absolute error for Case 2 at 7 = 0.2. G----O Upstream weighting versus l/N; X-X Quadratic finite elements versus 1/2N.

L.C. Young, Spatial approximationsfor simulating fluid displacements

27

The higher-order finite element methods are superior to linear elements or central differences for the miscible displacement problem. Figs. 17 through 20 indicate their potential for the immiscible displacement problem. Despite these results, higher-order methods are not pursued further in this paper due to the difficulty of determining a generalized artificial dispersion term. Further development of higher-order methods should be the subject of continued research. 4.2. Immiscible displacement in two dimensions In two dimensions, the problem studied by Yanosik and McCracken [2] is considered, since this is a difficult adverse mobility ratio displacement. The fractional flow and total mobility curves for this problem are plotted in Fig. 13. Traditionally, the mobility ratio for an immiscible displacement is evaluated using the average saturation behind the front to evaluate the driving fluid mobility. This saturation is given by the point where the tangent line intersects the top of the figure, i.e. S = 0.372. By this definition the mobility ratio for this problem is approximately four. If the mobility ratio is based on the saturation at the front, i.e. Sb = 0.287, the value is approximately two. It is difficult to develop data for a realistic immiscible displacement problem with a higher mobility ratio. The problem just described is considered using upstream weighting with five-point and nine-point differences and by the use of artificial dispersion with linear finite elements or central differences. Upstream weighting is normally applied to the individual equations for oil and water. Usually, the mobility of each phase is evaluated at the upstream location. When this approach is used with a five-point difference scheme, it is equivalent to solving. the following equations by central differences:

o=VkhVp+!AX&

(kl -$-$&)+fAY-$(kl

-j$$+

Wb)

This approach is equivalent to adding an artificial dispersion term to both the saturation equation (compare (32a) and (46a)) and the overall balance ((32b) and (46b)). If only the phase fractional flows were evaluated at the upstream location, the overall balance would not be altered. Equation (46) also suffers from the fact that the artificial dispersion terms are not invariant to a rotation of the coordinate system. For the miscible displacement problem, this was found to cause the grid orientation sensitivity. With the linear finite element method the following approximations to the water saturation equation are considered:

4 $$=V.kAf,Jp+{ V. goklVplV& V *Sol~lVfw * The oil saturation

equation

is treated

in an analogous

(474 W’W fashion.

In (47) either

artificial

dispersion term is invariant to grid rotation and similar in form to the physical dispersion term for the miscible displacement problem (see (16)). Equation (47a) will alter the overall balance in a manner analogous to (46b). Equation (47b) will not alter the overall balance since the fractional flow functions for oil and water sum to unity. In addition to the two different artificial dispersion terms, three methods of approximating the convection term are considered. In (47), the convection term is written in divergence form. Expanding the convection term by the chain rule, i.e., V - kAfJp

= f,J’ +khVp f khVp +Vfw ,

w-9

is equivalent to writing the equation in nondivergence form. The trapezoidal rule used with the linear finite elements is not sufficiently accurate to yield exact agreement for (48). For this reason the approximations will be different for the divergence and nondivergence form of the convection term, Since both of these are valid approximations, a linear combination of these approximations also is valid. For reasons which will become apparent, the average of the divergence and nondivergence approximations also is considered. We denote the average of these approximations as the median approximation. Appendix A lists the difference approximations for these various cases that were determined using linear finite elements and trapezoidal rule integration. For immiscible displacements in one dimension, the existence of a maximum principal was an important property for the satisfactory performance of methods when artificial dispersion like that in (47) is used. Therefore, the various approximations of (47) have been analyzed to determine sufficient conditions for the existence of a maximum principal. An example of this analysis is presented in Appendix A while the results are summarized in Table 3 for a uniform grid. The divergence approximation with artificial dispersion given by (47a) and the median approximation with artificial dispersion given by (47b) both possess a maximum principle when g, = iAX. The other approximations obey a maximum principle for go = ;AX only if the total mobility is constant. It was for this reason that the median approximation was developed. This approach also has the advantage that by using (47b), the pressure equation is not directly Table 3 Sufficient conditions uniform grid

for a miximum

principal

for equation

(47) and a

Convection term

Dispersion term

Minimum values of go

Divergence

V .

&AX

Nondivergence

V . gokIVp[Vh,

Divergence

V . golulVfw

Nondivergence

V . g&f]Vfw

Median

V *g0lulVfw

“The expressions jacent grid points.

gokIV,lVL

involving

$AX max($, dx

nlaX(hi,

21’ h i+l)

hi+Xi+l

AX maxb,

h iC1)

Ai + Ai+l

tAX Ai and A~+I apply for every pair of ad-

Fig. 21. Saturation contours at 0.3 Displaceable Hydrocarbon Fore Volumes d Injection predicted by upstream weighted nine-point. 20x 20 Diagonal grid; ---- 28 x 28 Parallel grid.

Fig, 22. Saturation contours at 0.3 Displaceable Hydrocarbon Pore VoIumes d Injection predicted by tinear finite elements or central differences with artificialdispersion. 20x 20 Diagonal grid; ---28 x 28 Parallel grid.

affected by the artificial dispersion term. For the results in the remainder of this paper a value of $AX is assumed for go, calculations with the divergence approximation use (47a), while calculations with the nondivergence or median approximations use (47b). Yanosik and McCracken [2] presented contours of water saturation for the upstream weighting methods. Their results showed that the five-point method exhibits a severe grid orientation sensitivity. Contours for the upstream weighted nine-point are presented in Fig. 21.

B 5 0.4 8 : d w 1 3 2 0.3 5 5 E z a.2

0.2

0,4

0.6

0.8

1.0

1.2

1.4

OISPLACEABLE HCPV INJECTED

Fig. 23. Predicted recovery for immiscible displacement problem 20 x 20 diagonal and 28 x Upstream weighted five-point, diagonal; --Upstream weighted five-point, parallel; nine-point; X X X Divergence five-point; 0 0 0 Nondivergence and median five-point.

grids. ---Upstream weighted

28 parallel

30

L.C. Young, Spatial approximations for simulating fluid displacements

Little grid orientation sensitivity is observed with this method. Contours using the median approximation are presented in Fig. 22. This method and the other artificial dispersion approximations (contours not shown) also display very little grid orientation sensitivity. Fig. 23 is a plot of recovery as a function of cumulative water injection. The upstream weighted five-point is the only method for which recovery is sensitive to grid orientation. For the methods which exhibit no grid orientation sensitivity there are only minor differences between the computed recoveries. These methods also have very little sensitivity to grid size. For example, the maximum difference in the fractional recovery predicted using a 20 x 20 diagonal grid and a 5 X 5 diagonal grid is less than 0.01 when either the upstream weighted nine-point or the median five-point is used. Except for the upstream weighted five-point, all of the methods give solutions to this problem which are sufficiently accurate for practical purposes.

5. Numerical approximations for multiphase multicomponent flow The preceding two sections considered miscible and immiscible flow separately. In this section, numerical approximations are considered for the more realistic case in which miscible and immiscible flows are occurring simultaneously in the same system. For example, in the following section a contact miscible tertiary displacement is considered. In that problem, miscible flow occurs between the solvent and oil, while immiscible flow occurs between the aqueous and nonaqueous phases. Numerical methods for these more general flow problems are developed using the insights gained in the preceding sections. The results presented in Section 3 suggested that for miscible displacements, physical dispersion should be incorporated into the approximating equations and for maximum efficiency upstream weighting should not be used. The results in Section 4 showed that for immiscible displacements upstream weighting or artificial dispersion must be used to obtain suitable results. For problems involving both types of flow, artificial dispersion of some type must be used for the immiscible aspects of the flow. Since we wish to resolve the effects of physical dispersion in the solution, it is highly desirable that the artificial dispersion not be incorporated in such a way that it tends to act with the physical dispersion to produce overly dispersive solutions. Equation (lb) is sufficiently complicated that an efficient means of incorporating artificial dispersion is difficult to determine. For the case with constant compressibility and ideal mixing, (13b), a suitable functional form can be deduced more readily. If for the moment the physical dispersion terms are neglected, experience from the preceding sections suggest the following approach for (13b): (49) For incompressible flow, or y = 1, (49) is directly analogous to (47b). An approach analogous to (47a) also could be used. If the last term in (49) is expanded by the chain rule, the following relationship is obtained:

L.C. Young, Spatial ~pproxjmatio~sfor simulating j&id displacements

31

If, in (50), go were replaced by the physical dispersivity, a, the first two terms on the right-hand side would be identical to the physical dispersion terms in (13b). This suggests that an appropriate functional form of artificial dispersion for (13b) is as follows:

$ (ySJ = -v *yjk + v *goylmi

$fJ

Equation (51) possesses the desired property that the artificial dispersion does not augment the physical dispersion. Having addressed the constant compressibility ideal mixing case, a suitable approach for the more general (lb) is now more obvious. Equation (51) generalizes for (lb) as follows:

4

$

(PoScJi

+ @gYi)

=

-v *( +

PokJi

V ’ (a -

+

Pgu*YI) + v *gOl

gO)&l

Ulv(PofoXi

uoIV& + V ’ (a -

+

Pgfgyi)

gO)PglUglVYt *

(52)

The finite element method with linear trial functions and trapezoidal rule integration is used to develop central difference approximations for (52). The approximations are analogous to those for (47) that are developed in Appendix A. The performance of this approach and the upstream weighted nine-point are considered in the following sections for problems involving both miscible and immiscible flows. For the artificial dispersion approach, all examples use go = fAX.

6. Contact miscible tertiary displaeeme~ts This problem is one involving both miscible and immiscible flow. The initial conditions correspond to those of a reservoir which has been thoroughly waterflooded such that water and residual oil remain in the reservoir. A solvent, which is first contact miscible with the oil, is injected to displace the residual oil. Injection of the solvent causes the formation of a tertiary oil bank which displaces the water ahead of it immiscibly. For this problem, the material balance for the water is the same as (la). However, since only one nonaqueous phase exists, the material balance on the nonaqueous components, (lb), reduces to

+

$

(Posoxi)

=

-V

’ (poU,Xi)

+

V ’ pea jUolV&

(53)

Table 4 Physical properties for contact miscible tertiary displacements Relative permeability

Viscosity

S, = 0.20, S, = 0.30, S, = o-35

water: pW = 0.7 cp solvent:

k,=“,5

( SW--SW)

J

l_S,,-SW,

/.&I=

0.1 ep

oil: p2 = 1.0 cp PO= [.%~j.Li~‘~+ XZ/.LF”~]-~

for i = 1,2. If we let component 1 represent the solvent, the initial and boundary conditions for the problem are as follows: (1) initially S, = S,,, x1 = 0; (2) at injection boundaries, the solvent flow rate equals the total flow rate; (3) at production boundaries, the dispersive flux is zero; (4) at all other boundaries, no flow occurs. Constant phase densities are assumed, while the other physical properties for the problem are given in Table 4. It is customary fur this type of displacement to base the injected and produced volumes on the original hydrocarbon pore volume (OHCPV), i.e., the hydrocarbon present at SW,. The oil relative permeability in Table 4 for this problem is the same as for the Case 1 immiscible displacement. The water relative permeability corresponds to a hysteresis curve, since k, = 0 at SW,> S,,. Fractional flow and total mobility curves have been calculated using pure solvent viscosity and pure oil viscosity. These functions are plotted in Fig. 24. The nobility ratio fur solvent displacing oil is unfavorable, white the nobility ratio for the displacement of water by the tertiary oil bank is favorable.

3.0

1.0 A (CP-’

)

0.5 0.3

0.1

Fig. 24. Fractional flow and total mobility curves for contact miscible tertiary viscosity, p= = p2 = 1.0 cp; --- Using soivent viscosity, pO = j.k, = 0.1 cp.

displacements.

-

Using oil

L.C. Young, Spatid approximationsfor simuluting &id dis~la~eme~~s

33

6.1. Contact miscible tertiary displacement in one dimension Fig. 25 is a plot of the oil saturation and solvent mass fraction profiles after 0.3 OHCPV injected for Pe = 100. At this point in time, only solvent and water are present near the inlet of the system. Ahead of the solvent is a transition zone to the tertiary oil bank where mobile oil and water exist. Ahead of the tertiary oil bank is residual oil and water. The tertiary oil bank has not reached the outlet, so only water is being produced. At later times water and solvent free oil are produced, while still later, solvent production begins followed by a continuous increase in the ratio of solvent to oil produced. The results in Fig. 25 were calculated using the median approximation. If N = 200 gives an accurate solution to this problem, then by comparing the results in Fig. 25, N = 25 also gives sufficient accuracy for practical purposes. Fig. 26 is a plot of the fraction of tertiary oil

X/L

25. Hydrocarbon saturation and solvent mass fraction at 0.3 Original HCPV injected for Pe = 100 using artificialdispersion.N = 200; 0, x N = 25.

Fig.

0.821034 0.0 0.02

0.04

0.06

0.08

0.10

0.12

(AX/L)O.'

Fig. 26. Effect of grid size on fraction tertiary oil recovered and produced solvent fraction after 0.7 Original HCPV injected. 0, x Calculated values.

34

L.C. Young, Spatial appro~imutio~s for simulating &id dis~luceme~rs

produced and the fraction of solvent present in the nonaqueous phase produced at 0.7 OHCPV of injection as a function of grid size. This time is just after solvent production has begun. It is difficult to accurately determine an asymptotic convergence rate for this problem, since the exact solution is not known, However, Fig. 26 shows that the method does converge as the grid is refined. For iV = 25 the errors from Fig. 26 are approximately 0.02 in fractional recovery and 0.01 in solvent fraction. These results show that this method is capable of giving accurate results. For a practical error tolerance, the number of grid points required is similar to that needed for the miscible (single phase) displacement problem studied in Section 3 (see Table 1). These results confirm that the functional form of the artificial dispersion in (52) does not augment the physical dispersion.

The same basic problem as that just described has been simulated in two dimensions using diagonal and parallel grids for a repeated five-spot pattern. Fig. 27 is a plot of cumulative oil recovery (as a fraction of the tertiary oil in place at the start of the simulation) versus the original hydrocarbon pore volumes of solvent injected. The calculated recovery for this problem is sensitive to the value of the Peclet number. An increase in the level of dispersion (decrease in Pe) tends to stabilize the displacement and thereby improve sweep efficiency and recovery. This same trend was observed for the single phase miscible displacement problem (see Fig. 9). By comparing Figs. 9 and 27, the sensitivity to dispersion level appears to be greater for the tertiary displacement. This is because for the tertiary case the recovery is plotted as a fraction of tertiary oil. For the two problems, the sensitivity to Pe is similar when recovery is compared on the basis of original oil in place. In either case the results tend to be quite sensitive to dispersion level. Fig. 28 shows the effect of grid size on recovery for the median and nondivergence methods with diagonal and parallel grids. For coarse grids, the results show a dependence on the method used, i.e., median or nondivergence, and the grid orientation. However, as the grids are refined, these differences diminish. It appears that in the limit as AX/L goes to zero there would be no effect of either the grid orientation or the method used. The finest grids shown in

ORIGINAL HCPV INJECTED

Fig. 27. Sensitivity of two-dimensional

contact miscible tertiary displacement

simulations

to dispersion level.

3s

J AX/L

Fig. 23. Convergence of artificial dispersion approach miscibte tertiary dispIacements in two dimensions. cr-0 Median five-point, diagonal grids; O---- 0 Median five-point, parallel grids; &--A Nondivergence five-point, diagonal grids; A ---A Nondivergence five-point, parallel grids.

for contact

Fig. 29. Effect of grid size on two-d~~e~s~onal contact miscible tertiary displacement calcuIations with diagonal grids. X-X Upstream weighted ninepoint; A-A Nondivergence five-point, a = go = :AX; W Median five-point, a = go = aAX.

Fig. 28 (30 x 30 diagonal. for Pe = 20 and 40 X 40 diagonal for Pe = 40) were used to plot the recovery curves in Fig. 27. For the Pe = 80 case in Fig. 28, a 40 x 40 diagonal grid was used with the median approximation. Based on the results in Fig. 28, the tertiary oil recoveries plotted in Fig. 27 should be accurate to within 0.005 for Pe = 20 and 40 and within 0.02 for Pe = 80. This problem also has been simulated using the upstream weighted nine-point with various grid sizes. Fig. 29 is a plot of cumulative recovery versus grid size. The effect of grid size on the solution is very similar to that observed for the single phase miscible displacement (see Fig. 11). Fig. 29 also contains results calculated using the median and nondivergence five-point methods for the ease when the dispersivity is proportional to the grid size, These results are qualitatively similar to those of the nine-point. As was the case for the single phase miscible displacement, methods which embody an effective dispersion level which is proportional to grid size do not converge as the spatial grid is refined. This occurs because for adverse mobility ratio displacements, the problem is not well posed in the absence of dispersion. When the physical ‘dispersion term is included in the appruximation~ a convergent solution can be obtained_

7. Displacement of a volatile oil component by a condensable gas component This problem

involves aspects of miscible and immiscible flow. Phase changes also occur

during the displacement. The two components in this problem are mutually soluble over a limited concentration range. If the gas component is added to the liquid, the fluid will remain single phase up to a certain bubble point composition. Once two phases are formed, the addition of more gas causes an increase in the fraction of the volume occupied by the gas phase. With the addition of still more of the gas component, a dew point composition is reached, and the mixture becomes a single gas phase. Hall and Geffen [lo] studied this problem to gain an understanding of how liquid propane is displaced by methane gas in the miscible slug displacement process. They studied this type of two-component displacement experimentally and compared the experimental results to an analytical solution. The analytical solution is for one dimension and neglects capillary pressure, dispersion, and the effect of pressure on the fluid properties. Hall and Geffen considered the nonideal mixing case; whereas here constant density is assumed. This problem is governed by (13b) and (15). The dimensionIess form of the equation for the gas component is

where 7 = tu/(&), 5 = X/L, Pe = L/a. In (54) the subscripts on S,f, co and c, have been omitted for convenience. These variables are defined in (lo), (12) and (9), respectively. Equation (54) reflects the assumptions listed above, with the exception that the physical dispersion terms have been included. The data considered for this problem are given in Table 5. For the bubble and dew point compositions of 0.25 and 0.75, respectively, the gas saturation and liquid and vapor compositions are plotted in Fig. 30. In the single phase regions the composition of the phase which is present corresponds to the overall composition, while in the two phase region, the composition of the liquid and vapor are constant and equal to the bubble and dew point compositions, respectively. In the single phase regions, the composition of the phase which is not present is undefined. This fact has important consequences with regard to simulating physical dispersion effects. The gas component fractional flow function and the total mobility are plotted in Fig. 31. From Fig. 31 it is readily apparent that these functions are not smooth at points where a phase transition occurs. Although the total mobility has no influence on the computed compositions in one dimension, it would have an impact on two-dimensional calculations (see Table 3). The lack of smoothness in these functions may degrade the performance of some numerical methods. In the absence of dispersion, the solution to this problem develops two shock fronts which Table 5 Physical properties for condensable

gas displacing voIatile oil

Relative ~~eability

Viscosity

3ubble and dew point composition

k, = (1 - s,/o.7y

/.& = [l - CO+ cO(lo)-‘“I-” cp /A0= [l - cp + c,(lo)-“4]-4 cp

Sbp = 0.25 s+ = 0.75

k,* = S:

0.1

0.2

Fig. 30. The dependence of phase compositions and gas saturation on the crverall volume fraction of gas component for two-component partially miscible disa placement.

travel at speeds ~~op~~t~u~al to the faster front conesponds to the point the sgower front corresponds to the shocks the profik is given by (35).

l.Q

0.8

s

s

Fig. 31. l?ractionaI Row and tatai mobility for partially miscible displacement.

functions

slupe~ of the two tangent lines drawn on Fig. 31. The where condensation to a single phase liquid occurs, whik point where complete vaporization occtfrs, Between the Fig. 32 contains a plot of the analytical solution. Also

Fig, 32. Composition profiles after 0.5 HCPV Injected for partially miscible displacement. N = 50 Upstream weighting. analytical solution; X-X N = 20 Upstream weighting; O---O

Zero dispersion

38

L.C. Young, Spatial approximations for simulating fluid displacements

plotted on Fig. 32 are numerical solutions calculated using upstream weighting with 20 and 50 grid intervals. The median and divergence approximations with (Y= go = $AX are equivalent to upstream weighting for (54). From Fig. 32 the numerical solutions appear to be converging to the zero dispersion analytical solution. Fig. 33 is a plot of the average absolute error in S as a function of grid size. The results in Fig. 33 confirm that the method is converging. The convergence rate is even somewhat higher (0.91) than it was for the immiscible displacement problems (0.75 and 0.77) that were studied in Section 4. For the reasons that have been discussed in the preceding sections, it is highly desirable to include the effects of physical dispersion in this problem. For example, from Fig. 31 the mobility ratios are approximately 1.2 across the condensation front and 6.2 across the vaporization front. From the results presented in Figs. 11 and 29 it is not likely that any method which neglects physical dispersion will converge for this problem in two dimensions. Numerical approximations for the dispersion terms require the calculation of differences in the phase compositions, c, and c,, between neighboring grid points. However, when a phase transition occurs the composition of the phase which disappears is undefined. Therefore, differences in composition of the disappearing phase cannot be directly calculated. For this particular problem, the phase compositions in the two phase region are known to be constant, so the phase composition gradients are zero. Knowledge of this fact could possibly be used to evaluate the dispersion terms. However, if the problem were generalized by allowing the fluid properties to vary with pressure or if more components were considered, the composition gradient in the two phase region would not be known a priori. It might be possible to extrapolate the composition gradient from known values at adjacent grid points or to neglect the composition differences across the grid interval where the phase transition occurs. Neither of these approaches have been attempted. Extrapolation of the gradient could easily encounter difficulties since there is no guarantee that the gradient can always be defined at the neighboring grid interval (this is especially true in two or three dimensions). Even if one of these approaches works, the lack of smoothness of the coefficients (see Fig. 31) may cause the solution to converge slowly.

l/N

Fig. 33. Average grid size.

absolute

error in gas component

overall volume

fraction

after 0.5 HCPV

Injected

as a function

of

The numerical difficulties just described are caused by the fact that the numerical approximation in effect ignores the phase transitions, The phase transitions should be treated as moving internal boundaries. Approaches which track the phase transition have been considered in other disciplines [IS], but apparently have not been attempted for petroleum reservoir problems. To the author’s knowledge, even the fact that these numerical di~~ulties exist has not been reported previously for this class of problems. Due to these difficulties two-dimensianal results for this problem are not reported. Based on the results in Figs. 11 and 29 it is a forgone conclusion that without the ability to simulate physical dispersion effects, two-dimensional calculations will fail to converge as the spatial grid is refined.

The final problem considered is a multiple-contact-miscible (MCM) displacement. The MCM displacement mechanism is operating in many pilot and commercial scale floods which employ CO, or hydrocarbon gases to displace oil or nitrogen to displace gas condensates and light oils at high pressure. Hutchiasun and Braun [I] defined the different MCM processes and the conditions under which they exist. The particular problem considered here is the displacement of a gas condensate by nitrogen. This problem was considered previously by Young and Stephenson [6]. No mobiIe water exists during the displacement. The fluid properties are described using a modified Redlich-Kwong equation of state with seventeen components (nitrogen plus sixteen components to describe the in-place gas). The fluid description predicts nonideal mixing of the components. Therefore, neither the physical problem nor upstream weighting possesses a maximum principal. The displacement consists of nitrogen displaying the in-place gas at 6000 psi. Only onedimensional flow is considered. Fig. 34 is a plot of the nitrogen overall mass fraction and hydrocarbon liquid saturation at 0.2 and 0.6HCPV of injection. This simulation employed upstream weighting with 50 grid intervals. In this case upstream weighting is equivalent to the

Fig. 34. Profiles of liquid phase saturation and nitrogen overall mass fraction after 0.2 and 0.6 HCPV Injected for nitrogen displacing a gas condensate at 6000 psi.

40

L.C. Young, Spatial approximations for simulating @id displacements

Fig. 35. Ternary diagram of pseudo component mass fractions for nitrogen displacing a gas condensate at 6000 psi. Phase envelope; ---Tie lines; . . . . . Overall composition after 0.6 HCPVI. 0 Original fluid composition; -

divergence approximation when LY= go = $AX. Under these conditions the median approximation differs slightly from upstream weighting because the artificial dispersion is of the form in (47b) rather than (47a). Fig. 35 is a pseudo ternary diagram plotted from the simulation results at 0.6 HCPVI. The critical tie line is used to define the displacement mechanism. Since the in-place fluid composition lies to the right of the critical tie line and the injected fluid composition is to the left, the displacement is multiple-contact-miscible with miscibility being generated by a vaporization process [l]. In this system, the nitrogen initially raises the dew point pressure of the in-place gas and condensation occurs. However, after repeated contacts the nitrogen strips the light and intermediate components from the condensed liquid. As a result the gas becomes enriched and takes on more of a solvent-like character. These effects are evident in Fig. 34. After 0.6 HCPVI, the average liquid saturation behind the displacement front is roughly half that at 0.2HCPVI. Analytical solutions have been developed for multicomponent displacements under the assumption of zero capillary pressure and zero dispersion. These analytical solutions predict a piston-type displacement for the condition present in this example [19,20]. Since the upstream weighting method should converge to the zero dispersion solution, the effect of grid size on computed results was investigated. Fig. 36 is a plot of cumulative liquid recovery versus HCPV injected for various grid sizes. The fluid produced in the simulation passes through a series of separators. The cumulative liquid from the separators is plotted in Fig. 36 as a fraction of the total liquid which would result by passing all of the original fluid-in-place through the separators. The results in Fig. 36 exhibit a sensitivity to grid size. Also plotted in Fig. 36 is the recovery curve for a piston displacement.

41

a

PISTONDlSPLACEMENT

iii 100 a aN=50 Y wN=%l o oN=lil

HYDROCARBON

AN

PORE VOLUMES INJECTED

Fig. 36. Effect of grid size on liquid recovery predicted with upstream weighting,

To further investigate the effect of grid size, the fractional recovery after 1.2HCPVI was subtracted from unity and plotted on Fig. 37. This was done in an effort to determine if this difference is, in effect, the error in the numerical approximation. From Fig, 37 it appears that the results are converging toward the zero dispersion (piston displacement) solution at a fairly slow rate. Experience with other problems suggests that the ‘convergence rate’ is highly dependent on the Auid property description, The zero dispersion solution is not the one of interest for this type of displacement. If this were true, there woufd be no reason to incur the expense of simulating numerous components with an equation of state; a two-component contact miscible system would be sufficient. The reason for performing the complex fluid property calculations is to obtain a more accurate representation of the displacement efficiency and the effect of the condensed liquid saturation on the fluid mobilities.

Fig. 37. Convergence displacement solution.

of liquid

recovery

calculated

with upstream

weighting

to the zero dispersion

piston

42

L.C. Young, Spatial approximations for simulating fluid displacements

In both laboratory and field scale displacements a residual heavy hydrocarbon phase is left behind the displacement front. Recoveries from laboratory systems often exhibit a sensitivity to the system length [21,22]. These observations suggest that physical dispersion or mechanical mixing plays an important role in these processes. It appears that in many laboratory systems, the mechanical mixing is due primarily to adverse viscosity ratios. Therefore, for these problems, the solution of interest is the one which includes the dispersion phenomenon. Unfortunately, the phase change which occurs during the displacement causes the same numerical difficulties that were described in the previous section. For the particular problem considered here, the liquid phase is immobile, so the liquid phase dispersion term is zero and can be evaluated. However, for the same displacement at lower pressures or for the displacement of an oil rather than a gas the difficulty of evaluating the dispersion term is again present. Improved numerical methods are needed to adequately solve this type of problem.

9. Summary and conclusions The major conclusions from this work are summarized in the following. (1) Resolution of physical dispersion effects in miscible displacement problems is accomplished more efficiently with central differences or conventional finite element methods than with either conventional or optimal upstream weighting methods. The following effects of dispersion have previously been identified: (a) Dispersion controls the length of mixing or transition zones between various fluid banks. (b) Dispersion increases sweep efficiency in adverse mobility ratio displacements. (c) For multiple-contact-miscible displacements, dispersion reduces displacement efficiency. Consequently, when these factors are important, the resolution of physical dispersion effects is needed. (2) If the physical dispersion term is included in the approximating equations, central differences and conventional finite element methods converge when applied to adverse mobility ratio miscible displacements in both one and two dimensions. If physical dispersion is neglected, the upstream weighted five-point method is sensitive to grid orientation and fails to converge for these problems in two dimensions. The upstream weighted nine-point method is not sensitive to grid orientation; however, it fails to converge for adverse mobility ratio displacements when physical dispersion is neglected. (3) Some form of artificial dispersion is required for most immiscible displacement problems of practical interest. For two-dimensional adverse mobility ratio displacements, a grid orientation sensitivity will occur unless the artificial dispersion term has a functional form which is invariant to a coordinate rotation. An artificial dispersion term of appropriate form has been developed for the five-point difference method or linear Lobatto-Galerkin method. This method yields results that are comparable to those of the upstream weighted nine-point method. (4) The quadratic Lobattc+Galerkin finite element method is more efficient than central difference (or linear Lobatto-Galerkin method) for problems involving only miscible flow.For immiscible displacements, the quadratic method can yield results that are more efficient than upstream weighting; however, an artificial dispersion term (without adjustable parameters) which applies to a wide range of problems is needed to make the method truly useful. This is an area where additional research is warranted.

(5) For problems which involve both miscible and immiscible flow, an artificial dispersion term is needed to complement the physical dispersion term. An appropriate artificial dispersion term has been developed for use with the linear Lobatto-Galerkin method (five-point central differences), The artificial dispersion term was constructed so that it would not act together with the physical dispersion term to produce overly dispersive results. convergence of this method has been demonstrated in one and two dimensions for problems without phase changes. (6) Numerical difficulties are encountered for problems with phase changes. Near the phase change, the numerical approximation of the physical dispersion term cannot be directly calculated. At the point where a phase transition occurs, the coefficients of the governing equations are not smooth. Better numerical techniques are needed for problems of this type. (7) When upstream weighting is used for one-dimensional multiple-contact-miscible displacements, the calculated results converge to the zero dispersion piston displacement solution as the grid is refined.

Appendix A. Difference approximations for (47) Difference approximations for the following equation have been derived using linear finite elements and trapezoidal rule integration: (b 2

= V - khf,Jp

+ V - g&lVfw .

(A.11

For convenience the subscript w will be omitted. difference approximations are as follows:

&Z!L+T iJ dt

Ti.j-l/2(Pij-Pi,j-l)+

+

Di-l&&j

-fi-l,j)+

+

m,j-1/2(hj

-&-1)

the various

Ti+1/2,jtPij_Pi+l,j)

i-1/2,j(pij-pi-1,j)+

+

With this simplified notation,

Ti,j+ll2(Pij -Pi,j+l) ~i~~,2,j~~j

+

-fi+l,j)

~i,j+li2($j_fi,j+l)=

0

*

(A*9

The pore volume is given by & = i&(AXi + AXi+l)(A k; Jr A Yj+l) where AXi = Xi - X;:-1 and

AYj= I$-

Yi_1.

The coefficients, T, of the convection term depend on whether the divergence, nondivergence or median approximation is selected. Each approximation can be represented by, e.g.

but (kAf)i+l/z,j is determined proximations

by a different

averaging

procedure.

The three different

ap-

are:

(1) Divergence

approximation,

(2) Nondivergence

approximation,

(3) Median approximations

The other coefficients, T, are determined dispersion term are as follows:

in an analogous fashion. The coefficients, D, of the

The other coefkients, Q are calculated in an analogous fashion. For the alternate form of artificial dispersion (47a), the f’s in (A.2) are replaced by L’s and in the expression for the coefficient, D, the A’s are omitted. A.I. Sufficient conditions for a maximum principal The approximation to the overall balance is determined equation fur the oil saturation. The resulting approximation convection term is treated,

by adding (A.2) to the analogous is the same, regardless of how the

+ TQ+~/z(p~~ -pi,i+1) = 0 3

64.3)

+

Tfj-l/2(&

-Phi-1)

L.C. Young, Spatial approximationsfor simulating fluid displacements

45

where E+1/2,

/ =

Ax+AYj+l 4AXi+1

(kijhij

+

ki+l,jhi+l.j)

.

The other coe~cients are determined in an analogous fashion. In (A.2), if the coefficients on (fij -fi-l,j), (fij -fi+l,j), (fij -f&i-i) and (fii -fi,j+l) are all positive for every grid point, then the approximation obeys a maximum principal. As an example, the coefficient of (fij -fi+l,i) is examined for the median approximation. To perform the analysis fij times (A.3) is subtracted from (A.2). We note that Tii1n.j - fijE’++l/2.j=

Ar;, i AY.,, i+t’ (~jjA~j+ ki+l,jhi+l,j)(fi+1,j_fi,i) G&AX

*

The following condition is needed on the coefficient of (fti - fi+l,j): Di+llZ.j

-

AY+AY.,, $+*’ (k&j + ki+l,jAi+l,j)(pij - pi+l,j) a 0 . 8AX

We note that I)i+l/z.i can be bounded Di+1/2,j

by

~&A~~~,yi’l (kij&j + ki+l,jAi+l,j)lpi+l,j-pijl r+l

,

so the following condition is sufficient for a positive coefficient:

It should be noted that for an irregular grid, a maximum principal can be maintained by using a different value of go in each element and making its value equal to the maximum of $AX, ;AY. For highly irregular grids, this choice of go can cause the dispersion terms to dominate convection in the direction of short grid spacing. This may impose stability restrictions if the dispersion terms are evaluated explicitly. The use of a different value of go in the two directions may be warranted in this case. This approach would cause the dispersion terms to depend on coordinate rotation; however, highly irregular grids normally occur in regions of the grid that are not of major interest.

Acknowledgment I thank John L. Yanosik encouragement.

and the late Howard

N. Hall for their support,

advice and

46

L.C. Young, Spatial approximations for simulating fluid displacements

References [l] C.A. Hutchinson, Jr. and P.H. Braun, Phase relations of miscible displacement in oil recovery, AIChE J. 7 (1961) 64-72. [2] J.L. Yanosik and T.A. McCracken, A nine-point, finite-difference reservoir simulator for realistic prediction of adverse mobility ratio displacements, Sot. Pet. Engrs. J. (1979) 253-262; Trans. AIME 257. [3] L.C. Young, A finite-element method for reservoir simulation, Sot. Pet. Engrs. J. (1981) 115-127. [4] T.K. Perkins and O.C. Johnston, A review of diffusion and dispersion in porous media, Sot. Pet. Engrs. J. (1963) 70-84. [5] J.W. Gardner, F.M. Orr and P.D. Patel, The effect of phase behavior on COZ flood displacement efficiency, J. Pet. Tech. (1981) 2067-2081. [6] L.C. Young and R.E. Stephenson, A generalized compositional approach for reservoir simulation, Sot. Pet. Engrs. J. (1983) 727-742. [7] J. Bear, Dynamics of Fluids in Porous Media (American Elsevier, New York, 1972) 613. [8] M. Delshad, D.J. MacAllister, G.A. Pope and B.A. Rouse, Multiphase dispersion and relative permeability experiments, Paper SPE 10201 presented at the 56th SPE Annual Technical Conference and Exhibition, San Antonio, TX, 1981. [9] S.J. Salter and K.K. Mohanty, Multiphase flow in porous media: I. Macroscopic observations and modeling, Paper SPE 11017 presented at 57th SPE Annual Technical Conference and Exhibition, New Orleans, LA, 1982. [lo] H.N. Hall and T.M. Geffen, A laboratory study of solvent flooding, Trans. AIME 210 (1957) 48-57. [ll] P.M. Gresho and R.L. Lee, Don’t suppress the wiggles-they’re telling you something, in: Finite Element Methods for Convection Dominated Flows, AMD Vol. 34 (ASME, New York, 1979). [12] A. Settari, H.S. Price and T. DuPont, Development and application of variational methods for simulation of miscible displacement in porous media, Sot. Pet. Engrs. J. (1977) 228-246; Trans. AIME 263. [13] H.J. Morel-Seytoux, Introduction to flow of immiscible liquids in porous media in: R.J.M. Dewiest, ed., Flow Through Porous Media (Academic Press, New York, 1969) 455-516. [14] S.E. Buckley and M.C. Leverett, Mechanism of fluid displacement in sands, Trans. AIME 146 (1942) 107-116. [15] H.J. Welge, A simplified method for computing oil recovery by gas or water drive, Trans. AIME 195 (1952) 91-98. [16] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Interscience, New York, 1967) 306308, 337. [17] P.L. Terwilliger, L.E. Wilsey, H.N. Hall, P.M. Bridges and R.A. Morse, An experimental and theoretical investigation of gravity drainage performance, Trans. AIME 192 (1951) 285-295. [18] D.R. Lynch, Unified approach to simulation on deforming elements with application to phase change problems, J. Comput. Phys. 47 (1982) 387-411. [19] H.J. Welge, E.F. Johnson, S.P. Ewing and F.H. Brinkman, The linear displacement of oil from porous media by enriched gas, J. Pet. Tech. (1961) 787-796; Trans. AIME 222. [20] F.G. Helfferich, Theory of multicomponent multiphase displacement in porous media, Sot. Pet. Engrs. J. (1981) 51-62. [21] J.J. Rathmell, F.I. Stalkup and R.C. Hassinger, A laboratory investigation of miscible displacement by carbon dioxide, Paper SPE 3483 presented at 46th SPE Annual Technical Conference and Exhibition, New Orleans, LA, 1971. [22] W.F. Yellig, Carbon dioxide displacement of a West Texas reservoir oil, SOC.Pet. Engrs. J. (1982) 805-815.