A third-order differencing scheme is described and compared with several alternative methods commonly used to discretize the convection term of the component conservation equations of compositional simulators. This new method can be used at higher cell Peclet numbers than those used with other common convective differencing schemes, resulting in a substantial saving in computation time for the same accuracy. Comparisons are made with analytical solutions and with 2D and 3D simulation results for single-point upstream weighting, two-point upstream weighting, and Chaudhari's method.
Because many EOR methods are slug processes, dispersive mixing and crossflow significantly affect performance. Proper numerical modeling and prediction of these effects are essential for successful design and evaluation of both laboratory corefloods and field-scale projects. In chemical flooding, a slug process where a small amount of oil-mobilizing chemical is injected into the formation, mixing and crossflow phenomena are especially important to performance and need to be modeled accurately. This mixing is characterized by dispersion in both the longitudinal and transverse directions. Interaction of dispersion and phase behavior may cause the chemical slug to lose its effectiveness through several mechanisms (e.g., phase trapping or changes in phase environment). In other cases, dispersion may enhance oil recovery efficiency by increasing the sweep efficiency (e.g., unfavorable mobility ratio and/or flow in layered reservoirs). Dispersion may also lead to early breakdown of the surfactant slug, depending on the relative position of the different banks that may develop during micellar/ polymer flooding.
In 3D simulation of tracer studies conducted before and during many EOR projects for reservoir characterization, one main objective is to determine the levels of physical dispersion and heterogeneity, which in practice are coupled on a reservoir scale. In numerical reservoir simulation, however, artificial numerical dispersion can further smear concentration fronts by increasing the level of dispersion, resulting in inaccurate predictions of recoveries and breakthrough times. In 1D simulation of corefloods where recovery mechanisms may be under investigation, the interaction of artificial numerical dispersion and phase behavior may result in erroneous conclusions. Therefore, accurate simulation of physical processes involved in EOR requires that numerical dispersion be essentially eliminated.
One would like to simulate physical dispersion accurately in a controllable manner through information entered into the simulator with no additional numerical dispersion, practical grid sizes, and results that are independent of grid orientation. The artificial numerical dispersion encountered in compositional reservoir simulation results primarily from the truncation error of the finitedifference operator used to discretize the first-order spatial derivative of the convection term. On the other hand, for convectiondominated flow, at high Peclet numbers, the hyperbolic nature of the equations can produce oscillatory solutions.
It is well known that lower-order methods, such as single-point up weighing, suffer from excessive numerical dispersion and grid orientation effects. For these reasons, other methods have been proposed, such as the method of characteristics and flux-correction schemes. Because of either the complexity of the solution algorithm or the computer requirements (both storage and run time), however, such methods have not been popular in compositional reservoir simulation of EOR processes.
A simpler second-order explicit method, proposed by Chaudhari and implemented in the U. of Texas chemical flooding simulator (UTCHEM), produces accurate results for chemical flooding applications when the cell Peclet number is kept below two; this is, unfortunately, a limiting factor in large-scale or convection-dominated applications. At low cell Peclet numbers, the two-point upstream-weighting method produces more accurate results than single-point upstream weighting; however, at higher cell Peclet numbers (convection-dominated flows), it also produces inaccurate results.
A third order convective differencing scheme based on Leonard's method is presented that is especially suited for convection-dominated flow problems and that results in a minimal increase m storage or computation time per timestep per gridblock compared with lower-order methods. This method, the "higher-order method," has been generalized to variable velocities in three dimensions with complex phase behavior and even allows for a variable number of phases locally. Comparison of the results of four differencing schemes (single-point upstream weighting, two-point ups weighting, Chaudhari's numerical dispersion control, and the higher-order method) used with UTCHEM for tracer flow, waterflooding, polymer flooding, and micellar/polymer flooding is presented.
The following differencing equations are suggested by Leonard: for u >0,
deltaC Ci+1-Ci-l Ci+l-3Ci+3Ci-l-Ci-2 u ------ =u ------------------------------- deltaX i 2delta Xi 6delta Xi
1 delta 4C ---------------delta Xi3+H).........................(1a) 12 delta X4
and for u less than 0,
delta C Ci+1-Ci-1 Ci+2-3Ci+1+3Ci-Ci-1 U------ =u --------------------------------- delta X i 2delta Xi 6 delta Xi
1 delta 4C +------------- delta Xi3+H)...........................(1b) 12 delta X4
These differencing equations can be thought of as being equivalent to central differencing with an upstream-weighted correction factor. They can be derived by Taylor series expansion for the first-order spatial derivative of concentration. This method is also equivalent to a weighted average of the two-point upstream-weighting method and central-differencing method with weighting coefficients of 1/3 and 2/3, respectively. The above equations are for constant velocity and constant grid sizes.
In compositional simulation, the differencing operator for the first-order spatial derivative of the convection term can generally be defined as
delta uC Ui+ 1/2 Ci+ 1/2 -Ui-1/2 Ci-1/2 -------- = ---------------------------------..............(2) delta X i delta Xi
The component and phase radices are not shown.