Vertex Approximate Gradient discretization preserving positivity for two-phase Darcy flows in heterogeneous porous media

Abstract

In this article, a new nodal discretization is proposed for two-phase Darcy flows in heterogeneous porous media. The scheme combines the Vertex Approximate Gradient (VAG) scheme for the approximation of the gradient fluxes with an Hybrid Upwind (HU) approximation of the mobility terms in the saturation equation. The discretization in space incorporates naturally nodal interface degrees of freedom (d.o.f.) allowing to capture the transmission conditions at the interface between different rock types for general capillary pressure curves. It is shown to guarantee the physical bounds for the saturation unknowns as well as a nonnegative lower bound on the capillary energy flux term. Numerical experiments on several test cases exhibit that the scheme is more robust compared with previous approaches allowing the simulation of 3D large Discrete Fracture Matrix (DFM) models. 1 discretizations like the Mixed Hybrid Finite Element (MHFE) method [34,2] or the Hybrid Finite Volume scheme [1]. The additional asset of nodal discretizations compared with cell centered discretizations is to include d.o.f. at rock type interfaces as long as the mesh is conforming with the heterogeneities of the porous medium. These nodal pressures and saturations unknowns will be used to enforce the Darcy normal flux continuity equations combined with the saturation jump condition [16,11,27]. In this work, we will consider the Vertex Approximate Gradient (VAG) discretization developed in [27,12,13] for two-phase Darcy flows in heterogeneous media. The VAG scheme is based on nodal d.o.f. like Control Volume Finite Element (CVFE) methods but it also includes the cell d.o.f. which are eliminated at the linear algebra level at each Newton iteration without any fill-in. These cell d.o.f. provide an additional flexibility in the choice of the control volumes in order to avoid mixing different rock types at nodal control volumes. It has been shown in [27] to be more accurate than usual CVFE discretizations.The time integration will be chosen implicit to avoid severe time step restrictions in high velocity regions such as fractures. It will be also fully coupled to account for the strong coupling between the pressure and saturation unknowns in the transmission conditions at different rock type interfaces. As noticed in [5], the pressure and saturation unknowns cannot be decoupled at such different rock type interfaces to preserve the stability of the discretization.The selected numerical method should also provide a robust nonlinear convergence of the Newton type solver to allow for large time steps, typically at the time scale of the matrix in DFM simulations. Let us refer to [41,37] for alternative linearization schemes to the Newton method which aim to the improvement of the nonlinear solver. In this work, the additional nonlinear solver robustness will be first achieved by extension of the Hybrid Upwind (HU) transport scheme to the VAG discretization framework. The HU transport scheme has been introduced in [24,28] as an alternative to the Phase Pot...

Nodal Discretization of Two-Phase Discrete Fracture Matrix Models

“…The first is to go back to face based discretizations allowing the elimination of the mf interface unknowns with a local nonlinear interface solver as in [1] using TPFA discretization on orthogonal meshes and in [2] using an HFV discretization. The second perspective is to use the more robust Hybrid Upwinding approximation of the mobilities to define the two-phase Darcy fluxes at mf interfaces as proposed in [6] for TPFA schemes and in [19] for the VAG discretization.…”

“…The main advantage of this framework, which applies to an arbitrary number of rock types, is to incorporate in the construction of the functions (19) the saturation jump condition at different rock type interfaces and to apply to general capillary pressure functions. In practice, we use τ = s nw for χ = {m} and the parametrization defined in [14] for χ = {m, f }.…”

“…As in Section 3.5, the definition of the primary and secondary unknowns at the d.o.f. located at the rock type interfaces is based on the parametrization of the capillary pressure graphs (19). To fix ideas, we assume the presence of 3 rock types RT = {m, f d , f b } where f d is a fracture drain rock type and f b is a fracture barrier rock type while m denote again the matrix rock type Let the fracture network Γ be partitioned into the networks Γ d of fractures acting as drains and the network Γ b of fractures acting as barriers.…”

Upstream mobility finite volumes for the Richards equation in heterogenous domains

“…Moreover, one can still make use of the parametrized cut-Newton method proposed in [5] to compute the solution to the nonlinear system corresponding to the scheme. This method appears to be very efficient, while it avoids the possibly difficult construction of compatible parametrizations at the interfaces as in [9][10][11]. An involved comparative study on the robustness of the Newton solver is presented in [6], where other strategies to capture the discontinuities related to rock changes are also addressed.…”

A hybrid-dimensional compositional two-phase flow model in fractured porous media with phase transitions and Fickian diffusion

“…In the case of two-phase Darcy flow, a natural extension of continuous pressure models is considered in [17,18,19,20,21,22,23,24] for fractures acting as drains. It is obtained by assuming that both phase pressures are continuous at mf interfaces.…”