Applying the finite--volume method for solving the convection--dispersion equation in radial coordinates.

Abstract

The finite-control volume technique is used to solve the convection-dispersion equation in radial coordinates. In the dispersion model presented, the dispersion coefficient is dependent on both velocity and diffusion coefficient. The finite-control volume scheme is shown to have negligible numerical dispersion. The code is verified by performing a convergence study to show that the numerical results are independent of the mesh and time-step sizes.

The solution allows analyzing the dispersion transport of a slug in a steady radial flow from an injection well fully penetrating a homogeneous reservoir of uniform thickness and finite arial extent. Results show that the dispersive mixing zone does not necessarily grow in proportion with the square root of time for all times, unlike for the linear displacement. Results reveal a typical minimum of the dispersive mixing zone with time for radial displacement. It appears that at some radial distance close to the wellbore, the dimensionless mixing zone decreases with the square root of time since there is an accumulation of the miscible slug in this particular space domain. However, when the dissolution of the slug starts taking effect particularly at large distances away from the injection point, the mixing zone begins to increase with increasing square root of time. These results indicate that the optimal slug size of miscible displacement is time-dependent.

At a fixed time value, the higher the value of dispersivity, the more mixing takes place, and hence the greater the concentration of the miscible fluid becomes. A higher dispersivity value aggravates the dissolution of the slug at the toe of the concentration profile, but gives rise to a larger concentration at the tail. The smaller the porosity value is, the larger the rate of change of concentration. This is expected since a smaller value of porosity acts like a choke on the porous medium and results on a greater degree of convective mixing than diffusion mixing.

Keywords: dispersion, porous media, finite-volume method.

Introduction

Dispersion is the mixing of a solvent with in-situ hydrocarbon caused by the combined effects of diffusion and convection. Convection is induced by variation in the velocity field caused by the tortuous flow path of the porous medium. Diffusion occurs as a consequence of a concentration gradient. Dispersion phenomena arise in a variety of applications such as heat flow in material, transport of pollutants in aquifers and lakes, and miscible displacement in hydrocarbon saturated porous media. [1] Dispersion is one of the most prominent dissipative effects in miscible displacements. It affects the amount of solvent mixing with the in-situ oil to promote miscibility under favorable conditions. Dispersion may be detrimental if it is convection dominated, or beneficial if it is diffusion dominated since it prevents fingering through the miscible zone.

When one fluid is miscibly displacing another in a porous medium and when the displacement is stable, the following convection-dispersion equation describes the overall transport and mixing of incompressible fluids flowing through a porous medium:

[phi] [partial derivative]C/[partial derivative]t + u [nabla]C = [phi][nabla]. (K[nabla]C) (1)

In this notation, C is the concentration which is a function of time (t) and position, u is the superficial velocity (Darcy velocity), and K is the dispersion coefficient. For the linear system with constant dispersion coefficient, various analytical solutions based on the error function are established. These solutions differ according to the boundary conditions imposed. Chase [2] presented a finite-element numerical solution for the physical problem discussed. However, Chase [2] concluded that, even with a finite-element method, significant numerical dispersion could still exist.

Numerical solutions to fluid-flow equations, usually nonlinear partial-differential equations, are generally obtained by making discrete approximations to derivatives. Whether finite-difference or finite-element methods are used, these approximations always introduce truncation errors that often can distort the accuracy and stability of the solution. The truncation error is often referred to as numerical dispersion because, to lowest order, it can be represented as a second spatial derivative term added to any true dispersion term in the problem. Artificial smearing as a result of numerical dispersion can render the numerical solution meaningless. [3-5]

