A New 3D Front-Tracking Approach to Modeling Displacement Processes in Complex Large-Scale Reservoirs

Abstract

Abstract A new 3D front-tracking model for processes with sharp fronts has been developed. This model has advantages over traditional finite-difference and finite-element methods, in that it distinctly recognizes fronts as discontinuities, in the process reducing both computational effort and numerical dispersion. The model treats the reservoir as having a number of displacement regions divided by sharp fronts. Fluid properties within regions are treated as uniform, but heterogeneous rock properties can be handled. Fronts are treated as piecewise planes. The pressure field is solved as a succession of steady- state calculations using standard finite-element methods. Fronts are moved according to explicitly calculated pressures. This model is well-suited to dealing with problems in which reservoir architecture is the primary control on process performance. Process physics, even those involved in complex processes such as in-situ combustion, can be handled by the model with sufficient accuracy to make it useful for field-scale studies, and/or those which require representation of small-scale heterogeneities, or those where a very large number of runs is required (e.g., for stochastic simulation). Introduction Most oil and gas recovery processes involve multiphase transport in porous media. Often these transport processes are convection-dominated and tend to develop sharp fronts. This is especially true when effective mobility ratios are low, such as water displacing gas or in thermal recovery processes. Reservoir simulation requires accurate and efficient models to simulate highly complex flow processes at large scales. Conventional finite-difference reservoir simulators have achieved great success in this regard, having become the dominant tools for reservoir simulation. However, they possess two well-known deficiencies in solving strongly convective processes:–Numerical dispersion and oscillations. For convection-dominated processes that develop sharp fronts, solutions around sharp fronts tend to be smeared or to oscillate. The effect can be reduced by grid refinement and other numerical techniques. - Grid-orientation effects. Solutions depend on the orientation and size of the grids. The effect can be reduced, but cannot be eliminated, by grid refinement. Grid-orientation effects are a result of coupling between the anisotropic numerical diffusion and the physical instability of the displacement front. To attain high level of accuracy for convection-dominated flow processes, conventional simulators require large numbers of nodes and significant computation times. This can make field-scale simulation very expensive. A number of techniques intended to overcome these deficiencies have been developed, such as higher-order difference schemes, adaptive local grid refinement, and moving finite-element methods. Front tracking is another technique that eliminates numerical dispersion and significantly reduces grid-orientation effects without appreciably increasing computational costs. There are a number of methods that can be called front-tracking methods. Among them are streamline and streamtube methods, method of characteristics, and shock-fitting methods. We reserve the term front tracking exclusively for the shock-fitting methods described below. In front-tracking methods, fluid interfaces (sharp fronts) are represented by mathematical discontinuities that separate the flow domain into a number of regions. A pressure equation (Darcy's law) is applied to each region. Different regions are coupled through front conditions. Fluid interfaces propagate with velocities that are obtained as part of the solution, or are derived from the solution. McBryan et al. developed a 2D front-tracking scheme using rectangular grids. Fronts are represented by piecewise linear element boundaries. Grids are locally adapted to match the fronts. Glimm et al. applied the same scheme to 2D reservoir simulation. The saturation equation is solved through the solution of interior regions and the propagation of interfaces. Interior regions are solved by a standard upwind scheme. P. 323^