One Approach to the Grid Orientation Problem in Reservoir Simulation

Consider a linear, incompressible, two-phase flow in the absence of capillary pressure. Assume that a saturation shock front exists as depicted schematically in Fig. 1, and that the mobility behind the front is extremely high. During any time step of a numerical simulation, the exact position of the shock front generally is not known, so we assume that in an average sense the position of the shock may be placed midway between the nodes. In any numerical reservoir simulator that contains an upstream formulation, the value of Ss on Fig. 1 would be representative of the saturation at which the interblock phase mobility is evaluated. This would result in the calculation of the pressure profile shown by the solid Line P. There is clearly some inconsistency, as the true pressure profile is more like the dashed Line P'. Therefore one possible source of the grid orientation effect is that the upstream finite-difference formulation produces far too small a pressure drop across a shock front, and that the directionally dependent pressure field truncation error is a significant fraction of this pressure drop.Equating fluid velocities on either side of the shock front gives (1) where lambda T represents the total fluid mobility. If the pressure at the shock front p s is eliminated from these equations, we obtain (2) whereby it can be seen that the total effective mobility is given by (3) The harmonic mobility given by Eq. 3 is closer to twice the smaller value of (lambda Tu, lambda Td) than to the upstream value, so it has the capability of producing larger pressure drops across a shock front.On Fig. 1, it can be seen that the fluid fraction flowing at the shock front is more similar to the upstream value. This suggests that for the individual phase mobilities an upstream formulation should be phase mobilities an upstream formulation should be used: (4) with a similar upstream expression for the oil phase.Inspection of Eq. 4 reveals that at unit mobility ratio it reduces to the standard single-point formulation. Also, in a one-dimensional incompressible system with fixed injection rate, upsilon T in Eq. 2 is constant. Under these circumstances, the saturation profiles calculated using the harmonic mobilities are profiles calculated using the harmonic mobilities are identical to the standard single-point upstream, although the pressure profiles are different. This similarity to single-point upstream under many circumstances and the fact that the results of Todd et al. already show marked grid orientation effect at unit mobility ratio for single-point upstream suggest that the single-point harmonic formulation should be extended to two-point harmonic. The simplest way to achieve this is to replace the value of lambda wu in Eq. 4 by the two-point value. The overshoot constraints we have used are that the extrapolated value must not be less than min (lambda wu, lambda wd) nor greater than max (lambda wu, lambda wd). We shall refer to this scheme as "two-point harmonic." The net result is very simple. The normal upstream mobilities lambda in any numerical reservoir simulator are to be modified by the term 2 lambda Td/lambda Tu + lambda Td). P. 160