For miscible displacement in radial geometry, a number of authors interested in the stability of displacements in porous media solved the convection dispersion equation analytically and numerically with various degrees of simplicity. Tan and Homsy [6] employed a quasi-steady-state analysis. They showed that the unfavorable viscosity gradient resulted in an algebraically growing perturbation, rather than the exponential behavior known from rectilinear displacements. Pankiewitz and Meiburg [7] extended this analysis to fluid combinations that give rise to non-monotonic viscosity profiles, which they found to be destabilizing. Riaz and Meiburg [8] reported on the stability of axial and helical perturbations in three-dimensional displacements. Yortsos [9] incorporated the effects of equilibrium adsorption and showed the existence of a mathematical transformation that related radial flows to rectilinear ones. Tang and Peaceman [5] reported both analytic and a numerical solutions for the convection-dispersion problem for the case when the dispersion coefficient was only dependent on velocity but did not depend on diffusion. Their numerical solution which used a time-implicit, centered-in-space difference scheme gave results that compared almost exactly with their analytic solution. The emphasis of their study was on the type of boundary condition used at the wellbore, though. They reported that the concentration profiles were insensitive to whether the constant concentration, or the flux boundary condition was used.

The purpose of this research is to apply the finite-volume technique [10] (Patankar technique) to obtain a numerical solution for the radial convection-dispersion equation. The work of Tang and Peaceman [5] is extended by including the effects of both velocity and diffusion in the dispersion coefficient. An analytic solution for this problem becomes very difficult to obtain because a number of coefficients of the convection-dispersion partial differential equation are variable. The application of Patankar technique for solving this problem associated with miscible displacement becomes unique since it eliminates the effects of numerical dispersion. The new formulation allows to study the effects of petrographic properties like porosity, tortuosity, dispersivity on the dispersion phenomenon in general.

Model Formulation

For the linear system with constant dispersion coefficient, various analytical solutions, based on the error function, are established for Equation (1). These solutions differ according to the boundary conditions imposed. [1,11] In radial coordinates, Equation (1) is transformed into the following mathematical statement:

[phi][partial derivative]C/[partial derivative]t + u [partial derivative]C/[partial derivative]r = [phi]/r [partial derivative]/[partial derivative]r([rK.sub.r] [partial derivative]C/[partial derivative]r) (2)

subject to the initial condition

C = 0, for t = 0 for all radii (3)

With the far boundary condition

C = 0, when r [right arrow] [infinity] at all time (4)

And the inlet flux boundary condition:

2[pi][Hr.sup.*.sub.w](uC - [alpha]u [partial derivative]C/[partial derivative]r) = [iC.sub.J] (5)

where H is the reservoir thickness, [r.sup.*.sub.w] is the wellbore radius, [alpha] is the dispersivity, and "i" is the injection rate. In this analysis, we consider only the continuous injection case:

[C.sub.J] = 1 (6)

For this simple geometry, the Darcy velocity at any radius is given by the injection rate divided by the cross sectional area available for flow as follows:

u = i/2[pi]Hr (7)

In this analysis, the dimensionless variables are defined as:

[t.sub.D] = [integral]idt/[V.sub.p] (8)

[C.sub.D] = C - [C.sub.I]/[C.sub.J] - [C.sub.I] (9)

[r.sub.D] = r/[r.sub.e] (10)

In this notation [C.sub.I] is the initial concentration and [C.sub.I] is the injected concentration. Since r.sup.*.sub.e] >> [r.sup.*.sub.w], the following approximation for the pore volume is used:

[V.sub.P] [approximately equal to] [pi][r.sup.*2.sub.e] H [phi] (11)

Taking the derivatives of Equations (8) and (9), and substituting expressions for the derivatives (dt) and (dr) into Equation (2), and using an approximation for the pore volume given by Equation (11), we obtain the following:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

The above formulation identifies a non-dimensional group which may be termed as the Local Peclet number given by

[NL.sub.pe] = ru/[phi][K.sub.r] (13)

Further development of the right-hand-side derivative of Equation (12) leads to the following non-dimensional convection-diffusion equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

Equation (14) identifies a second non-dimensional term [N.sub.R] given by

[N.sub.R] = 1/[r.sub.D][K.sub.r] [partial derivative][K.sub.r]/[partial derivative][r.sub.D] (15)

[N.sub.R] is a non-dimensional rate of change of radial dispersion coefficient with respect to radius. Both [NL.sub.pe] and [N.sub.R] are a function of radius. An average Peclet number, which is defined as the ratio of convective to dispersive transport, is therefore deduced from Equation (13), and is given by the following expression:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

