ogst
Accueil Tous les numéros Volume 73 (2018) Oil & Gas Science and Technology - Rev. IFP Energies nouvelles, 73 (2018) 78 Full HTML
Numerical methods and HPC
Open Access
Numéro
Oil & Gas Science and Technology - Rev. IFP Energies nouvelles
Volume 73, 2018
Numerical methods and HPC
Numéro d'article 78
Nombre de pages 18
DOI https://doi.org/10.2516/ogst/2018067
Publié en ligne 12 décembre 2018

© C. Cancès, published by IFP Energies nouvelles, 2018

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

1.1 Two-phase porous media flows

Incompressible two-phase porous media flows are often modeled by the following set of equations. In the absence of source terms, the volume of the phase α ∈ {n, w} (n and w stand for non-wetting and wetting respectively) is locally conserved along time as a consequence of(1)where ϕ denotes the porosity of the rock (supposed to be constant w.r.t. time), where s α denotes the saturation of the phase α and v α denotes the filtration speed of the phase α. It is classically assumed that the phase filtration speeds obey the generalized Darcy law(2)

The intrinsic permeability tensor of the porous medium is a definite positive and symmetric tensor field whereas the relative permeability k r,α of the phase α is an increasing function of the saturation satisfying ,k r,α (0) = 0 (we neglect the residual saturations). The variations of the viscosities μ α and of the densities ρ α are neglected, and g = ∇(g· x ) denotes the gravity vector. The phase pressures p  = (p n, p w) are the main unknowns of the system together with the saturations s  = (s n, s w). Two algebraic relations are imposed to close the problem. The first one is a constraint coming from the fact that the whole pore volume is saturated by the fluid:(3)

The second one links the capillary pressure to the non-wetting phase saturation in a monotone way:(4)where π is a maximal monotone graph from [0,1] to that may also depend on the space variable in the case where the geological environment is made of several different rocks (see for instance [1–5]).

Once complemented by no-flux conditions on the boundary of the porous domain Ω, the model (1)–(4) has a very particular variational structure. As depicted in [6] (see also [7–9]), this problem can be reinterpreted as the generalized gradient flow [10] of the energy(5)with

This energy is made of the capillary energy, of the gravitational potential energy, and of a contribution related to the constraint (3):

The capillary potential is a convex function that is finite in [0,1], infinite outside of [0,1], and satisfies(6)whereis the subdifferential of Π at s. The functional is convex and its subdifferential ∂E( s ) is made of the couples h = (h n, h w) = (p n − ρ n g⋅ x , p w − ρ w g⋅ x ) such that the relation (4) holds.

We will not go deep into details on the description of the gradient flow structure since it is not the purpose of this paper. We refer to [7, 11] for a more complete discussion on the Wasserstein gradient flow interpretation of porous media flows and to [7, 12] for a general presentation on gradient flows in metric spaces. Let us only stress important points that will motivate our discussion on numerical schemes. First, the gradient flow structure implies that ℰ( s ) is a decreasing function of time and that the dynamics aims at maximizing this decay. This yields an energy/dissipation of the form(7)for all t ≥ 0. The energy/dissipation relation (7) can be obtained by multiplying formally the equation (1) by h α and by summing over α ∈ {n,w}. As a consequence, stable steady states s ∞ are local minimizers of ℰ for which each phase is hydrostatic on its support, i.e.,(8)

In order to compute in a accurate way the long-time behavior of the system (this is very important in the context of basin modeling), the numerical scheme has to be designed to make the discrete counterpart of the energy ℰ( s ) decay along time and, as much as possible, to be exact on the equilibrium (8) (see [13]). In the hydrostatic zones, the equilibrium consists in a balance between the diffusion and the convection. Thus we will avoid operator splitting and discretize the convection and the diffusion simultaneously to recover this equilibrium.

1.2 A simplified model problem

Our purpose can already be illustrated on the single scalar equation(9)

Π being a convex and coercive internal energy functional as previously and Ψ being a smooth external (possibly gravitational) potential. The mobility function η is nondecreasing on and satisfies η(0) = 0 and η(s) > 0 if s > 0. Here, we used the notation for the domain of Π and we assume that 0 ∈ D(Π). If Π is defined on the whole , i.e., , then we assume moreover that Π is superlinear:

The solutions to (9) corresponding to nonnegative initial data remain nonnegative along time. This property might be destroyed by the numerical approximation, so that we need to extend the definitions of the functions Π, π, and η for negative values of s when this is possible. In the case where ∂Π(0) contains a finite value, i.e., if(10)then Π can be artificially extended on (−∞,0) ∪ D(Π) into a convex function (still denoted by Π) with a single valued subdifferential at , i.e., ∂Π(0) = π°(0) = π(0). This can be done for instance by prescribing(11)henceif s < 0. The function η can for instance be extended by setting η(s) = −s if s < 0.

