Abstract
A Discrete Duality Finite Volume (DDFV) method to solve on unstructured meshes the flow problems in anisotropic nonhomogeneous porous media with full Neumann boundary conditions is proposed in the present work. We start with the derivation of the discrete problem. A result of existence and uniqueness of a solution for that problem is given thanks to the properties of its associated matrix combined with adequate assumptions on data. Their theoretical properties, namely, stability and error estimates (in discrete energy norms and -norm), are investigated. Numerical test is provided.
1. Introduction and the Model Problem
Efficient schemes are required for addressing flow problems in geologically complex media. The most important criteria of efficiency are mass conservation in grid blocks, accurate approximation of Darcy velocity, capability for dealing with anisotropic flow on unstructured grids and diverse heterogeneities (relevant to absolute permeability, porosity, etc.), and robust and easy implementation. Schemes well known in the literature for meeting many of the previous criteria are the following: Mixed Finite Element methods (see, e.g., [1, 2]), Control-Volume Finite Element methods (see, e.g., [3, 4]), Mimetic Finite Difference methods (see, e.g., [5, 6] and references therein), Cell-Centered Finite Volume methods (see, e.g., [7–11] and certain references therein; see also [12–14]), Multipoint Flux Approximation (see, e.g., [15–18] and some contributions to convergence analysis of MPFA O-scheme like [19]) and Discrete Duality Finite Volume methods (DDFV methods for short).
The DDFV methods come in two formulations. The first formulation is based on interface flux computations for primary and dual meshes, accounting with the interface flux continuity (see, e.g., [20, 21]) and the second formulation of DDFV is based on pressure gradient reconstructions over a diamond grid (see [22–24]). Note that this second formulation attracted the attention of some mathematicians as Andreianov, Boyer, and Hubert who have greatly contributed to its mathematical development. Indeed, key ideas involved in the pressure gradient reconstruction approach have been generalized by these authors to nonlinear operators of Leray-Lions type. Motivated by the possibility of increasing the order of convergence of the pressure gradient reconstruction method for nonlinear operators, Boyer and Hubert have proposed in [25] the so-called modified DDFV.
Beyond flow problems, we find the applications of DDFV methods in many areas: the numerical modeling of the surface erosion occurring at a fluid/soil interface undergoing a flow process in [26], the discretisation of partial differential equation appearing in image processing in [27], the assessment of nuclear waste repository safety in the context of simulating flow, transport in porous media in [28], and so on.
For presenting our analysis of DDFV method, let us consider the 2D diffusion problem consisting of finding a function which satisfies the following partial differential equation associated with nonhomogeneous Neumann boundary conditions:where is a given open polygonal domain, is its boundary, and are two given functions (defined, resp., in and , at least almost everywhere in Lebesgue measure sense), and is the unit normal vector to outward to . The permeability , with , may be a full matrix depending on solely spatial variables and obeying the following conditions.(i)The primary mesh is such that the discontinuity of lies on mesh interfaces.(ii)Symmetry is(iii)Uniform ellipticity and boundedness arewhere denotes the Euclidian norm in and where are the components of satisfyingwhere denotes the subspace of made of piecewise constant functions defined in .
One also assume that and satisfy the following conditions.(i)Compatibility condition is(ii)Null average condition is(iii)Let us emphasize that the novelty of this work is that Neumann boundary conditions are imposed on the whole boundary, which causes additional difficulties in the analysis.
2. A Finite Volume Formulation of the Model Problem
We briefly present finite volume formulation of the model problem (1)-(2) on unstructured meshes into the same spirit as [20, 25, 29–31]. We assume that the diffusion matrix is a piecewise constant function over .This assumption is currently used at least in industrial problems (e.g., reservoir simulation problems). Indeed, a subsurface area is made up of a collection of various geologic formations that may be characterized at intermediate scales by averaged full permeability tensors over grid blocks of the primary grid: for more details on this topic, see [32, 33].
2.1. Formulation of the Discrete Problem
First of all, notice that the method under consideration is analyzed in this work for general polygonal domains covered with unstructured matching primary meshes made up of arbitrary convex polygons (see Figure 1). Let us consider some definitions needed in what follows.
Definition 1. A mesh defined on is compatible with the discontinuities of the permeability tensor if these discontinuities are located along the mesh interfaces.
Definition 2. One defines an edge-point as any point located over an edge and different from the extremities of that edge.
Definition 3. Two edge-points and are named “neighboring edge-points” if they share the same vertex in the sense that and belong to two different edges that intersect in .
Let us recall our main objectives in this work:(i)Compute at the cell-points (to be defined later) and at the interior vertices from the mesh the values of the unknown function as a solution (expected unique) of a linear system.(ii)Analyze the stability and the convergence of this solution in some discrete energy norm similar to the one introduced in [31].
In the context of unstructured primary meshes (including square primary meshes with cell-points chosen different from cell-centers), the definition of a discrete energy norm requires that the cell-points lie inside certain perimeters. In this connection, the main steps for defining the cell-points are as follows:(a)Choose arbitrarily a unique point (different from a vertex) on each edge of the mesh . This process generates a finite family of edge-points.(b)Join every pair of neighboring edge-points by a dotted straight line. By this way, we generate an auxiliary mesh denoted .(c)Fix arbitrarily a unique point inside each intersection of a primary cell and an auxiliary cell. These points define a finite family of cell-points.Figure 2 illustrates the location of the edge-points and cell-points within both primary and auxiliary meshes.
Remark 4. Note that, in the 3D framework, for a given primary mesh the corresponding auxiliary mesh is generated very easily as follows. Any primary cell involves a finite number of faces and to each face is assigned one and only one face point lying necessarily on the boundary of . Therefore, one could associate with these face points the smallest convex polygon containing all of them. This is the process leading to generation of the auxiliary mesh associated with the primary mesh under consideration.
We should emphasize that we need to take cell-points inside auxiliary cells (see Figure 2) in view of achieving the following goal which is to define a vector space for the piecewise constant solution (named “weak approximate” solution in the sequel) and to equip it with a discrete energy norm.
Remark 5. Note that there exists a DDFV approach based on pressure gradient reconstructions that can address nonlinear elliptic problems: see [25, 34] for more details. Nevertheless in these works, the choice of edge-points depends strongly on that of cell-points as each edge-point coincides (by Definition 2) with the intersection of a primary edge and the straight line joining two cell-points located in both sides of that edge.
From the boundary-value problem theory (see, e.g., [35]), system (1)-(2) possesses a unique solution in under assumptions (3)–(7). We have assumed that the diffusion coefficient is a piecewise constant function over . The discontinuities of naturally divide into a finite number of convex subdomains . We now make the additional assumption that the restriction over of the exact solution to system (1)-(2), denoted by , satisfies
Let us now focus on a finite volume formulation of problem (1)-(2) in terms of a linear system which should be derived from the elimination of auxiliary unknowns, namely, edge-point pressures, in flux balance equations over primary cells and also dual cells (to be introduced later). This linear system should involve the real numbers and as discrete unknowns which are expected to be reasonable approximations of (cell-point pressures) and (interior vertex pressures), respectively, where and and where stands for dual mesh. We now give a description of the procedure leading to the linear discrete system.
Let be a primary cell, where is the corresponding cell-point. We integrate the two sides of the mass balance equation (1) in . Applying Ostrogradsky’s theorem to the integral in the left-hand side leads to computing the flux across the boundary of . Thanks to a suitable quadrature formula this computation yields a relation involving edge-point pressures. Due to the flux continuity over the mesh interfaces, the edge-point pressures can be eliminated from the previous relation.
As an illustration of this technique for computing the fluxes across primary cell boundaries, we consider the edge associated with the primary cell (see Figure 3).
(a)
(b)
Let be the diffusion tensor of the primary cell Denoting as the unit normal vector to exterior to , the flux expression over the edge viewed as part of the boundary of is given bywhere is the truncation error and where we introduce the following necessary ingredients:where (resp., ) is the angle defined by the vectors and (resp., and ) and where (resp., ) denotes the unit normal vector to (resp., to ) exterior to the half-plane from containing the point and bordered by the straight line (resp., ). Notice that , . Moreover if the primary mesh () is regular in the sense of Definition 6, we have and therefore , where is a certain real number not depending on .
In what follows, we will need the following notations. We denote by the set of all edge-points (from the primary mesh of course), the subset of made up of internal edge-points, that is, edge-points lying on , the subset of made up of edge-points lying on edges included in the boundary of , and (for ) the subset of made up of edge-points lying on the boundary of the primary cell .
We now introduce the notion of regular primary mesh that should play a central role in the sequel.
Definition 6. The set defines a regular mesh system if the following conditions are fulfilled.
There exist and , both of them mesh independent, such thatwhere stands for the set of primary cells and where is the distance between a cell-point or a vertex from and an edge-point from .
Proposition 7. Assume that (i) the primary mesh is compatible with the discontinuities of the permeability tensor (see Definition 1), (ii) the set is a regular mesh system, and (iii) relation (8) is honored.
Then, there exists a strictly positive number mesh independent such that
Using the previous notations and thanks to the consistency of the flux approximation across cell edges (see Proposition 7), one reasonably can approximate the flux balance within any primary cell as follows:where , , and are such that and where . It is clear that the number of discrete unknowns and is greater than the number of equations in system (14). It is then natural to complete this system with discrete equations obtained from mass flux balance over dual cells; see Figure 4 for the definition of the dual mesh.
In what follows, we will need to use the notion of pseudoedge associated with dual cells. Let us define now this notion that is illustrated in Figure 5.
Definition 8. Let and be two cell-points from the primary mesh (i.e., ) such that the corresponding primary cells and are adjacent, and consider (recall that , for , is the set of edge-points lying in the boundary of the primary cell ). The line defines a pseudoedge denoted by and “centered” on , with extremities and .
We will say that a pseudoedge is associated with a dual cell if it is lying in the boundary of .
Remark 9. The boundary of each dual cell is a union of a finite number of pseudoedges (see Figure 4).
Let us now look for discrete flux balance equations over dual cells. Integrating the two sides of the balance equation (1) in a dual cell , applying Ostrogradsky’s theorem for the left-hand side, and exploiting Remark 9 lead towhere stands for the outward unit normal vector to the boundary of and where is a pseudoedge associated with the dual cell . Recall that is the set of edge-points lying in the boundary of the dual cell .
Let us look for a flux approximation across the pseudoedge viewed as part of the boundary of . Denoting by the exact flux over , it can be expressed by the relation
By using the same process, the computation of the flux across iswhere (recall that denotes the angle defined by the vectors and ).
Similarly, the computation of the flux across leads towhere We should now formulate a global estimate for the truncation errors associated with the flux approximation over the pseudoedge.
Proposition 10. Under the same assumptions as those of Proposition 7, there exists a positive number mesh independent such that
The above inequality shows the consistency of the flux approximation across the pseudoedge . We summarize what precedes aswhere we have setWe deduce from (22) that an approximate flux balance equation over any dual cell can be formulated as follows:where , , and are such that and and where , a boundary-vertex are such that and .
Systems (14) and (24) naturally suggest defining a finite volume formulation of the diffusion problem (1)-(2) as follows.
Find and such that
Remark 11. It is useful to note that
2.2. Existence of Solutions to System (25)-(26) Conditions for Uniqueness
(i) Matrix Properties of the DDFV Problem (25)-(26). It is easily seen that the symmetry of the matrix associated with the linear system (25)-(26) essentially follows from the symmetry of the diffusion coefficient (see assumption (3)). We assume that all the cell-points and all the vertices (with respect to the primary mesh) are numbered. Therefore one can identify and with two disjoint subsets of the set of positive integers. To fix ideas, let us setThen denotes the total number of cell-points and stands for the total number of vertices. On the other hand, define the subvectors and byand set wherewhere is a subvector with components defined by the right-hand side of (25) and is a subvector with components defined by the right-hand side of (26).
Remark 12. The matrix satisfies the following properties:(1);(2);(3).
Let us introduce two vectors of named and and defined by
Proposition 13 (characterization of Kernel space of ). (i) The matrix is singular.(ii) Moreover, let denote the subset of defined as follows: (named Kernel space of in the sequel); then we have
Sketch for the Proof(i)The singularity of is an immediate consequence of Remark 12.(ii)Define the space byand endow it with the seminorm defined asThen, find the above identification of the Kernel space of using the following Lemma.
Lemma 14. The matrix is positive; that is,
Proof. Let , where and .
It follows from what precedes thatDefinewhere are two adjacent primary cells sharing the edge as interface containing the edge-point . Then, relation (37) becomeswhere denotes the subset of made of primary cells adjacent to the domain boundary . Let us prove that the homogenized symmetric matrix is positive definite; that is,Settingit is easy to check thatwhere we have setwhich are strictly positive numbers.
Since the diffusion matrix is symmetric and positive definite (see assumptions (3)-(4)), Cauchy-Schwarz’s inequality for the inner product associated with ensures thatas either and or and are not collinear. Therefore, and thus is a symmetric and positive definite matrix.
It follows from what precedes that the matrix possesses strictly positive eigenvalues. Let be its least eigenvalue. So we haveRemarking that there exist two real numbers and strictly positive depending exclusively on the geological structure of the domain such thatthe following result is easily seen.
Lemma 15.Moreover Thanks to inequality (45) and to Lemma the proof of Lemma 14 ends.
Lemma 16. The matrix satisfies the relationwhere is a strictly positive number mesh independent and where is a seminorm defined on by relation (35).
Proof. It is conducted with the same arguments as those developed in the previous proof, except that one should go much farther by proving that there exists mesh independent such thatThe eigenvalues of the symmetric positive definite matrix satisfy the so-called characteristic equation associated with ; that is,The least eigenvalue of denoted bywhere , is a strictly positive number. One can easily deduce thatwhere is a strictly positive number. We should bound the quantities , , and by mesh independent strictly positive numbers. Let us start first with We consider a change of coordinates by moving from the initial Cartesian coordinates to a local one, namely, , where is the edge-point on the interface between the cells and and where is a vector orthogonal to and oriented such that the basis change matrix is a rotation. Denoting the permeability tensor of the cell by in the initial Cartesian coordinates and by in the local coordinates we haveSimilarly, we get for the cell Then it is easy to check thatthat is,where denotes the determinant. Similarly, we have for the cell It follows from what precedes thatwhere and are given, respectively, by relations (43).
On one hand, we can deduce from (4) and Definition 6 thatOn the other hand, we remark thatwhere the set (introduced in (8)) depends exclusively on the lithologic structure of the medium . Then we deduce from (59), (60), and (61) thatRemarking thatand exploiting again (4) and Definition 6 lead to the following inequality:Thanks again to (4) and Definition 6 one can easily check thatLemma 16 follows from the combination of (53), (62), (64), and (65).
Proposition 17 (discrete compatibility condition). The right-hand side of the discrete system (25)-(26) satisfies the following discrete compatibility condition:
Proof. First of all, note that the double summation in the previous proposition is equal to zero thanks to Remark 11. Let us consider two vectors of called and and defined byAccording to the matrix form of the discrete system (see relation (31)), we haveIt follows from Remark 12 thatThereforeSimilarly we haveThis ends the proof of Proposition 17.
An immediate consequence of the previous result is the following existence result for a solution to the discrete problem (25)-(26).
Proposition 18 (existence result). (1) The linear system (25)-(26) possesses an infinite number of solutions.
(2) More precisely, if is a solution to the discrete system (25)-(26), then is the set of all solutions for this system.
The question to know how to get the physical solution of the model problem has a natural answer; that is, one should follow the same way as the continuous problem analysis. Since the dimension of is , required are two (linearly independent) discrete versions of the null average condition (7). This naturally leads to the following result.
Proposition 19 (uniqueness result). Under the following (null average) conditions,the discrete problem (25)-(26) possesses a unique solution.
Proof. Let us consider the spaceOne knows from Lemma 14 thatThen it follows that, for all in the subspace of made up of such that conditions (72) are fulfilled, one has .
3. Stability and Error Estimates
We deal here with a theoretical analysis of the solution for the discrete system (25), (26), and (72). Recall that the existence and uniqueness of that solution (under explicit conditions) is proven in the previous section. We assume in what follows that the primary mesh is regular in the sense of Definition 6.
3.1. Preliminaries and a Stability Result
Let us consider the auxiliary mesh introduced in the previous section (see Figure 4). Note that each mesh cell of contains either one cell-point or one corner point and only one which should be lying inside or on the boundary of . In the sequel, a node is a cell-point or a corner point with respect to the primary mesh.
Definition 20. An auxiliary mesh cell is degenerate if the corresponding node is lying on the boundary of
In the sequel, we will say sometimes “auxiliary cell” instead of “auxiliary mesh cell.” We denote by ) the subspace of made of functions which satisfies conditions (72). This space is obviously nonempty as there is the null function. Moreover, the solution of the discrete system (25), (26), and (72) could clearly be identified with one (and only one) function from ). In the sequel we denote by such a function named “cellwise constant (approximate) solution” of the diffusion problem (1)-(2). Let us endow with the following discrete seminorm.
For all , definewhere “s.t.” means “such that.” As this seminorm is actually a discrete energy norm for the space ), we denote it by as far as functions from ) are concerned. A norm on the space is defined by the mapping We focus here on the proof of the stability of the solution to the system of (25), (26), and (72). For this purpose, we need to introduce as in [31] a preliminary result, namely, a discrete version of Poincaré inequality which is based upon the following ingredients. Consider the following spacesand the linear operators and defined as follows: with and , where and are, respectively, subspaces of and made of functions satisfying conditions (72)-(i) and (72)-(ii), respectively.
Let us equip the function spaces and with the following discrete energy norms:These mappings are only seminorms, respectively, for and that one can transform into norms for these spaces as follows:where we have set and . The following results are key ingredients for proving the stability of the discrete solution in the sense of the discrete energy norm defined previously on the space .
Lemma 21 (discrete versions of Poincaré inequality). We have the following inequalities:where represents diverse strictly positive numbers mesh independent.
Note that the right-hand side of inequalities (81) and (82) is, respectively, norms for and . Since is bounded, the previous lemma permits seeing that these norms are equivalent to standard ones, namely, (80). For proving the preceding lemma one can use the same arguments as in [36].
Lemma 22 (a key-result). Let be a regular mesh system (defined on) in the sense of Definition 6 and denote by the corresponding dual mesh. For , denote by the value of in the control volume . Let be a discrete trace function defined a.e. (for the ()-Lebesgue measure) by on , for all , where denotes the set of mesh elements adjacent to the domain boundary . ThenSimilarly, we havewhere stands for diverse positive numbers mesh independent.
Proof. See, for instance, [36].
Remark 23. It is more than useful to note that
Let us give now one of the main results of this section.
Proposition 24 (stability result). The cellwise constant approximate solution of the diffusion problem (1)-(2) satisfies the following inequality:where is a strictly positive real number not depending on the spatial discretization.
Proof. Multiplying (25) by and (26) by and summing over and , respectively, lead toRecall that