Numerical Solution of 2-D Diffusion Problems Using Discrete Duality Finite Volume Method on General Boundary Conditions

Hubert Donfack¹, Aubin Kinfack Jeutsa²

¹Faculty of Science, The University of Bamenda, Bamenda, Cameroon.
²Higher Technical Teachers’ Training College, The University of Buea, Kumba, Cameroon.

DOI: 10.4236/jamp.2022.106135                

Abstract

In this paper, we present a discrete duality finite volume (DDFV) method for 2-D flow problems in nonhomogeneous anisotropic porous media under diverse boundary conditions. We use the discrete gradient defined in diamond cells to compute the fluxes. We focus on the case of Dirichlet, full Neumann and periodic boundary conditions. Taking into account the periodicity is the main new ingredient with respect to our recent works. We explain the procedures step by step, for numerical solutions. We develop a matlab code for algebraic equations. Numerical tests were provided to confirm our theoretical results.

Keywords

Flow Problems, Nonhomogeneous Anisotropic Media, Finite Volumes, Test Simulations

1. Introduction

The discrete ion finite volume is one of the new generation finite volumes very popular today in Geoscience Engineering. The works from [1] [2] [3] [4] [5] have been the first to propose the DDFV method as used today for anisotropic flow problems. The formulation in terms of explicit discrete duality has been introduced in [5] [6]. The aim of this paper is to reformulate approach our DDFV approach by using the discrete gradient operators defined on the diamond meshes in order to show that, it is well adapted to very general meshes including the case of non-conformal locally refined meshes. We focus on general grids a survey of DDFV formulation for anisotropic flow problems under: 1) Dirichlet; 2) Full Neumann and 3) periodic boundary conditions. With regard to the mathematical analysis, we will refer to [5] [6] for details. An example of such a problem reads as follows: Find a function $\phi$ defined in $\bar{\Omega }$, that satisfies the following partial differential equation:

$-div\left(Kgrad\left(\phi \right)\right)=f\text{in}\Omega =\left]\mathrm{0,1}\right[\times \left]\mathrm{0,1}\right[\text{with}\text{appropriate}\text{boundary}\text{conditions}$ (1)

where f is a given function and $\bar{\Omega }$ is a closure of $\Omega$ that is a given open polygonal domain (not necessarily convex). On the other hand, $K=K\left(x\right)$, with $x={\left(x_1,x_2\right)}^{\text{T}}\in \Omega$ is a piecewise constant function in $\Omega$ and is symmetric, i.e.

$K_{ij}\left(x\right)=K_{ij}\left(x\right)\forall 1\le i;j\le 2$ (2)

$\exists {\delta }_{\pm }\in {ℝ}_{+}^{\star }\text{such}\text{that}\forall \xi \in {ℝ}^2,{\delta }_{-}{\Vert \xi \Vert }^2\le {\left(\xi \right)}^{\text{T}}K\left(x\right)\xi \le {\delta }_{+}{\Vert \xi \Vert }^2$ (3)

in $\Omega$, where $\Vert \Vert$ and ${\left(\mathrm{.}\right)}^{\text{T}}$ denote respectively the euclidian norm and the transposition operator in ${ℝ}^2$.

2. DDFV Formulation for Model Problem (1) on General Boundary Conditions

We first recall Gauss’ divergence theorem which is very useful in this section.

Theorem 1 (Gauss’ divergence theorem) Let’s S be a closed surface bounding a soled D, oriented outwards. Let $\stackrel{\to }{F}$ be a vector field with continuous partial derivatives then

${\displaystyle {\iint }_S\stackrel{\to }{F}\text{d}\stackrel{\to }{s}}={\displaystyle {\iiint }_D\nabla \cdot \stackrel{\to }{F}\text{d}v}$ (4)

Let us focus on a DDFV formulation of the problem (1) in terms of a linear system which involves ${\left\{{\phi }_P\right\}}_{P\in \mathcal{P}}$ and ${\left\{{\phi }_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}$ as discrete unknowns expected to be close approximations of ${\left\{{\bar{\phi }}_P\right\}}_{P\in \mathcal{P}}$ (cell-point pressures) and

${\left\{{\bar{\phi }}_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}$ (vertex pressures) respectively, where ${\bar{\phi }}_p=\bar{\phi }\left(x_1^P,x_2^P\right)$ and

${\bar{\phi }}_{D^{\star }}=\bar{\phi }\left(x_{\text{1}}^{D^{\star }}\mathrm{,}{x}_{\text{2}}^{D^{\star }}\right)$ and where $\mathcal{D}$ represents the dual mesh. The DDFV theory

exposed in this work is inspired from the one developed in [6]. However we use the discrete gradient operators to simplify the heaviness concerns the DDFV discretization of balance equation for boundary cells from both primal and dual meshes.

2.1. Domain Discretizations: Some Definitions and Notations

For getting the DDFV formulation of the continuous problem (1), we first introduce a primary mesh $\mathcal{P}$ over $\Omega$ possessing the following properties: 1) $\mathcal{P}$ divides $\bar{\Omega }$ int o a finite number of convex mesh elements. The mesh element vertices are named the primary mesh vertices. 2) The discontinuity points of $K=K\left(x\right)$ are located in the mesh element interfaces. 3) For any edge of $\mathcal{P}$ arbitrarily exhibited one point and named edge-point in what follows, denote by the set of edge-points, ${\mathcal{E}}^{int}$ the set of internal edge-point and ${\mathcal{E}}^{ext}$ the set of boundary edge-points.

Joining the edge points in suitable way defines an auxiliary mesh $\mathcal{A}$ necessary for locating a family of cellpoints in mesh element (with one point per mesh element). The auxiliary mesh plays a key role in the theoretical analysis of our DDFV solution see [6] [7] for details. Note the (Figure 1) below illustrate the process of discretization. Note that degenerate dual cells are involved along the boundary for Neumann and periodic boundary conditions.

Given two adjacent cellpoints P and L sharing $\sigma =\left[A^{\star }{B}^{\star }\right]$ as common interface, for each $I\in \left[A^{\star }{B}^{\star }\right]$, can we introduce the quadrangle $\left(P\mathrm{,}{A}^{\star }\mathrm{,}L\mathrm{,}{B}^{\star }\right)$ diamond cell ${\mathfrak{D}}_{\sigma {\sigma }^{\star }}$ whose diagonals are $\sigma =\left[A^{\star }{B}^{\star }\right]$ and ${\sigma }^{\star }=\left[PL\right]$, as shown in (Figure 2) below. Notice that the diamond cells are the union of four disjoint convex triangles. Furthermore, if $I\in \left[A^{\star }{B}^{\star }\right]\cap \partial \Omega$ then the quadrilateral ${\mathfrak{D}}_{\sigma {\sigma }^{\star }}$ degenerate in a single triangle. The set of diamond cells is denoted by $D$ and we have

$\bar{\Omega }=\underset{\mathfrak{D}\in D}{\cup }\bar{\mathfrak{D}}$

Remark 1 Note that the set $\mathcal{E}$ will sometimes be identified with the set of diamond meshes since there is a trivial bijection between the two sets.

• ${\alpha }_D$ angle between ${\tau }_{A^{\star }{B}^{\star }}$ and ${\tau }_{PL}$

• ${\alpha }_{B^{\star }IL}$ angle between ${\tau }_{A^{\star }{B}^{\star }}$ and ${\sigma }_L$

• ${\alpha }_{B^{\star }IL}=\frac{\Pi }{2}+{\theta }_h^{L\mathrm{,}I}$

Figure 1. Domain discretizations for Dirichlet and Neumann boundary conditions (left) and for periodic boundary (right). The including is in black full lines, the corresponding auxiliary mesh in black dotted lines and the associated dual mesh in red discontinuous lines.

Figure 2. Examples of diamond cells. Left: A normal diamond cell. Right: A degenerate diamond cell.

• ${\xi }_{\left[PL\right]}^{A^{\star }}$ the unit normal to $\left[PL\right]$ exterior to triangle $PL{A}^{\star }$

• ${\xi }_{\left[A^{\star }{B}^{\star }\right]}^L=-{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P$

• $m_{\sigma }$ the length of $\sigma$

• $m_{{\sigma }^{\star }}$ the length of ${\sigma }^{\star }$

• $m_{\mathfrak{D}}$ the measure of the diamond cell

• $D_{A^{\star }}=\left\{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D\backslash {\mathfrak{D}}_{\sigma {\sigma }^{\star }}\cap C_{A^{\star }}\ne \varnothing \right\}$

• $D_P=\left\{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D\backslash {\mathfrak{D}}_{\sigma {\sigma }^{\star }}\cap C_P\ne \varnothing \right\}$

• ${\mathfrak{Q}}_{L{A}^{\star }}$ quarter diamond defines by the triangle ( $L{A}^{\star }I$ )

• ${\tilde{\mathfrak{D}}}^L$ half diamond defines by the triangle ( $L{A}^{\star }{B}^{\star }$ )

• ${\mathfrak{D}}_{\sigma {\sigma }^{\star }}={\tilde{\mathfrak{D}}}^L\cup {\tilde{\mathfrak{D}}}^P\mathrm{<mo>\; =\; </mo>}{\mathfrak{Q}}_{L{A}^{\star }}\cup {\mathfrak{Q}}_{L{B}^{\star }}\cup {\mathfrak{Q}}_{P{A}^{\star }}\cup {\mathfrak{Q}}_{P{B}^{\star }}$

An elementary geometry in triangles (PLI), (PID) and (LID) with reference (Figure 3) permit us to write:

$(\begin{array}{l}{h}_{IL}{\xi }_{\left[LI\right]}^{A^{\star }}+h_{IP}{\xi }_{\left[PI\right]}^{A^{\star }}-h_{PL}{\xi }_{\left[PL\right]}^{A^{\star }}={\stackrel{\to }{0}}_{{ℝ}^2}\\ h_{DI}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P+h_{IP}{\xi }_{\left[PI\right]}^{A^{\star }}-h_{PD}{\xi }_{\left[PL\right]}^{A^{\star }}={\stackrel{\to }{0}}_{{ℝ}^2}\\ h_{IL}{\xi }_{\left[LI\right]}^{A^{\star }}-h_{DI}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P-h_{DL}{\xi }_{\left[PL\right]}^{A^{\star }}={\stackrel{\to }{0}}_{{ℝ}^2}\end{array}$ (5)

this implies,

${\xi }_{\left[LI\right]}^{A^{\star }}=\frac{h_{DI}}{h_{IL}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P+\frac{h_{DL}}{h_{IL}}{\xi }_{\left[PL\right]}^{A^{\star }}$ (6)

${\xi }_{\left[PI\right]}^{A^{\star }}=-\frac{h_{DI}}{h_{IP}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P+\frac{h_{DP}}{h_{IP}}{\xi }_{\left[PL\right]}^{A^{\star }}$ (7)

Main assumptions:

( ${\mathcal{A}}_1$ ) We assume that the primary mesh $\mathcal{P}$ is compatible with the discontinuities of the permeability tensor K defined in $\Omega$.

( ${\mathcal{A}}_2$ ) The permeability discontinuities divide $\Omega$ into a finite number of convex polygonal subsets denoted by ${\left\{{\Omega }_s\right\}}_{s\in S}$.

We suppose that the restriction over ${\bar{\Omega }}_s$ of the exact solution $\phi$ to the model problem (1), denoted by ${\phi |}_{{\Omega }_s}$, satisfies the following property:

${\phi |}_{{\Omega }_s}\in C^2\left({\bar{\Omega }}_s\right)\forall s\in S$ (8)

Figure 3. Diamond cells and notations.

Let us recall below the diverse boundary conditions that we are going to investigate in this paper:

1) Dirichlet boundary conditions

$\phi =0\text{on}\Gamma =\partial \Omega$ (9)

2) Full Neumann boundary conditions

$-Kgrad\left(\phi \right)\cdot \eta =g\text{on}\Gamma =\partial \Omega$ (10)

where $\eta$ is a outward unit normal vector and g is a given function

3) Periodic boundary conditions

$\begin{array}{l}\begin{array}{l}\phi \left(\mathrm{.,0}\right)=\phi \left(\mathrm{.,1}\right)\\ K\left(\mathrm{.,0}\right)grad\phi \left(\mathrm{.,0}\right)\cdot n_{{\Gamma }_{So}}+K\left(\mathrm{.,1}\right)grad\phi \left(\mathrm{.,1}\right)\cdot n_{{\Gamma }_{No}}=0\end{array}\}\text{in}\left[\mathrm{0,1}\right]\\ \begin{array}{l}\phi \left(0,.\right)=\phi \left(\mathrm{.,1}\right)\\ K\left(\mathrm{.,0}\right)grad\phi \left(\mathrm{0,.}\right)\cdot n_{{\Gamma }_{We}}+K\left(\mathrm{1,.}\right)grad\phi \left(\mathrm{1,.}\right)\cdot n_{{\Gamma }_{Es}}=0\end{array}\}\text{in}\left[\mathrm{0,1}\right]\end{array}$ (11)