The framework of (9) is already interesting since it includes Richards equation for which s is the water saturation and Ψ = −ρ g· x . The capillary energy function Π is then defined on as the antiderivative of the capillary pressure function. In usual models (see [14], pp. 343–345), D(Π) = [0,1] and condition (10) does not hold. This also includes a linear Fokker-Planck equation (but written in a nonlinear form) when η(s) = s and π(s) = log(s) or models for crowd motions [15]. We are interested in the computation of non-negative solutions corresponding to initial data with(12)

There is still a generalized gradient flow structure corresponding to (9), the energy being given by(13)

The counterpart to (7) is obtained by multiplying formally (9) by p + Ψ and writes(14)

Similarly to what was discussed in Section 1.1, one aims to build numerical schemes for (9) that are accurate in the long-time limit. The steady states for (9) are now given by(15)

One of the difficulties arising in the resolution of the problem (9) is that the variable p and s are related by the maximal monotone graph π that may have vertical or horizontal parts. Similar difficulties occur for instance in the context of reactive flows in porous media, cf. [16]. For Richards equation, the graph π has vertical parts as depicted on Figure 1. The solution (s, p) to the problem (9) can be deduced from the knowledge of p, but not from s. But the choice of p as a primary unknown in a naive numerical method can yield severe difficulties. Typically mass conservation can be lost, and severe troubles can be encountered in the convergence of the iterative procedure to compute the solution to the nonlinear system arising from the numerical scheme. These difficulties motivated the development of several strategies to optimize the convergence properties of the iterative procedures. An popular approach consist in making used of robust fixed point procedures with linear convergence speed rather than with Newton’s method (see for instance [17–22]). There were also important efforts carried out to fix the difficulties of Newton’s method [23–25]. Comparisons between the fixed point and the Newton’s strategies are presented for instance in [26, 27] (see also [28]). In [29], the authors combine a Picard-type fixed point strategy with Newton’s method (i.e., they perform a few fixed points iterations before running Newton’s algorithm). An alternative approach would consist in keeping both s and p as unknowns together with the additional equation p ∈ π(s) that is often rephrased as a complementary constraint and then solving the problem with a non-smooth Newton method (see for instance [30–32]). Another classical solution consists in partitioning Ω at each time t into a part Ω s (t) where s is chosen as a primary variable and a part Ω p (t) where p is chosen as a primary variable. Switching from one variable to another is then compulsory [33] in this case. This can be seen as a particular case of the approach proposed recently in [34, 35], which consists in parametrizing the graph π in a suitable way. More precisely, the graph can be described by two functions and such that(16) cf. Figure 1. We assume moreover that the parametrization is non-degenerate in the sense that is a strictly increasing function.

thumbnail Fig. 1

The maximal monotone graph π depicted on top is parametrized as with some monotone functions (red) and (blue) depicted in the bottom picture. Here, and satisfy .

Then the knowledge of the unphysical variable w allows to reconstruct both s and p, so that on can use w as a primary variable in the computations. There are an infinity of possible choices for the parametrization of π. One example is the resolvents of π and π −1, i.e., and . The idea in [34] is to take advantage of this flexibility and to choose the parametrization in order to optimize the convergence of the Newton’s method.

Before discretizing the nonlinear transient problems (1)–(4) or (9), we need to introduce some material concerning the discretization of elliptic equations.

1.3 A diffusion building block for numerical approximation

The goal of this section is to introduce a rather general framework that can fit to many numerical methods to solve the elliptic equation(17)subject to no-flux boundary conditions. There is a huge literature about the numerical resolution of the above problem. Besides classical conforming and non-conforming finite elements (see for instance [36, 37]), let us mention some approaches based on mixed finite elements [38, 39], Multi-Point Flux Approximation (MPFA) finite volumes [40–43], or Discontinuous Galerkin method [44–46]. The list of aforementioned methods and related publications is of course far from being exhaustive, and some other numerical methods will be discussed in what follows.

Let be the set of geometrical entities to which correspond the unknowns (i.e., the cells for cell centered methods, the edges for hybrid methods, or the vertices for vertex centered methods), then we are interested in numerical methods for solving (17) that reduce into a linear system of the form(18)where are such that ∑ K m K  = |Ω|. The compatibility condition on the right-hand side reduces to ∑ K m K f K  = 0. We denote by the vector unknown and by the matrix corresponding to the system (18), i.e., where b K  = m K f K . We require the transmissivity coefficients a KL (and thus the matrix ) to be symmetric(19)

In what follows, denotes the set of the couples such that the transmissivity a KL is different from 0. Thanks to the symmetry property (19), the scheme (18) is equivalent to: ,(20)

Many numerical schemes can write as (18) and (19). This is for instance the case of the classical Two-Point Flux Approximation (TPFA) finite volume scheme [47–49], where(21)as soon as the cells K and L share an edge σ KL . In (21), x K (resp. x L ) denotes the center of the cell K, whereas the coefficient m K in (18) is the Lebesgue measure of the cell K. The TPFA scheme requires strong assumptions for its consistency. The permeability tensor must be isotropic – in (21) is the harmonic mean of the permeabilities and associated to the cells K and L –, and the edge σ KL must be orthogonal to the straight line [ x K , x L ] joining the cell centers K and L.