Similarly, an average non-dimensional rate of change of radial dispersion coefficient with respect to radius ([[bar.N].sub.R]) is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

Substituting the expression of NR from Equation (15) into Equation (17), and using the approximation ([r.sup.*2.sub.e] - [r.sup.*2.sub.w]) [approximately equal to] [r.sup.*2.sub.e] in the limit when [r.sup.*.sub.e] is large compared to [r.sup.*.sub.w], we obtain the following expression:

[[bar.N].sub.R] = ln [K.sub.re]/[K.sub.rw] (18)

where [K.sub.re] and [K.sub.rw] are the dispersion coefficients at [r.sup.*.sub.e] and [r.sup.*.sub.w], respectively. In general, the radial dispersion coefficient [K.sub.r] is defined by

[K.sub.r] = [D.sub.m]/[phi][F.sup.*] + [alpha]u (19)

In this notation, [D.sub.m] is the molecular diffusion coefficient, [F.sup.*] is the resistivity formation factor, and [alpha] is the dispersivity. We use Pirson's model [12] that relates tortuosity to the resistivity formation factor and rock porosity to rewrite Equation (19).

Adapting the Finite-Volume Methodology

In order to describe the finite-volume method in a classical way, we first use the general form of the diffusion-convection equation. Later, we apply this method on the governing equation of this manuscript. The generic conservation equation [10] for a general quantity [eta] is written as follows:

[partial derivative]([rho][eta])/[partial derivative]/[partial derivative][x.sub.i] ([rho][u.sub.i][eta])-[partial derivative]/[partial derivative][x.sub.i]([GAMMA] [partial derivative][eta]/[partial derivative][x.sub.i]) = S' (20)

Conservation of mass, momentum, and energy can be formulated in the form of the above equation, for which S' stands for the source term and [GAMMA] stands for the diffusion coefficient.

In this method, the first step is the integration of the transport equation for quantity [eta] over a three-dimensional control volume V without any approximation. For simplicity, the method is illustrated by taking the steady-state case of Equation (20):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

Next, Gauss's divergence theorem is applied to the left-hand-side of Equation (21) as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

In this notation, [n.sub.i] is the unit vector normal to the surface element dA. This theorem transforms volume integrals of divergence terms into surface integrals of fluxes all around the control volume. Applying Gauss' theorem leads to the following integrated conservation equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (23)

In general, this equation cannot be solved analytically and a numerical solution will require that we need to transform it into an algebraic equation. This transformation is called the discretization. It requires approximations to evaluate integrals and to perform interpolations. Next, the procedure is illustrated using cylindrical two-dimensional grid represented for the governing equation of this convection dispersion model (Equation (2)) which can be rearranged as follows:

[partial derivative]C/[partial derivative]t + u/[phi] [partial derivative]C/[partial derivative]r = 1/r [partial derivative]/[partial derivative]r (r[K.sub.r] [partial derivative]C/[partial derivative]r) (24)

From Equation (7), it is clear that

[partial derivative]/[partial derivative]r (r u/[phi]) = 0 (25)

By using Equation (25), Equation (24) can be re-arranged as follows:

