Discontinuous and Mixed Finite Elements for Two-Phase Incompressible Flow

Abstract

Summary The simulation of multiphase flow presents several difficulties, includingthe occurrence of sharp moving fronts when convection is dominating,the need for a good approximation of velocities to calculate the convective terms of the equation, andflow singularities around wells. To handle the first difficulty, we propose a Godunov-type higher-order scheme based on a piecewise linear approximation of the saturation associated with a multidimensional slope limiter. With respect to the second, the pressure equation is approximated by means of a mixed-hybrid formulation equivalent to the classic mixed formulation but yielding a positive-definite linear system. To solve the third difficulty, we introduce macroelements around wells. Numerical experiments illustrate the capabilities of the method. Introduction In a previous paper,1 we proposed a finite-element method (FEM) for 2D diphasic simulations. In addition to the global pressure formulation,1,2 the basic tools for this method were a higher-order Godunov-type3 discontinuous FEM to approximate the saturation equation with little numerical diffusion and a mixed finite-element approximation of the pressure equation to reduce grid-orientation effects through better coupling of the pressure and saturation equa-tions. In this paper, we present new developments added to this method concerning the approximation of the saturation equation, the resolution of the pressure equation, and the representation of wells. The discontinuous FEM consists of a discontinuous, piecewise, linear or bilinear approximation of the saturation that allows us to design a higher-order scheme. Inside each cell, the water-conservation equation is multiplied by test functions and integrated by parts. Numerical fluxes approximating water flow rates on the edges are calculated at Gauss integration points on the edges by use of Godunov equations solving a 1D Riemann problem and involving the two traces of the saturation at the Gauss points. Extending ideas developed by Van Leer,4 we added a multidimensional slope limiter that greatly enhances the scheme's stability and that prevents overshoots and oscillations2 to the discontinuous FEM (see Ref. 5 for a similar approach). For the pressure equation, the usual mixed finite-element formulation2,6,7 yields a linear system that is difficult to solve because it is not positive-definite. To overcome this drawback, we use the equivalent mixed-hybrid formulation.2,8 In addition to the pressure inside the cells and the total flow rates across the edges, this formulation introduces as new unknowns the pressure values on the edges. The original degrees of freedom can be eliminated to obtain a linear system in the new unknowns that is symmetric, positive-definite, and easy to solve. Finally, well models called macroelements9 are presented. A union of cells around a well is divided into several sectors inside which the flow is assumed to be radial and is calculated through a simplified 1D simulator in radial coordinates. For each well, a finite number of 1D radial problems (the macroelement) links the finite-element variables to the well variables, allowing for direct calculation of the bottomhole pressure (BHP), for example, and for precise calculation of water breakthrough. With such a method, one can simulate in particular nonradial flow around producing wells. All these features are incorporated in the BIDIMIX simulator for tw0-phase incompressible flow developed jointly by INRIA, Inst. Français du Pétrole (IFP), and Elf Aquitaine. Performance of the simulator is demonstrated on four different numerical tests: a comparison with a first-order finite-difference method on a simulation of a quarter of a five-spot pattern, an edgedrive coning problem, an imbibition problem, and a field-scale problem. Saturation Equation Ref. 1 describes the saturation equation in detail. The water-conservation law isEquation 1 where the water volumetric flow vector, uw, can be written asEquation 2 Here, the usual terms have been rearranged in such a way that each physical effect is clearly identified: r=-?(x)Pcmax(x)?[a(Swr)] is the capillary diffusion term; b0(Swr) ut, is the standard transport term, with the fractional flow, b0, and the total water/oil flow vector, ut; b1(Swr)q1 describes the differential action of gravity on the two fluids, with q1 proportional to ?D(x); and b2(Swr)q2 describes the differential action of capillary pressure heterogeneity on the two fluids with q2 proportional to ?Pcmax(x). The last three terms will be approximated in the same way; therefore, it is convenient to introduce the transport termEquation 3 Typical boundary conditions for these equations areEquation 4 on the outer boundary, ? l, of the field,Equation 5 (given water injection rate) on the boundary, ?e, of an injection well; andEquation 6 (capillary effects neglected) on the boundary, ?s, of a production well. The initial condition isEquation 7 Discontinuous Finite-Element Approximation of Saturation. The reduced saturation, S wr, is approximated by Swra?M1, which is linear on every Triangle (bilinear on every Parallelogram) K of the finite-element mesh covering the field domain (see Fig. 1). By use of forward differencing in time, the saturation is obtained at each timestep through a two-step calculation. The first step is a finite-element calculation, and the second step limits the slope of the saturation calculated in the first step. Discontinuous Finite-Element Approximation of Saturation. The reduced saturation, S wr, is approximated by Swra?M1, which is linear on every Triangle (bilinear on every Parallelogram) K of the finite-element mesh covering the field domain (see Fig. 1). By use of forward differencing in time, the saturation is obtained at each timestep through a two-step calculation. The first step is a finite-element calculation, and the second step limits the slope of the saturation calculated in the first step.

A new treatment of capillarity to improve the stability of IMPES two-phase flow formulation

“…The classical fractional flow formulation used in [9][10][11][12]26] becomes inconsistent. In this paper, we follow the twophase flow formulation [32,38,48].…”

“…The methods that need the gradient of saturation in spacial dimension can not be well used to the case of different capillary pressure functions for multiple rocktypes, because of the discontinuity of saturation across rock interface. The classical fractional flow formulation used in [9][10][11][12]26] may not be suitable for highly heterogeneous media because of its inconsistency [32]. A two-phase flow formulation has been proposed in [32,33] that is based on the conception that the wettingphase pressure is always continuous as long as none of the phases is immobile, and can correctly describe discontinuities in saturation due to continuity of capillary pressure and discontinuities of different capillary pressure functions on the interface of regions.…”

Numerical simulation and homogenization of two-phase flow in heterogeneous porous media

“…A numerical code [4] has been written to solve problem (20) and compute the homogenized tensor obtained in Section 4.1, for any heterogeneous distributions in a given cell. The results of the approximate solution of Equation (20) are used below in several test problems. 1.…”

Compositional Modeling by the Combined Discontinuous Galerkin and Mixed Methods

“…Chavent et al (1990) used these methods to solve two-phase incompressible, immiscible fluid flow problems. The combined methods were used to solve convectiondiffusion equations by Siegel et al (1997).…”