Another classical example of scheme satisfying (18) is the classical conformal finite elements with mass lumping, which can be seen as a particular box scheme [50]. In this case, the unknowns are located at the vertices of a simplicial mesh of Ω. Denoting by ϕ K the Lagrange basis function associated to the vertex , the transmissivities a KL are defined by(22)and the lumped mass associated to the vertex K is Note that in this context, the transmissivities a KL may be negative (for instance in the case where does not satisfy the Delaunay condition if ).

The Hybrid-Mixed-Mimetic (HMM) methods [51] containing Hybrid Finite Volumes (HFV, see [52]), Mixed Finite Volumes [53], and Mimetic Finite Differences [54, 55] also enter the framework (18)–(19), as well as Discrete Duality Finite Volumes (DDFV) [56]. We refer to Droniou’s review paper [57] and to the book [58] for a more complete overview of the methods entering our framework.

In the case where the transmissivities are nonnegative(23)the scheme (18) is monotone and it fulfills the maximum principle. The property (23) is lost as soon as the mesh and the anisotropy tensor do not fulfill restrictive conditions. But for many methods, the transmissivities remain blockwise positive. This means that there exist geometrical entities and positive definite symmetric matrices , , such that the scheme (20) rewrites(24)

The vector δ M u represents the inner variations of u inside . Let us illustrate formula (24) on some classical methods. First, for the TPFA scheme, since a KL  > 0, then one can choose , as the set of the edges, , and where M is the diamond cell corresponding to the edge (K, L). The case conformal finite elements is more interesting. In this case, denotes the set of the simplices, whereas is equal to the space dimension d. The vertices of the simplex are denoted by , andwhereas the matrix is defined by(25)

A crucial property of the local diffusion matrices is that the condition number of can be bounded by a quantity depending only on the condition number of and on the regularity (in Ciarlet’s sense [36], see also [37]) of the simplex M, i.e.,(26)

The HFV (or SUSHI [52]) also enter the framework (24). In this case, is the set of the control volumes that can be quite general (non-convex cells having various number of edges). To each cell M, there is one cell unknown and edge unknowns where denotes the number of edges of the element M, and the inner variation vector δ M v is given by

The matrix is built thanks to a formula similar to (25) but with an approximate gradient that is piecewise constant on the half diamonds. Here again, the conditioning of the matrix only depends on the one of and on the regularity of the mesh. The Vertex Approximate Gradient (VAG) scheme [59, 60] has a very similar structure and also enters our framework (see [61]). The main difference with SUSHI is that the unknowns are located at the vertices of the cells instead of the edges. A last example concerns the DDFV method that, at least in 2D, also enters the framework (24) as highlighted in [62, 63]. In this context, is the set of the diamond cells and once again, the local matrices have bounded conditioning numbers if has a bounded conditioning number and if the mesh fulfills some weak regularity assumption.

Before coming back to our nonlinear parabolic problem, there is a last point that we need to point out here. The coercivity of the numerical method (20) amounts to claiming that the matrix is a symmetric definite positive matrix, the lowest eigenvalue of being bounded away from 0 uniformly w.r.t. the mesh size:where

Combining this information with (26), it was proven in [64] that(27)for some C depending only on the mesh regularity and on the anisotropy of the permeability tensor. A similar property was derived for the SUSHI scheme in [52]. The abstract framework of polytopal toolboxes of Droniou et al. [58] also allows to derive inequalities of this type.

1.4 About the paper content

The (physical) energy stability of numerical methods appears to be a secondary point for a large part of the mathematical literature concerning porous media flows with capillary diffusion. This originates probably from the fact that the mathematical analysis for such problems often relies on the use of the so-called Kirchhoff transform and global pressure (see for instance [65–68]) which may mask the considerations regarding the physical energy. The question of the convergence of numerical methods is more often addressed. For degenerate parabolic scalar equations, we refer for instance to [69–79], whereas for multiphase porous media flows, we refer to [80–85]. Here again, the list is of course far from being exhaustive.

In the remaining of this paper, we will focus on the construction of schemes for (1)–(4) or more simply for (9) that make the discrete counterpart of the energy decay with time. The goal is to design schemes that capture in an accurate way the long-time behavior of the continuous model, and in particular the equilibriums (8) and (15). Two strategies will be discussed. The first one is detailed in Section 2 and consists in working with the formulation (20) and to use an appropriate upwinding of the mobilities. The second one, which is presented in Section 3, rather uses the formulation (24) and aims at preserving a similar structure even in the presence of a degenerate mobility.

For the sake of simplicity, we do not discuss here about possible source terms or other boundary conditions. Adding these terms to the problem is possible but the energy of the system would no longer decrease in general. However, under reasonable assumptions, the problem remains energy stable. We refer to [86] for a detailed treatment of source terms and inhomogeneous boundary conditions in the context mass-lumped finite elements.

2 Upstream mobility schemes

Before addressing the two-phase flow problem (1)–(4), let us first focus on the simpler scalar problem (9).

