Numerical Calculation of Immiscible Displacement By a Moving Reference Point Method
Authors: Rachford H.H. 
Published: December 2007

Abstract
Numerical solutions of immiscible flow problems in which dispersive effects of capillarity are dominated by convection require excessively fine grid spacing with attendant high computing costs. The use of coarser spacing reduces cost but often produces oscillation or undue dispersion associated with displacement fronts.
A numerical formulation is proposed here which should be applicable to two - dimensional flow problems. It is in part analogous to an approach previously tested for miscible systems. The convective transport is approximated using a change of variables to yield a coordinate system moving approximately with the local characteristic velocity. The capillarity- induced dispersive terms in the differential system describing the process are approximated with respect to a fixed coordinate system by the usual implicit formulation.
One-dimensional tests of the procedure yielded results in which the saturation profiles tended smoothly to the zero-capillary pressure solution as the ratio of viscous to capillary forces was successively increased in a sequence of calculations. This contrasted favorably with solutions by other numerical procedures which would require attendant grid refinements to approach the zero capillary pressure results.

INTRODUCTION
Numerical solution of displacement problems has until recently relied on applying methods developed primarily for transient heat-flow problems. Such problems are classified as parabolic in type, and where the heat transport is purely by diffusion their solutions are characterized by a high degree of smoothness. It is not surprising, therefore, that for approximating these solutions available finite difference methods are quite adequate. In flow problems the transport is partly by diffusion, partly by convection or flow.
Although the problem remains of parabolic type because the dispersive effects of capillary forces or diffusion play some role in every displacement, at high flow rates the problem is dominated by convection, and solutions tend toward those of equations of the hyperbolic type. Solutions of hyperbolic problems are characterized by the translation of fronts, or discontinuities, that may progressively increase in sharpness. Numerical methods for treating parabolic problems become less and less satisfactory as displacement rates increase and the role of dispersion due to concentration or capillary pressure gradients becomes small relative to transport due to flow.
In computation the difficulty manifests itself as an error associated with the grid size chosen.1-6 In summary, if the heat-flow type approximations are to include the terms arising due to convection, one of several choices may be made:an upstream (to the direction of flow) approximation for the convection terms may be used;a centered-in-distance (CID) approximation may be used; ora recently developed approximation based on the theory of oscillation matrices may be chosen.6
The last appears to have significant promise for one-dimensional flow problems; its extendibility to two or three dimensions is an open question. In either of the first two approaches, a suitably small ratio of v?x/D must be maintained, where v is the velocity, ?x is the grid spacing and D the effective dispersivity in the direction of flow.
In the first choice, the approximation of the convective part is only first-order correct and errors introduced appear as a numerically induced dispersivity of magnitude proportional to v?x. In the CID choice, the approximation can be second - order correct, but the difference formulation fails to satisfy the maximum principle unless a condition on v?x/D is met. Practically, this means that for high flow rates oscillatory solutions may result in the neighborhood of a front unless exceedingly small grid intervals are taken. While the procedure proposed by Stone and Brian4 permits a less severe limitation to be placed on this ratio, ultimately the flow rates increase relative to the dispersivity the oscillation obtains.