where ${\Gamma }_{No}$ (northern boundary), ${\Gamma }_{So}$ (southern boundary), ${\Gamma }_{Es}$ (eastern boundary) and ${\Gamma }_{We}$ (western boundary) define a partition of the domain boundary denoted by $\Gamma$, and $n_{{\Gamma }_{No}}\mathrm{,}{n}_{{\Gamma }_{So}}\mathrm{,}{n}_{{\Gamma }_{Es}}\mathrm{,}{n}_{{\Gamma }_{We}}$ the corresponding outward unit normal vectors

Due to periodicity conditions on the flux (see Equations (11), the source-term f should satisfy the following compatibility condition:

${\displaystyle {\int }_{\Omega }}f\text{d}x=0$ (12)

Note that the set $\mathcal{V}$ of primal mesh vertices contains the set:

made of four corner-points. To get advantage of the periodicity setting of the problem, the boundary-vertices are distributed along the domain boundary in such a way that the orthogonal projection of a boundary-vertex M (different from a corner-point) on the opposite side of $\Gamma$ is also a boundary-vertex (different from a corner-point) denoted by $M^{\perp }$. Figure 1 (right) illustrates this fact that one can express in other words by:

2.2. The DDFV Balance Equation in an Internal Cell

For an internal cell the way for getting the DDFV balance equation is the same for Dirichlet, Neumann and periodic boundary conditions. Let $C_P$ be an internal primary cell where P is the corresponding cellpoint. We start, with introducing a discrete gradient operator ${\nabla }^{\mathcal{T}}$ defined to be constant on each quarter diamond cell ${\mathfrak{Q}}_{L{A}^{\star }}$ associated to the control volume $C_L$

${\nabla }_{{\mathfrak{Q}}_{L{A}^{\star }}}^{\mathcal{T}}{\phi }^{\mathcal{T}}=\frac{1}{\cos \left({\theta }_h^{L,I}\right)}\left[\frac{{\phi }_I-{\phi }_{A^{\star }}}{h_{I{A}^{\star }}}{\xi }_{\left[IL\right]}^{A^{\star }}+\frac{{\phi }_I-{\phi }_L}{h_{IL}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\right]$ (13)

Then integrating two sides of the balance Equation (1) in $C_L$ and apply Gauss’s divergence theorem to the left-hand side. The discrete balance equation follows from flux computations over the boundary of $C_L$ by using a suitable gradient operator ${\nabla }^{\tau }$ defined above. Thus we have:

$\begin{array}{c}{F}_{\left[A^{\star }{B}^{\star }\right]}^L=-{\displaystyle {\int }_{\left[A^{\star }I\right]}}{K}^{L}grad\left(\phi \right){\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\text{d}\gamma -{\displaystyle {\int }_{\left[I{B}^{\star }\right]}}{K}^{L}grad\left(\phi \right){\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\text{d}\gamma \\ =-h_{I{A}^{\star }}\langle K^L{\nabla }_{{\mathfrak{Q}}_{L{A}^{\star }}}^{\tau }\phi \mathrm{,}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\rangle -h_{I{A}^{\star }}\langle K^L{\nabla }_{{\mathfrak{Q}}_{L{A}^{\star }}}^{\tau }\phi \mathrm{,}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\rangle +T_{\left[A^{\star }{B}^{\star }\right]}^L\\ =b_h\left(K^L\right)\left({\phi }_{A^{\star }}-{\phi }_{B^{\star }}\right)+\frac{h_{A^{\star }{B}^{\star }}}{h_{IL}}{a}_h\left(K^L\right)\left({\phi }_L-{\phi }_I\right)+T_{\left[A^{\star }{B}^{\star }\right]}^L\end{array}$ (14)

where

$a_h\left(K^L\right)=\frac{\langle K^L{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\mathrm{,}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\rangle }{\cos \left({\theta }_h^{L,I}\right)}$

$b_h\left(K^L\right)=\frac{\langle K^L{\xi }_{\left[IL\right]}^{A^{\star }}\mathrm{,}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L\rangle }{\cos \left({\theta }_h^{L,I}\right)}$

${\theta }_h^{L\mathrm{,}I}=meas\left(\widehat{{\sigma }_L\mathrm{,}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^L}\right)$

$h=\max \left\{size\left(\mathcal{P}\right),size\left(\mathcal{D}\right)\right\}$, Figure 3 gives an illustration of previous comment. Note that $0\le {\theta }_h^{P,I},{\theta }_h^{L,I}<\frac{\pi }{2}$ and therefore $0<\cos \left({\theta }_h^{P,I}\right),\cos \left({\theta }_h^{L,I}\right)\le 1$.

Definition 1 The system $\left(\mathcal{P}\mathrm{<mo>\; ;\; </mo>}\mathcal{D}\right)$ defines an eligible system of meshes if the following conditions are fulfilled:

1) There exists $\theta \in \left]0\mathrm{,}\frac{\pi }{2}\right]$, mesh independent, such that:

$0\le {\theta }_h^{P\mathrm{,}I}\le \frac{\pi }{2}-\theta \forall P\in \mathcal{P}\forall I\in {\mathcal{E}}^P$ (15)

2) There exists $0<\varpi \le 1$, mesh independent, such that:

$\forall P\in \mathcal{P}\forall I\in {\mathcal{E}}^P\varpi h\le h_{PI}\mathrm{,}{h}_{A^{\mathrm{<mo>\; *\; </mo>}}{B}^{\mathrm{<mo>\; *\; </mo>}}}\le h$ (16)

where $h_{PI}=\Vert \stackrel{\leftarrow }{PI}\Vert$ and $h_{A^{\star }{B}^{\star }}=\Vert \stackrel{\leftarrow }{A^{\star }{B}^{\star }}\Vert$.

Following ideas we have exposed in [6], and the fact that approximate fluxes $F_{\left[A^{\star }\mathrm{<mo>\; ;\; </mo>}{B}^{\star }\right]}^P$ and $F_{\left[A^{\star }\mathrm{<mo>\; ;\; </mo>}{B}^{\star }\right]}^L$ meet the principle of flux continuity over the interface

between $C_P$ and $C_L$ if and only if the approximate edge point pressure ${\phi }_I$ satisfies to the following relation:

$\begin{array}{c}{\phi }_I=\frac{b_h\left(K^P\right)+b_h\left(K^L\right)}{a_h\left(K^P\right)\frac{h_{A^{\star }{B}^{\star }}}{h_{PI}}+a_h\left(K^L\right)\frac{h_{A^{\star }{B}^{\star }}}{h_{IL}}}\{\left({\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right)+a_h\left(K^L\right)\frac{h_{A^{\star }{B}^{\star }}}{h_{IL}}{\phi }_L\\ +a_h\left(K^P\right)\frac{h_{A^{\star }{B}^{\star }}}{h_{PI}}{\phi }_P\}\end{array}$ (17)

Proposition 2 Under the only assumption (8) we have:

$\begin{array}{c}{F}_{\left[A^{\ast }{B}^{\ast }\right]}^P\approx \left[\frac{a_h\left(K^P\right)a_h\left(K^L\right)h_{A^{\ast }{B}^{\ast }}}{a_h\left(K^P\right)h_{IL}+a_h\left(K^L\right)h_{PI}}\right]\left[{\phi }_P-{\phi }_L\right]\\ +\left[\frac{a_h\left(K^L\right)b_h\left(K^P\right)h_{PI}-a_h\left(K^P\right)b_h\left(K^L\right)h_{IL}}{a_h\left(K^L\right)h_{PI}+a_h\left(K^P\right)h_{IL}}\right]\left[{\phi }_{B^{\ast }}-{\phi }_{A^{\ast }}\right]\end{array}$ (18)

In addition, if the system of meshes $\left(\mathcal{P}\mathrm{<mo>\; ;\; </mo>}\mathcal{D}\right)$ is eligible in the sense of Definition 1, the truncation error $T_{\left[A^{\ast }{B}^{\ast }\right]}^P$ (also denoted by $T^{P\mathrm{,}I}$ ) associated with this flux approximation satisfies the following inequality

$\left|T^{P\mathrm{,}I}\right|\equiv \left|T_{\left[A^{\ast }{B}^{\ast }\right]}^P\right|\le C{h}^2$ (19)

where $A^{\star }=A^{\star }\left(P\mathrm{,}I\right)$, $B^{\star }=B^{\star }\left(P\mathrm{,}I\right)$, $L=L\left(P\mathrm{,}I\right)$ and where C is a mesh independent positive number. $\square$

Recall that only the case ${\mathcal{E}}^P\cap {\mathcal{E}}^{ext}=\varnothing$, with $P\in \mathcal{P}$, is concerned by the previous result. Summing the two sides of relation (18) on $I\in {\mathcal{E}}^P$ leads to the following result.

Proposition 3 Let us suppose that the system of meshes $\left(\mathcal{P}\mathrm{<mo>\; ;\; </mo>}\mathcal{D}\right)$ is eligible in the sense of Definition 1. Under the assumption (8), the discrete balance equation in any primary cell $C_P$, with ${\mathcal{E}}^P\cap {\mathcal{E}}^{ext}=\varnothing$, reads as:

$\begin{array}{l}{\displaystyle \underset{I\in {\mathcal{E}}^P}{\sum }}\{T^{P\mathrm{,}I}+\left[\frac{a_h\left(K^P\right)a_h\left(K^L\right)h_{A^{\star }{B}^{\star }}}{a_h\left(K^P\right)h_{IL}+a_h\left(K^L\right)h_{IP}}\right]\left[{\phi }_P-{\phi }_L\right]\\ +\left[\frac{a_h\left(K^P\right)b_h\left(K^L\right)h_{IL}-b_h\left(K^P\right)a_h\left(K^L\right)h_{IP}}{a_h\left(K^P\right)h_{IL}+a_h\left(K^L\right)h_{IP}}\right]\left[{\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right]\}={\displaystyle {\int }_{C_P}}f\text{d}x\end{array}$ (20)

where $A^{\star }=A^{\star }\left(P\mathrm{,}I\right)$, $B^{\star }=B^{\star }\left(P\mathrm{,}I\right)$ and $L=L\left(P\mathrm{,}I\right)$ (see Section 2 devoted mainly to Notations).

Moreover, the truncation error $T^{P\mathrm{,}I}$ satisfies the inequality (19).

It is clear that the number of discrete unknowns ${\left\{{\phi }_P\right\}}_{P\in \mathcal{P}}$ and ${\left\{{\phi }_{A^{\star }}\right\}}_{A^{\star }\in \mathcal{D}}$ is greater than the number of discrete balance equations given by the system of Equation (20) valid for all $P\in \mathcal{P}$. We naturally should close this system with discrete equations obtained from mass balance equations over dual cells. It is our purpose now to look for discrete balance equations over dual cells $C_{A^{\star }}$. So, we integrate the two sides of (1) in $C_{A^{\star }}$. The left-hand side is the flow exchanged between $C_{A^{\star }}$ and outside of this cell, whereas the right-hand side is the term-source contribution (for a fixed time period) Let us look or a flux approximation across the pseudo-edge $\left[PIL\right]$ viewed as part of the boundary of $C_{B^{\star }}$. Denoted by $F_{\left[PIL\right]}^{B^{\star }}$, it can be expressed by the relation

$F_{\left[PIL\right]}^{B^{\star }}=-{\displaystyle {\int }_{\left[PI\right]}}{K}^{P}grad\left(\phi \right){\xi }_{\left[PI\right]}^{B^{\star }}\text{d}\gamma -{\displaystyle {\int }_{\left[IL\right]}}{K}^{L}grad\left(\phi \right){\xi }_{\left[IL\right]}^{B^{\star }}\text{d}\gamma$ (21)

By using the same process, the approximate flux across the pseudo-edge is

$F_{\left[PIL\right]}^{B^{\star }}=-h_{PI}\langle K^P{\nabla }_{{\mathfrak{Q}}_{P{B}^{\star }}}^{\tau }\phi \mathrm{,}{\xi }_{\left[PI\right]}^{B^{\star }}\rangle -h_{IL}\langle K^L{\nabla }_{{\mathfrak{Q}}_{L{B}^{\star }}}^{\tau }\phi \mathrm{,}{\xi }_{\left[IL\right]}^{B^{\star }}\rangle +T^{B^{\star }\mathrm{,}I}$ (22)

Summing the two sides of the Equation (22) on $I\in {\mathcal{E}}^{B^{\star }}$ leads to what follows.

Proposition 4 Under the assumption (8), the flux balance equation for any interior dual cell $C_{B^{\star }}$ reads,:

$\underset{I\in {\mathcal{E}}^{B^{\star }}}{{\displaystyle \sum }}{F}_{\left[PIL\right]}^{B^{\star }}={\displaystyle {\int }_{C_{B^{\star }}}}f\left(x\right)\text{d}x\forall B^{\star }\in C_{B^{\star }}$ (23)

In addition, if the system of meshes $\left(\mathcal{P}\mathrm{<mo>\; ;\; </mo>}\mathcal{D}\right)$ is eligible in the sense of Definition 1, the truncation error $T^{B^{\star }\mathrm{,}I}$ obeys the following inequality:

$\left|T^{B^{\star }\mathrm{,}I}\right|\le C{h}^2$ (24)

where C is a mesh independent positive number.

2.3. The DDFV Balance Equation in Boundary Cells

2.3.1. Case of Dirichlet Boundary Conditions

As shown in (Figure 1), there is no need of degenerate dual cell in the domain discretization for Dirichlet boundary conditions. Thus the set of boundary cells reduces to the set of primary cells adjacent to the boundary. Note that for $f\in L^2\left(\Omega \right)$ one easily proves that Equations (1) and (9) get a unique variational

solution in the Sobolev space $H_0^1\left(\Omega \right)$. Letting $\underset{I\in {\mathcal{E}}^{int}\cap {\mathcal{E}}^P}{{\displaystyle \sum }}{F}_{\left[A^{\star }{B}^{\star }\right]}^P$ represent the left of the balance equation in primary cell $C_P$ adjacent to $\Gamma$ reads as:

$\underset{I\in {\mathcal{E}}^{int}\cap {\mathcal{E}}^P}{{\displaystyle \sum }}{F}_{\left[A^{\star }{B}^{\star }\right]}^P+\underset{I\in {\mathcal{E}}^{ext}\cap {\mathcal{E}}^P}{{\displaystyle \sum }}{a}_h\left(K^P\right)\frac{h_{A^{\star }{B}^{\star }}}{h_{PI}}{\phi }_P={\displaystyle {\int }_{C_P}}f\text{d}x$ (25)

Because the degenerate discrete gradient reads as:

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{ext}}^{\tau }{\phi }^{\tau }=\frac{1}{\sin \left({\alpha }_D\right)}\left(\frac{-{\phi }_P}{h_{LP}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P\right)$ (26)

2.3.2. Case of Full Neumann Boundary Conditions

In the Sobolev space $H^1\left(\Omega \right)$ the existence and uniqueness of a variational solution to (1) - (10) are ensured if:

$f\in L^2\left(\Omega \right)g\in L^2\left(\Gamma \right){\displaystyle {\int }_{\Omega }}f\text{d}x-{\displaystyle {\int }_{\Gamma }}g\text{d}\gamma =0$ (27)

Let ${\Gamma }_P$ denote the boundary of any primary cell $C_P$, then the DDFV balance equation in a primary cell adjacent to domain boundary by using a suitable external discrete gradient operator reads as:

$\underset{I\in {\mathcal{E}}^{ext}\cap {\mathcal{E}}^P}{{\displaystyle \sum }}\left[\mathrm{...}\right]={\displaystyle {\int }_{C_P}}f\text{d}x-{\displaystyle {\int }_{\Gamma \cap {\Gamma }_P}}g\text{d}\gamma =0$ (28)

Note that the degenerate discrete gradient reads as:

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{ext}}^{\tau }{\phi }^{\tau }=\frac{1}{\sin \left({\alpha }_D\right)}\left(\frac{{\phi }_{B^{\star }}-{\phi }_{A^{\star }}}{h_{B^{\star }{A}^{\star }}}{\xi }_{\left[IP\right]}^{B^{\star }}+\frac{{\phi }_I-{\phi }_P}{h_{PI}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P\right)$ (29)

And ${\phi }_I$ is given by Equation (17) with $K^L=0$, finally the degenerate discrete gradient reads as:

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{ext}}^{\tau }{\phi }^{\tau }=\frac{{\phi }_{B^{\star }}-{\phi }_{A^{\star }}}{h_{B^{\star }{A}^{\star }}\sin \left({\alpha }_D\right)}\left({\xi }_{\left[IP\right]}^{B^{\star }}+\frac{b_h\left(K^P\right)}{a_h\left(K^P\right)}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P\right)$ (30)

It is easy to check that the balance equation in a dual cell (adjacent or not) to the domain boundary reads as:

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{ext}}^{\tau }{\phi }^{\tau }=\frac{{\phi }_{B^{\star }}-{\phi }_{A^{\star }}}{h_{B^{\star }{A}^{\star }}\sin \left({\alpha }_D\right)}\left({\xi }_{\left[IP\right]}^{B^{\star }}+\frac{b_h\left(K^P\right)}{a_h\left(K^P\right)}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P\right)$

2.3.3. Case of Periodic Boundary Conditions

Note that the periodicity condition (11) implies the periodicity of the discrete solution. Clearly the source f and the tensor K should be extended to ${ℝ}^2$ by periodicity. So there exists a periodic partition of whole space ${ℝ}^2$ into control volumes. Our strategy to obtain discrete balance equations associated with degenerated dual cells is to consider a DDFV formulation in whole space ${ℝ}^2$ (without distinguishing between boundary and interior control volumes). We extend the domain $\Omega$ by introducing the fictitious domain around the boundary of the initial domain $\Omega$, the fictitious points (cell points and vertex points) are defined so as to have their corresponding in $\Omega$, due to periodicity. To doing so, let’s compute the flux across the recomposed dual cells $C_{E^{\star }{E}^{\star \perp }}$ and $C_{\pi }$, we define the fictitious domain ${\Omega }^F$. Note that the original domain $\Omega$ is embedded inside a fictitious domain ${\Omega }^F$ such that ${\bar{\Omega }}^F=\bar{\Omega }\cup \Sigma \cup \partial {\Omega }^F$ The computational domain $\Omega$ is meshed with respect to the periodicity setting of the problems. The boundary-vertices ( $A^{\star }\in \partial {\Omega }^F$ ) are distributed along the domain boundary $\partial {\Omega }^F$ in such a way that the orthogonal projection of a boundary-vertex on the opposite side of $\partial {\Omega }^F$ is a latest vertex-point inside the original domain $\Omega$ denoted by $A^{\star \perp }$ which means that ${\left({\phi }_{A^{\star \perp }}\right)}_{A^{\star \perp }\in \partial {\Omega }^F}={\left({\phi }_{A^{\star }}\right)}_{A^{\star }\in \Omega }$. The cell point pressure L are distributed inside the domain $\Sigma$ in such a way that the orthogonal projection of a boundary dual mesh on the opposite of $\partial \Omega$ is a least cell point inside the original domain $\Omega$ denoted by $L^{\perp }$ which means that ${\left({\phi }_L\right)}_{L\in \Sigma }={\left({\phi }_{L^{\perp }}\right)}_{L^{\perp }\in \Omega }$. As illustrated in (Figure 4) in the particular case of the uniform rectangular mesh. The strategy allows to have all the dual meshes inside the fictitious domain which greatly simplifies the writing of discrete balance equations. It’s important to note that the introduction of the fictitious domain also makes possible to deal efficiently the Neunman boundary conditions, for that the permeability tensor must be null in $\Sigma$ domain. Then we use the same process to compute the flux across the boundary cells by defining the appropriate discrete gradient operators. It follows that three types of discrete gradient operators on the diamond mesh are necessary for the computation of the flux whether for the internal or external meshes. Note that this domain extension technique was used in [8] for numerical treatment of initial-Boundary value problems with Mixed Boundary Conditions.

Definition 2 A diamond mesh is said to be fictitious or recomposed if at least one of its vertices is in $\Sigma$ domain. Let us denote $D=D^{int}\cup D^{Fict}$ where

Figure 4. A schematic illustration of the fictitious domain delimited by $\partial {\Omega }^F-\partial \Omega$ and recomposed dual cells around the domain $\partial \Omega$ painted in blue.

$D^{int}$ is the set of internal diamond cells and $D^{Fict}$ is the set of fictitious diamond cells. It follows that the discrete gradients are written as follows

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{int}}^{\tau }{\phi }^{\tau }=\frac{1}{\sin \left({\alpha }_D\right)}\left(\frac{{\phi }_{B^{\star }}-{\phi }_{A^{\star }}}{h_{B^{\star }{A}^{\star }}}{\xi }_{\left[PL\right]}^{B^{\star }}+\frac{{\phi }_L-{\phi }_P}{h_{LP}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P\right)\mathrm{<mo>\; :\; </mo>}\forall {\sigma }^{\star }\in {\mathcal{E}}^{int}$ (31)

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{Fict}}^{\tau }{\phi }^{\tau }=\frac{1}{\sin \left({\alpha }_D\right)}\left(\frac{{\phi }_{B^{\star }}-{\phi }_{A^{\star \perp }}}{h_{B^{\star }{A}^{\star \perp }}}{\xi }_{\left[P^{\perp }{L}^{\perp }\right]}^{B^{\star }}+\frac{{\phi }_{L^{\perp }}-{\phi }_{P^{\perp }}}{h_{P^{\perp }{L}^{\perp }}}{\xi }_{\left[A^{\star \perp }{B}^{\star }\right]}^{P^{\perp }}\right)\mathrm{<mo>\; :\; </mo>}\forall \sigma \mathrm{,}{\sigma }^{\star }\in \Sigma$ (32)