2.1 The scheme for the simplified problem (9)

The goal is to tune the numerical method (18) to approximate the solution to (9) while preserving a good energetic behavior at the discrete level. Concerning the time discretization, we stick to the backward Euler scheme. Let be parametrizations of the graph π as in (16), then the numerical scheme writes: (28)where we have used the notation for the neighbors of K. The time step τ n is not necessarily uniform and can depend on n ≥ 1. The initial saturation is assumed to be given and once has been computed for n ≥ 1, one sets

Following [64] and[87], the mobility is chosen thank to an upwinding that takes the sign of a KL into account:(29)

Note that , so that one can write the scheme (28) under the equivalent form(30)

In the above formula and as in (20), v is an arbitrary element of .

2.2 Some properties of the scheme

As a consequence of the choice (29) for the mobility , the method preserves the non-negativity. More precisely, this means that if is a solution of the scheme (28), then so does its projection on , where is such that and (cf. Fig. 1). So without loss of generality, we can assume that(31)

This property still holds for transmissivities a KL violating the monotonicity condition (23). Note that the inequality (31) is of interest only if .

Let us now turn to the energy estimate. Choosing in the weak formulation (30) and using the convexity inequalityfor all and for all p ∈ ∂Π(s) = π(s), we obtain the following discrete counterpart of (14):(32)where(33)

If the monotonicity condition (23) holds (like for instance for TPFA finite volume schemes), the inequality (32) is enough to ensure that the second term is non-negative, thus the discrete energy is diminishing along time, i.e.,

Moreover, summing over n provides a control on the numerical dissipation (34)

This estimate together with (31) is enough to prove the existence of one solution ( w n ) to the scheme for all n ≥ 1. Thanks to the monotonicity of the scheme and the non-degeneracy of the parametrization , the solution is unique (see [34]). Finally, the convergence of the scheme can be proved if the approximation of the diffusion operator is consistent. This requires some classical conditions on the mesh and on the anisotropy tensor for the TPFA scheme [49].

The case where the condition (23) is not satisfied is more intricate. Indeed, the estimate (32) still hold but the fact that 0 is a lower bound in (34) is no longer true in general. deducing a control on the energy and on its dissipation from (32) is not straightforward. Indeed, the second term does not have an obvious sign and we are not able to claim that the dissipation is positive along time. However, it is possible to bound from below the dissipation rate. To do so, let us first remark that the definition (29) of the upstream mobilities implies that(35)whatever the non-decreasing function . In the above relations, denotes the interval with extremal points and and we have used the monotonicity of η. In the case where is bounded (this is natural in the porous media flow context since is), we can take advantage of this formulation with , of the inequality (27), on the Lipschitz regularity of the external potential Ψ to show (cf. [88], Sect. 3.2) that(36)

A similar inequality can be obtained even in the case where is not bounded under suitable assumptions on the nonlinearities, see for instance ([61], Sect. 3.1). Inequality (36) allows to show that the discrete energy grows at most linearly with time, as well as the existence of (at least) one solution to the scheme thanks to a topological degree argument. Moreover, it allows one to show the convergence of the scheme towards the unique solution to the continuous problem if the mobility function η is bounded and if π −1 is a function (i.e., when there are no hyperbolic degeneracy in the problem (9)) when the mesh size and the time step τ n tend to 0. We refer to [88] for the details concerning the convergence analysis of the scheme (28).

The upstream mobility scheme can be tested by other quantities of the form thanks to (35). For s ∈ [0, 1], we denote by the largest interval of [−∞,+∞] such that for all . Then given a non-decreasing onto function , there corresponds a convex and coercive function such that . In particular, there holds(37)for all .

We can then take benefits of the characterization (35) of the upstream mobility to get enhanced regularity estimates. For instance, choosing(38)one gets that(39)provided Ψ is regular enough and the initial entropy is finite:

Details for estimate (39) are provided in Appendix A.1. Such an estimate is the corner-stone in the study [89] where the convergence of an upstream mobility scheme was established for the Dupuit approximation of multiphase porous media flows. It was also used in [8], or in [90] where a degenerate Cahn-Hilliard system with a very similar mathematical structure to (1)–(4) (see [91, 92]) was considered.

2.3 Convergence of the scheme

Let us illustrate the convergence of the scheme (28) when the mesh size and the time step tend to 0. As a model problem, we consider the very simple convection diffusion equation(40)where Ω = (0, 1)2, x = (x 1, x 2) t , and t f  = 0.05. The permeability tensor is assumed to be diagonal of the form

The equation boils down into the very simple linear equation(41)where e 2 = (0, 1) t . We choose the initial condition and the no-flux boundary condition in accordance with the exact solution(42)with α = π 2 + 1/4. Note that s ini ( x ) = s ex ( x ,0) vanishes when x 2 = 1.

The solution is computed thanks to a Control Volume Finite Element scheme [87, 88]. This scheme requires conformal triangulations of Ω. We use successively refined Delaunay grids from the FVCA5 benchmark [93]. In the isotropic case (corresponding to κ = 1), the scheme is monotone, i.e., condition (23) is fulfilled. This is no longer true in the anisotropic case κ = 20.