[partial derivative]C/[partial derivative]t + 1/r [partial derivative]/[partial derivative]r (r u/[phi] C) = 1/r [partial derivative]/[partial derivative]r (r [k.sub.r] [partial derivative]C/[partial derivative]r (27)

The general transport equation, for planar, axisymmetric, radial flow where the fluid properties depend on the radius r only, is of the form:

[partial derivative]([rho]C)/[partial derivative]t + 1/r [partial derivative]/[partial derivative]r (r [rho] [u.sub.r]C) = 1/r [partial derivative]/[partial derivative]r(r [k.sub.r] [partial derivative]C/[partial derivative]) (2)

where [u.sub.r] is the intersticial velocity and [rho] is the fluid density. It is clear that

Equation (26) is similar to Equation (27) with [u.sub.r] = u/[phi] and [rho] = 1. Using Equation (7), we can check that the continuity equation [partial derivative][rho]/[partial derivative]t + 1/r [partial derivative]/[partial derivative]r (r [u.sub.r]) = 0 is automatically satisfied. We re-write Equation (26) as follows:

[partial derivative]C/[partial derivative]t + 1/r [partial derivative]/[partial derivative]r(r u/[phi] C - r [k.sub.r] [partial derivative]C/[partial derivative]r) = 0 (28)

We denote [J.sub.r] = r u/[phi] C - r [k.sub.r] [partial derivative]C/[partial derivative]r. The variable [J.sub.r] denotes the radial flux per unit

of circumferential angle. It includes both the convective and diffusive fluxes. Using [J.sub.r], Equation (28) can be re-written as follows:

[partial derivative]C/[partial derivative]t + 1/r [partial derivative][J.sub.r]/[partial derivative]r = 0 (29)

Since the flow is planar, axisymmetric and radial, and the fluid properties depend only on the radius r, we need only to consider a radial sector from r = [r.sup.*.sub.w] to r = [r.sup.*.sub.e] and spanning a small angle [DELTA][theta]. A set of equally spaced N grid points are selected along the mid-line of the sector. Grid points 1 and N correspond to the nodes at r = [r.sup.*.sub.w] (the inner radius) and r = [r.sup.*.sub.e] (the outer radius), respectively. Grid points i such that N 1 < i < N correspond to the internal nodes. Control volumes are generated around each grid point. The first and last control volumes are the boundary control volumes. The remaining N-2 control volumes are the internal ones. A representative internal control volume, with center grid point P, is shown in Figure 1. The grid points E and W are the east and west neighbors, respectively, of the grid point P, in the radial direction. The grid points N and S are the north and south neighbors, respectively, of the grid point P, in the circumferential direction. The letters e, w, s and n represent the control volume faces. The control volume is assumed to be of unit thickness in the z-direction (out of the page).

To obtain the discrete equations of the concentration C, we will integrate Equation (29) over the control volume shown in Figure 1 and over the time interval from time t to time t + [DELTA]t, with [DELTA]t being the integration time step. The values of C at time t, i.e. the old values of C, will be tagged with a superscript 0. The values of C at time t + [DELTA]t, i.e. the new values of C, will be tagged with a superscript 1. Integrating Equation (29), we obtain the following expression:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (30)

To carry out the integration in Equation (30), we assume that at any time t, the concentration C at P, [C.sub.P], prevails throughout the control volume. Therefore, Equation (30) is re-arranged as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

[FIGURE 1 OMITTED]

where [DELTA]r = [r.sub.e] - [r.sub.w] and [DELTA][theta] = [[theta].sub.n] - [[theta].sub.s]. In conformity with the fully implicit scheme, we assume that the value [C.sup.1.sub.P] prevails throughout the control volume over the entire time step [DELTA]t. Similar assumptions are also made for the neighbor control volumes. Equation (31) can hence be re-written as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

where [J.sup.1.sub.e] = [J.sup.1.sub.r] (r = [r.sub.e]) and [J.sup.1.sub.w] = [J.sup.1.sub.r] (r = [r.sub.w]) are the total fluxes of the concentration per unit of circumferential angle, at time t + [DELTA]t, at the east and west faces of the control volume, respectively. Substituting Equation (31) and Equation (32) into Equation (30) and dividing by [DELTA][theta] [DELTA]t, we obtain the following:

([C.sup.1.sub.P] - [C.sup.0.sub.P) [r.sub.P] [DELTA]r/[DELTA]t + [J.sup.1.sub.e] - [J.sup.1.sub.w] 0 (33)

From here on, we omit the superscript 1, i.e. variables without superscript refer to the corresponding values at the new time t + [DELTA]t. Equation (33) can be re-written as

([C.sub.P] - [C.sup.0.sub.P]) [r.sub.P] [DELTA]r/[DELTA]t + [J.sub.e] - [J.sub.w] = 0 (34)

Applying the continuity equation, in the case of incompressible, planar, axisymmetric radial flow, for our control volume, we obtain the following:

[r.sub.e] [DELTA][theta] [u.sub.e]/[phi] = [r.sub.w] [DELTA][theta] [u.sub.w]/[phi] (35)

where [u.sub.e] = u(r = [r.sub.e]) and [u.sub.w] = u(r = [r.sub.w]). We define F = r v/[phi] and D = [k.sub.r] r/[DELTA]r . F indicates the strength of convection and D is a diffusion conductance term. [F.sub.e] and [D.sub.e] are the values of F and D at the east face of the control volume, respectively. [F.sub.w] and [D.sub.w] are the corresponding values at the west face of the control volume. Equation (35) can be re-written as

[F.sub.e] = [F.sub.w] (36)

Using Equation (36), Equation (34) can be re-written as

([C.sub.P] - [C.sup.0.sub.P]) [r.sub.P] [DELTA]r/[DELTA]t + ([J.sub.e] - [F.sub.e] [C.sub.P]) - ([J.sub.w] - [F.sub.w] [C.sub.P]) = 0 (37)

The Patankar Peclet number [??] is given ad [??] = F/D. Patankar [10] introduced [J.sup.*.sub.e] = [J.sub.e]/[D.sub.e]

and [J.sup.*.sub.w] = [J.sub.w]/[D.sub.w]. Using these latter expressions, Equation (37) can be re-written as

follows:

([C.sub.P] - [C.sup.0.sub.P]) [r.sub.P] [DELTA]r/[DELTA]t + ([J.sup.*.sub.e] - [[??].sub.e] [C.sub.P])[D.sub.e] - ([J.sup.*.sub.w] - [[??].sub.w] [C.sub.P])[D.sub.w] = 0 (38)

Patankar [10] proved that for the uniaxial case with constant conductance, at the face between two consecutive grid points i and i+1, the following relationship:

[J.sup.*] = B x [C.sub.i] - A x [C.sub.i+l] (39)

where A and B are dimensionless coefficients that are functions of the Peclet number [??]. The same result can be approximated fairly accurately to our case of radial flow provided that [DELTA]r/r is kept small enough. Patankar [10], also, proved that

B = A + [??] (40)

and

A([??]) = A ([absolute value of [??]])+ Max (-[??], 0) (41)

and

B([??]) = A ([absolute value of [??]])+ Max ([??], 0) (42)

Applying Equation (39) and Equation (40) to the east face of our control volume, i.e. nodes i and i+1 corresponding to the P and E grid points, respectively, we obtain the following:

[J.sup.*.sub.e] - [[??].sub.e] [C.sub.P] = (A([[??].sub.e])+[[??].sub.e]) [C.sub.P] - A([[??].sub.e]) [C.sub.E] - [[??].sub.e] [C.sub.P] = A([[??].sub.e])([C.sub.P] - [C.sub.E]) (43)

Similarly, Equations (39) and (40) are applied to the west face of the control volume, i.e. nodes i-1 and i corresponding to the W and P grid points, respectively. The following expression is obtained:

[J.sup.*.sub.e] - [[??].sub.e] [C.sub.P] = (A([[??].sub.e])+[[??].sub.e]) [C.sub.P] - A([[??].sub.e]) [C.sub.E] - [[??].sub.e] [C.sub.P] = A([[??].sub.e])([C.sub.P] - [C.sub.E]) (43)

By substituting Equation (43) and Equation (44) into Equation (38), Equation (38) is rearranged as follows:

([C.sub.P] - [C.sup.0.sub.P]) [r.sub.P] [DELTA]r/[DELTA]t + A([??].sub.e])[D.sub.e]([C.sub.P] - [C.sub.E]) - B([[??].sub.w] [D.sub.w] ([C.sub.w] - [C.sub.P]) = 0 (45)

Using Equations (41), (42), and the definition of the Peclet number, we obtain:

A([[??].sub.e])[D.sub.e] = (A([[absolute value of [??]].sub.e])+ Max(-[[??].sub.e],0))[D.sub.e] = A([[absolute value of [??]].sub.e])[D.sub.e] + Max(-[F.sub.e],0) (46)

and

B([[??].sub.w])[D.sub.w] = (A([[absolute value of [??]].sub.w])+ Max(-[[??].sub.w],0))[D.sub.w] = A([[absolute value of [??]].sub.w])[D.sub.w] + Max([F.sub.w],0) (47)

Using Equations (46) and (47), Equation (45) can be re-arranged in the following general format:

[a.sub.P] [C.sub.P] = [a.sub.E] [C.sub.E] + [a.sub.W] [C.sub.W] + b (48)

where

[a.sub.E] = (A([[absolute value of [??]].sub.e])[D.sub.e] + Max(-[F.sub.e],0) (49a)

[a.sub.W] = (A([[absolute value of [??]].sub.w])[D.sub.w] + Max([F.sub.w],0) (49b)

[a.sup.0.sub.P] = [r.sub.P] [[DELTA]r/[DELTA]t (49c)

b = [C.sup.0.sub.P] [a.sup.0.sub.P] (49e)

[a.sub.P] = [a.sub.E] + [a.sub.W] + [a.sup.0.sub.P] (49e)

The choice of A([absolute value of [??]]) depends on the scheme to be used. [10] Our results are obtained using the Exponential (exact) scheme, given by:

A([[absolute value of [??]]=[absolute value of [??]]/exp[absolute value of [??]]-1 (50)

A total of N-2 linear, algebraic equations of type Equation (48) are generated for the N-2 internal nodes. The remaining 2 equations are obtained from the boundary conditions. The boundary condition at r = [r.sup.*.sub.w] is given by Equation (5) and implies the following:

(1 + [alpha]/dr)[C.sub.1] - a/dr [C.sub.2] = 1 (51a)

and the boundary condition at r = [r.sup.*.sub.e] gives

[C.sub.N] = 0 (51b)

Equations (48), (51a) and (51b) constitute a system of N linear, algebraic equations and N unknown C values corresponding to the concentrations at the present time (the old values of C are assumed known). This system of equations is solved using the TriDiagonal Matrix Algorithm (Patakar, 1980). [10]

The procedure described above allows the calculation of the concentration solution at time t + [DELTA]t given the concentration solution at time t . Starting from the known concentration values at time t = 0, i.e. [C.sub.i] = 0, for 1 [less than or equal to] i [less than or equal to] N, the concentration solution is obtained at time t = [DELTA]t. Next, the concentration solution is computed at t = 2 [DELTA]t, using the concentration solution at time t = [DELTA]t. The procedure can be repeated till the desired final time is reached.

[FIGURE 2 OMITTED]

Discussion

Profiles of concentration of injected miscible fluid are presented as a function of dimensionless time and radius in Figures 2 through 10. The space domain has been divided into 600 increments whereas the time step was updated every one second. A time step less than one second, and a discretization of the space domain in more than 600 increments did not change the concentration profiles.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

Figure 2 shows the profiles of concentration versus radius for early, intermediate, and late times. At an early time, the dimensionless concentration at the wellbore face is much less than unity. At an intermediate time this concentration increases progressively to a value of one. At late times, the displacement front is well established and keeps moving to the reservoir external boundary with a degree of smearing that increases with increasing radius. This front dissipation is more visible in Figure 3. At a fixed time, closer radii to the wellbore face have a sharper front than those radii farther away in the direction of the reservoir boundary.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

In line with expectation, the rate of change of concentration with time decreases as the dimensionless radius value increases (Figure 3). Figure 4 illustrates the effect of dispersivity ([??]) on the concentration profile at a fixed radius value. At a fixed time value, the higher the value of dispersivity, the more mixing takes place, and hence the greater the concentration of the miscible fluid becomes. Figure 5 reveals the effect of this mixing on the displacement front. A higher dispersivity value aggravates the dissolution of the slug at the toe of the concentration profile, but gives rise to a larger concentration at the tail.

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

Figure 6 displays the effect of porosity on the concentration profile. The smaller the porosity value is, the larger the rate of change of concentration. This is expected since a smaller value of porosity acts like a choke on the porous medium and results therefore on a greater degree of convective mixing than diffusion mixing. This concept is even reasserted in Figure 7 that shows a snapshot of the displacement front behavior. The smaller the value of porosity is, the steeper the displacement front becomes. Figures 8 and 9 display the effect of volumetric flow rate on the concentration history and profile, respectively. The higher the injection rate value is, the faster the buildup of the displacement front is. This observation indicates that the influence of velocity induced dispersion is more important for large rate values for this homogeneous isotropic porous medium, in line with expectation. The effect of tortuosity on the front history is displayed in Figure 10..

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

Low tortuosity values can have a slight stabilizing influence on the displacement. At a relatively low injection rate, the tortuosity value did not appear to influence the rate of change of concentration, though

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

Figures 11 through 14 display plots of the dimensionless mixing zone as a function of the square root of time. The plots vary the dispersivity, injection rate, porosity, and the tortuosity values, respectively. These Figures show that the dispersive mixing zone does not necessarily grow in proportion with the square root of time for all times, unlike the linear displacement. These results reveal a typical minimum of the dispersive mixing zone with time for radial displacement. It appears that at some radial distance close to the wellbbore, the dimensionless mixing zone decreases with the square root of time since there is an accumulation of the miscible slug in this particular space domain.

[FIGURE 15 OMITTED]

[FIGURE 16 OMITTED]

However, when the dissolution of the slug starts taking effect, particularly at large distances away from the injection point, the mixing zone begins to increase with increasing square root of time. This result may be inferred from a plot of dispersive mixing zone as a function of the average Peclet number (Figure 15). At an early time, the dispersive mixing zone displays a minimum locus. However, at late times there is an apparent monotonic increase of the dispersive mixing zone. The latter result is consistent with an increasing growth rate of instability with increasing Peclet number. These results indicate that the optimal slug size of miscible displacement is time-dependent.

[FIGURE 16b OMITTED]

[FIGURE 16c OMITTED]

[FIGURE 16d OMITTED]

Figures 16a through 16d display a comparison between the solution of the model presented in this study and Tang and Peaceman [5] solution. At relatively high injection rate (0.1 [m.sup.3]/sec to 0.01 [m.sup.3]/sec) there is a reasonable agreement between the two solutions. However, as the injection rate decreases (from 0.01 [m.sup.3]/sec to 0.0001 [m.sup.3]/sec) the disparity between the two solutions increases. This is not a surprising result since the dispersion coefficient, in Tang and Peaceman model, is only a function of velocity and does not account for diffusion effects; Whereas, in the model presented in this study, the dispersion coefficient accounts for diffusion as well. At relatively low flow rates, the mixing between the insitu crude and solvent become diffusion dominated, a fact that aggravates the disparity between our solution and Tang and Peaceman solution.

Conclusions

The finite-control volume technique is used for optimizing the accuracy of the numerical solution to the convection-dispersion equation in radial coordinates. This technique is favored because it is known to eliminate the detrimental effects of numerical dispersion. The convection-dispersion equation arises in the simulation of Enhanced Oil Recovery (EOR) processes where front smearing resulting from numerical dispersion can distort the behavior of the miscible displacement process. This front smearing can, for instance, lead to false predictions of breakthrough times. Our solution should allow to better estimate the optimum slug size for miscible displacement projects in hydrocarbon reservoirs. This in return allows a better estimate of the total amount of solvent needed to maintain miscibility conditions at the displacement front in the bulk of the reservoir, leading to better reservoir management and to a more adequate economic forecast of EOR processes.

The optimal slug size appears to be time-dependent. The profile of dimensionless mixing zone versus the square root of dimensionless time displays a minimum locus. This result appears to be unique for radial displacement.

Acknowledgments

This work was supported by Kuwait University Research Grant No. [EP01/05].

References

[1] Lake, L.W., 1989, Enhanced Oil Recovery, Prentice Hall: Englewood Cliffs, New Jersey.

[2] Chase, C.A., 1971, "Finite-Element Analysis of Single Well Backflow Tracer Test in a Homogeneous Reservoir," SPE Annual Meeting, New Orleans, October 1971; Paper 3485.

[3] Peaceman, D.W., 1977, Fundamentals of Numerical Reservoir Simulation, Elsevier Scientific Publishing Company, New York.

[4] Settari, A., and Price, H.S., 1977, "Development and Application of Variational Methods for Simulation of Miscible Displacement in Porous Media," SPE Journal, June issue, p. 228.

[5] Tang, D.H.E., and Peaceman, D.W., 1987, "New Analytical and Numerical Solutions for the Radial Convection-Dispersion Problem," SPE Reservoir Engineering, August Issue, p. 343.

[6] Tan, C.T., and Homsy, G.M., 1987, "Stability of Miscible Displacement in Porous Media: Radial Source Flow," Phys. Fluids, 30, p. 1239.

[7] Pankiewitz, C., and Meiburg, E., 1999, "Miscible Porous Media Displacements in the Quarter Five-Spot Configuration. Part 3. Non-Monotonic Viscosity Profiles," Journal of Fluid Mechanics, 388, p.171.

[8] Riaz, A., and Meiburg, E. Radial Source Flows in Porous Media: Linear Stability Analysis of Axial and Helical Perturbations in Miscible Displacements. Phys. Fluids 2003, 15, 938.

[9] Yortsos, Y.C. Stability of Displacement Processes in Porous Media in Radial Flow Geometries. Phys. Fluids 1987, 30, 2928.

[10] Patankar, S.V. Numerical Heat Transfer and Fluid Flow, Taylor and Francis, 1980.

[11] Peters, E .J, Gharbi, R., and Afzal, N. A Look at Dispersion in Porous Media Through Computed Tomography Imaging. Journal of Petroleum Science and Engineering 1996, 15, 23.

[12] Garrouch, A.A., Ali, L., and Qasem, F. Using Diffusion and Electrical Measurements to Assess Tortuosity of Porous Media. Industrial Engineering and Chemistry Research Journal 2001, 40, 4363.

A. A. Garrouch

Petroleum Engineering Department, Kuwait University, P.O. Box 5969, Safat 13060, Kuwait

M. M. Al-Dousari

Petroleum Engineering Department, Kuwait University, P.O. Box 5969, Safat 13060, Kuwait

Nomenclature

A                 dimensionless coefficient
B                 dimensionless coefficient
C                 concentration
[C.sub.D]         dimensionless concentration
[C.sub.I]         initial concentration
[C.sub.J]         injected concentration
[C.sub.P]         concentration C at point P
dA                surface element
[D.sub.e]         diffusion conductance value at the east face of the
                  control volume
[D.sub.m]         molecular diffusion coefficient
[F.sup.*]         resistivity formation factor
[F.sub.e]         convection strength value at the east face of the
                  control volume
H                 reservoir thickness
[J.sub.r]         radial flux per unit of circumferential angle
i                 volumetric injection rate
K                 dispersion coefficient
[K.sub.r]         radial dispersion coefficient
[K.sub.re]        dispersion coefficients at [r.sup.*.sub.e]
[K.sub.rw]        dispersion coefficients at [r.sup.*.sub.w]
N                 number of grid points
[n.sub.i]         unit vector normal to the surface element dA
[NL.sub.pe]       local Peclet number
[N.sub.pe]        average Peclet number
[N.sub.R]         local non-dimensional rate of change of radial
                  dispersion coefficient with respect to radius
[[bar.N].sub.R]   average non-dimensional rate of change of radial
                  dispersion coefficient with respect to radius
P                 Patankar definition of Peclet number
r                 radius
[r.sub.D]         dimensionless radius
[r.sub.P]         radius at P

[r.sup.*.sub.e]   drainage radius

[r.sup.*.sub.w]   wellbore radius
S'                source term
t                 time
[t.sub.D]         dimensionless time
u                 superficial velocity
[u.sub.e]         superficial velocity value at the east face of the
                  control volume
[u.sub.r]         interstitial velocity
[u.sub.w]         superficial velocity value at the west face of the
                  control volume V three-dimensional control volume
[V.sub.P]         pore volume

Greeck            Symbols

[??]              dispersivity
[phi]             porosity
[??]              general quantity
[GAMMA]           diffusion coefficient
[??]              angle
[rho]             fluid density

Subscripts

e, w, n, s        east, west, north and south face of the control
                  volume
p                 at point P.