${\nabla }_{{\mathcal{D}}_{\sigma {\sigma }^{\star }}^{Fict}}^{\tau }{\phi }^{\tau }=\frac{1}{\sin \left({\alpha }_D\right)}\left(\frac{{\phi }_{B^{\star }}-{\phi }_{A^{\star }}}{h_{B^{\star }{A}^{\star }}}{\xi }_{\left[PL\right]}^{B^{\star }}+\frac{{\phi }_L-{\phi }_P}{h_{P{L}^{\perp }}}{\xi }_{\left[A^{\star }{B}^{\star }\right]}^P\right)\mathrm{<mo>\; :\; </mo>}\forall {\sigma }^{\star }\in \partial \Omega$ (33)

Remark 2 Note that the set $\mathcal{E}$ primary edges will sometimes be identified with the set $D$ of diamond meshes since there is a trivial bijection between the two sets i.e. $\forall {\sigma }^{\star }\in \mathcal{P}$, $\exists \sigma \in \mathcal{D}$ such that ${\mathcal{D}}_{\sigma {\sigma }^{\star }}^{int}\in D$.

2.3.4. DDFV Scheme of the Model Problem

To obtain a DDFV scheme of the problem (1), we rewrite these fluxes across the primary edges and pseudo-edges calculated previously in the base $\left({\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\right)$ taking in consideration the vectorial relations (6) and (7). By neglecting truncation errors in the Equations (19)-(24), using equally remark (2) and definition (2) we get the following discrete scheme of the problem (1).

Proposition 5 The DDFV scheme associated of the model problem (1)

We can define a DDFV of problem as find ${\bar{\phi }}^{\tau }=\left({\left({\bar{\phi }}_P\right)}_{P\in \mathcal{P}}\mathrm{,}{\left({\bar{\phi }}_{A^{\star }}\right)}_{A^{\star }\in \mathcal{D}}\right)$ such that

$(\begin{array}{ll}-\underset{\mathfrak{D}\in D_P}{{\displaystyle \sum }}{m}_{\sigma }\langle K_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^{\tau }{\bar{\phi }}^{\tau }\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle =meas\left(C_P\right)f_P\hfill & \forall C_P\in \mathcal{P}\hfill \\ -\underset{\mathfrak{D}\in D_{A^{\star }}}{{\displaystyle \sum }}{m}_{{\sigma }^{\star }}\langle K_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^{\tau }{\bar{\phi }}^{\tau }\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle =meas\left(C_{A^{\star }}\right)f_{A^{\star }}\hfill & \forall C_{A^{\star }}\in \mathcal{D}\hfill \end{array}$ (34)

Which can be written taking in consideration remark (2)as follows

$(\begin{array}{ll}-\underset{\sigma \in D_P}{{\displaystyle \sum }}{m}_{\sigma }\langle K_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^{\tau }{\bar{\phi }}^{\tau }\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle =meas\left(C_P\right)f_P\hfill & \forall C_P\in \mathcal{P}\hfill \\ -\underset{{\sigma }^{\star }\in D_{A^{\star }}}{{\displaystyle \sum }}{m}_{{\sigma }^{\star }}\langle K_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^{\tau }{\bar{\phi }}^{\tau }\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle =meas\left(C_{A^{\star }}\right)f_{A^{\star }}\hfill & \forall C_{A^{\star }}\in \mathcal{D}\hfill \end{array}$ (35)

where $f_P$ and $f_{A^{\star }}$ denote the mean value of f over $C_P$ and $C_{A^{\star }}$ respectively,the equivalent diffusion tensor ${\bar{K}}_{\mathfrak{D}}$ satisfies

$\langle {\bar{K}}_{\mathfrak{D}}{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle =\frac{sin\left({\alpha }_D\right)\times m_{\sigma }\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle \langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle }{\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle \cos \left({\theta }_h^{L\mathrm{,}I}\right)h_{IL}+\langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle \cos \left({\theta }_h^{P\mathrm{,}I}\right)h_{IP}}$ (36)

$\langle {\bar{K}}_{\mathfrak{D}}{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle =\frac{\left(\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle \langle K^L{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle h_{DL}+\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle \langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle h_{DP}\right)}{\sin {\left({\alpha }_D\right)}^{-1}\left(\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle \cos \left({\theta }_h^{L\mathrm{,}I}\right)h_{IL}+\langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle \cos \left({\theta }_h^{P\mathrm{,}I}\right)h_{IP}\right)}$ (37)

$\begin{array}{c}\langle {\bar{K}}_{\mathfrak{D}}{\xi }_{\sigma \mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle =\frac{\sin \left({\alpha }_D\right)}{m_{\sigma }}\{\left(\frac{h_{DL}^2}{h_{IL}}\frac{\langle K^L{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle }{\cos \left({\theta }_h^{L,I}\right)}+\frac{h_{DP}^2}{h_{IP}}\frac{\langle K^P{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle }{\cos \left({\theta }_h^{P,I}\right)}\right)\begin{array}{c}\begin{array}{l}\\ \end{array}\\ \\ \\ \end{array}\\ -\frac{1}{h_{IL}{h}_{IP}}{\left(\frac{h_{IL}{h}_{DP}\frac{\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle }{\cos \left({\theta }_h^{P,I}\right)}-h_{IP}{h}_{DL}\frac{\langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle }{\cos \left({\theta }_h^{L,I}\right)}}{\frac{\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle }{\cos \left({\theta }_h^{L,I}\right)}{h}_{IL}+\frac{\langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle }{\cos \left({\theta }_h^{P,I}\right)}{h}_{IP}}\right)}^2\}\end{array}$ (38)

Remark 3 When the points P, I and L are aligned, which means $D=I$, then we have ${\theta }_h^{L,I}={\theta }_h^{P,I}=\mp \left(\frac{\pi }{2}-{\alpha }_D\right)$, $m_{\sigma }=m_{{\sigma }_L}+m_{{\sigma }_P}$,

$\sin \left({\alpha }_D\right)=\cos \left({\theta }_h^{P,I}\right)=\cos \left({\theta }_h^{L,I}\right)$, $m_{{\sigma }_L}=h_{IL}$, $m_{{\sigma }_P}=h_{IP}$, in this case the matrix ${\bar{K}}_{\mathfrak{D}}$ is then defined by:

$\langle {\bar{K}}_{\mathfrak{D}}{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle =\frac{\left(m_{{\sigma }_L}+m_{{\sigma }_P}\right)\langle K^P{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle \langle K^L{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle }{\langle K^P{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_L}+\langle K^L{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_P}}$ (39)

$\langle {\bar{K}}_{\mathfrak{D}}{\xi }_{\sigma .}^P,{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle =\frac{\langle K^P{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle \langle K^L{\xi }_{{\sigma }^{\star }.}^{A^{\star }},{\xi }_{\sigma .}^P\rangle m_{{\sigma }_L}+\langle K^P{\xi }_{\sigma .}^P,{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle \langle K^L{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_P}}{\langle K^P{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_L}+\langle K^L{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_P}}$ (40)

$\begin{array}{c}\langle {\bar{K}}_{\mathfrak{D}}{\xi }_{\sigma .}^{A^{\star }},{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle =\frac{m_{{\sigma }_L}\langle K^L{\xi }_{{\sigma }^{\star }.}^{A^{\star }},{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle +m_{{\sigma }_P}\langle K^P{\xi }_{{\sigma }^{\star }.}^{A^{\star }},{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle }{m_{{\sigma }_L}+m_{{\sigma }_P}}\\ -\frac{m_{{\sigma }_L}\times m_{{\sigma }_P}{\left(\langle K^P{\xi }_{\sigma .}^P,{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle -\langle K^L{\xi }_{\sigma .}^P,{\xi }_{{\sigma }^{\star }.}^{A^{\star }}\rangle \right)}^2}{\langle K^P{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_L}+\langle K^L{\xi }_{\sigma .}^P,{\xi }_{\sigma .}^P\rangle m_{{\sigma }_P}}\end{array}$ (41)

We recognize in (39) the weighted harmonic mean-value of $\langle K^P{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle$ and $\langle K^L{\xi }_{\sigma \mathrm{.}}^P\mathrm{,}{\xi }_{\sigma \mathrm{.}}^P\rangle$ and in the first term of (41) the weighted arithemetic mean-value of $\langle K^L{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle$ and $\langle K^P{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\mathrm{,}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}\rangle$. In this particular case this scheme was already proposed in [9], we also recognize in particular case the median-dual mesh based on the primal centers and the midpoint of the edges proposed in [4].

Proposition 6 Discrete integration by parts formula associated to the model problem (1) For all ${\bar{\phi }}^{\tau },{\bar{\varphi }}^{\tau }\in {ℝ}^{\tau }$ we have

$\begin{array}{l}\underset{\mathfrak{D}\in D}{{\displaystyle \sum }}meas\left(\mathfrak{D}\right)\langle K_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^{\tau }{\bar{\phi }}^{\tau }\mathrm{,}{\nabla }_{\mathfrak{D}}^{\tau }{\bar{\varphi }}^{\tau }\rangle \\ =\underset{C_P\in \mathcal{P}}{{\displaystyle \sum }}meas\left(C_P\right){\varphi }_P{f}_P+\underset{C_{A^{\star }}\in \mathcal{D}}{{\displaystyle \sum }}meas\left(C_{A^{\star }}\right){\varphi }_{A^{\star }}{f}_{A^{\star }}\end{array}$

Proof. Since $m_{{\sigma }^{\star }}{\xi }_{{\sigma }^{\star }\mathrm{.}}^{A^{\star }}=h_{IL}{\xi }_{\left[IL\right]\mathrm{.}}^{A^{\star }}+h_{IP}{\xi }_{\left[IP\right]\mathrm{.}}^{A^{\star }}$ the proof is similar to the one exposed in [4].

3. Existence of Discrete Solutions of the Model Problem

3.1. Existence and Uniqueness of DDFV Scheme (Case of Dirichlet Boundary Conditions)

Let us check the existence and uniqueness of the discrete solution.

Proposition 7 The matrix associated with the linear system (34)-(31)-(26) is symmetric and positive definite

Proof. It is easily seen that the symmetry of the matrix M associated with the system (26), (31) and (34) essentially follows from the symmetry of the diffusion coefficient K. We should now prove the positive definiteness of that matrix.

Development of the quadratic expression:

$\begin{array}{c}\langle M{\phi }^{\mathcal{T}}\mathrm{,}{\phi }^{\mathcal{T}}\rangle =2\underset{\mathfrak{D}\in D}{{\displaystyle \sum }}meas\left(\mathfrak{D}\right)\langle {\bar{K}}_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^{\mathcal{T}}{\phi }^{\mathcal{T}}\mathrm{,}{\nabla }_{\mathfrak{D}}^{\mathcal{T}}{\phi }^{\mathcal{T}}\rangle \\ =2\underset{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D}{{\displaystyle \sum }}meas\left({\mathfrak{D}}_{\sigma {\sigma }^{\star }}\right)\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\nabla }_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}^{\mathcal{T}}{\phi }^{\mathcal{T}}\mathrm{,}{\nabla }_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}^{\mathcal{T}}{\phi }^{\mathcal{T}}\rangle \\ =\underset{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D}{{\displaystyle \sum }}\{\frac{mes\left({\sigma }^{\star }\right)}{mes\left(\sigma \right)\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle {\left({\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right)}^2\\ +\frac{2\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle }{\sin \left({\alpha }_D\right)}\left({\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right)\left({\phi }_L-{\phi }_P\right)\\ +\frac{mes\left(\sigma \right)}{mes\left({\sigma }^{\star }\right)\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^P\rangle {\left({\phi }_L-{\phi }_P\right)}^2\}\end{array}$

This quadratic form is elliptic if and only if:

$\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle \langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^P\rangle -{\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle }^2\succ 0$ (42)

In effect we have:

$\begin{array}{c}{\Delta }^{{\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}}=\frac{mes\left(\sigma \right)}{mes\left({\sigma }^{\star }\right)\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^P\rangle \times \frac{mes\left({\sigma }^{\star }\right)}{mes\left(\sigma \right)\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle \\ -{\left[\frac{1}{\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle \right]}^2\end{array}$

$\begin{array}{c}={\left(\frac{1}{\sin \left({\alpha }_D\right)}\right)}^2\left[\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle \langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^P\rangle -{\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle }^2\right]\\ =N_1\times {\Delta }^{K^P}+N_2\times {\Delta }^{K^L}\end{array}$

where $N_1$ and $N_2$ are strictly positive numbers defined as

$\begin{array}{c}{N}_1={\left[\frac{h_{PD}\langle K^L{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle }{\langle K^P{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle \cos \left({\theta }_h^{L,I}\right)h_{IL}+\langle K^L{x}_s^P,x_s^P\rangle \cos \left(q_h^{P,I}\right)h_{IP}}\right]}^2\\ +\frac{\langle K^L{\xi }_{\sigma }^{A^{\star }},{\xi }_{\sigma }^{A^{\star }}\rangle \langle K^L{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle h_{PD}^2\times \cos \left({\theta }_h^{L,I}\right)h_{IL}}{h_{IP}\cos \left({\theta }_h^{P,I}\right){\left(\langle K^P{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle \cos \left({\theta }_h^{L,I}\right)h_{IL}+\langle K^L{x}_s^P,x_s^P\rangle \cos \left(q_h^{P,I}\right)h_{IP}\right)}^2}\end{array}$

$\begin{array}{c}{N}_2={\left[\frac{h_{LD}\langle K^P{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle }{\langle K^P{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle \cos \left({\theta }_h^{L,I}\right)h_{IL}+\langle K^L{x}_s^P,x_s^P\rangle \cos \left(q_h^{P,I}\right)h_{IP}}\right]}^2\\ +\frac{\langle K^L{\xi }_{\sigma }^{A^{\star }},{\xi }_{\sigma }^{A^{\star }}\rangle \langle K^L{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle h_{LD}^2\times \cos \left({\theta }_h^{P,I}\right)h_{IP}}{h_{IL}\cos \left({\theta }_h^{L,I}\right){\left(\langle K^P{\xi }_{\sigma }^P,{\xi }_{\sigma }^P\rangle \cos \left({\theta }_h^{L,I}\right)h_{IL}+\langle K^L{x}_s^P,x_s^P\rangle \cos \left(q_h^{P,I}\right)h_{IP}\right)}^2}\end{array}$

Since the diffusion matrix K is symmetric and positive definite, the Cauchy-Schwartz inequality for the inner product associated with K ensures that

${\Delta }^{K^P}=\langle K^P{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle \langle K^P{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^P\rangle -{\langle K^P{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle }^2\succ 0$

and

${\Delta }^{K^L}=\langle K^L{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle \langle K^L{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^P\rangle -{\langle K^L{\xi }_{\sigma }^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle }^2\succ 0$

as either ${\xi }_{\sigma }^P$ and ${\xi }_{\sigma }^{A^{\star }}$ are not collinear. Therefore ${\Delta }^{{\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}}\succ 0$, thus the matrix ${\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}$ are symmetric and positive definite matrices. Under the assumption 1 the matrix ${\bar{K}}_{\mathfrak{D}}$ possesses strictly positive eigenvalues. Let ${\delta }_{\min }^{\mathfrak{D}}$ be its least eigenvalue. so we have

${\delta }_{\min }^{\mathfrak{D}}\underset{\mathfrak{D}\in D}{{\displaystyle \sum }}meas\left(\mathfrak{D}\right){\Vert {\nabla }_{\mathfrak{D}}^{\tau }{\bar{\phi }}^{\tau }\Vert }^2\le \frac{1}{2}\langle M{\bar{\phi }}^{\tau },{\bar{\phi }}^{\tau }\rangle$ (43)

the equality holds in Equation (43) if and only if ${\phi }^{\mathcal{T}}=0$. Thus, the positive definiteness of the matrix M is proven.

Proposition 8 The linear system (26), (31) and (34) possesses a unique solution.

Proof. It follows from the previous proposition.

3.2. Existence and Uniqueness of DDFV Scheme (Case of Periodic Boundary Conditions)

Remark 4 The proof of existence and uniqueness of discrete solutions with Neumann conditions is very similar with that of periodic boundary conditions. For Neumann case we refer to [10] for more details

Let us start this subsection with some preliminary remarks and results. First of all, we assume that the cell-points and the vertex-points from the primary grid are numbered. The numbering is performed in such a way that any boundary-vertex (different from a corner-point) and its orthogonal projection on the opposite side get the same number. On the other hand, the four corner-points $O_1^{\mathrm{<mo>\; *\; </mo>}}$, $O_2^{\mathrm{<mo>\; *\; </mo>}}$, $O_3^{\mathrm{<mo>\; *\; </mo>}}$ and $O_4^{\mathrm{<mo>\; *\; </mo>}}$ are given the same number. This way of numbering accounts with the periodicity setting of the discrete problem.

Remark 5 Note that to obtain a square system, due to the periodicity we will only compute the fictitious diamond cells located at the northern and eastern boundary or western and southern boundary

Proposition 9 Let $M_h$ be the total number of discrete unknowns according to the previous numbering of cell-points and vertex-points. Set that: ${\mathcal{D}}_{\ast }={\mathcal{D}}^{int}\cup \mathcal{D}\left({\Gamma }_{So}\right)\cup \mathcal{D}\left({\Gamma }_{We}\right)\cup \left\{D_{O_1^{\star }}\right\}$, where ${\mathcal{D}}^{int}$ is the set that consists of dual cells associated with interior vertexpoints and where $\mathcal{D}\left({\Gamma }_{So}\right)$ and $\mathcal{D}\left({\Gamma }_{We}\right)$ respectively denote the sets of dual cells associated with boundary vertexpoints different from cornerpoints and located in the $x_1$ -axis and the $x_2$ -axis. Of course $D_{O_1^{\star }}$ is the dual cell associated with $O_1^{\star }$ . The two following discrete problems are equivalent:

(DP₁): Find $\left({\left\{{\bar{\phi }}_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{{\bar{\phi }}_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}\right)\in {ℝ}^{\mathcal{P}}\times {ℝ}^{\mathcal{D}}$ such that Equations (31),(32),(33)and (34)are satisfied

And

(DP₂): Find $\left({\left\{{\bar{\phi }}_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{{\bar{\phi }}_{D^{\star }}\right\}}_{D^{\star }\in {\mathcal{D}}_{\mathrm{<mo>\; *\; </mo>}}}\right)\in {ℝ}^{M_h}\equiv {ℝ}^{\mathcal{P}}\times {ℝ}^{{\mathcal{D}}_{\mathrm{<mo>\; *\; </mo>}}}$ such that Equations (31) and (32) and (33)and (34)are satisfied

Where, for a given set of mesh elements $\mathcal{K}$, one has set:

${ℝ}^{\mathcal{K}}=\left\{{\left\{v_K\right\}}_{K\in \mathcal{K}}\mathrm{<mo>\; ;\; </mo>}\forall K\in \mathcal{K},v_K\in ℝ\right\}$. (44)

For the sake of clarity of the exposition, $P\in \mathcal{P}$ will be used either for denoting a cellpoint or its associated number, idem for $D^{\star }\in \mathcal{V}\equiv \mathcal{D}$ concerning primal vertex-points.

3.2.1. Some Useful Vector Spaces and Preliminary Results

In view to algebraically address the linear discrete system (31) - (34), we introduce some adequate vector spaces. Recall that $\mathcal{A}$ denotes the auxiliary mesh i.e. an intermediate mesh between the primal mesh $\mathcal{P}$ and the dual mesh $\mathcal{D}$. The main feature of $\mathcal{A}$ is that each auxiliary cell involves either one only cellpoint or one only vertexpoint. Consequently, we have the following relation between $\mathcal{A}$, $\mathcal{P}$ and $\mathcal{D}$ :

$\mathcal{A}={\left\{{\Omega }_P\right\}}_{P\in \mathcal{P}}\cup {\left\{{\Omega }_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}$ (45)

where ${\Omega }_P$ and ${\Omega }_{D^{\star }}$ denote the two kinds of auxiliary cells emerging from the definition of the mesh $\mathcal{A}$. These cells will play a key role for the definition of DDFV solutions to the system of Equations (1) - (10) or (1) - (11). Let us introduce the following vector spaces:

${ℝ}^{\mathcal{A}}=\left\{\left({\left\{v_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{v_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}\right)\mathrm{<mo>\; ;\; </mo>}\forall P\in \mathcal{P}\mathrm{,}\forall D^{\star }\in \mathcal{D}\mathrm{,}{v}_P\mathrm{,}{v}_{D^{\star }}\in ℝ\right\}$ (46)

$\begin{array}{l}{ℝ}_{Perio}^{\mathcal{A}}=\{\left(\left\{v_P\right\}\mathrm{<mo>\; ;\; </mo>}\left\{v_{D^{\star }}\right\}\right)\in {ℝ}^{\mathcal{A}}\mathrm{<mo>\; ;\; </mo>}\forall D^{\star }\in \mathcal{D}\left({\Gamma }_{So}\right)\cup \mathcal{D}\left({\Gamma }_{We}\right),v_{D^{\star }}=v_{D^{\star \perp }}\\ \text{and}{v}_{O_1^{\star }}=v_{O_2^{\star }}=v_{O_3^{\star }}=v_{O_4^{\star }}\}\end{array}$ (47)

Recall that $\mathcal{D}\left({\Gamma }_{So}\right)$ and $\mathcal{D}\left({\Gamma }_{We}\right)$ respectively denote the sets of boundary vertex-points different from corner-points and located in the $x_1$ -axis and the $x_2$ -axis. In the sequel, we will sometimes do the following identification:

${ℝ}^{\mathcal{A}}\equiv {ℝ}^{\mathcal{P}}\times {ℝ}^{\mathcal{D}}\mathrm{.}$

Note that ${ℝ}_{Perio}^{\mathcal{A}}$ is an $M_h$ -dimensional subspace of ${ℝ}^{\mathcal{A}}$. So, the following identification will be sometimes considered:

${ℝ}_{Perio}^{\mathcal{A}}\equiv {ℝ}^{M_h}\equiv {ℝ}^{\mathcal{P}}\times {ℝ}^{{\mathcal{D}}_{\ast }}\mathrm{.}$

Proposition 10 The matrix ${\mathbb{M}}_h$ associated with the discrete system (DP₂) is singular.

Proof. The proof is elementary. Indeed it suffices to remark from the discrete system (31) - (34) that

${\displaystyle \underset{1\le j\le M_h}{\sum }}{\left({\mathbb{M}}_h\right)}_{ij}=0\forall 1\le i\le M_h\mathrm{.}$

Thus, the kernel $Ker\left({\mathbb{M}}_h\right)$ of ${\mathbb{M}}_h$ defined by

$Ker\left({\mathbb{M}}_h\right)=\left\{V_h\in {ℝ}^{M_h}\mathrm{<mo>\; ;\; </mo>}{\mathbb{M}}_h{V}_h=0\right\}$

involves a nonzero vector, namely ${\mathbb{I}}_h=\left(\mathrm{1,1,}\cdots \mathrm{,1}\right)\in {ℝ}^{M_h}$, and therefore ${\mathbb{M}}_h$ is singular.

We know from the previous proposition that the discrete problem (DP₂) could get either no solution or an infinite number of solutions. Indeed, existence of solutions to this problem depends on whether the right-hand side to (DP₂) is in the orthogonal of the kernel of ${\mathbb{M}}_h$. Our purpose now is to give a characterization of $Ker\left({\mathbb{M}}_h\right)$, the kernel of ${\mathbb{M}}_h$. Before that, we recall a result we need for proving the characterization of $Ker\left({\mathbb{M}}_h\right)$.

Proposition 11 (Characterization of $Ker\left({\mathbb{M}}_h\right)$ )

• $Ker\left({\mathbb{M}}_h\right)=\left\{V_h\in {ℝ}^{M_h}\equiv {ℝ}_{Perio}^{\mathcal{A}}\mathrm{<mo>\; ;\; </mo>}{\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0\right\}$ ;

• ${\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0\iff V_h={\sigma }_{\mathcal{P}}{\mathbb{I}}_{\mathcal{P}}+{\sigma }_{{\mathcal{D}}_{\ast }}{\mathbb{I}}_{{\mathcal{D}}_{\ast }}\forall V_h\in {ℝ}^{M_h}\equiv {ℝ}_{Perio}^{\mathcal{A}}$,

where we have set:

${\left({\mathbb{I}}_{\mathcal{P}}\right)}_i=(\begin{array}{cc}1& \text{if}i\in \mathcal{P}\\ 0& \text{otherwise}\end{array}\text{and}{\left({\mathbb{I}}_{{\mathcal{D}}_{\ast }}\right)}_i=(\begin{array}{cc}1& \text{if}i\in {\mathcal{D}}_{\ast }\\ 0& \text{otherwise}\text{.}\end{array}$

and where ${\sigma }_{\mathcal{P}}$ and ${\sigma }_{{\mathcal{D}}_{\ast }}$ are given real numbers.

Proof. The inclusion of $Ker\left({\mathbb{M}}_h\right)$ in $\left\{V_h\in {ℝ}^{M_h}\mathrm{<mo>\; ;\; </mo>}{\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0\right\}$ is trivial. Let us concentrate on the proof of the inclusion of $\left\{V_h\in {ℝ}^{M_h}\mathrm{<mo>\; ;\; </mo>}{\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0\right\}$ in $Ker\left({\mathbb{M}}_h\right)$. For this purpose, let $V_h=\left({\left\{V_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{V_{D^{\star }}\right\}}_{D^{\star }\in {\mathcal{D}}_{\ast }}\right)\in {ℝ}^{M_h}\equiv {ℝ}^{\mathcal{P}\cup {\mathcal{D}}_{\ast }}$ be a vector such that

${\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0.$

where one has set

${\mathcal{D}}_{\ast }={\mathcal{D}}^{int}\cup \mathcal{D}\left({\Gamma }_{So}\right)\cup \mathcal{D}\left({\Gamma }_{We}\right)\cup \left\{D_{O_1^{\star }}\right\}\mathrm{.}$

Therefore, it follows from the integration by parts formula that:

$\begin{array}{c}\langle {\mathbb{M}}_h{\phi }^h,j^h\rangle =2\underset{\mathfrak{D}\in D_{\star }}{{\displaystyle \sum }}meas\left(\mathfrak{D}\right)\langle {\bar{K}}_{\mathfrak{D}}{\nabla }_{\mathfrak{D}}^h{\phi }^h\mathrm{,}{\nabla }_{\mathfrak{D}}^h{\phi }^h\rangle \\ =2\underset{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D_{\star }}{{\displaystyle \sum }}meas\left({\mathfrak{D}}_{\sigma {\sigma }^{\star }}\right)\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\nabla }_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}^h{\phi }^h\mathrm{,}{\nabla }_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}^h{\phi }^h\rangle \\ =\underset{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D_{\star }}{{\displaystyle \sum }}\{\frac{meas\left({\sigma }^{\star }\right)}{meas\left(\sigma \right)\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{\sigma }^{A^{\star }}\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle {\left({\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right)}^2\\ +\frac{2\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{{\sigma }^{\star }}^P\mathrm{,}{\xi }_{\sigma }^{A^{\star }}\rangle }{\sin \left({\alpha }_D\right)}\left({\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right)\left({\phi }_L-{\phi }_P\right)\\ +\frac{meas\left(\sigma \right)}{meas\left({\sigma }^{\star }\right)\sin \left({\alpha }_D\right)}\langle {\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}{\xi }_{{\sigma }^{\star }}^P\mathrm{,}{\xi }_{{\sigma }^{\star }}^P\rangle {\left({\phi }_L-{\phi }_P\right)}^2\}\end{array}$ (48)

According to [6], there exists ${\lambda }_{\min }$ a mesh-independent, strictly positive real number that minimizes both of least eigenvalues of ${\bar{K}}_{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}}$. We then deduce from Equation (48) that

${\lambda }_{\min }\underset{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in D_{\star }}{{\displaystyle \sum }}\left\{{\left({\phi }_{B^{\star }}-{\phi }_{A^{\star }}\right)}^2+{\left({\phi }_L-{\phi }_P\right)}^2\right\}\le \langle {\mathbb{M}}_h{\phi }^h,j^h\rangle$ (49)

It is obvious from the previous inequality that the space $\left\{V_h\in {ℝ}^{M_h}\mathrm{<mo>\; ;\; </mo>}{\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0\right\}$ is included in the one spanned by the vectors ${\mathbb{I}}_{\mathcal{P}}$ and ${\mathbb{I}}_{{\mathcal{D}}_{\ast }}$ defined by

${\left({\mathbb{I}}_{\mathcal{P}}\right)}_i=(\begin{array}{cc}1& \text{if}i\in \mathcal{P}\text{i}\text{.e}\text{.}i\text{is}\text{a}\text{cellpoint}\text{number}\\ 0& \text{otherwise}\end{array}$

and

${\left({\mathbb{I}}_{{\mathcal{D}}_{\ast }}\right)}_i=(\begin{array}{cc}1& \text{if}i\in {\mathcal{D}}_{\ast }\\ 0& \text{otherwise}\text{.}\end{array}$

Remarking that ${\mathbb{I}}_{\mathcal{P}}$ and ${\mathbb{I}}_{{\mathcal{D}}_{\ast }}$ are two (linearly independent) eigenvectors with zero as corresponding eigenvalue with respect to the matrix ${\mathbb{M}}_h$, the proof of Proposition 11 is completed.

Remark 6 Note that the dimension of the space $Ker\left({\mathbb{M}}_h\right)$ is equal to 2. This information plays a key-role for uniqueness conditions investigated in the next section.

3.2.2. Existence and Conditions for Uniqueness of a Solution to (DP₂)

Let us start with an existence result for the discrete problem (DP₂) which involves the system of Equations (31)-(34).

Proposition 12 (Existence)

Under the assumption (12), the discrete system (DP₂) possesses an infinite number of solutions.

Proof. Just remark that the condition (12) ensures that the right-hand side of the system of Equations (34) is a vector of ${ℝ}^{M_h}$ orthogonal to $Ker\left({\mathbb{M}}_h\right)$.

Recall that for the continuous problem (1), the assumptions (2) - (12) ensure the existence of a family of variational solutions i.e. set of functions $\phi$ living in the space

$\mathcal{H}=\left\{v\in H^1\left(\Omega \right);v\left(.,0\right)=v\left(.,b\right)\text{in}\left[0,a\right]\text{and}v\left(0,.\right)=v\left(a,.\right)\text{in}\left[0,b\right]\right\}\mathrm{,}$

and such that

${\displaystyle {\int }_{\Omega }}K\left(x\right)grad\phi \left(x\right)\cdot gradv\left(x\right)\text{d}x={\displaystyle {\int }_{\Omega }}f\left(x\right)v\left(x\right)\text{d}x\forall v\in \mathcal{H}\mathrm{.}$

The uniqueness of a solution to this continuous problem is got from the subspace of $\mathcal{H}$ made up of functions v satisfying the following condition:

${\displaystyle {\int }_{\Omega }}v\left(x\right)\text{d}x=0.$ (50)

This condition makes obvious the necessity of associating a discrete function $v_h$ (defined almost everywhere in $\Omega$ ) with any vector $\left({\left\{v_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{v_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}\right)$ from ${ℝ}_{perio}^{\mathcal{A}}$.

• Basic discrete function spaces. Denote by $E\left({ℝ}_{perio}^{\mathcal{A}}\right)$ the space of such discrete functions $v_h$ and let ${\mathbb{I}}_S$ be the characteristic function of a subset $S$ of $\Omega$ i.e. ${\mathbb{I}}_S\left(x\right)=1$ if $x\in S$ and 0 otherwise. We define the space $E\left({ℝ}_{perio}^{\mathcal{A}}\right)$ as:

$\begin{array}{l}E\left({ℝ}_{Perio}^{\mathcal{A}}\right)=\{v\mathrm{<mo>\; :\; </mo>}\Omega \to ℝ\mathrm{<mo>\; ;\; </mo>}\exists \left({\left\{v_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{v_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}\right)\in {ℝ}_{Perio}^{\mathcal{A}}\text{such}\text{that}\begin{array}{c}\\ \end{array}\\ v\left(x\right)=\underset{P\in \mathcal{P}}{{\displaystyle \sum }}{v}_P{\mathbb{I}}_{{\Omega }_P}\left(x\right)+\underset{D^{\star }\in \mathcal{D}}{{\displaystyle \sum }}{v}_{D^{\star }}{\mathbb{I}}_{{\Omega }_{D^{\star }}}\left(x\right)\text{a}\text{.e}\mathrm{.}\text{in}\Omega \}\end{array}$ (51)

where ${\Omega }_P$ and ${\Omega }_{D^{\star }}$ are two generic notations of the two kinds of auxiliary cells from $\mathcal{A}$. Let us introduce the following mappings (viewed as projections in some sense):

$v_h\in E\left({ℝ}_{Perio}^{\mathcal{A}}\right)\mapsto {\Pi }^{\mathcal{P}}\left[v_h\right]\left(x\right)\equiv v_h^{\mathcal{P}}\left(x\right)=\underset{P\in \mathcal{P}}{{\displaystyle \sum }}{v}_P{\mathbb{I}}_{C_P}\left(x\right)\in E\left({ℝ}^{\mathcal{P}}\right)$ (52)

$v_h\in E\left({ℝ}_{Perio}^{\mathcal{A}}\right)\mapsto {\Pi }^{\mathcal{D}}\left[v_h\right]\left(x\right)\equiv v_h^{\mathcal{D}}\left(x\right)=\underset{D^{\star }\in \mathcal{D}}{{\displaystyle \sum }}{v}_{D^{\star }}{\mathbb{I}}_{C_{D^{\star }}}\left(x\right)\in E\left({ℝ}^{\mathcal{D}}\right)$ (53)

where $E\left({ℝ}^{\mathcal{P}}\right)$ and $E\left({ℝ}^{\mathcal{D}}\right)$ are function spaces respectively associated with vector spaces ${ℝ}^{\mathcal{P}}$ and ${ℝ}^{\mathcal{D}}$ (see relation (44) for the definition of these vector spaces).

• Looking for uniqueness conditions. Since the dimension of $Ker\left({\mathbb{M}}_h\right)$ is equal to 2, we should look for two discrete analogues of (50) (namely quadrature formulas), linearly independent, that should ensure uniqueness for a solution to the discrete problem (DP₂). These discrete analogues are defined as:

${\displaystyle {\int }_{\Omega }}{v}_h^{\mathcal{P}}\left(x\right)\text{d}x=0\text{and}{\displaystyle {\int }_{\Omega }}{v}_h^{\mathcal{D}}\left(x\right)\text{d}x=0$

In other words,

$\underset{P\in \mathcal{P}}{{\displaystyle \sum }}mes\left(C_P\right)v_P=0\text{and}\underset{D^{\star }\in \mathcal{D}}{{\displaystyle \sum }}mes\left(C_{D^{\star }}\right)v_{D^{\star }}=0$

Proposition 13 Define two quadrature expressions ${\mathbb{Q}}^{\mathcal{P}}\left(\mathrm{.}\right)$ and ${\mathbb{Q}}^{\mathcal{D}}\left(\mathrm{.}\right)$ on the discrete function space $E\left({ℝ}^{\mathcal{A}}\right)$ as follows: for all

$v_h=\left(\left\{v_P\right\}\mathrm{<mo>\; ;\; </mo>}\left\{v_{D^{\star }}\right\}\right)\in E(\; ℝ\; A\; )$

${\mathbb{Q}}^{\mathcal{P}}\left(v_h\right)=\underset{P\in \mathcal{P}}{{\displaystyle \sum }}mes\left(C_P\right)v_P\text{and}{\mathbb{Q}}^{\mathcal{D}}\left(v_h\right)=\underset{D^{\star }\in \mathcal{D}}{{\displaystyle \sum }}mes\left(C_{D^{\star }}\right)v_{D^{\star }}\mathrm{.}$

Then, ${\mathbb{Q}}^{\mathcal{P}}\left(\mathrm{.}\right)$ and ${\mathbb{Q}}^{\mathcal{D}}\left(\mathrm{.}\right)$ are linearly independent linear forms on $E\left({ℝ}^{\mathcal{A}}\right)$. $\square$

Proof. Let ${\lambda }_{\mathcal{P}}$ and ${\lambda }_{\mathcal{D}}$ be two real numbers such that ${\lambda }_{\mathcal{P}}{\mathbb{Q}}^{\mathcal{P}}\left(\mathrm{.}\right)+{\lambda }_{\mathcal{D}}{\mathbb{Q}}^{\mathcal{D}}\left(\mathrm{.}\right)=0$ on $E\left({ℝ}^{\mathcal{A}}\right)$, that is, ${\lambda }_{\mathcal{P}}{\mathbb{Q}}^{\mathcal{P}}\left(v_h\right)+{\lambda }_{\mathcal{D}}{\mathbb{Q}}^{\mathcal{D}}\left(v_h\right)=0\forall v_h\in E\left({ℝ}^{\mathcal{A}}\right)\mathrm{.}$

Taking $v_h$ to be successively $\underset{P\in \mathcal{P}}{{\displaystyle \sum }}{\mathbb{I}}_{{\Omega }_P}$ and $\underset{D^{\star }\in \mathcal{D}}{{\displaystyle \sum }}{\mathbb{I}}_{{\Omega }_{D^{\star }}}$, one easily see that ${\lambda }_{\mathcal{P}}={\lambda }_{\mathcal{Q}}=0$. This proves the proposition.

We have gathered all the ingredients for existence and uniqueness of a solution to the discrete problem (DP₂).

Proposition 14 (Existence and Uniqueness) The problem that consists in finding $\left({\left\{{\bar{\phi }}_P\right\}}_{P\in \mathcal{P}}\mathrm{<mo>\; ;\; </mo>}{\left\{{\bar{\phi }}_{D^{\star }}\right\}}_{D^{\star }\in \mathcal{D}}\right)\in {ℝ}_{Perio}^{\mathcal{A}}$ such that the Equations (31) - (34) are satisfied possesses a unique solution if condition (12) is fulfilled and if

$\underset{P\in \mathcal{P}}{{\displaystyle \sum }}mes\left(C_P\right){\bar{\phi }}_P=\underset{D^{\star }\in \mathcal{D}}{{\displaystyle \sum }}mes\left(C_{D^{\star }}\right){\bar{\phi }}_{D^{\star }}=0.$ (54)

Proof. We already know from Proposition 12 that under the condition (12), the discrete system (31)-(34) possesses an infinite number of solutions in ${ℝ}_{Perio}^{\mathcal{A}}$. Moreover, in the proof of Proposition 11, we have shown the fact that

${\left(V_h\right)}^t{\mathbb{M}}_h{V}_h\ge 0\forall V_h\in {ℝ}_{Perio}^{\mathcal{A}}$ (55)

So, to end the proof of Proposition 14, we just need to prove the positive definiteness of ${\mathbb{M}}_h$ over the subspace ${ℝ}_{Perio\mathrm{,0}}^{\mathcal{A}}$ of ${ℝ}_{Perio}^{\mathcal{A}}$ made up of vectors $V_h$ satisfying the conditions (54). For this purpose, it suffices to show the assertion:

${\left(V_h\right)}^t{\mathbb{M}}_h{V}_h=0\Rightarrow V_h=0\forall V_h\in {ℝ}_{Perio\mathrm{,0}}^{\mathcal{A}}\mathrm{.}$

The second part of Proposition 11 lets this assertion be trivial.

4. A Short Survey of the 2D Implementation of DDFV Scheme on Uniform Mesh

4.1. DDFV Formulation on the Internal Meshes

Let us consider the model problem (1), in two dimensional region $\Omega =\left(\text{0,1}\right)\times \left(\text{0,1}\right)$ (where the permeability D is used instead of K). Before starting the implementation step, let us introduce the controle volume (Figure 5) and the main notations. We consider a primal mesh composed of rectangular cells

$K_{i;j}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]\times \left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right]$, $\forall i;j\in \left\{1,2,\cdots ,N\right\}$. For the sake of simplicity,

Figure 5. Controle volume of DDFV scheme up (primal) and down (dual).

we will assume that the mesh is uniform, and we enforce

$x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}=y_{i+\frac{1}{2}}-y_{i-\frac{1}{2}}=h$ for all $i;j\in \left\{1,2,\cdots ,n\right\}$ where $h=\frac{1}{N}\succ 0$ is fixed. we denoted by $\left(x_{i+\frac{\text{1}}{\text{2}}}\mathrm{,}{y}_{j+\frac{\text{1}}{\text{2}}}\right)$ the vertices of the primal mesh and by $\left(x_i\mathrm{,}{y}_j\right)$ the center of the primal cell $K_{i;j}$. Around each vertex of the primal mesh $\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)$, we construct a dual cell $K_{i+\frac{1}{2};j+\frac{1}{2}}=\left[x_i,x_{i+1}\right]\times \left[y_j,y_{j+1}\right]$. The set of the dual cells ${\left(K_{i+\frac{1}{2};j+\frac{1}{2}}\right)}_{i,j\in \left\{1,2,\cdots ,N-1\right\}}$ constitutes a second mesh wish will call dual mesh. The diamond meshes are defined by the vertices $\left(x_i,y_{j+1}\right),\left(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}}\right),\left(x_i,y_j\right),\left(x_{i+\frac{\text{1}}{\text{2}}}\mathrm{,}{y}_{j+\frac{\text{1}}{\text{2}}}\right)$ or $\left(x_{i-\text{1}}\mathrm{,}{y}_j\right)\mathrm{,}\left(x_{i+\frac{\text{1}}{\text{2}}}\mathrm{,}{y}_{j+\frac{\text{1}}{\text{2}}}\right)\mathrm{,}\left(x_i\mathrm{,}{y}_j\right)\mathrm{,}\left(x_{i+\frac{\text{1}}{\text{2}}}\mathrm{,}{y}_{j-\frac{\text{1}}{\text{2}}}\right)$ and denoted respectively by ${\mathcal{D}}_{i+\frac{1}{2}\mathrm{,}j}$ or ${\mathcal{D}}_{i\mathrm{,}j+\frac{1}{2}}$. As you can see in Figure 6, the centers of the dual cells are the vertices of the primal mesh and conversely. The set $\Delta h={\left\{\left(x_i,y_{j+1}\right),\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)\right\}}_{i,j\in \left\{0,1,\cdots ,N\right\}}$ made of finite number of nodes is called the computational grid. Recall that:

$x_i=\frac{x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}}}{2},y_j=\frac{y_{j+\frac{1}{2}}+y_{j-\frac{1}{2}}}{2},\forall i,j\in \left\{1,2,\cdots ,N\right\}$

We also adopt the following conventions: $x_0=x_{\frac{1}{2}}=y_0=y_{\frac{1}{2}}=0$ and.

$x_{N+1}=x_{N+\frac{1}{2}}=y_{N+1}=y_{N+\frac{1}{2}}=1$ We are looking for values ${\phi }_{i,j},{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}$ which approximate $\phi \left(x_i\mathrm{,}{y}_j\right)$ and $\phi \left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)$ respectevely. By rewriting Formula

(34) in this particular case, the DDFV formulation of the model problem (1) on the internal meshes is written as follows

Figure 6. The natural ordering nodes for $N=4$ : Dirichlet boundary conditions. (LEFT), for Neumann boundary conditions (middle) and for periodic boundary (Right).

$\begin{array}{l}\frac{2D_{22}^{i,j}{D}_{22}^{i,j+1}}{D_{22}^{i,j}+D_{22}^{i,j+1}}\left({\phi }_{i,j}-{\phi }_{i,j+1}\right)+\frac{D_{22}^{i,j}{D}_{21}^{i,j+1}+D_{21}^{i,j}{D}_{22}^{i,j+1}}{D_{22}^{i,j}+D_{22}^{i,j+1}}\left({\phi }_{i-\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right)\\ +\frac{2D_{22}^{i,j}{D}_{22}^{i,j-1}}{D_{22}^{i,j}+D_{22}^{i,j-1}}\left({\phi }_{i,j}-{\phi }_{i,j-1}\right)+\frac{D_{22}^{i,j}{D}_{21}^{i,j-1}+D_{21}^{i,j}{D}_{22}^{i,j-1}}{D_{22}^{i,j}+D_{22}^{i,j-1}}\left({\phi }_{i+\frac{1}{2},j-\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j-\frac{1}{2}}\right)\\ +\frac{2D_{11}^{i,j}{D}_{11}^{i+1,j}}{D_{11}^{i,j}+D_{11}^{i+1,j}}\left({\phi }_{i,j}-{\phi }_{i+1,j}\right)+\frac{D_{11}^{i,j}{D}_{12}^{i+1,j}+D_{12}^{i,j}{D}_{11}^{i+1,j}}{D_{11}^{i,j}+D_{11}^{i+1,j}}\left({\phi }_{i+\frac{1}{2},j-\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right)\\ +\frac{2D_{11}^{i,j}{D}_{11}^{i-1,j}}{D_{11}^{i,j}+D_{11}^{i-1,j}}\left({\phi }_{i,j}-{\phi }_{i-1,j}\right)+\frac{D_{11}^{i,j}{D}_{12}^{i-1,j}+D_{12}^{i,j}{D}_{11}^{i-1,j}}{D_{11}^{i,j}+D_{11}^{i+1,j}}\left({\phi }_{i-\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j-\frac{1}{2}}\right)\\ ={\displaystyle {\int }_{K_{i,j}}}f\text{d}x\forall 1\le i,j\le N\end{array}$ (56)

and

$\begin{array}{l}\left[\frac{\left(D_{12}^{i+1,j+1}-D_{12}^{i+1,j+1}\right)\left(D_{21}^{i,j+1}-D_{21}^{i+1,j+1}\right)}{2\left(D_{11}^{i,j+1}+D_{11}^{i+1,j+1}\right)}+\frac{D_{22}^{i,j+1}+D_{22}^{i+1,j+1}}{2}\right]\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{3}{2}}\right)\\ +\left[\frac{D_{11}^{i,j+1}{D}_{21}^{i+1,j+1}+D_{11}^{i+1,j+1}{D}_{21}^{i,j+1}}{D_{11}^{i,j+1}+D_{11}^{i+1,j+1}}\right]\left({\phi }_{i,j+1}-{\phi }_{i+1,j+1}\right)\\ +\left[\frac{\left(D_{12}^{i+1,j}-D_{12}^{i,j}\right)\left(D_{21}^{i,j}-D_{21}^{i+1,j}\right)}{2\left(D_{11}^{i,j}+D_{11}^{i+1,j}\right)}+\frac{D_{22}^{i,j+1}+D_{22}^{i+1,j}}{2}\right]\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}\right)\\ +\left[\frac{D_{11}^{i+1,j}{D}_{21}^{i,j}+D_{11}^{i,j}{D}_{21}^{i+1,j}}{D_{11}^{i,j+1}+D_{11}^{i+1,j+1}}\right]\left({\phi }_{i,j+1}-{\phi }_{i,j}\right)\\ +\left[\frac{\left(D_{21}^{i+1,j+1}-D_{21}^{i+1,j+1}\right)\left(D_{12}^{i+1,j}-D_{12}^{i+1,j+1}\right)}{2\left(D_{22}^{i+1,j}+D_{22}^{i+1,j+1}\right)}+\frac{D_{11}^{i,j+1}+D_{11}^{i+1,j+1}}{2}\right]\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{3}{2},j+\frac{1}{2}}\right)\end{array}$

$\begin{array}{l}+\left[\frac{D_{22}^{i+1,j}{D}_{12}^{i+1,j+1}+D_{22}^{i+1,j+1}{D}_{12}^{i+1,j}}{D_{22}^{i+1,j}+D_{22}^{i+1,j+1}}\right]\left({\phi }_{i+1,j}-{\phi }_{i+1,j+1}\right)\\ +\left[\frac{\left(D_{21}^{i,j+1}-D_{21}^{i,j}\right)\left(D_{12}^{i,j}-D_{12}^{i,j+1}\right)}{2\left(D_{22}^{i,j}+D_{22}^{i,j+1}\right)}+\frac{D_{11}^{i,j}+D_{11}^{i,j+1}}{2}\right]\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}\right)\\ +\left[\frac{D_{22}^{i,j+1}{D}_{12}^{i,j}+D_{22}^{i,j}{D}_{12}^{i,j+1}}{D_{22}^{i,j}+D_{22}^{i,j+1}}\right]\left({\phi }_{i,j+1}-{\phi }_{i,j}\right)\\ ={\displaystyle {\int }_{K_{i+\frac{1}{2},j+\frac{1}{2}}}}f\text{d}x\forall 1\le i,j\le N-1\end{array}$ (57)

The DDFV method can be viewed as a double classical finite volume acting on the primal meshes and on the dual meshes. When the medium is isotropic systems (56) and (57) are decoupled and therefore independent. In the homogeneous and nonhomogeneous anisotropic case there is a connectivity between the

unknown ${\phi }_{i,j}$ and ${\phi }_{i+\frac{1}{2},j+\frac{1}{2}}$. For every fixed primary mesh

$K_{i;j}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]\times \left[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}\right]$ or dual mesh $K_{i+\frac{1}{2};j+\frac{1}{2}}=\left[x_i,x_{i+1}\right]\times \left[y_j,y_{j+1}\right]$

the corresponding Equation (56) or (57) involves nine unknown nodal values. For that reason in this particular case DDFV scheme is so called the nine point scheme.

4.2. Matrix Form of the Linear System (56) and (57) and Node Numbering

If we adopt the lexicographic (Figure 6) order according to which nodes (correspondingly, the unknown components) are numbered by proceeding the primary mesh to the dual mesh, from left to right and from the bottom to the top. By so doing we obtain a symmetric linear system whose matrix form is:

$\left[\begin{array}{cc}A& B\\ B^{\text{T}}& C\end{array}\right]\left[\begin{array}{c}{\Phi }_{cc}\\ {\Phi }_{vc}\end{array}\right]=\left[\begin{array}{c}{F}_{cc}\\ F_{vc}\end{array}\right]$ (58)

where $F_{cc}$ is a subvector with $card\left(\mathcal{P}\right)$ components, $F_{cv}$ is a subvector with $card\left(\mathcal{D}\right)$, A is a $card\left(\mathcal{P}\right)\times card\left(\mathcal{P}\right)$ symmetric matrix, C is a $card\left(\mathcal{D}\right)\times card\left(\mathcal{D}\right)$ symmetric matrix, B is a $card\left(\mathcal{P}\right)\times card\left(\mathcal{D}\right)$ matrix and $B^{\text{T}}$ it its associated transpose matrix. Note that the components of $F_{cc}$ and $F_{vc}$ depend on f.

4.3. The Mesh Data Structure and Connectivity between General Nodal and Its Neighbors

Global indices generally used for building the full system of equations over the computational domain while local indices are employed to define the local stencil for an element information that is useful during the discretization process. In this case system local indices are used interchangeably as global indices $\left(i\mathrm{,}j\right)$, $\left(i+\frac{\text{1}}{\text{2}}\mathrm{,}j+\frac{\text{1}}{\text{2}}\right)\leftarrow \text{local}\text{index}$, $\mathcal{N}\leftarrow \text{global}\text{index}$ because they can be readily translated to global indices and vice versa, as illustrated in Tables 1-3. As you can see on the equation systems (57) and (56), the unknowns located at the vertex of the primal cells have the half-integer indices, which complicates the task in terms of Matlab coding, to break obstacle we modify the local indices so as to have only integer indices, the numbering of the local indices remains unchanged. Prior to assembling the nodal equations in matrix form, each node needs to be assigned a unique number (global index), this because the solution to the system shown by Equation (58) is a vector ${\left[{\Phi }_{cc}\mathrm{,}{\Phi }_{vc}\right]}^{\text{T}}$ i.e. it is a 1D column matrix, not a 2D matrix. Essentially this implies combining the new i and j indices into a

Table 1. Connectivity between indices and coordinates of the associated node case of Dirichlet boundary conditions $h=0.5$.

Table 2. Connectivity between indices and coordinates of the associated nodes case of Neumann boundary conditions $h=0.5$.

Table 3. Connectivity between indices and coordinates of the associated nodes case of periodic boundary conditions $h=0.5$.

single index, which in fact is aforementioned unique number assigned to the node. But this case does not fulfill these requirements because the unknowns located on the Eastern ${\Gamma }_E$ and Western ${\Gamma }_W$ boundary are equal, similarly for those on the Northern ${\Gamma }_N$ and Southern ${\Gamma }_S$ boundary. This requires special treatment afterwards. We start by showing the full matrix assembly strategy. This requires writing the function named IndexGlobalnode which takes as input the global index and the mesh size h and as output the associated local indices, then the function Globalnode Index which takes as input the local index, the mesh size h and returns global index. Let us set a general nodal point P associated to the global index k nodal located inside the domain, which can be a vertex or a center of the primary mesh. This general nodal point P as you can see in Figure 7 is defined by its neighbors nodes on the West, East, North, South, North-West, North-East, South-west and South-East are defined as follow; W, E,

Figure 7. Grid showing both node numbers and general nodal point K with its neighbors in the case of periodic boundary conditions. (left) Connectivity node and its neighbors in the primal cell: $K\le N^2$ : right Connectivity node and its neighbors in the dual cell: $K\ge N^2+1$.

N, S, NE, NW, SW and SE. The discretized Equations (56) and (57) have found to take the following general form (59). So we have all gradients to assemble the matrix in Dirichlet and Neumann boundary case. Let’s now examine the periodic boundary case.

$C_P{\phi }_P+{\displaystyle \underset{n}{\sum }}{C}_{nb}{\phi }_{nb}=\langle f_P\rangle$ (59)

where ${\displaystyle \sum }$ indicates summation overall neighboring nodes (nb), $C_{nb}$ are the neighboring coefficients $C_W,C_W,C_S,C_N,C_{NE},C_{NW},C_{SE},C_{SW}$ and ${\phi }_{nb}$ are the approximate values of pressure at the neighboring nodes. Table 2 shows the connectivity between the indices and the coordinates of the associated nodes in the case of Neumann boundary conditions. We can see that the boundary nodes are unknowns. However for the periodic case, we select from Table 3 only the nodes belonging to $\bar{\Omega }\backslash \left({\Gamma }_N\cup {\Gamma }_E\right)$ of discretization domain, keeping in mind that to get the others by periodicity. Matlab offers several ways (meshgrid, ndgrid) to generate the list of nodes and the corresponding coordinates. The challenge remains to find the connectivity between the neighboring nodes constituting the nine point stencil. Figure 7 illustrates how we proceed to convert Equations (56) and (57) in the form (59).

5. Numerical Results

In this section, we confirm with some numerical tests the theoretical results we have proven in the previous section. For each test-case, a uniform rectangular mesh is used with different levels of refinement materialized by successive decreasing values assigned to the mesh size $h=\frac{1}{N}$. In the following test-cases we have taken: $h=\frac{1}{4},\frac{1}{8},\frac{1}{16},\frac{1}{32}$ and $\frac{1}{64}$.

5.1. Notations

In the following tables, ndu denotes the number of discrete unknowns. Recall that ${\left|\mathrm{.}\right|}_{\mathcal{A}}$ stands for the discrete $H_0^1$ - norm and let ${\Vert \mathrm{.}\Vert }_{\mathrm{1,}\mathcal{A}}$ denote the discrete $H^1$ -norm defined as:

${\Vert v\Vert }_{\mathrm{1,}\mathcal{A}}={\left[{\left|v\right|}_{\mathcal{A}}^2+{\Vert v\Vert }_{L^2\left(\Omega \right)}^2\right]}^{\frac{1}{2}}$ (60)

When the exact solution is available, the relative discrete $L^2$ -norm of the error for the exact potential is defined as:

$er{L}^2={\left(\frac{\underset{P\in \mathcal{P}}{{\displaystyle \sum }}meas\left(C_P\right){\left[\phi \left(x^p\right)-{\bar{\phi }}_P\right]}^2}{\underset{P\in \mathcal{P}}{{\displaystyle \sum }}meas\left(C_P\right){\left[\phi \left(x^p\right)\right]}^2}\right)}^{\frac{1}{2}}$

Defined by analogy, $er-{\Vert \mathrm{.}\Vert }_{\mathrm{1,}\mathcal{A}}$ is the relative discrete $H^1$ -norm of the error for the exact potential. For a given mesh, since different levels i of refinement may be considered, we denote by $er{L}^2\left(i\right)$ and $er-{\Vert \mathrm{.}\Vert }_{\mathrm{1,}\mathcal{A}}\left(i\right)$ the corresponding relative discrete $L^2$ -norm and relative discrete $H^1$ -norm of the exact potential. Let us set for any integer i (with $i\ge 2$ ):

$ra{L}^2\left(i\right)=-2\frac{ln\left[er{L}^2\left(i\right)\right]-ln\left[er{L}^2\left(i-1\right)\right]}{ln\left[ndu\left(i\right)\right]-ln\left[ndu\left(i-1\right)\right]}$

We define $ra-{\Vert \mathrm{.}\Vert }_{\mathrm{1,}\mathcal{A}}\left(i\right)$, for any integer $i\ge 2$, with the same relation as for $ra{L}^2\left(i\right)$, except that in this relation $er{L}^2\left(i\right)$ is replaced with $er-{\Vert \mathrm{.}\Vert }_{\mathrm{1,}\mathcal{A}}\left(i\right)$. On the other hand, $erflm$ stands for $L^{\infty }$ -norm of the error on the exact mean value of the flux across the mesh edges. So we have:

$erflm=\underset{\sigma \in \mathcal{E}}{\max }\frac{1}{meas\left(\sigma \right)}\left|Q_{\sigma }-{\bar{Q}}_{\sigma }\right|$

where $Q_{\sigma }$ and ${\bar{Q}}_{\sigma }$ are respectively the exact and the approximate flux across $\sigma$ which is either a primal edge or a dual edge. The symbol $er{L}^{\infty }$ denotes the pressure error for $L^{\infty }$ -norm.

Let $ocv\left[\times \times \times \right]$ denote the error order of convergence to zero for the norm $\left[\times \times \times \right]$ which may be taken to be one of the following norms ${\Vert \mathrm{.}\Vert }_{L^2}$, ${\Vert \mathrm{.}\Vert }_{\mathrm{1,}\mathcal{A}}$ and ${\Vert \mathrm{.}\Vert }_{L^{\infty }}$. The first two norms are used for pressure error estimates while the last one serves for the flux error estimate. The quantity $ocv\left[\times \times \times \right]$ is defined as:

$ocv\left[\times \times \times \right]=\frac{ln\left(er\left[\times \times \times \right]\left(imax\right)-ln\left(er\left[\times \times \times \right]\left(imax-1\right)\right)\right)}{ln\left(h_{imax}\right)-ln\left(h_{imax-1}\right)}$

where $imax$ is the maximum level of refinement of a given primal mesh.

At last, we denote by $slope\left[\times \times \times \right]$ the error order of convergence to zero for the norm $\left[\times \times \times \right]$ when computed with the least-square method. The quantity $slope\left[\times \times \times \right]$ obeys the relation:

$ln\left[er\left[\times \times \times \right]\left(i\right)\right]=slope\left[\times \times \times \right]ln\left[h\left(i\right)\right]+\nu \mathrm{,}$

where $\nu$ is a real number.

5.2. Test Problems

We consider two groups of test problems. In the first group, the medium is homogeneous and so is spatially periodic. In the second group, the medium is taken to be heterogeneous and spatially periodic.

5.2.1. Group I

We consider in this group two cases: a homogeneous isotropic medium and a homogeneous anisotropic medium.

Data from Test problem 1: The medium $\Omega =\left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[\times \left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[$ is associated with the following diffusion coefficient matrix.

$K\left(x\mathrm{,}y\right)=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ (61)

The exact solution to Equations (1) such that:

${\displaystyle {\int }_{\Omega }}u\left(x\right)\text{d}x=0$ (62)

is

$u\left(x,y\right)=\sin \left(2\pi x\right)\sin \left(2\pi y\right).$

Note that it is easy to determine the corresponding (source term) function f and check that this function satisfies the compatibility condition (12). The same remark remains valid for all the test problems analyzed in this section.

According to Table 4, DDFV computations of the approximate solution to Test problem number one display a quadratic convergence for $L^2$ -norm and discrete $H^1$ -norm concerning the pressure. The same rate of convergence is observed for ${\Vert \mathrm{.}\Vert }_{\infty }$ -norm concerning the interface fluxes. The quadratic convergence for $L^2$ -norm and discrete $H^1$ -norm numerically obtained in the homogeneous diffusion analyzed here is not in contradiction with our theoretical results. Indeed, the order one of convergence is based upon much weaker assumptions on the diffusion coefficient which is supposed to be piecewise constant. Note that similar results have been obtained for Dirichlet and Neumann boundary conditions by diverse authors in FVCA5 Benchmark [11].

Table 4. Convergence rate of flux and pressure errors for $L^{\infty }$ -norm, $L^2$ -norm and discrete $H^1$ -norm in Test-problem number 1.

Data from Test-problem 2: Let $\Omega =\left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[\times \left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[$ be a square with the following diffusion coefficient:

$K\left(x\mathrm{,}y\right)=\left[\begin{array}{cc}1& 1/2\\ 1/2& 1\end{array}\right]$ (63)

The exact solution to (1) - (11) such that

${\displaystyle {\int }_{\Omega }}u\left(x\mathrm{,}y\right)\text{d}x=0$ (64)

is

$u\left(x,y\right)=\sin \left(2\pi x\right)\cos \left(2\pi y\right).$

Results from Table 5 confirm the comment we have developed about the homogeneous flow in $\Omega$ exposed in Table 1. The result trends do not change even if one considers homogeneous media with contrasting diffusion coefficients like

$K\left(x\right)=\left[\begin{array}{cc}1& 10\\ 10& 1000\end{array}\right]$ (65)

5.2.2. Group II

Now we consider a nonhomogeneous isotropic and anisotropic porous domain $\Omega$. Computation results are presented in (Table 6) below.

Data for Test-problem 3: We have taken $\Omega =\left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[\times \left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[$ associated with the following diffusion coefficient:

$K\left(x,y\right)=\left[\begin{array}{cc}2& \sin \left(\pi x\right)\sin \left(\pi y\right)\\ \sin \left(\pi x\right)\sin \left(\pi y\right)& 1\end{array}\right]$

The exact solution to the system (1) - (11) such that

${\displaystyle {\int }_{\Omega }}u\left(x\mathrm{,}y\right)\text{d}x=0$ (66)

is what follows:

$u\left(x,y\right)=2\sin \left[\pi \left(x+y\right)\right]\cos \left[\pi \left(x+y\right)\right].$

Table 5. Convergence rate of flux and pressure errors for $L^{\infty }$ -norm, $L^2$ -norm and discrete $H^1$ -norm in Test-problem number 2.

Table 6. Convergence rate of pressure error for $L^{\infty }$ -norm, $L^2$ -norm and for discrete $H^1$ -norm in Test-problems 3 and 4.

Data for Test-problem 4: Let $\Omega$ be the square $\left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[\times \left]\mathrm{0\; <mo>\; ;\; </mo>1}\right[$ associated with the following diffusion coefficient:

$K\left(x\right)=(\begin{array}{l}\left[\begin{array}{cc}1000& 0\\ 0& 1000\end{array}\right]\text{if}x\in \left[\frac{1}{4};\frac{3}{4}\right]\times \left[\frac{1}{4};\frac{3}{4}\right]\\ \left[\begin{array}{cc}750& 0\\ 0& 2000\end{array}\right]\text{otherwise}.\end{array}$ (67)

The exact solution to the system (1) - (11) such that

${\displaystyle {\int }_{\Omega }}u\left(x\mathrm{,}y\right)\text{d}x=0$ (68)

is what follows:

$u\left(x\mathrm{,}y\right)=2\sin \left[\pi \left(x+y\right)\right]\cos \left[\pi \left(x+y\right)\right]$.

5.2.3. Concluding Remarks

The numerical simulations in comparison with the exact solutions show the accuracy of the method as you can observe it in Figure 8 and Figure 9.

The numerical experiments were performed on uniform square meshes and have shown that:

1) For isotropic homogeneous media, one gets a quadratic convergence of the approximate pressure for $L^{\infty }$ -norm, $L^2$ -norm and discrete $H^1$ -norm as well. The same convergence rate is observed concerning the flux for the vector max-norm. These results are in accordance with those published in the Finite Volume literature (see test-problem number 1).

2) For anisotropic homogeneous media, one can see that the rate of convergence remains globally the same, except for the $L^{\infty }$ -norm that displays a linear convergence (see test-problem number 2).

3) For spatially varying diffusion coefficients D, the cell mean value of D is taken to be the cell-center diffusion coefficients. Approximations of pressure are performed with the order one of convergence for $L^{\infty }$ -norm and discrete $H^1$ -norm

Figure 8. A numerical solution test-problem Group I $h=\frac{1}{2^6}$ LEFT: Approximation solution; Right: Exact solution.

Figure 9. A numerical solution test-problem Group II $h=\frac{1}{2^6}$ LEFT: Approximation solution; Right: Exact solution.

while a 1.50 order of convergence is got for $L^2$ -norm (see test-problem number 3).

4) For piecewise constant diffusion coefficients D, the same results as for spatially varying diffusion coefficients are obtained except for the $L^{\infty }$ -norm that led to a 0.5 order of convergence (see test-problem number 4).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Njifenjou, A. and Moukouop-Nguena, I. (2001)) Traitement des anisotropies de perméabilité en simulation d’écoulement en milieu poreux par les volumes finis. Systèmes Informatiques pour la Gestion de l’Environnement, Douala, 22-25 October 2001, 245-259.
[2]	Njifenjou, A. and Moukouop-Nguena, I. (2006) A Finite Volume Approximation for Second Order Elliptic Problems with a Full Matrix on Quadrilateral Grids. International Journal on Finite, 3, 1-30.
[3]	Domelevo, K. and Omnès, P. (2005) A Finite Volume Method for the Laplace Equation on Almost Arbitrary Two-Dimensional Grids. ESAIM: Mathematical Modelling and Numerical Analysis, 39, 1203-1249. https://doi.org/10.1051/m2an:2005047
[4]	Hermeline, F. (2003) Approximation of Diffusion Operators with Discontinuous Tensor Coefficients on Distorted Meshes. Computer Methods in Applied Mechanics and Engineering, 192, 1939-1959. https://doi.org/10.1016/S0045-7825(02)00644-8
[5]	Njifenjou, A. and Kinfack, A.J. (2008) Convergence Analysis of an MPFA Method for Flow Problems in Anisotropic Heterogeneous Porous Media. International Journal of Pure and Applied Mathematics, 5, 35-86.
[6]	Njifenjou, A., Donfack, H. and Moukouop-Nguena, I. (2013) Analysis on General Meshes of a Discrete Duality Finite Volume Method for Subsurface Flow Problems. Computational Geosciences, 17, 391-415. https://doi.org/10.1007/s10596-012-9339-6
[7]	Donfack, H., Njifenjou, A. and Jeutsa, A.K. (2019) A Discrete Duality Finite Volume Method for Flow Problems with Prescribed Periodic Boundary Conditions. Development and Mathematical Analysis of Finite Volumes for Complex Physics, 50 p.
[8]	Abdullah, A.N. and Omar, B.H. (2018) Numerical Treatment of Initial-Boundary Value Problems with Mixed Boundary Conditions. American Journal of Computational Mathematics, 8, 153-174. https://doi.org/10.4236/ajcm.2018.82012
[9]	Hermeline, F. (2007) A Finite Volume Method for Approximating 3D Diffusion Operators on General Meshes. Computer Methods in Applied Mechanics and Engineering, 196, 2497-2526. https://doi.org/10.1016/j.cma.2007.01.005
[10]	Kinfack, A., Njifenjou, A. and Nganhou, J. (2012) Convergence Analysis Unstructured Meshes of a DDFV Method for Flow Problems with Full Neumann Boundary Conditions. Journal of Applied Mathematics, 2016, Article ID: 5891064. https://doi.org/10.1155/2016/5891064
[11]	Herbin, R. and Hubert, F. (2008) Benchmark on Discretization Schemes for Anisotropic Diffusion Problems on General Grids. In: Eymard, R. and Herard, J.M. Eds., Finite Volume for Complex Applications V, John Wiley & Sons, France, 659-692.