We also compute the solution to our scheme with an anisotropy ratio κ = 20. The numerical solutions are compared with those computed with the following linear scheme with centered convection: (43)

Here, Ψ K  =  x K· e 2. The scheme (43) is second order accurate in space. In the isotropic case and for fine enough grids (leading to small enough Peclet numbers), the scheme turns out to be monotone. This is no longer the case in the anisotropic case κ = 20 even though the scheme remains second order accurate in space.

In order to make sure that the error corresponding to the time discretization remains small when compared to the error corresponding to the space discretization, we divide the time step by 4 when the mesh size is divided by 2. We present on Figure 2 the L 2 (Ω × (0, t f)) error corresponding to the numerical solutions produced by the schemes (28) and (43) for Delaunay triangular grids from [93].

thumbnail Fig. 2

L 2 (Ω × (0, t f)) error as a function of the mesh size for the scheme (28) in the isotropic case κ = 1 (solid red) and anisotropic case κ = 20 (solid blue) corresponding to upstream mobilities (29). Comparison with the solution to the linear scheme (43) in the isotropic case (dashed red) and anisotropic case (dashed blue).

The method (28) preserves positivity whatever the anisotropy ratio. This is not the case of the linear scheme (43) that produces undershoots in the anisotropic case κ = 20. As expected because of the upwinding procedure, the scheme (28) is first order accurate in space, i.e.,(44)where u h denotes the piecewise constant reconstruction on the dual mesh and where h denotes the mesh size. But it appears on Figure 2 that the constant C appearing in (44) strongly depends on the anisotropy ratio.

To sum up, the method (28) enjoys a very strong stability when used in the non-monotone context, but this is mainly due to the excessive numerical diffusion. Even though the method is still first order accurate w.r.t. space, the constant strongly depends on the anisotropy ratio. This makes the method inefficient in the case of large anisotropy ratios or of poor quality meshes. However, the upstream mobility finite volume scheme remains a very robust method for solving complex but isotropic problems, like for instance some degenerate Cahn-Hilliard systems [90] or multiphase porous media flows [8, 89, 94].

2.4 Long-time behavior of the scheme

It is natural to wonder if the numerical scheme is able to reproduce in an accurate way the long-time behavior (15) of the continuous problem. This question is indeed of broad interest for porous media flows, in particular in the context where the time scales are long like for instance in basin modeling or for nuclear waste repository management. Therefore, the study of the long-time behavior of Finite Volume schemes for convection-diffusion problems has been the purpose of several contributions (see, e.g., [13, 62, 95–97]).

We go back to the two-phase flow model presented in Section 1.1 that we discretize on a Delaunay mesh thanks to a classical fully implicit upstream mobility finite volume scheme (see for instance [94] or [8]). In particular, the monotonicity relation (23) is fulfilled and the discrete saturations remain bounded between 0 and 1.

At the continuous level, the relation (8) prescribing the long-time limit boils down to the following alternative:(45)for some . The parameter γ is determined by the conservation of mass

The steady state (45) can be discretized directly into(46)with γ' fixed so that(47)

In the above relation, is the discrete porosity. In particular, in the classical case where is a (single-valued) function, the following alternative holds:for all

Let be the discrete counterpart of the energy (5), i.e.,(48)then the energy is decreasing along time, i.e.,

One can show (cf. Appendix A.2) that is a minimizer of under the constraint (47). Therefore the relative energy is non-negative. If the capillary pressure π is an increasing function on (0, 1), then the relative energy vanishes if and only if . This quantity can be used to illustrate the convergence of towards as n tends to ∞. In the case depicted on Figure 3, we set Ω = (−1/2,1/2)2, , k r,α (s) = s, μ n = 10, μ w = 1, ρ n = 0.87, ρ w = 1, g = (0, −9.81) T , ϕ ≡ 1, , and

thumbnail Fig. 3

Evolution of the relative energy as a function of time. We observe the exponential convergence towards until the machine precision is reached.

We see on Figure 3 that the relative energy converges exponentially fast towards 0 until it reaches the machine precision. This shows both that the energy is effectively dissipated along time and that given by (45) and (46) is a steady solution to the scheme.

3 Schemes with local positive dissipation tensors

The upstream mobility numerical schemes presented in the previous section enjoy very nice properties but they are merely first order accurate in space and lack robustness w.r.t. the anisotropy ratio (or to the mesh regularity) as illustrated in Figure 2. This motivates the development of alternative numerical methods with the following specifications.

  1. No Kirchhoff transform: the scheme should only involve quantities with a clear physical meaning.

  2. Nonlinear stability: the scheme must fulfill a discrete counterpart of the energy/energy dissipation relation (14).

  3. Convergence: when the mesh size tends to 0 and under mild regularity assumptions, the approximate solution produced by the numerical scheme must converge towards a solution to the PDE (9).

  4. Second order accuracy: if the solution to (9) is smooth enough (say C 2), the error between the exact and the approximate solution must behave as h 2 where h stands for the mesh size.

  5. Robustness: the scheme should allow general grids and general anisotropy tensors. The accuracy should not be excessively impacted by the anisotropy ratio or the mesh regularity.

