Implementation of a Flux-Continuous Finite-Difference Method for Stratigraphic, Hexahedron Grids

Abstract

Summary In this paper we present a 3D flux-continuous finite-difference formulation designed for flow simulation of models with nonorthogonal hexahedron grids with general tensor permeability. Our development follows that of Aavatsmark et al.,1 but we do not operate in transformed space. The new 27-point discretization formula has been implemented in a finite-difference reservoir simulator. This stencil has many desirable properties, including collapsing into a consistent form in two dimensions. We demonstrate that there are many practical situations in which neglecting the influence of nonorthogonality and general tensors results in first-order errors in flow predictions. A rigorous implementation for this 27-point difference operator as a control-volume finite-difference method determines the upwinding of convection terms associated with multiphase computations. Results and issues associated with implementation of this operator in a conventional finite-difference reservoir simulator are discussed. As an alternative to directly solving the linear matrix associated with the 27-point stencil of the flux-continuous operator, we examine iterative methods that split the matrix into a 7-point stencil part and a remainder. The 7-point stencil part is solved by a direct or iterative method, with the remainder part updated from the previous timestep or iteration. This split operator may permit retention of the linear solver for the standard 7-point formulation while retaining nonorthogonal grid and tensor information. Introduction Reservoir characterization models continue to grow in geologic and geometric complexity. Highly detailed descriptions of reservoir properties and complex stratigraphic and structural features are being included into real-field applications. This continued expansion in the complexity of reservoir models is stretching the capabilities of conventional, or brick-like, grids. Nonorthogonal grids, on the other hand, can describe geometrically complex surfaces and volumes accurately. While tetrahedron grids are the most flexible grids in three dimensions, their application in reservoir flow simulation is still limited. A number of important prerequisites are necessary in order to benefit from the flexibility of general unstructured grids. These include rigorous discretization schemes, improved solution algorithms, and more powerful computational and visualization resources. Stratigraphic grids (S-grids) offer an alternative, and they are common in flow simulation practice. Based on hexahedra, S-grids can accommodate complex geologic and geometric features while retaining many of the desirable features of structured grids. For example, S-grids aligned with layer-boundary surfaces are generally structured, but nonorthogonal. To translate the benefits of S-grids to more accurate flow predictions in geometrically complex reservoirs, a rigorous treatment is needed. Discretization of the flow equations for domains with nonorthogonal S-grids yields full tensors. Full tensors can also result from considerations related to the scaleup of detailed description of reservoir properties (e.g., permeability). The scaled-up permeability is a full tensor even if the fine-scale description is isotropic2 and the grid is orthogonal. Furthermore, if the simulation grid is not aligned with the principal directions of permeability, the permeability of the simulation grid becomes a full tensor with nonzero off-diagonal terms.3 For some systems, the cross terms of the permeability tensor can be quite important. These include formations containing complex crossbedding, dipping layers not aligned with the coordinate system,2 or extensive fracturing. 4,5 In these cases the simulation model will contain full tensor permeabilities with significant off-diagonal terms, which must be accommodated by the simulator if reservoir flow performance is to be predicted accurately. Edwards and Rogers6 and Aavatsmark et al.7 independently proposed a flux-continuous finite-difference formulation for two-dimensional (2D) flow. Crumpton et al.8 developed a similar formula based on flux continuity at cell boundaries. Lee et al.3 derived an explicit form of the flux-continuous finite-difference method for 2D Cartesian grids with tensor permeability. They also investigated real-field flow problems. Aavatsmark et al.1 recently extended a similar development to three-dimensional (3D) flow. In their development, Aavatsmark et al. translated the problem from real-space to orthogonal Cartesian grid using a generalized coordinate transform. The computations of the discretization scheme are then carried out in transformed space. Following the development of Aavatsmark et al.,1 we derive a 3D, flux-continuous, finite-difference operator designed to handle nonorthogonal grids with heterogeneous full-tensor permeability. Unlike Aavatsmark et al.,1 however, we do not use a coordinate transform. For general 3D hexahedron grids with heterogeneous permeability, the coordinate transformation does not necessarily yield a simpler derivation of the flux-continuous finite-difference formula. Our scheme produces a 27-point stencil for 3D hexahedron grids. The flux continuity conditions among the 27 neighboring cells entails solving a 24×24 matrix. Because the transmissibilities are computed only once in the initialization step of reservoir simulation, the relatively complex transmissibility computation of the flux-continuous method should have a minimal impact on the overall numerical performance.

Accurate local upscaling with variable compact multipoint transmissibility calculations

“…For example, nonorthogonal grids lead to full-tensor effects in the pressure equation that governs flow even when the underlying permeability field is isotropic. Many authors have investigated MPFA approaches for solving single-scale flow problems with full-tensor permeability that is assumed constant over grid cells [3,17,31]. Naturally, MPFA methods lead to larger stencils than TPFA.…”

Global variable compact multipoint methods for accurate upscaling with full-tensor effects

“…Full-tensor anisotropy can also appear from grid nonorthogonality effects. It is generally accepted that multipoint flux approximations (MPFA) are required to accurately represent full-tensor effects in finitevolume flow simulators [2,10,23]. These methods express the flux between two adjacent grid blocks not only in terms of the pressure in those grid blocks, as in twopoint flux approximations (TPFA), but also in terms of pressures in a number of other grid blocks near the face.…”

Nonlinear two-point flux approximation for modeling full-tensor effects in subsurface flow simulations

“…Such techniques, investigated extensively by a number of researchers (e.g., [2,3,10,17]), represent the flux from block i to block j in terms of not only the pressures in blocks i and j, but also the pressures in a number of other nearby blocks. Specifically, for a structured, logically rectangular grid in two dimensions, the MPFA o-method (see [2]) leads to a nine-point stencil for the discrete pressure equation; in three dimensions, a 27-point stencil results.…”