3.1 Presentation of the methodology

Assume that the building block (18) is second order accurate, then a natural scheme to fulfill the specifications (i) and (iv) of the previous list for problem (9) is(49)with a centered choice for the mobility:

Multiplying the scheme by and summing over leads once again to(50)

Already in the case developed in Section 2 where was chosen thanks to upwinding, the sign of the second term of the left-hand side was unclear. It was however possible to get a sufficient control to claim that the energy was growing at most linearly, cf. (36). This conclusion does not hold any longer in general for a centered choice of the mobilities and no control on the energy can be deduced from (50). Hence the specification (ii) is not satisfied by scheme (49).

In order to correct this, we propose a scheme based on the formalism (24), that is(51)

The above relation must hold for any . The matrix is called the local dissipation tensor. It must incorporate the local diffusion but also the mobilities . In order to ensure the dissipation property for the scheme, we want to be symmetric semi-definite positive. In [61], we proposed to choose(52)with(53)for the particular choice(54)

The fact that the energy is diminishing along time is obtained by choosing and by applying a simple convexity inequality, leading to(55)the second term being non-negative since is semi-definite positive.

Additionally, the method is globally mass conservative, i.e.,(56)where . This estimate is obtained by choosing v  = 1 in (50). Let us stress that in the particular cases of the schemes studied in [61, 62, 86], the scheme is also locally conservative since fluxes can be constructed.

In general, the quantity appearing in (52) is chosen as a convex combination of . In order to assess that the scheme (50)–(52) converges, the numerical mobilities have to satisfy an additional coercivity condition, that is(57)for some uniform α > 0. The condition (56) is clearly satisfied by the choice (53) with α = 1/2, but it prohibits the choice that would have been quite natural in the context of the SUSHI [52] or VAG [59] schemes.

The implementation of the scheme (50)–(53) can appear to be too involved. An easy way to simplify it is to choose independent on i, i.e.,

The matrix then reduces to and commutes with , leading to the simpler formula . In view of the constraint (56), a natural choice for η M ( v ) is(58)

This choice was successfully used in [86] for a method based on conformal finite elements with mass lumping (here ), and in [62, 63] for a nonlinear DDFV method (here d = 2 and ).

3.2 Conditional positivity preservation and convergence w.r.t. the grid

The scheme (50) amounts at each time step to a system of nonlinear equations of the form(59)

The functions η, and are uniformly continuous on , then is also uniformly continuous. In order to ensure the existence of a solution w n to the system, we need some bounds on w n (that might depend on the mesh and the time step).

In the case where (10) holds, the estimate (55) (in particular the dissipation term) provides a sufficient bound on w n in order to apply a topological degree argument and to claim that there exists (at least) one solution to the scheme (51). In the more intricate situation where (10) is no longer satisfied, corresponding to the situation where(60)(with a slight abuse of notation), then one can show that(61)for some ζ depending on the time step τ n and on the mesh. This is done for instance in [61] of Lemma 3.7 in the context of the VAG scheme, or in [63] of Lemma 3.5 for a DDFV scheme. In both case, the proof strongly relies on the coercivity assumption (57).

As a consequence of (61), the scheme (51) preserves the positivity as soon as (60) holds. This property is unfortunately lost in general when (10) is satisfied. To illustrate this fact, we show results of [61] where the solution of the porous medium equation(62)is rewritten under the form(63)for three different choices of nonlinearities, that are

  1. η(s) = 1 and π(s) = |s|s (recall that we need to extend the nonlinearities for negative saturations if (10) holds);

  2. η(s) = 2|s| and π(s) = s;

  3. η(s) = 2s 2 and π(s) = log(s).

The case (a) does not enter our framework since η(0) = 1 ≠ 0, but it is interesting since it corresponds to the most natural approach to solve (62). The condition (10) holds in cases (a) and (b), whereas (60) holds in case (c). Therefore, the positivity of the solutions should be guaranteed only in this last case. To illustrate this fact, let us choose Ω = {(x, y) ∈ (0, 1)2}, and let us approximate the exact solution to (62) defined by(64)thanks to the VAG scheme [61] on successively refined triangular meshes from the FVCA5 benchmark [93]. The problem is here complemented with Dirichlet boundary conditions.

We observe in Tables 1–3 that second order convergence is destroyed for all the three schemes because of the lack of regularity of the exact solution. As expected, the discrete solution corresponding to the choice (c) remains positive because condition (60) is verified. This is no longer the case for the choices (a) and (b) and undershoots are observed. But the choice (b) appears to be both cheaper and more accurate than the choice (a), whereas the amplitude of the undershoots is reduced.

Table 1

Choice (a) of mobility and pressure functions, convergence towards (64).

Table 2

Choice (b) of mobility and pressure functions, convergence towards (63).

Table 3

Choice (c) of mobility and pressure functions, convergence towards (63).

Let us illustrate again the ability of the approach. We consider the linear and isotropic convection diffusion equation

that we rewrite under the nonlinear form(65)

We aim to approximate the exact solution (42) thanks to the nonlinear DDFV method proposed in [62, 63]. We compute the approximate solution corresponding to the sequence of Kershaw meshes from [93], see Figure 4. The numerical results are presented in Table 4.

thumbnail Fig. 4

The Kershaw meshes are highly deformed topologically cartesian grids. The coarsest mesh of the family is depicted here.

Table 4

Approximation of (42) with a nonlinear DDFV scheme [62, 63] on the Kershaw mesh family from [93], final time T = 0.25. M is the mesh index, τ is the time step, errgs and errs respectively stand for the L 2(Ω × (0, T)) error on ∇ s and the L ∞((0, T); L 2(Ø)) error on s, whereas ordgs and ords are the corresponding convergence orders. N max and N mean are the maximal and mean numbers of Newton iterations, and s min minimal value of the approximation of s during the whole simulation.

As expected, since π(s) = log(s) fulfills (60), the solution remains positive although the grid is very irregular. The convergence at order 2 on u is not affected by the poor mesh regularity. We observe a super-convergence for the gradient.

Finally, let us show that the methodology presented in this section is often much more robust w.r.t. anisotropy than the methodology presented in Section 2. To this end, we stick to the test-case of the linear Fokker-Planck equation written in a nonlinear form described in Section 2.3. We discretize it with the energy stable finite element scheme with mass lumping proposed in [86]. In this scheme, the mobilities η M ( v ) are chosen according to formula (58). The final time is set to t f = 0.25. The L 2 ((0, t f) × Ω) error as a function of the mesh size is plotted on Figure 5. As expected, the method is of order 2 whatever the anisotropy ratio. But it is worth noticing that the accuracy is almost not affected by the anisotropy ratio .

thumbnail Fig. 5

Convergence towards the analytical solution (42) of the Fokker-Planck equation (41) for different values of the anisotropy ratio .

3.3 About the long-time behavior of the scheme

We aim now to illustrate the long-time behavior of the scheme. One discretizes the equation (65) complemented with no-flux boundary conditions following the methodology of [61]. In particular, we take , hence we denote by s n  =  w n the vector of the unknown saturations. The scheme (51) then rewrites(66)for any . The mobilities are discretized thanks to formula (54). The mass is conserved along time, i.e.,(67)

For any ρ > 0, the vector of defined by(68)is a steady solution to the scheme (66) since for all . Then the expected long-time limit as n → ∞ is s ∞ where ρ has been tuned so that

To illustrate this fact, we plot on Figure 6 the evolution of the relative energy as a function of the discrete time nτ for the sequence of successively refined triangular meshes used in Tables 1–3 and for the Kershaw mesh depicted on Figure 4. The relative energy is proved to be decreasing along time and vanishes if and only if s n  =  s ∞.

thumbnail Fig. 6

Plot of the log of the relative energy (in log scale) as a function of time.

We show in Figure 6 that the relative energy converges exponentially fast towards 0, providing a discrete counterpart to the relationthat is deduced from (42). We observe on Figure 6 that the scheme preserves exactly (up to machine precision) the long-time behavior of the equation, i.e., the long-time limit s ∞ is a discretization through (67) of the exact long-time behavior of (64). Figure 6 suggests that the convergence speed is sensitive to mesh regularity (the convergence is slightly too fast on the Kershaw mesh) and to mesh size (the convergence is too slow on coarse triangular meshes). Note that the relative energy is defined only for non-negative s n . For a behavior like the one depicted on Figure 6, the scheme must (at least) preserve positivity and admit s ∞ defined by (67) as a steady state.

3.4 Application to two-phase porous media flows

We show now preliminary results obtained thanks to the methodology of this section to simulate two-phase porous media flows. We consider an anisotropic quarter five spot problem with injection in a top-left corner and extraction in the bottom-right corner. We neglect gravity, so that the equations boil down into(69) (70)

In the above system, we eliminated the unknowns p n and s w. They can be deduced from p w and s n by the relations s w = 1−s n and p n = p w + π(s n). We fix the nonlinearities η n, η w, π, f n, and f w as in [98]:and

The pure wetting phase is injected in ω inj = (0, 0.05) × (0.95, 1) near the top-left corner, while the fluid occupying the area ω prod = (0.98, 1) × (0, 0.02) near the bottom-right corner is extracted (see Fig. 7). Defining the fractional flow functions bythen the source terms q n and q w are defined by

thumbnail Fig. 7

Schematic representation of the injection and production wells.

The permeability tensor is anisotropic and we take a constant saturation profile ϕ ≡ 1.

The system (69)–(70) is complemented with no-flux boundary conditions and with the initial condition s n (·,0) ≡ 0.9.

For the discretization, we use conforming finite elements with mass lumping on a structured triangulation made of 5000 triangles and a constant time step τ = 10−3. For the phase mobilities, we make a choice of type (58). As explained in [86], the method is locally conservative: one can build equilibrated fluxes in the Raviart-Thomas-Nedelec space such thatwhere denotes the approximation of the saturation s α at time .

We plot snapshots of the saturation sn on Figure 8.

thumbnail Fig. 8

Approximate saturation s n,h at time t = 0.2 (top left), t = 0.5 (top right), t = 1.5 (bottom left) and t = 2.5 (bottom right).

Acknowledgments

The author warmly thanks the anonymous referees for their precious feedback. He overall wants to thank his numerous collaborators on this topics, namely Ahmed Ait Hammou Oulhaj, Konstantin Brenner, Claire Chainais-Hillairet, Thomas O. Gallouët, Cindy Guichard, Stella Krell, Maxime Laborde, Léonard Monsaingeon, Flore Nabet, and Martin Vohralík. They actively contributed to the realisation of this work. This research was supported by the French National Research Agency through project GEOPOR (ANR-13-JS01-0007-01) and Labex CEMPI (ANR-11-LABX-0007-01).

References

A Some technical details

A.1 About Estimate (39)

We now give some details on the derivation of estimate (39). Note that the derivation of estimate (36) relies on similar arguments.

Let be defined by (38), then multiplying (28) by and summing over leads to(71)where

The term A n can be underestimated thanks to the convexity inequality (37):(72)

Thanks to the characterization (35) of the upstream mobility, the term B n can be underestimated bywhatever the mean value between and . Choosingwhich lies between and because of the definition (38) of , one gets thatwhereis the quantity that we want to bound from above, and

The quantitywhich is an approximation of at x K , is supposed to be bounded from below by some quantity M depending only on the regularity of the mesh. Since the method preserves the positivity of the saturations , one obtains that(73)thanks to the global conservativity of the scheme (28). As a consequence, one gets that

One concludes by noticing that Φ is bounded from below, hence

A.2 Relation (45) as an optimality condition

Before justifying why (45) can be seen as an optimality condition for the discrete energy under the mass contraint (47), let us first notice that sincethe discrete counterpartof the energy (5) can be rewritten as

The second term in the above equality does not depend on s n , hence it can be omitted in (48) and we can write the discrete energy as a function of only, i.e., .

A second preliminary remark is the following: if ϕ K m K  = 0 for some , then has no influence on the energy . Its value can be fixed arbitrarily, for instance by (45). We assume for simplicity that ϕ K m K  > 0 for all even though this assumption can be easily bypassed.

Let us now go to the constrained optimization problemwhere, setting , we denoted

Note that s n necessarily belongs to otherwise would be infinite. The problem is equivalent to the saddle-point problem

We can swap the and the and optimize w.r.t. s n , so that we get the optimality conditionwhich is equivalent to (45).

All Tables

Table 1

Choice (a) of mobility and pressure functions, convergence towards (64).

In the text
Table 2

Choice (b) of mobility and pressure functions, convergence towards (63).

In the text
Table 3

Choice (c) of mobility and pressure functions, convergence towards (63).

In the text
Table 4

Approximation of (42) with a nonlinear DDFV scheme [62, 63] on the Kershaw mesh family from [93], final time T = 0.25. M is the mesh index, τ is the time step, errgs and errs respectively stand for the L 2(Ω × (0, T)) error on ∇ s and the L ∞((0, T); L 2(Ø)) error on s, whereas ordgs and ords are the corresponding convergence orders. N max and N mean are the maximal and mean numbers of Newton iterations, and s min minimal value of the approximation of s during the whole simulation.

In the text

All Figures

thumbnail Fig. 1

The maximal monotone graph π depicted on top is parametrized as with some monotone functions (red) and (blue) depicted in the bottom picture. Here, and satisfy .
In the text
thumbnail Fig. 2

L 2 (Ω × (0, t f)) error as a function of the mesh size for the scheme (28) in the isotropic case κ = 1 (solid red) and anisotropic case κ = 20 (solid blue) corresponding to upstream mobilities (29). Comparison with the solution to the linear scheme (43) in the isotropic case (dashed red) and anisotropic case (dashed blue).
In the text
thumbnail Fig. 3

Evolution of the relative energy as a function of time. We observe the exponential convergence towards until the machine precision is reached.
In the text
thumbnail Fig. 4

The Kershaw meshes are highly deformed topologically cartesian grids. The coarsest mesh of the family is depicted here.
In the text
thumbnail Fig. 5

Convergence towards the analytical solution (42) of the Fokker-Planck equation (41) for different values of the anisotropy ratio .
In the text
thumbnail Fig. 6

Plot of the log of the relative energy (in log scale) as a function of time.
In the text
thumbnail Fig. 7

Schematic representation of the injection and production wells.
In the text
thumbnail Fig. 8

Approximate saturation s n,h at time t = 0.2 (top left), t = 0.5 (top right), t = 1.5 (bottom left) and t = 2.5 (bottom right).
In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.