Discrete Duality Finite Volume for Anisotropic Diffusion Problems with Prescribed Robin Boundary Conditions

Hubert Donfack

Department of Mathematics and Computer Sciences, The Bamenda University, Bamenda, Cameroon.

DOI: 10.4236/jamp.2025.1311225                

Abstract

This paper presents and analyzes a Discrete Duality Finite Volume (DDFV) method to solve 2D diffusion problems under prescribed Robin boundary conditions. The derivation of a symmetric discrete problem is established. The existence and uniqueness of a solution to this discrete problem are shown via the positive definiteness of its associated matrix. We show that the discrete scheme meets the Neumann problem when the parameter $\alpha \to 0$ (and, in a sense, when $\alpha \to \infty$ the Dirichlet problem). This work is a continuation of our work regarding the development of DDFV methods. The main innovation here is taking into account Robin’s boundary conditions. We provide a few steps of Matlab implementation and numerical tests to confirm the effectiveness of the method.

Keywords

Finite Volume Schemes, Discontinuous Coefficients, Diffusion Problems, Homogeneous Porous Media

1. Introduction and the Model Problem

Numerical methods play a crucial role in approximating solutions to complex mathematical problems that cannot be solved analytically. The discrete duality finite volume (DDFV) method is one of the new generation finite volumes, very popular today in Geoscience Engineering. The work from [1]-[5] has 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 work is to formulate and conduct a theoretical analysis of the flux-based DDFV method on general grids for flow problems in polygonal anisotropic nonhomogeneous media under Prescribed Robin Boundary conditions, which is the combination of Dirichlet and Neumann boundary conditions.

It is important to note that the DDFV methods come from two formulations. The first, known as the flux-based DDFV, focuses on interface flux computations for primary and dual meshes, ensuring interface flux continuity. This approach was pioneered by [4] and [5]. The second formulation, known as the gradient-based DDFV, is based on pressure gradient reconstructions over a diamond grid, first introduced in [3]. This second formulation has been further developed by mathematicians such as Andreianov, Boyer, and Hubert, who have extended key ideas to nonlinear operators of the Leray-Lions type (see [7]). Motivated by the potential for increasing the order of convergence of the gradient-based method for nonlinear operators, Boyer and Hubert proposed the modified DDFV in [8]. Here, “flux” refers to the outward normal component of the Darcy velocity on cell boundaries.

Donfack and Jeutsa [9] presented a DDFV method for 2-D flow problems in nonhomogeneous anisotropic porous media under diverse boundary conditions. They focus on the case of Dirichlet, full Neumann and periodic boundary conditions, taking into account the periodicity. They use the discrete gradient defined in diamond cells to compute the fluxes. Matlab code was also developed for algebraic equations.

Martin et al. [10] proposed a DDFV method for non-overlapping optimized Schwarz method with Robin transmission conditions to solve anisotropic diffusion problems, taking into account the anisotropic diffusion across sub-domain interfaces. They prove convergence using energy estimates for general decomposition, including cross points and fully anisotropic diffusion. Their analysis reveals that primal and dual meshes might be coupled using different optimized Robin parameters in the optimized Schwarz methods.

Njifenjou et al. [6] presented and analyzed, on unstructured grids, a discrete duality finite volume method for 2D flow problems in nonhomogeneous anisotropic porous media. For presenting our DDFV formulations, let us consider the 2D diffusion problem: Find a function $\varphi$ defined in $\overline{\Omega }$ , that satisfies the following partial differential equation:

$-div\left(A\nabla \varphi \right)+c\varphi =f\text{in}\Omega$ (1.1)

$\alpha \varphi +\left(n\cdot A\nabla \varphi \right)=g\text{on}\Gamma =\partial \Omega$ (1.2)

where $f,g,\alpha >0$ and $c$ are given functions, $\text{Γ}$ is the boundary of Ω, $n$ is the outward unit normal to the boundary and $\overline{\Omega }$ is a closure of Ω that is a given open polygonal domain (not necessarily convex). On the other hand, $A=A\left(x\right)$ , with $x={\left(x_1,x_2\right)}^t\in \Omega$ is a piecewise constant in Ω and is symmetric, that is,

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

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

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

The existence and uniqueness results of Equations (1.1) and (1.2) are well documented in the literature and summarized by the following propositions. See, for example, reference [11] for more details.

Proposition 1.1. Assume that $f\in L^2\left(\Omega \right)$ , $g\in L^2\left(\Gamma \right)$ , $c\in L^{\infty }\left(\Omega \right)$ , $\alpha \in L^{\infty }\left(\Gamma \right)$ and $V=H^1\left(\Omega \right)$ . Let $v\in V$ , then

$a\left(\varphi ,v\right)=l\left(v\right)$ (1.5)

where

$a\left(\varphi ,v\right)={\displaystyle {\int }_{\Omega }\left(A\nabla \varphi ,\nabla v+c\varphi v\right)\text{d}x}+{\displaystyle {\int }_{\Gamma }\alpha {\gamma }_0\left(\varphi \right){\gamma }_0\left(v\right)\text{d}\Gamma }$

$l\left(v\right)={\displaystyle {\int }_{\Omega }fv\text{d}x}+{\displaystyle \int g{\gamma }_0\left(v\right)\text{d}\Gamma }$

with $\varphi \in V$ , defines a variational formulation for problem model (1.1) and (1.2).

Proposition 1.2. Under the assumptions (1.3) and (1.4),$f\in L^2\left(\Omega \right)$ , $c\in L^{\infty }\left(\Omega \right)$ , $\alpha \in L^{\infty }\left(\Gamma \right)$ , $g\in L^{\infty }\left(\Gamma \right)$ , $\alpha \ge 0$ . Then the problem (1.5) has one and only one solution.

Proof. It follows from the Lax-Milgram theorem □

For sake of simplicity in the presentation, we will consider the functions $\alpha$ and $c$ as positive constants.

2. The Discrete Problem: Existence and Uniqueness

We start with exhibiting the discrete problem. Then, after we focus on the existence and uniqueness of a discrete solution.

2.1. Discretization Domain

Figure 1. Example of primary matching unstructured mesh.

The DDFV theory exposed here is inspired by the one developed in [6]. Therefore, our exposition is based upon a primal mesh, an auxiliary mesh and the associated dual mesh clearly defined in what follows. The spatial domain Ω is split into a family of convex open polygons with matching interfaces. This family defines what is called primal mesh denoted by . Figure 1 below illustrates an example of a primal mesh made up of convex polygons. Following an original idea from [6] one can introduce an auxiliary mesh denoted by $\mathcal{A}$ . This mesh plays a key role in either the construction of our DDFV scheme or the stability analysis and error estimates developed later. For deriving an auxiliary mesh from a given primal mesh, one starts with arbitrary fixing one point per internal edge. Let denote by ${ℰ}^{int}$ the set of these edge-points. For a good understanding of the way the auxiliary mesh should be built, one should know the concept of neighboring edge-points introduced in [6]. Two edge-points are neighboring if they share a primal cell and their corresponding edges intersect at the same vertex of this primal cell. The auxiliary mesh $\mathcal{A}$ is generated by joining with a straight line all neighboring edge-points (see Figure 2). Note that this auxiliary mesh permits us to fix perimeters in which should necessary be located a finite family of primal cellpoints of ensuring the regularity of the gridding. There is a trivial bijective map between the family of cellpoints and the family of primal cells. So, to any cellpoint $P$ one may associate a unique primal cell denoted by $C_P$ , and vice versa. It is then clear that both families could be indexed by . Joining any cellpoint with any edgepoint sharing the same primal cell leads is called (in the sequel) a dual mesh. (see Figure 3)

Figure 2. A primary mesh and the associated auxiliary mesh (dooted lines), including edgepoints and cellpoints in black and blue colors.

Figure 3. Combination of a primary mesh and the associated dual mesh (red discontinuous lines) including the auxiliary mesh.

Main Assumption

• We assume that the primary mesh is compatible with the discontinuities of the permeability tensor $A$ defined in Ω.

• The Permeability discontinuities divide Ω into a finite number of convex polygonal subsets denoted by ${\left\{{\Omega }_s\right\}}_{s\in S}$ .

We suppose that the restriction over ${\overline{\Omega }}_s$ of the exact solution $\varphi$ to Equations (1.1)-(1.2), denoted by ${\varphi |}_{{\Omega }_s}$ , satisfies the following property:

${\varphi |}_{{\Omega }_s}\in C^2\left({\overline{\Omega }}_s\right)\forall s\in S$(2.1)

2.2. DDFV Formulation of the Model Problem

We shall now focus on a DDFV formulation of the problem (1.1)-(1.2) which involves as discrete unknowns and ${\left\{{\varphi }_{D^{*}}\right\}}_{D^{*}\in \mathcal{D}}$ expected to be close approximations of (cell-point pressures) and ${\left\{{\overline{\varphi }}_{D^{*}}\right\}}_{D^{*}\in \mathcal{D}}$ (vertex-point pressures) respectively where ${\varphi }_P=\varphi \left(x_1^P,x_2^P\right)$ and ${\varphi }_{D^{\star }}=\varphi \left(x_1^{D^{\star }},x_2^{D^{\star }}\right)$ and where $\mathcal{D}$ represent the dual mesh. As developed in [6], we are inspired to present the DDFV theory here.

2.2.1. DDFV Flux Approximation in Primal and Dual Cells

Consider $C_P$ as a primary cell, with $P$ representing the corresponding cell point. We integrate both sides of the balance Equation (1.1) within $C_P$ . The application of Ostrogradsky’s theorem to the integral on the left-hand side of this equation involves computing the flux across the boundary $C_P$ .

${\displaystyle \sum _{\sigma \in \partial C_P}\underset{\text{Flux}}{\underbrace{-{\displaystyle {\int }_{\partial C_P}A\nabla \varphi \cdot {\eta }_{\partial C_P}\text{d}\gamma }}}}+c{\displaystyle {\int }_{C_P}\varphi \text{d}x}={\displaystyle {\int }_{C_P}f\text{d}x}$ (2.2)

where ${\eta }_{\partial C_P}$ is the unit normal to the boundary of $\partial C_P$

Using a suitable quadrature formula for approximating this flux, a discrete balance equation is derived. To illustrate our ideas, we consider an edge $\sigma =\left[A^{\text{*}}{B}^{\text{*}}\right]$ associated with the primary cell $C_P$ see Figure 4.

Figure 4. Two molecules for a DDFV computation of the flux across the edge $\sigma =\left[A^{\text{*}}{B}^{\text{*}}\right]$ . Left: $\left[A^{\text{*}}{B}^{\text{*}}\right]$ lies inside Ω, Right: $\left[A^{\text{*}}{B}^{\text{*}}\right]$ is part of the boundary of Ω.

Let $A^P$ be the absolute permeability tensor of the cell $C_P$ . Denoted by ${\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P$ the unit normal vector to $\left[A^{\text{*}}{B}^{\text{*}}\right]$ exterior to $C_P$ , the flux expression over the edge $\left[A^{\text{*}}{B}^{\text{*}}\right]$ viewed as part of the boundary of $C_P$ is given by

$F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P=-{\displaystyle {\int }_{A^{\text{*}}I}\nabla \varphi \cdot \left(A^P{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\right)\text{d}\gamma }-{\displaystyle {\int }_{I{B}^{\text{*}}}\nabla \varphi \cdot \left(A^P{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\right)\text{d}\gamma }$ (2.3)

Before starting with the flux computations across the subedges $\left[A^{\text{*}}I\right]$ and $\left[I{B}^{\text{*}}\right]$ , let us set

Then, it is easily seen that the following identity holds:

$A^P{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P=a_h\left(A^P\right){\sigma }_P+b_h\left(A^P\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}$ (2.4)

where the real numbers $a_h\left(A^P\right)$ and $b_h\left(A^p\right)$ are given by the relations

$a_h\left(A^P\right)=\frac{{\left({\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\right)}^t{A}^P{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P}{\text{cos}\left({\theta }_h^{P,I}\right)}$

$b_h\left(A^P\right)=\frac{{\left({\xi }_{\left[PI\right]}^{A^{\text{*}}}\right)}^t{A}^P{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P}{\text{cos}\left({\theta }_h^{P,I}\right)}$

${\theta }_h^{P,I}=mes\left({\sigma }_P,\widehat{{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P}\right)$

Using assumption (2.1), the decomposition (2.4) and the Taylor expansion series the flux through the edge $\sigma =\left[A^{\star }{B}^{\star }\right]$ is written as follows

$\begin{array}{c}{F}_{\left[A^{\star }{B}^{\star }\right]}^P=b_h\left(A^P\right)\left({\varphi }_{A^{\star }}-{\varphi }_I\right)+a_h\left(A^P\right)\frac{h_{I{A}^{\star }}}{h_{PI}}\left({\varphi }_P-{\varphi }_I\right)+T_{\left[A^{\star }I\right]}^P\\ +b_h\left(A^P\right)\left({\varphi }_I-{\varphi }_{B^{\star }}\right)+a_h\left(A^P\right)\frac{h_{I{B}^{\star }}}{h_{PI}}\left({\varphi }_P-{\varphi }_I\right)+T_{\left[I{B}^{\star }\right]}^P\\ =b_h\left(A^P\right)\left({\varphi }_{A^{\star }}-{\varphi }_{B^{\star }}\right)+a_h\left(A^P\right)\frac{h_{A^{\star }{B}^{\star }}}{h_{PI}}\left({\varphi }_P-{\varphi }_I\right)+T_{\left[A^{\star }{B}^{\star }\right]}^P\end{array}$ (2.5)

where $h_{I{A}^{\text{*}}}=\left|\overrightarrow{I{A}^{\text{*}}}\right|,h_{I{B}^{\text{*}}}=\left|\overrightarrow{I{B}^{\text{*}}}\right|,h_{PI}=\left|\overrightarrow{PI}\right|$ and the truncation errors $T_{\left[A^{\text{*}}I\right]}^P, T_{\left[I{B}^{\text{*}}\right]}^P$ are given by

$\begin{array}{l}{T}_{\left[A^{\text{*}}I\right]}^P=\frac{a_h\left(A^P\right)}{2}\left[h_{I{A}^{\text{*}}}^2{\left({\sigma }_P\right)}^t{\varphi }^{″}\left(M\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}-h_{I{A}^{\text{*}}}{h}_{IP}{\sigma }_P^t{\varphi }^{″}\left(Q\right){\sigma }_P\right]\\ T_{\left[I{B}^{\text{*}}\right]}^P=\frac{a_h\left(A^P\right)}{2}\left[h_{I{B}^{\text{*}}}^2{\left({\sigma }_P\right)}^t{\varphi }^{″}\left(N\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}-h_{I{A}^{\text{*}}}{h}_{IP}{\sigma }_P^t{\varphi }^{″}\left(Q\right){\sigma }_P\right]\\ T_{I\left[A^{\text{*}}{B}^{\text{*}}\right]}=T_{\left[A^{\text{*}}I\right]}^P+T_{\left[I{B}^{\text{*}}\right]}^P\end{array}$ (2.6)

with $M\in \left[A^{\text{*}},I\right],Q\in \left[P,I\right],N\in \left[B^{\text{*}},I\right]$ and ${\varphi }^{″}(.)$ being the Hessian of $\varphi$ . Note that ${\varphi }^{″}(.)$ exists, thanks to Equation (2.1).

For estimating the truncation error $T_{I\left[A^{\text{*}}{B}^{\text{*}}\right]}$ , we start with some useful notations. First of all, recall that $ℰ$ is the set made of edge-points (from the primary mesh of course). We denote by ${ℰ}^{int}$ the subset of $ℰ$ made up of edge-points lying on Ω, ${ℰ}^{ext}$ the subset of $ℰ$ made up of edge-points lying on the boundary of Ω, and ${ℰ}^P$ (for ) the subset of $ℰ$ made up of edge-points lying on the boundary of the primary cell $C_P$ .

Remark 2.1. Note that the set $ℰ$ will sometimes be identified with the set of primary edges since there is a trivial bijection between the two sets. In the same order of ideas, the set of cellpoints and the set of vertices will sometimes be identified with the set of primary cells and the set $\mathcal{D}$ of dual cells, respectively. At last, primary mesh and primal mesh mean the same thing in this work [6].

Definition 2.2. The system defines a regular mesh if the following condition is fulfilled: There exists $\theta \in \left]0,\frac{\pi }{2}\right[$ , not depending on $h,$ such that

(2.7)

Ingredients are gathered for estimating truncation errors. In this connection, it is easily seen that the following result holds:

Proposition 2.3. Assume that the mesh system defines a regular mesh system in the sense of the previous definition and that there exists $0<\varpi <1$ , mesh independent, such that

(2.8)

where $A_I^{\text{*}},B_I^{\text{*}}\in \mathcal{D}$ are extremities of the only edge (from the cell $C_P$ ) involving the edge-point $I$ .

Under the assumptions (2.1) and Equation (2.7), there exists a strictly positive number $C$ , mesh independent, such that

$\left|T_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\right|\le C{h}^2$ (2.9)

The relation (2.9) is a consistency property and so naturally allow us to approximate $F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P$ as follows:

$F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\approx b_h\left(A^P\right)\left({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}}\right)+a_h\left(A^P\right)\frac{h_{A^{\text{*}}{B}^{\text{*}}}}{h_{PI}}\left({\varphi }_P-{\varphi }_I\right)$ (2.10)

When $\left[A^{\text{*}}{B}^{\text{*}}\right]$ is part of the domain boundary, the edge pressure ${\varphi }_I$ is determined by the Robin conditions. it is worth noting that in this case, both points $A^{\text{*}}$ and $B^{\text{*}}$ must be situated on the boundary as well. However, if $\left[A^{\text{*}}{B}^{\text{*}}\right]$ serve as an interface between cell $C_P$ and some primary cell denoted by $C_L$ , the edge pressure ${\varphi }_I$ is an auxiliary unknown. But, thanks to the principle of flux continuity, one can approximate it with a linear function of ${\varphi }_P,{\varphi }_L,{\varphi }_{A^{\text{*}}}$ , and ${\varphi }_{B^{\text{*}}}$ . For investigating the above-mentioned linear function, we compute the flux across $\left[A^{\text{*}}{B}^{\text{*}}\right]$ viewed as part of the boundary of $C_L$ . For this purpose, we set

$A^L{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^L={\widehat{a}}_h\left(A^L\right){\sigma }_L+{\widehat{b}}_h\left(A^L\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}$(2.11)

where the real numbers ${\widehat{a}}_h\left(A^L\right)$ and ${\widehat{b}}_h\left(A^L\right)$ are given by the relations

${\widehat{a}}_h\left(A^L\right)=\frac{{\left({\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\right)}^t{A}^L{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P}{\text{cos}\left({\theta }_h^{L,I}\right)}$

${\widehat{b}}_h\left(A^L\right)=\frac{{\left({\xi }_{\left[IL\right]}^{A^{\text{*}}}\right)}^t{A}^L{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P}{\text{cos}\left({\theta }_h^{P,I}\right)}$

${\theta }_h^{L,I}=mes\left({\sigma }_L,\widehat{{\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P}\right)$

where ${\xi }_{\left[IL\right]}^{A^{\text{*}}}$ denotes the unit normal vector to $\left[IL\right]$ exterior to the triangle $\left(LI{A}^{\text{*}}\right)$ . Performing flux computation over the interelement $\left[A^{\text{*}}{B}^{\text{*}}\right]$ , viewed as part of the boundary of $C_L$ (see Figure 4) leads to

$F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^L={\widehat{b}}_h\left(A^L\right)\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)+{\widehat{a}}_h\left(A^L\right)\frac{h_{A^{\text{*}}{B}^{\text{*}}}}{h_{IL}}\left({\varphi }_L-{\varphi }_I\right)+T_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^L$(2.12)

where for a fixed $Q\in \left[IL\right]$ and fixed $M,N\in \left[A^{\text{*}}{B}^{\text{*}}\right]$ , we have set

$\begin{array}{c}{T}_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^L=\frac{{\widehat{a}}_h\left(A^L\right)}{2}[h_{I{A}^{\text{*}}}^2{\left({\sigma }_L\right)}^t{\varphi }^{″}\left(M\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}+h_{I{B}^{\text{*}}}^2{\left({\sigma }_L\right)}^t{\varphi }^{″}\left(N\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}\\ -h_{A^{\text{*}}{B}^{\text{*}}}{h}_{IP}{\sigma }_L^t{\varphi }^{″}\left(Q\right){\sigma }_L]\end{array}$(2.13)

Thus, this flux can be approximated with the expression

$F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^L\approx {\widehat{b}}_h\left(A^L\right)\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)+{\widehat{a}}_h\left(A^L\right)\frac{h_{A^{\text{*}}{B}^{\text{*}}}}{h_{IL}}\left({\varphi }_L-{\varphi }_I\right)$(2.14)

The approximate fluxes $F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P$ and $F_{\left[A^{\text{*}}{B}^{\text{*}}\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 ${\varphi }_I$ satisfies to the following relation:

$\begin{array}{c}{\varphi }_I=\frac{h_{IP}{h}_{IL}}{h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)}\{\frac{a_h\left(A^P\right)}{h_{IP}}{\varphi }_P+\frac{{\widehat{a}}_h\left(A^L\right)}{h_{IL}}{\varphi }_L\\ +\frac{1}{h_{A^{\text{*}}{B}^{\text{*}}}}\left(b_h\left(A^P\right)-\widehat{b}\left(A^L\right)\right)\left({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}}\right)\}\end{array}$(2.15)

for all $I\in {ℰ}^{int}$ .

This is a consistent approximation for ${\varphi }_I$ in the sense that the corresponding truncation error vanishes when $h$ goes to zero. So, replacing ${\varphi }_I$ in Equation (2.14) by its approximate value yields the following conservative scheme:

$\begin{array}{c}{F}_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P\approx \left[\frac{a_h\left(A^P\right){\widehat{b}}_h\left(A^L\right)h_{IL}+h_{IP}{\widehat{a}}_h\left(A^L\right)b_h\left(A^P\right)}{h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)}\right]\left({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}}\right)\\ +\left[\frac{a_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{A^{\text{*}}{B}^{\text{*}}}}{h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)}\right]\left({\varphi }_P-{\varphi }_L\right)\end{array}$(2.16)

Case where $I\in {ℰ}^P\cap \partial \Omega$ : How can we approximate ${\varphi }_I$ ?

Recall that: $-\alpha \varphi -\left(\eta \cdot A▽\varphi \right)=-g$ on $\Gamma =\partial \Omega$ .

$\begin{array}{l}-{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}\alpha \varphi \text{d}\gamma }-{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}A▽\varphi \cdot {\eta }_{\partial \text{Ω}}\text{d}\gamma }=-{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}g\text{d}\gamma }\\ F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P=\alpha h_{\left[A^{\text{*}}{B}^{\text{*}}\right]}{\varphi }_I-{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}g\text{d}\gamma }\end{array}$ (2.17)

substituting the value of $F_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P$ from Equation (2.17) and simplifying, we have the following

$\begin{array}{c}{\varphi }_I=\frac{1}{\left(\alpha h_{IP}+a_h\left(A^P\right)\right)h_{A^{\text{*}}{B}^{\text{*}}}}\left\{h_{IP}{b}_h\left(A^P\right)\left({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}}\right)+a_h\left(A^P\right)h_{A^{\text{*}}{B}^{\text{*}}}{\varphi }_P\right\}\\ +\frac{h_{IP}}{\left(\alpha h_{IP}+a_h\left(A^P\right)\right)h_{A^{\text{*}}{B}^{\text{*}}}}{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}g\text{d}\gamma }\end{array}$ (2.18)

substituting ${\varphi }_I$ as given in Equation (2.18) into (2.17), we have the following flux approximation for boundary primal cell

$\begin{array}{c}{F}_{\left[A^{\star }{B}^{\star }\right]}^P=\alpha h_{A^{\star }{B}^{\star }}{\varphi }_I-{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}\gamma }\\ ={\displaystyle \sum _{I\in {ℰ}^P\cap {ℰ}^{ext}}[\frac{\alpha h_{IP}{b}_h\left(A^P\right)}{\alpha h_{IP}+a_h\left(A^P\right)}\left({\varphi }_{A^{\star }}-{\varphi }_{B^{\star }}\right)+\frac{\alpha a_h\left(A^P\right)h_{A^{\star }{B}^{\star }}}{\alpha h_{IP}+a_h\left(A^P\right)}{\varphi }_P]}\\ +{\displaystyle \sum _{I\in {ℰ}^P\cap {ℰ}^{ext}}\frac{-a_h\left(A^P\right)}{\alpha h_{IP}+a_h\left(A^P\right)}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}\gamma }}\end{array}$ (2.19)

Using the previous notations and thanks to the consistency of the previous DDFV flux approximations, the approximate flux balance equation within $C_P$ is

(2.20)

where $A^{\text{*}}$ and $B^{\text{*}}$ are extremities of the edge containing $I$ and lying on the boundary of the primary cell $C_P$ and where is such that ${\overline{C}}_P\cap {\overline{C}}_L=\left[A^{\text{*}}{B}^{\text{*}}\right]$.

2.2.2. DDFV Flux Approximation in Dual Cells

It is clear that the number of discrete unknowns and ${\left\{{\varphi }_{P^{\text{*}}}\right\}}_{P^{\text{*}}\in \mathcal{D}}$ is greater than the number of discrete balance equations given by the system of Equation (2.20) valid for all . 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_{B^{\text{*}}}$ . For carrying out our technique, we need to introduce the notion of pseudo-edge associated with dual cells see Figure 5.

Figure 5. Illustration of a dual cell (blue dotted line) with its four pseudo-edges that intersect primal edges at the red point interface of $ℰ$ .

Definition 2.4. Let $P$ and $L$ be two cellpoints from the primary mesh (that is, ) such that the corresponding primary cells $C_P$ and $C_L$ are adjacent, and consider $I\in {ℰ}^P\cap {ℰ}^L$ (recall that ${ℰ}^E$ , for , is the set of edge-points from $ℰ$ lying on the boundary of cell $C_E$ ). The line $\left[PI\right]\cup \left[IL\right]$ defines a pseudo-edge, denoted by $PIL$ , whose extremities are $P$ and $L$ . We will say that a pseudo-edge is associated with a dual cell $C_{B^{*}}$ if it is part of the boundary of $C_{B^{*}}$

Remark 2.5. Note that the boundary of any dual cell is a union of finite number of pseudo-edges.

Now, we will examine discrete balance equations across dual cells. This process involves the following steps: We begin by integrating both sides of Equation (1.1) within a dual cell $C_{B^{\text{*}}}$ as depicted in Figure 3. The application of Ostrogradski’s theorem and the utilization of Remark (2.5) result in:

${\displaystyle \sum _{I\in {ℰ}^{B^{\star }}}-{\displaystyle {\int }_{\left[PIL\right]}A\nabla \varphi \cdot {\eta }_{B^{\star }}\text{d}\gamma }}+c{\displaystyle {\int }_{C^{B^{\star }}}\varphi \cdot \text{d}x}={\displaystyle {\int }_{C^{B^{\star }}}f\cdot \text{d}x}$ (2.21)

where ${\eta }_{B^{\text{*}}}$ is the outward unit normal vector of the boundary of $C_{B^{\text{*}}}$ and $\left[PIL\right]$ is a pseudo-edge associated with the dual cell $C_{B^{\text{*}}}$ . Recall that ${\epsilon }^{B^{\text{*}}}$ is the set of edge points lying in the boundary of the dual cell $C_{B^{\text{*}}}$ . We will seek a flux approximation along the pseudo-edge $\left[PIL\right]$ , considering it as a segment of the boundary of $C_{B^{\text{*}}}$ . Denoting the exact flux over $\left[PIL\right]$ as $F_{\left[PIL\right]}^{B^{\text{*}}}$ , it can be expressed by the following relation:

$\begin{array}{c}{F}_{\left[PIL\right]}^{B^{\text{*}}}=-{\displaystyle {\int }_{\left[PI\right]}\nabla \varphi \cdot \left(A^P{\xi }_{\left[PI\right]}^{B^{\text{*}}}\right)\text{d}\gamma }-{\displaystyle {\int }_{\left[IL\right]}\nabla \varphi \cdot \left(A^L{\xi }_{IL}^{B^{\text{*}}}\right)\text{d}\gamma }\\ =F_{\left[PI\right]}^{B^{\text{*}}}+F_{\left[IL\right]}^{B^{\text{*}}}\end{array}$ (2.22)

Initially, our attention should be directed towards computing the flux across $\left[PI\right]$ . Therefore, let’s define the subsequent decomposition of $A^P{\xi }_{\left[PI\right]}^{B^{\text{*}}}$ :

$A^P{\xi }_{\left[PI\right]}^{B^{\text{*}}}=c_h\left(A^P\right){\sigma }_P-d_h\left(A^P\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}$ (2.23)

It is clear that

$c_h\left(A^P\right)=b_h\left(A^P\right)$

$d_h\left(A^P\right)=\frac{{\left({\xi }_{\left[PI\right]}^{B^{\text{*}}}\right)}^t{A}^P{\xi }_{\left[PI\right]}^{B^{\text{*}}}}{\text{cos}\left({\theta }_h^{P,I}\right)}$

Therefore

$\begin{array}{c}{F}_{\left[PI\right]}^{B^{\text{*}}}=-{\displaystyle {\int }_{\left[PI\right]}\nabla \varphi \cdot \left(A^P{\xi }_{\left[PI\right]}^{B^{\text{*}}}\right)\text{d}\gamma }\\ =b_h\left(A^P\right){\displaystyle {\int }_{\left[PI\right]}\nabla \varphi \cdot {\sigma }_P\text{d}\gamma }+d_h{\displaystyle {\int }_{\left[PI\right]}\nabla \varphi \cdot {\tau }_{A^{\text{*}}{B}^{\text{*}}}\text{d}\gamma }\\ =b_h\left(A^P\right)\left({\varphi }_P-{\varphi }_I\right)+\frac{h_{PI}}{h_{A^{\text{*}}{B}^{\text{*}}}}{d}_h\left(A^P\right)\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)+T_{\left[PI\right]}^{B^{\text{*}}}\end{array}$ (2.24)

where

$\begin{array}{c}{T}_{\left[PI\right]}^{B^{\text{*}}}=\frac{d_h\left(A^P\right)}{2}[h_{PI}^2{\left({\tau }_h\right)}^t{\varphi }^{″}\left(Q_{PI}\right){\sigma }_P+\frac{h_{PI}{h}_{I{B}^{\text{*}}}^2}{h_{A^{\text{*}}{B}^{\text{*}}}}{\left({\tau }_h\right)}^t{\varphi }^{″}\left(I{B}^{\text{*}}\right){\tau }_h\\ -\frac{h_{PI}{h}_{I{A}^{\text{*}}}^2}{h_{A^{\text{*}}{B}^{\text{*}}}}{\left({\tau }_h\right)}^t{\varphi }^{″}\left(I{A}^{\text{*}}\right){\tau }_h]\end{array}$ (2.25)

with ${\tau }_h={\tau }_{A^{\text{*}}{B}^{\text{*}}},Q_{PI}\in \left[PI\right],Q_{I{A}^{\text{*}}}\in \left[I{A}^{\text{*}}\right]$ and $Q_{I{B}^{\text{*}}}\in \left[I{B}^{\text{*}}\right]$ .

Neglecting $T_{\left[PI\right]}^{B^{\text{*}}}$ and exploiting Equation (2.15) yields the following approximation:

$\begin{array}{l}{F}_{\left[PI\right]}^{B^{\text{*}}}\approx \frac{h_{IP}\left[d_h\left(A^P\right)\left\{h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right\}+b_h\left(A^P\right)h_{IL}\left\{{\widehat{b}}_h\left(A^L\right)-b_h\left(A^P\right)\right\}\right]}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]{\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)}^{-1}}\\ \frac{{\widehat{a}}_h\left(A^L\right)b_h\left(A^P\right)h_{IP}}{\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_P-{\varphi }_L\right)\end{array}$(2.26)

Likewise, let’s concentrate on computing the flux across $\left[IL\right]$ . To achieve this, we establish the following setting:

$A^L{\xi }_{\left[IL\right]}^{B^{\text{*}}}=-{\widehat{c}}_h\left(A^L\right){\sigma }_L-{\widehat{d}}_h\left(A^L\right){\tau }_{A^{\text{*}}{B}^{\text{*}}}$(2.27)

where

${\widehat{c}}_h\left(A^L\right)={\widehat{b}}_h\left(A^L\right)$

${\widehat{d}}_h\left(A^L\right)=\frac{{\left({\xi }_{\left[IL\right]}^{B^{\text{*}}}\right)}^t{A}^L{\xi }_{\left[IL\right]}^{B^{\text{*}}}}{\text{cos}\left({\theta }_h^{L,I}\right)}$

Following the same procedure as for $\left[PI\right]$ , we obtain

$\begin{array}{c}{F}_{\left[IL\right]}^L={\widehat{c}}_h\left(A^L\right){\displaystyle {\int }_{\left[IL\right]}\nabla \varphi \cdot {\sigma }_L\text{d}\gamma }+{\widehat{d}}_h\left(A^L\right){\displaystyle {\int }_{\left[IL\right]}\nabla \varphi \cdot {\tau }_{A^{\text{*}}{B}^{\text{*}}}\text{d}\gamma }\\ ={\widehat{b}}_h\left(A^L\right)\left({\varphi }_I-{\varphi }_L\right)+\frac{h_{IL}}{h_{A^{\text{*}}{B}^{\text{*}}}}{\widehat{d}}_h\left(A^L\right)\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)+T_{\left[IL\right]}^{B^{\text{*}}}\end{array}$(2.28)

where

$\begin{array}{c}{T}_{\left[IL\right]}^{B^{\text{*}}}=\frac{d_h\left(A^L\right)}{2}[h_{IL}^2{\left({\tau }_h\right)}^t{\varphi }^{″}\left(Q_{IL}\right){\sigma }_L+\frac{h_{IL}{h}_{I{B}^{\text{*}}}^2}{h_{A^{\text{*}}{B}^{\text{*}}}}{\left({\tau }_h\right)}^t{\varphi }^{″}\left(I{B}^{\text{*}}\right){\tau }_h\\ -\frac{h_{IL}{h}_{I{A}^{\text{*}}}^2}{h_{A^{\text{*}}{B}^{\text{*}}}}{\left({\tau }_h\right)}^t{\varphi }^{″}\left(I{A}^{\text{*}}\right){\tau }_h]\end{array}$ (2.29)

Neglecting the error term and substituting the value of ${\varphi }_I$ from Equation (2.28) yield

$\begin{array}{c}{F}_{\left[IL\right]}^{B^{\text{*}}}\approx \frac{h_{IL}\left[{\widehat{d}}_h\left(A^L\right)\left\{h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right\}+{\widehat{b}}_h\left(A^L\right)h_{PI}\left\{b_h\left(A^P\right)-{\widehat{b}}_h\left(A^L\right)\right\}\right]}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]{\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)}^{-1}}\\ +\frac{{\widehat{b}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}}{\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_P-{\varphi }_L\right)\end{array}$(2.30)

Proposition 2.6. Under the same assumptions as those of Proposition 2.3, there exists a strictly positive number $C$ , mesh independent, such that

$\left|T_{\left[PI\right]}^{B^{\text{*}}}\right|+\left|T_{\left[IL\right]}^{B^{\text{*}}}\right|\le C{h}^2$ (2.31)

We summarize the flux approximation across the pseudo-edge $\left[PIL\right]$ in the following way:

$\begin{array}{c}{F}_{\left[PIL\right]}^{B^{\text{*}}}\approx \frac{{\widehat{b}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}+b_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{IP}}{h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)}\left({\varphi }_P-{\varphi }_L\right)\\ +\frac{w_h\left(P,L,I\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]{\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)}^{-1}}\end{array}$(2.32)

where we have set

$\begin{array}{c}{w}_h\left(P,L,I\right)=\left[a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{PI}\right]\times \left[d_h\left(A^P\right)h_{PI}+{\widehat{d}}_h\left(A^L\right)h_{IL}\right]\\ +h_{PI}{h}_{IL}\left[{\widehat{b}}_h\left(A^L\right)-b_h\left(A^P\right)\right]\times \left[b_h\left(A^P\right)-\widehat{b}\left(A^L\right)\right]\end{array}$(2.33)

We have to compute the flux around the degenerated dual cells i.e. with $I\in {ℰ}^{B^{\star }}\cap \partial \Omega$ , we substitute ${\varphi }_I$ as given in Equation (2.15) into Equation (2.28). That is

$\begin{array}{l}{F}_{PI}^{B^{\text{*}}}=\left[b_h\left(A^P\right)\left({\varphi }_P-{\varphi }_I\right)+\frac{h_{PI}}{h_{A^{\text{*}}{B}^{\text{*}}}}{d}_h\left(A^P\right)\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)\right]\\ ={\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}\left[\frac{\alpha b_h\left(A^P\right)h_{IP}}{\alpha h_{IP}+a_h\left(A^P\right)}{\varphi }_P+\frac{h_{IP}\left(\alpha d_h\left(A^P\right)h_{IP}+\left(a_h\left(A^P\right)d_h\left(A^P\right)-b_h^2\left(A^P\right)\right)\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left(\alpha h_{IP}+a_h\left(A^P\right)\right){({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}})}^{-1}}\right]\\ -{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}\frac{b_h\left(A^P\right)h_{IP}}{\left(\alpha h_{IP}+a_h\left(A^P\right)\right)h_{A^{\text{*}}{B}^{\text{*}}}}{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}g\text{d}\gamma }\end{array}$ (2.34)

We deduce from the previous developments that the discrete balance equation in any dual cell $C_{B^{\text{*}}}$ reads

$\begin{array}{l}{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{int}}\{\frac{{\widehat{b}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}+b_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{PI}}{\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_P-{\varphi }_L\right)\\ +\frac{w_h\left(P,L,I\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)\}+c{\displaystyle \sum _{B^{\star }\in C_{B^{\star }}}meas\left(C_{B^{\star }}\right){\varphi }_{B^{\star }}}\\ {\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}\left[\frac{\alpha b_h\left(A^P\right)h_{IP}}{\alpha h_{IP}+a_h\left(A^P\right)}{\varphi }_P+\frac{h_{IP}\left(\alpha d_h\left(A^P\right)h_{IP}+\left(a_h\left(A^P\right)d_h\left(A^P\right)-b_h^2\left(A^P\right)\right)\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left(\alpha h_{IP}+a_h\left(A^P\right)\right){({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}})}^{-1}}\right]\\ =meas\left(C_{B^{\text{*}}}\right)f_{B^{\text{*}}}+{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}\left[\frac{b_h\left(A^P\right)h_{IP}}{\left(\alpha h_{IP}+a_h\left(A^P\right)\right)h_{A^{\text{*}}{B}^{\text{*}}}}\right]{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}g\text{d}\gamma }\\ +{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}{\displaystyle {\int }_{C_{B^{\text{*}}}\cap \partial \Omega }g\text{d}\gamma }\forall B^{\text{*}}\in \mathcal{D}\text{and}{B}^{\text{*}}\in D\end{array}$(2.35)

2.2.3. DDFV Scheme of the Model Problem (1.1)-(1.2)

From Equations (2.34) and (2.35), we define a DDFV formulation of the model problem (1.1)-(1.2) as: Find and ${\left\{{\varphi }_{D^{\text{*}}}\right\}}_{D^{\text{*}}\in \mathcal{D}}$ such that

(2.36)

$\begin{array}{l}{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{int}}\{\frac{{\widehat{b}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}+b_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{PI}}{\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_P-{\varphi }_L\right)\\ +\frac{w_h\left(P,L,I\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)\}+c{\displaystyle \sum _{B^{\star }\in C_{B^{\star }}}meas\left(C_{B^{\star }}\right){\varphi }_{B^{\star }}}\\ {\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}\left[\frac{\alpha b_h\left(A^P\right)h_{IP}}{\alpha h_{IP}+a_h\left(A^P\right)}{\varphi }_P+\frac{h_{IP}\left(\alpha d_h\left(A^P\right)h_{IP}+\left(a_h\left(A^P\right)d_h\left(A^P\right)-b_h^2\left(A^P\right)\right)\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left(\alpha h_{IP}+a_h\left(A^P\right)\right){({\varphi }_{A^{\text{*}}}-{\varphi }_{B^{\text{*}}})}^{-1}}\right]\\ =meas\left(C_{B^{\text{*}}}\right)f_{B^{\text{*}}}+{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}\left[\frac{b_h\left(A^P\right)h_{IP}}{\left(\alpha h_{IP}+a_h\left(A^P\right)\right)h_{A^{\text{*}}{B}^{\text{*}}}}\right]{\displaystyle {\int }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}g\text{d}\gamma }\\ +{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}{\displaystyle {\int }_{C_{B^{\text{*}}}\cap \partial \text{Ω}}g\text{d}\gamma }\forall B^{\text{*}}\in \mathcal{D}\end{array}$(2.37)

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

Remark 2.7. If $\alpha =0$ we obtain the balance equation of diffusion reaction problem with prescribed Neumann bounbary conditions which is written as follows

(2.38)

$\begin{array}{l}{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{int}}\{\frac{{\widehat{b}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}+b_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{PI}}{\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_P-{\varphi }_L\right)\\ +\frac{w_h\left(P,L,I\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)\}+c{\displaystyle \sum _{B^{\star }\in C_{B^{\star }}}meas\left(C_{B^{\star }}\right){\varphi }_{B^{\star }}}\end{array}$

$\begin{array}{l}{\displaystyle \sum }_{I\in {ℰ}^{B^{\star }}\cap {ℰ}^{ext}}\frac{h_{IP}\left(a_h\left(A^P\right)d_h\left(A^P\right)-b_h^2\left(A^P\right)\right)}{h_{A^{\star }{B}^{\star }}}\left({\varphi }_{A^{\star }}-{\varphi }_{B^{\star }}\right)\\ =meas\left(C_{B^{\text{*}}}\right)f_{B^{\text{*}}}+{\displaystyle \sum }_{I\in {ℰ}^{B^{\star }}\cap {ℰ}^{ext}}\frac{h_{IP}{b}_h\left(A^P\right)}{h_{A^{\star }{B}^{\star }}{a}_h\left(A^P\right)}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}\gamma }\\ +{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{ext}}{\displaystyle {\int }_{C_{B^{\text{*}}}\cap \partial \text{Ω}}g\text{d}\gamma }\forall B^{\text{*}}\in \mathcal{D}\end{array}$ (2.39)

Remark 2.8. When $\alpha \to \infty$ and $g=0$ we recognize the balance equation of diffusion reaction problem with homogeneous Dirichlet bounbary conditions which is written as follows

(2.40)

$\begin{array}{l}{\displaystyle \sum }_{I\in {ℰ}^{B^{\text{*}}}\cap {ℰ}^{int}}\{\frac{{\widehat{b}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}+b_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{PI}}{\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_P-{\varphi }_L\right)\\ +\frac{w_h\left(P,L,I\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left[h_{IL}{a}_h\left(A^P\right)+h_{IP}{\widehat{a}}_h\left(A^L\right)\right]}\left({\varphi }_{B^{\text{*}}}-{\varphi }_{A^{\text{*}}}\right)\}\\ +c{\displaystyle \sum _{B^{\star }\in C_{B^{\star }}}meas\left(C_{B^{\star }}\right){\varphi }_{B^{\star }}}=meas\left(C_{B^{\text{*}}}\right)f_{B^{\text{*}}}\forall B^{\text{*}}\in \mathcal{D}\end{array}$(2.41)

3. Existence and Uniqueness for a Discrete Scheme

Assume that all the cellpoints and all the vertices (with respect to the primary mesh) are numbered. On the other hand and $Card\left(\mathcal{D}\right)$ denote respectively the total number of cellpoints and vertices. A preliminary result is:

Proposition 3.1. Under the assumption (1.3)-(1.4), the matrix associated with the linear system (2.36)-(2.37) is symmetric and positive definite.

Proposition 3.2. Under the assumption (1.3)-(1.4), the linear system (2.36)-(2.37) possesses a unique solution.

It is sufficient to prove Proposition 3.1 as it implies Proposition 3.2.

Proof. It is easily seen that the symmetry of the matrix $M$ associated with the linear system (2.36)-(2.37) essentially follows from the symmetry of the diffusion coefficient $A$ . We now should prove the positive definiteness of that matrix.

Let us set:

(3.1)

$M=\left[\begin{array}{ll}\text{A}\hfill & \text{B}\hfill \\ {\text{B}}^t\hfill & \text{D}\hfill \end{array}\right]$ (3.2)

where

$\text{M}\left[\begin{array}{l}{\varphi }_{cc}\hfill \\ {\varphi }_{vc}\hfill \end{array}\right]=\left[\begin{array}{l}{\text{F}}_{cc}\hfill \\ {\text{F}}_{vc}\hfill \end{array}\right]$ (3.3)

where $\text{A}$ is a symmetric matrix, $\text{D}$ is a $Card\left(\mathcal{D}\right)\times Card\left(\mathcal{D}\right)$ symmetric matrix, $\text{B}$ is a matrix and ${\text{B}}^t$ is the transpose of $\text{B}$ . In this context, ${\text{F}}_{cc}$ represents a subvector with a number of components equal to derive from the right-hand side of Equation (2.36). Similarly, ${\text{F}}_{vc}$ is a subvector with $Card\left(\mathcal{D}\right)$ components obtained from the right-hand side of Equation (2.37). We start by developing the quadratic expression

$\left[{\Phi }_{cc},{\Phi }_{vc}\right]M\left[\begin{array}{l}{\Phi }_{cc}\hfill \\ {\Phi }_{vc}\hfill \end{array}\right]$

of the system (2.36)-(2.37). It is clear that

(3.4)

where $D$ is a set diamond cells (see Figure 6) and ${\mathfrak{D}}_{\sigma {\sigma }^{\star }}$ a diamond cell whose diagonals are $\sigma =\left[A^{\text{*}}{B}^{\text{*}}\right]$ and ${\sigma }^{\text{*}}=\left[PL\right]$ . This diamond cell is widely used to establish the stability result and error estimates. Set

$A_{11}^{PL}=\frac{a_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{A^{\text{*}}{B}^{\text{*}}}}{a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^P\right)h_{PI}}$

$A_{22}^{PL}=\frac{w_h\left(P,L,I\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left[a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{IP}\right]}$

$A_{12}^{PL}=A_{21}^{PL}=\frac{{\widehat{c}}_h\left(A^L\right)a_h\left(A^P\right)h_{IL}+c_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{PI}}{a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{IP}}$

$A_{11}^{PI}=\frac{\alpha a_h\left(A^P\right)h_{A^{\text{*}}{B}^{\text{*}}}}{\alpha h_{IP}+a_h\left(A^P\right)}$

$A_{22}^{PI}=\frac{h_{IP}\left(\alpha d_h\left(A^P\right)h_{IP}+\left(a_h\left(A^P\right)d_h\left(A^P\right)-b_h^2\left(A^P\right)\right)\right)}{h_{A^{\text{*}}{B}^{\text{*}}}\left(\alpha h_{IP}+a_h\left(A^P\right)\right)}$

$A_{21}^{PI}=\frac{2\alpha h_{IP}{b}_h\left(A^P\right)}{\alpha h_{IP}+a_h\left(A^P\right)}$

where are two adjacent primary cells sharing the edge $\left[A^{\text{*}}{B}^{\text{*}}\right]$ as interface containing the edge-point $I$ . Let us prove that the homogenized symmetric permeability tensor $A^{PL}$ is positive definite, that is $A_{11}^{PL}{A}_{22}^{PL}-{\left(A_{21}^{PL}\right)}^2>0$ . Setting

${△}^{PL}=A_{11}^{PL}{A}_{22}^{PL}-{\left(A_{21}^{PL}\right)}^2,$ (3.5)

we have that

$\begin{array}{l}{\Delta }^{PL}=N_1\left[a_h\left(A^P\right)d_h\left(A^P\right)-{\left(b_h\left(A^P\right)\right)}^2\right]+N_2\left[{\widehat{a}}_h\left(A^L\right){\widehat{d}}_h\left(A^L\right)-{\left(\widehat{b}\left(A^L\right)\right)}^2\right]\hfill \end{array}$ (3.6)

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

$\begin{array}{l}{N}_1={\left(\frac{{\widehat{a}}_h\left(A^L\right)h_{IL}}{a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{IP}}\right)}^2+\frac{a_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{IP}{h}_{IL}}{{\left[a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{IP}\right]}^2}\\ \text{and}\\ N_2={\left(\frac{a_h\left(A^P\right)h_{IP}}{a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{IP}}\right)}^2+\frac{a_h\left(A^P\right){\widehat{a}}_h\left(A^L\right)h_{IP}{h}_{IL}}{{\left[a_h\left(A^P\right)h_{IL}+{\widehat{a}}_h\left(A^L\right)h_{IP}\right]}^2}\end{array}$(3.7)

Given that the diffusion coefficient $A$ adhere to assumptions (1.3)-(1.4), symmetric and positive definite, Cauchy-Schwarz’s inequality applied to the inner product associated with $A$ guarantees that

$\begin{array}{l}{a}_h\left(A^P\right)d_h\left(A^P\right)-{\left(b_h\left(A^P\right)\right)}^2>0\\ \text{and}\\ {\widehat{a}}_h\left(A^L\right){\widehat{d}}_h\left(A^L\right)-{\left(\widehat{b}\left(A^L\right)\right)}^2>0\end{array}$(3.8)

Since either ${\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P$ and ${\xi }_{\left[PI\right]}^{B^{\text{*}}}$ or ${\xi }_{\left[A^{\text{*}}{B}^{\text{*}}\right]}^P$ and ${\xi }_{\left[IL\right]}^{B^{\text{*}}}$ are not collinear, it follows that ${\Delta }^{PL}>0$ . Consequently, $A^{PL}$ forms a symmetric positive definite matrix. For the second quadratic form, set

${\Delta }^{PI}=\left[\begin{array}{ll}{A}_{11}^{PI}\hfill & A_{21}^{PI}\hfill \\ A_{21}^{PI}\hfill & A_{22}^{PI}\hfill \end{array}\right]$ (3.9)

The second quadratic form is positive if

$\text{det}\left({\Delta }^{PI}\right)=A_{11}^{PI}{A}_{22}^{PI}-{\left(A_{21}^{PI}\right)}^2>0$ (3.10)

We now show that the bilinear form is quadratic and elliptic. That is

$\text{det}\left({\Delta }^{PI}\right)=\left({\alpha }^2{h}_{IP}^2+\alpha a_h\left(A^P\right)h_{IP}\right)\left(a_h\left(A^P\right)d_h\left(A^P\right)-b_h^2\left(A^P\right)\right)\ge 0$ (3.11)

The Cauchy-Schwarz inequality for the inner product associated with $A$ ensure that relation (3.11) is positive.

It follows from the previous step that the matrices $A^{PL}$ and $A^{PI}$ have respectively ${\lambda }_{min}^{PL}$ and ${\lambda }_{min}^{PI}$ as strictly positive eigenvalues. So, we have

(3.12)

The equality holds in Equation (3.12) if and only if ${\Phi }_{cc}=0$ and ${\Phi }_{vc}=0$ . Thus, the positive definiteness of the matrix $M$ is proven.

4. DDFV Formulation of the Problem in Uniform Grid

In this section, we derive our discrete problem via the discrete balance equations over the primary and dual cell interfaces. We consider the domain $\Omega =\left]0,1\right[\times \left]0,1\right[$ covered with a primary mesh whose grid size is $h=\frac{1}{N}$ , where $N$ is a given strictly positive integer. We denote by $K_{ij}$ the primary grid block given by

$K_{ij}=\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]$ where

$x_{i+\frac{1}{2}}=x_{i-\frac{1}{2}}+h$ , $y_{j+\frac{1}{2}}=y_{j-\frac{1}{2}}+h$ for $i,j=1,2,3,\cdots ,N-1$ .

We also consider a second grid block $K_{i+\frac{1}{2},j+\frac{1}{2}}$ centered in $\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)$ and called dual grid block and defined by $K_{i+\frac{1}{2},j+\frac{1}{2}}=\left[x_{i-1},x_{i+1}\right]\times \left[y_{j-1},y_{j+1}\right]$ where $x_{i+1}=x_i+h$ , $y_{j+1}=y_j+h$ for $i,j=1,2,3,\cdots ,N-1$ For any cell along the boundary of the domain, we adopt the convention $x_{\frac{1}{2}}=y_{\frac{1}{2}}=0$ , $x_{N+\frac{1}{2}}=y_{N+\frac{1}{2}}=1$ and $x_0=y_0=0$ , $x_{N+1}=y_{N+1}=1$ . We now focus our attention on finding for the unknowns ${\left({\phi }_{i,j}\right)}_{1\le i,j\le N}$ and ${\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right)}_{0\le i,j\le N}$ which are reasonable approximate solutions of Equation (1.1)-(1.2), where ${\phi }_{i,j}=\phi \left(x_i,y_j\right)$ and ${\phi }_{i+\frac{1}{2},j+\frac{1}{2}}=\phi \left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)$ . Figure 7 is a computational grid,

where the local coordinates are given by the integer and fractional indices. We give now a description of the procedure leading to the linear discrete system. Integrating the balance Equation (1.1)-(1.2) in the grid-block $K_{i,j}$ and $K_{i+\frac{1}{2},j+\frac{1}{2}}$ respectively, applying Ostrogradski’s theorem leads to compute fluxes through the boundary of $K_{i,j}$ and $K_{i+\frac{1}{2},j+\frac{1}{2}}$ as illustrated by Figure 7. The integration in primary cell $K_{ij}$ is a direct result of taking the interface flux values over $i\pm \frac{1}{2},j$ and $i,j\pm \frac{1}{2}$ as you can see in the following formulas.

Figure 7. 2D computational grid in uniform mesh with primal cells (black lines) and dual cells (black discontinuous lines).

$\begin{array}{l}{\displaystyle {\int }_{K_{i,j}}-div\left(A^{ij}\nabla \phi \right)\text{d}x}={\displaystyle \sum _{\sigma \in \partial K_{i,j}}{\displaystyle {\int }_{\sigma }-A^{ij}\nabla \phi \cdot {\eta }_{ij}^{\sigma }\text{d}\sigma }}\\ =F_{i+\frac{1}{2},j}^{K_{i,j}}+F_{i-\frac{1}{2},j}^{K_{i,j}}+F_{i,j+\frac{1}{2}}^{K_{i,j}}+F_{i,j-\frac{1}{2}}^{K_{i,j}}\forall 1\le i,j\le N\end{array}$ (4.1)

Likewise the integration in primary cell $K_{i+\frac{1}{2},j+\frac{1}{2}}$ is a direct result of taking the flux values over $i+\frac{1}{2},j+1$ , $i+\frac{1}{2},j$ , $i+1,j+\frac{1}{2}$ and $i,j+\frac{1}{2}$ as you can see in the following formulas.

$\begin{array}{l}{\displaystyle {\int }_{K_{i+\frac{1}{2},j+\frac{1}{2}}}-div\left(A^{i+\frac{1}{2},j+\frac{1}{2}}\nabla \phi \right)\text{d}x}={\displaystyle \sum _{{\sigma }^{\star }\in \partial K_{i+\frac{1}{2},j+\frac{1}{2}}}{\displaystyle {\int }_{{\sigma }^{\star }}-A^{i+\frac{1}{2},j+\frac{1}{2}}\nabla \phi \cdot {\eta }_{i+\frac{1}{2},j+\frac{1}{2}}^{{\sigma }^{\star }}\text{d}{\sigma }^{\star }}}\\ =F_{i+\frac{1}{2},j+1}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+F_{i+\frac{1}{2},j}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+F_{i+1,j+\frac{1}{2}}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+F_{i,j+\frac{1}{2}}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}\forall 0\le i,j\le N\end{array}$ (4.2)

4.1. Computation of Fluxes

Let us consider the segment $\left[A^{\star }{B}^{\star }\right]$ delimited by points $x_{i-\frac{1}{2},j+\frac{1}{2}}$ and $x_{i+\frac{1}{2},j+\frac{1}{2}}$ (see Figure 7), le flux the edge $\left[A^{\star }{B}^{\star }\right]$ seen as element of the primal cell $K_{ij}$ and denoted by $F_{i,j+\frac{1}{2}}^{K_{i,j}}$ is calculated as follows.

$\begin{array}{c}{F}_{i,j+\frac{1}{2}}^{K_{i,j}}={\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}-A^{i,j}\nabla \phi \cdot {\eta }_{ij}^{i,j+\frac{1}{2}}\text{d}s}\\ ={\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}\langle -\left[\begin{array}{ll}{A}_{11}^{i,j}\hfill & A_{12}^{i,j}\hfill \\ A_{21}^{i,j}\hfill & A_{22}^{i,j}\hfill \end{array}\right]\left[\begin{array}{l}\frac{\partial \phi }{\partial x}\hfill \\ \frac{\partial \phi }{\partial y}\hfill \end{array}\right],\left(\begin{array}{l}0\hfill \\ 1\hfill \end{array}\right)\rangle \text{d}x}\\ =-{\displaystyle {\int }_{\left[A^{\star }I\right]}\left(A_{21}^{i,j}\frac{\partial \phi }{\partial x}+A_{22}^{i,j}\frac{\partial \phi }{\partial y}\right)\text{d}x}-{\displaystyle {\int }_{\left[I{B}^{\star }\right]}\left(A_{21}^{i,j}\frac{\partial \phi }{\partial x}+A_{22}^{i,j}\frac{\partial \phi }{\partial y}\right)\text{d}x}\end{array}$

where ${\eta }_{ij}^{i,j+\frac{1}{2}}$ is the unit normal vector exterior from $K_{i,j}$ and $I$ the midpoint or interface point of segment $\left[A^{\star }{B}^{\star }\right]$ . We now use suitable gradient operators to compute the flux $F_{i,j+1/2}^{K_{ij}}$ as follows

$\begin{array}{c}{F}_{i,j+\frac{1}{2}}^{K_{i,j}}=-\frac{h}{2}\left(\frac{{\phi }_{i,j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}}{\frac{h}{2}}\right)A_{21}^{i,j}-\frac{h}{2}\left(\frac{{\phi }_{i,j+\frac{1}{2}}-{\phi }_{i,j}}{\frac{h}{2}}\right)A_{22}^{i,j}\\ -\frac{h}{2}\left(\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i,j+\frac{1}{2}}}{\frac{h}{2}}\right)A_{21}^{i,j}-\frac{h}{2}\left(\frac{{\phi }_{i,j+\frac{1}{2}}-{\phi }_{i,j}}{\frac{h}{2}}\right)A_{22}^{i,j}+T_{i,j+\frac{1}{2}}^{K_{i,j}}\\ =-A_{21}^{i,j}\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}\right)+2A_{22}^{i,j}\left({\phi }_{i,j}-{\phi }_{i,j+\frac{1}{2}}\right)+T_{i,j+\frac{1}{2}}^{K_{i,j}}\end{array}$ (4.3)

$T_{i,j+\frac{1}{2}}^{K_{i,j}}$ represents the truncation error so that $\left|T_{i,j+\frac{1}{2}}^{K_{i,j}}\right|\le Ch$ , where C is a mesh independent positive constant that depends exclusively on Ω. Considering the principle of flux continuity across the interface joining two adjacent cells $K_{i,j}$ and $K_{i,j+1}$ we have that,

$F_{i,j+\frac{1}{2}}^{K_{i,j}}=F_{i,j+\frac{1}{2}}^{K_{i,j+1}}$

where, $F_{i,j+1/2}^{K_{i,j+1}}$ represents the flux over $K_{ij+1}$ . We therefore obtain that the edge midpoint potential ${\phi }_{i,j+1/2}$ is given by

${\phi }_{i,j+\frac{1}{2}}=\frac{1}{A_{22}^{i,j}+A_{22}^{i,j+1}}\left\{\frac{A_{21}^{i,j+1}-A_{21}^{i,j}}{2}\left({\phi }_{i-\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right)+A_{22}^{i,j}{\phi }_{i,j}+A_{22}^{i,j+1}{\phi }_{i,j+1}\right\}$ (4.4)

Substituting the edge midpoint potential ${\phi }_{i,j+1/2}$ in the flux equation $F_{i,j+\frac{1}{2}}^{K_{i,j}}$ we obtain that

$F_{i,j+\frac{1}{2}}^{K_{i,j}}=\frac{2A_{22}^{i,j}{A}_{22}^{i,j+1}}{A_{22}^{i,j}+A_{22}^{i,j+1}}\left[{\phi }_{i,j}-{\phi }_{i,j+1}\right]+\frac{A_{22}^{i,j}{A}_{21}^{i,j+1}+A_{22}^{i,j+1}{A}_{21}^{i,j}}{A_{22}^{i,j}+A_{22}^{i,j+1}}\left[{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right]$

If we perform similar calculations on the neighboring primal cells of $K_{i,j}$ that is on the cells, $K_{i+1,j}$ , $K_{i-1,j}$ , $K_{i,j+1}$ . We obtain the following balance equation in internal primary cell.

$\begin{array}{l}\frac{2A_{22}^{i,j}{A}_{22}^{i,j+1}}{A_{22}^{i,j}+A_{22}^{i,j+1}}\left[{\phi }_{i,j}-{\phi }_{i,j+1}\right]+\frac{A_{22}^{i,j}{A}_{21}^{i,j+1}+A_{22}^{i,j+1}{A}_{21}^{i,j}}{A_{22}^{i,j}+A_{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{2A_{22}^{i,j}{A}_{22}^{i,j-1}}{A_{22}^{i,j}+A_{22}^{i,j-1}}\left[{\phi }_{i,j}-{\phi }_{i,j-1}\right]+\frac{A_{22}^{i,j}{A}_{21}^{i,j-1}+A_{22}^{i,j-1}{A}_{21}^{i,j}}{A_{22}^{i,j}+A_{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{2A_{11}^{i,j}{A}_{11}^{i-1,j}}{A_{11}^{i,j}+A_{11}^{i-1,j}}\left[{\phi }_{i,j}-{\phi }_{i-1,j}\right]+\frac{A_{11}^{i,j}{A}_{12}^{i-1,j}+A_{11}^{i-1,j}{A}_{12}^{i,j}}{A_{11}^{i,j}+A_{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{2A_{11}^{i,j}{A}_{11}^{i+1,j}}{A_{11}^{i,j}+A_{11}^{i+1,j}}\left[{\phi }_{i,j}-{\phi }_{i+1,j}\right]+\frac{A_{11}^{i,j}{A}_{21}^{i+1,j}+A_{11}^{i+1,j}{A}_{21}^{i,j}}{A_{11}^{i,j}+A_{11}^{i+1,j}}\left[{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right]\\ +meas\left(K_{i,j}\right)c_{ij}{\phi }_{i,j}=meas\left(K_{i,j}\right)f_{i,j}\forall 2\le i,j\le N\end{array}$ (4.5)

We develop in order to close system 4.1, the balance equations over the interface of the interior dual cell $K_{i+\frac{1}{2},j+\frac{1}{2}}$ (see Figure 7 for clarity), Following the same procedure as shown above when computing the balance equation for the primal cells using the dual cell $K_{i+\frac{1}{2},j+\frac{1}{2}}$ and the neighboring cells $K_{i+\frac{1}{2},j-\frac{1}{2}}$ ,$K_{i+\frac{1}{2},j+\frac{3}{2}}$ , $K_{i+\frac{3}{2},j+\frac{1}{2}}$ and $K_{i-\frac{1}{2},j+\frac{1}{2}}$ , we obtain the balance equation.

$\begin{array}{l}\frac{A_{11}^{i,j+1}{A}_{21}^{i+1,j+1}+A_{11}^{i+1,j+1}{A}_{21}^{i,j+1}}{A_{11}^{i,j+1}+A_{11}^{i+1,j+1}}\left[{\phi }_{i,j+1}-{\phi }_{i+1,j+1}\right]+\frac{A_{11}^{i+1,j}{A}_{21}^{i,j}+A_{11}^{i,j}{A}_{21}^{i+1,j}}{A_{11}^{i,j}+A_{11}^{i+1,j}}\left[{\phi }_{i+1,j}-{\phi }_{i,j}\right]\\ +\frac{A_{22}^{i+1,j}{A}_{12}^{i+1,j+1}+A_{22}^{i+1,j+1}{A}_{21}^{i+1,j}}{A_{22}^{i+1,j}+A_{22}^{i+1,j+1}}\left[{\phi }_{i+1,j}-{\phi }_{i+1,j+1}\right]+\frac{A_{22}^{i,j+1}{A}_{12}^{i,j}+A_{22}^{i,j}{A}_{12}^{i,j+1}}{A_{22}^{i,j}+A_{22}^{i,j+1}}\left[{\phi }_{i,j+1}-{\phi }_{i,j}\right]\\ +\left(\frac{A_{22}^{i,j+1}+A_{22}^{i+1,j+1}}{2}-\frac{{\left(A_{12}^{i+1,j+1}-A_{12}^{i,j+1}\right)}^2}{2\left(A_{11}^{i,j+1}+A_{11}^{i+1,j+1}\right)}\right)\left[{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j+\frac{3}{2}}\right]\\ +\left(\frac{A_{22}^{i,j}+A_{22}^{i+1,j}}{2}-\frac{{\left(A_{12}^{i+1,j}-A_{12}^{i,j}\right)}^2}{2\left(A_{11}^{i,j}+A_{11}^{i+1,j}\right)}\right)\left[{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}\right]\\ +\left(\frac{A_{11}^{i+1,j+1}+A_{11}^{i+1,j}}{2}-\frac{{\left(A_{21}^{i+1,j+1}-A_{21}^{i+1,j}\right)}^2}{2\left(A_{11}^{i+1,j+1}+A_{11}^{i+1,j}\right)}\right)\left[{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{3}{2},j+\frac{1}{2}}\right]\\ +\left(\frac{A_{11}^{i,j}+A_{11}^{i,j+1}}{2}-\frac{{\left(A_{21}^{i,j+1}-A_{21}^{i,j}\right)}^2}{2\left(A_{22}^{i,j+1}+A_{22}^{i,j}\right)}\right)\left[{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}\right]\\ +meas\left(K_{i+\frac{1}{2},j+\frac{1}{2}}\right)c_{i+\frac{1}{2},j+\frac{1}{2}}{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\\ =meas\left(K_{i+\frac{1}{2},j+\frac{1}{2}}\right)f_{i+\frac{1}{2},j+\frac{1}{2}}\forall 1\le i,j\le N-1\end{array}$ (4.6)

4.2. DDFV Flux Approximation along the Boundary Cells

It clearly results in Equations (4.1) and (4.2) that the case $j=N$ requires special treatment for example, due to the fact that the associated dual and primal cells meet the Robin boundary conditions. These are the fluxes $F_{i,N+\frac{1}{2}}^{K_{i,N}}$ , $F_{i+\frac{1}{2},N}^{K_{i,N}}$ , $F_{i,N+\frac{1}{2}}^{K_{i+\frac{1}{2},N+\frac{1}{2}}}$ and $F_{i+1,N+\frac{1}{2}}^{K_{i+\frac{1}{2},N+\frac{1}{2}}}$ as illustrated in Figure 8.

Figure 8. 2D computational grid along the boundary cells (black lines) and dual cells (hatted in black discontinuous lines).

How can we approximate fluxes $F_{i,N+\frac{1}{2}}^{K_{i,N}}$ and $F_{i,N+\frac{1}{2}}^{K_{i+\frac{1}{2},N+\frac{1}{2}}}$ ?. Recall from equation (1.2) that:

$-\alpha \varphi -\left(n\cdot A\nabla \varphi \right)=-g$

It follows by integration of both sides of the previous equation that:

$F_{i,N+\frac{1}{2}}^{K_{i,N}}=\alpha h{\phi }_{i,N+\frac{1}{2}}-{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s}$ (4.7)

where

$F_{i,N+\frac{1}{2}}^{K_{i,N}}=2A_{22}^{i,N}\left({\phi }_{i,N}-{\phi }_{i,N+\frac{1}{2}}\right)+A_{12}^{i,N}\left({\phi }_{i-\frac{1}{2},N+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},N+\frac{1}{2}}\right)$ (4.8)

By combining the two previous equations, we obtain edgepoint pressure on the boundary ${\phi }_{i,N+\frac{1}{2}}$ , that is:

$\begin{array}{c}{\phi }_{i,N+\frac{1}{2}}=\frac{2A_{22}^{i,N}}{\alpha h+2A_{22}^{i,N}}{\phi }_{i,N}+\frac{A_{12}^{i,N}}{\alpha h+2A_{22}^{i,N}}\left({\phi }_{i-\frac{1}{2},N+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},N+\frac{1}{2}}\right)\\ +\frac{1}{\alpha h+2A_{22}^{i,N}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s}\end{array}$

Therefore the fluxes $F_{i,N+\frac{1}{2}}^{K_{i,N}}$ and $F_{i,N+\frac{1}{2}}^{K_{i+\frac{1}{2},N+\frac{1}{2}}}$ are given by the following formulas:

$\begin{array}{l}{F}_{i,N+\frac{1}{2}}^{K_{i,N}}=\frac{2A_{22}^{i,N}\alpha h}{\alpha h+2A_{22}^{i,N}}{\phi }_{i,N}+\frac{2A_{12}^{i,N}\alpha h}{\alpha h+2A_{22}^{i,N}}\left({\phi }_{i-\frac{1}{2},N+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},N+\frac{1}{2}}\right)\\ -\frac{2A_{22}^{i,N}}{\alpha h+2A_{22}^{i,N}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s},\text{for}i=1,2,3,\cdots ,N\end{array}$ (4.9)

$\begin{array}{l}{F}_{i,N+\frac{1}{2}}^{K_{i+\frac{1}{2},N+\frac{1}{2}}}=\left[\frac{A_{11}^{i,N}}{2}-\frac{{\left(A_{12}^{i,N}\right)}^2}{\alpha h+2A_{22}^{i,N}}\right]\left({\phi }_{i+\frac{1}{2},N+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},N+\frac{1}{2}}\right)-\frac{2A_{12}^{i,N}\alpha h}{\alpha h+2A_{22}^{i,N}}{\phi }_{i,N}\\ +\frac{2A_{12}^{i,N}}{\alpha h+2A_{22}^{i,N}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s},\text{for}i=1,2,3,\cdots ,N\end{array}$ (4.10)

We can follow the same procedure to derive the other boundary fluxes over degenerated dual and primary cells. Then we have the following

$\begin{array}{l}{F}_{i,\frac{1}{2}}^{K_{i,1}}=\frac{2A_{22}^{i,1}\alpha h}{\alpha h+2A_{22}^{i,1}}{\phi }_{i,1}+\frac{2A_{21}^{i,1}\alpha h}{\alpha h+2A_{22}^{i,1}}\left({\phi }_{i+\frac{1}{2},\frac{1}{2}}-{\phi }_{i-\frac{1}{2},\frac{1}{2}}\right)\\ -\frac{2A_{22}^{i,1}}{\alpha h+2A_{22}^{i,N}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s}\text{for}i=1,2,3,\cdots ,N\end{array}$ (4.11)

$\begin{array}{l}{F}_{i,\frac{1}{2}}^{K_{i-\frac{1}{2},\frac{1}{2}}}=\left[\frac{A_{11}^{i,1}}{2}-\frac{{\left(A_{12}^{i,1}\right)}^2}{\alpha h+2A_{22}^{i,1}}\right]\left({\phi }_{i-\frac{1}{2},\frac{1}{2}}-{\phi }_{i+\frac{1}{2},\frac{1}{2}}\right)+\frac{2A_{12}^{i,1}\alpha h}{\alpha h+2A_{22}^{i,1}}{\phi }_{i,1}\\ +\frac{2A_{12}^{i,1}}{\alpha h+2A_{22}^{i,1}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s},\text{for}i=1,2,3,\cdots ,N\end{array}$ (4.12)

$\begin{array}{l}{F}_{N+\frac{1}{2},j}^{K_{N,j}}=\frac{2A_{12}^{N,j}\alpha h}{\alpha h+2A_{11}^{i,1}}\left({\phi }_{N+\frac{1}{2},j-\frac{1}{2}}-{\phi }_{N+\frac{1}{2},j+\frac{1}{2}}\right)+\frac{2A_{11}^{N,j}\alpha h}{\alpha h+2A_{11}^{N,j}}{\phi }_{N,j}\\ -\frac{2A_{22}^{N,j}}{\alpha h+2A_{11}^{N,j}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s}\text{for}j=1,2,3,\cdots ,N\end{array}$ (4.13)

$\begin{array}{l}{F}_{i,\frac{1}{2}}^{K_{N+\frac{1}{2},j+\frac{1}{2}}}=\left[\frac{A_{22}^{N,j}}{2}-\frac{{\left(A_{12}^{N,j}\right)}^2}{\alpha h+2A_{11}^{N,j}}\right]\left({\phi }_{N+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{N+\frac{1}{2},j-\frac{1}{2}}\right)-\frac{2A_{21}^{N,j}\alpha h}{\alpha h+2A_{11}^{N,j}}{\phi }_{N,j}\\ +\frac{2A_{11}^{N,j}}{\alpha h+2A_{11}^{N,j}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s},\text{for}j=1,2,3,\cdots ,N\end{array}$ (4.14)

$\begin{array}{l}{F}_{\frac{1}{2},j}^{K_{1,j}}=\frac{2A_{11}^{1,j}\alpha h}{\alpha h+2A_{11}^{1,j}}{\phi }_{1,j}+\frac{2A_{12}^{1,j}\alpha h}{\alpha h+2A_{11}^{1,j}}\left({\phi }_{\frac{1}{2},j+\frac{1}{2}}-{\phi }_{\frac{1}{2},j-\frac{1}{2}}\right)\\ -\frac{2A_{11}^{1,j}}{\alpha h+2A_{11}^{1,j}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s}\text{for}j=1,2,3,\cdots ,N\end{array}$ (4.15)

$\begin{array}{l}{F}_{i,\frac{1}{2}}^{K_{\frac{1}{2},j-\frac{1}{2}}}=\left[\frac{A_{22}^{1,j}}{2}-\frac{{\left(A_{12}^{1,j}\right)}^2}{\alpha h+2A_{11}^{1,j}}\right]\left({\phi }_{\frac{1}{2},j-\frac{1}{2}}-{\phi }_{\frac{1}{2},j+\frac{1}{2}}\right)-\frac{2A_{21}^{1,j}\alpha h}{\alpha h+2A_{11}^{1,j}}{\phi }_{1,j}\\ +\frac{2A_{21}^{1,j}}{\alpha h+2A_{11}^{N,j}}{\displaystyle {\int }_{\left[A^{\star }{B}^{\star }\right]}g\text{d}s},\text{for}j=1,2,3,\cdots ,N\end{array}$ (4.16)

Now all the ingredients are gathered for define our DDFV algebraic linear system assciated with the governing equations (1.1) and (1.2) in uniform mesh

Proposition 4.1. The DDFV scheme associated of the governing equations (1.1) and (1.2)

We define our DDFV of the problem as find ${\Phi }^h=\left({\left({\phi }_{i,j}\right)}_{1\le i,j\le N},{\left({\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\right)}_{0\le i,j\le N}\right)$ such that

$\begin{array}{l}{F}_{i+\frac{1}{2},j}^{K_{ij}}+F_{i-\frac{1}{2},j}^{K_{ij}}+F_{i,j+\frac{1}{2}}^{K_{ij}}+F_{i,j-\frac{1}{2}}^{K_{ij}}+meas\left(K_{ij}\right)c_{i,j}{\phi }_{ij}\\ =meas\left(K_{ij}\right)f_{ij}\text{for}i,j\in \left\{1,2,3,\cdots ,N\right\}\end{array}$ (4.17)

$\begin{array}{l}{F}_{i+\frac{1}{2},j+1}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+F_{i+\frac{1}{2},j}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+F_{i+1,j+\frac{1}{2}}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+F_{i,j+\frac{1}{2}}^{K_{i+\frac{1}{2},j+\frac{1}{2}}}+meas\left(K_{i+\frac{1}{2},j+\frac{1}{2}}\right)c_{i+\frac{1}{2},j+\frac{1}{2}}{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}\\ =meas\left(K_{i+\frac{1}{2},j+\frac{1}{2}}\right)f_{i+\frac{1}{2},j+\frac{1}{2}}\text{for}i,j\in \left\{0,1,2,3,\cdots ,N\right\}\end{array}$ (4.18)

where,

$F_{0,j+\frac{1}{2}}^{K_{\frac{1}{2},j+\frac{1}{2}}}=F_{i+\frac{1}{2},N+1}^{K_{i+\frac{1}{2},N+\frac{1}{2}}}=F_{i+\frac{1}{2},0}^{K_{i+\frac{1}{2},\frac{1}{2}}}=0\text{for}i,j\in \left\{0,1,2,3,\cdots ,N\right\}$

Which can be written for $g=0$ as

$\begin{array}{l}{\alpha }_{i,j}{\phi }_{i,j}+{\alpha }_{i,j+1}{\phi }_{i,j+1}+{\alpha }_{i,j-1}{\phi }_{i,j-1}+{\alpha }_{i-1,j}{\phi }_{i-1,j}+{\alpha }_{i+1,j}{\phi }_{i+1,j}\\ +{\alpha }_{i+\frac{1}{2},j+\frac{1}{2}}{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}+{\alpha }_{i-\frac{1}{2},j+\frac{1}{2}}{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}+{\alpha }_{i+\frac{1}{2},j-\frac{1}{2}}{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}+{\alpha }_{i-\frac{1}{2},j-\frac{1}{2}}{\phi }_{i-\frac{1}{2},j-\frac{1}{2}}\\ =meas\left(K_{i,j}\right)f_{i,j}\text{for}i,j\in \left\{1,2,3,\cdots ,N\right\}\end{array}$ (4.19)

$\begin{array}{l}{\beta }_{i+\frac{1}{2},j+\frac{1}{2}}{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}+{\beta }_{i+\frac{1}{2},j+\frac{3}{2}}{\phi }_{i+\frac{1}{2},j+\frac{3}{2}}+{\beta }_{i+\frac{1}{2},j-\frac{1}{2}}{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}+{\beta }_{i-\frac{1}{2},j+\frac{1}{2}}{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}\\ +{\beta }_{i+1,j+1}{\phi }_{i+1,j+1}+{\beta }_{i,j+1}{\phi }_{i,j+1}+{\beta }_{i+1,j}{\phi }_{i+1,j}+{\beta }_{i,j-1}{\phi }_{i,j-1}+{\beta }_{i,j}{\phi }_{i,j}\\ =meas\left(K_{i+\frac{1}{2},j+\frac{1}{2}}\right)f_{i+\frac{1}{2},j+\frac{1}{2}}\text{for}i,j\in \left\{0,1,\cdots ,N\right\}\end{array}$ (4.20)

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 an isotropic system, 4.19 and 4.20 are decoupled and therefore independent. In the homogeneous and nonhomogeneous anisotropic case there is a connectivity between the unknown ${\phi }_{i+\frac{1}{2},j+\frac{1}{2}}$ and ${\phi }_{i,j}$ . It is clear that the total number of unknows is $N^2+{\left(N+1\right)}^2$ . Depending on the numbering process (see Figure 9), the general form of the associated linear system is as follows

Figure 9. Node numbering (mesh size $N=4$ ) from primal cell points, left to right, bottom to top, to dual cell points left to right and bottom to top.

$\left(\begin{array}{l}A B^t\hfill \\ B C\hfill \end{array}\right)\left(\begin{array}{l}{\Phi }_{cc}\hfill \\ {\Phi }_{cv}\hfill \end{array}\right)=\left(\begin{array}{l}{F}_p\hfill \\ F_v\hfill \end{array}\right)$

With,

${\Phi }_{cc}={\left\{{\phi }_{\left(j-1\right)N+i}\right\}}_{1\le i,j\le N}$

${\Phi }_{cv}={\left\{{\phi }_{j\left(N+1\right)+i+1+N^2}\right\}}_{0\le i,j\le N}$

$F_{cc}={\left\{F_{\left(j-1\right)N+i}\right\}}_{1\le i,j\le N}$

$F_{cv}={\left\{F_{j\left(N+1\right)+i+1+N^2}\right\}}_{0\le i,j\le N}$

where,

$A$ is a $N^2\times N^2$ symmetric and positive define matrix, $B$ is a $N^2\times {\left(N+1\right)}^2$ , $C$ is a ${\left(N+1\right)}^2\times {\left(N+1\right)}^2$ symmetric and positive define matrix and $B^t$ the transpose of $B$ .

4.3. Implementation Steps

In this part, we briefly present the basis data structure of our problem, we discuss what input data are required in general, how matrix equation is assembled, how the boundary conditions are applied in the matrix equation. The key of implementation is to establish a bijective relationship between the nodes of the grid and the pairs of integer numbers $\left(i,j\right)$ for $i,j=1,2,3,\cdots ,N$ such that given the node indices $\left(i,j\right)$ and his eigth neighbors $\left(i+1,j\right)$ , $\left(i-1,j\right)$

$\left(i,j-1\right)$ , $\left(i,j+1\right)$ , $\left(i-\frac{1}{2},j-\frac{1}{2}\right)$ , $\left(i+\frac{1}{2},j-\frac{1}{2}\right)$ , $\left(i+\frac{1}{2},j+\frac{1}{2}\right)$ , $\left(i-\frac{1}{2},j+\frac{1}{2}\right)$ because they can readly translated to global indices and vice versa

as illustrated in Figure 9 and Figure 7. 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 O, as you can see in Figure 10, is defined by its neighbors nodes on the West, East, North, South, North-West, North-East, South-west and South-East are defined as follows: W, E, N, S, NE, NW, SW and SE.

The discretized Equations (4.19) and (4.20) have been found to take the following general form

$C_P{\phi }_P+{\displaystyle \sum }_n{C}_{nb}{\phi }_{nb}=\langle f_P\rangle$ (15)

where $\sum$ indicates summation overall neighboring nodes ($nb$ ), $C_{nb}$ are the neighboring coefficients $C_E,C_W,C_S,C_N,C_{NE},C_{NW},C_{SE},C_{SW}$ and ${\phi }_{nb}$ are the

Figure 10. Nine point stencils.

approximate values of pressure at the neighboring nodes. 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. The snippet Matlab codes below shows how we proceed to convert Equations (4.19) and (4.20) in the last form i.e. that is conversion from 2D to 1D

5. Some Numerical Simulations

This section is firstly dedicated to numerically validating the analysis of the behavior of the solution proposed in the abstract of this paper, and secondly to showing the numerical test. We consider two test problems on the square domain $\Omega =\left]0,1\right[\times \left]0,1\right[$ . For each test problem a uniform rectangular mesh is used with different levels of refinement materialized by successive decreasing values assigned to the mesh size $h=\frac{1}{2^i}$ ($i=2,3,4,5$ ), we select a diffusion tensor $D$ and the exact solution $\phi \left(x,y\right)$ . Then after we calculate the corresponding right hand side. The error and convergence orders are performed in the following norms.

5.1. Error Norms Used and Convergence Orders

When the analytical solution is available and the mesh is refined, we set $H_0^1\left(\Omega \right)$ -norm or energy norm

$\begin{array}{c}Er{H}_0^1=\{{\displaystyle \sum _{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in {\mathfrak{D}}^{int}}\left\{{\left({\phi }_P-{\phi }_L\right)}^2+{\left({\phi }_{A^{\star }}-{\phi }_{B^{\star }}\right)}^2\right\}}\\ +{{\displaystyle \sum _{{\mathfrak{D}}_{\sigma {\sigma }^{\star }}\in {\mathfrak{D}}^{ext}}\left\{{\left({\phi }_P\right)}^2+{\left({\phi }_{A^{\star }}-{\phi }_{B^{\star }}\right)}^2\right\}}\}}^{\frac{1}{2}}\end{array}$ (5.1)

$L^2$ -norm or energy norm

(5.2)

and

(5.3)

which is the relative discrete $L^2$ -norm of the error for the exact potential. Defined by analogy, $ergradL2$ is the relative discrete $L^2$ -norm of the error on the exact potential gradient. We denote by $erL2\left(i\right)$ (resp. $ergradL2\left(i\right)$ ) the relative discrete $L^2$ -norm (resp. $ergradL2$ ) of the error on the exact potential (resp. exact potential gradient) corresponding to a level $i$ (integer $\ge 2$ ) of refinement for a given mesh. Let us set for all integers $i\ge 2$

$ratioL2\left(i\right)=-2\frac{\ln \left[erL2\left(i\right)\right]-\ln \left[erL2\left(i-1\right)\right]}{\ln \left[nunkw\left(i\right)\right]-\ln \left[nunkw\left(i-1\right)\right]}$ (5.4)

We define $ratiogradL2\left(i\right)$ , for all integers $i\ge 2$ , with the same relation as for $ratioL2\left(i\right)$ , except that $erL2$ is replaced by $ergradL2$ .

We denote by $ocvL2$ (resp. $ocvgradL2$ ) the order of convergence of the approximate potential (gradient of potential) in $L^2$ -norm. $ocvL2$ is given by the formula

$ocvL=\frac{\ln \left[erL2\left(imax\right)\right]-\ln \left[erL2\left(imax-1\right)\right]}{\ln \left[h\left(imax\right)\right]-\ln \left[h\left(imax-1\right)\right]}$ (5.5)

where $imax$ is the maximum level of refinement of a given mesh and $h\left(imax\right)$ the corresponding mesh size. $ocvgradL2$ is defined by the same formula as for $ocvL2$ except that $erL2$ is replaced by $ergradL2$ .

5.2. Numerical Results

Test 1: Heterogeneous rotating anisotropy. We consider a case: Where $\alpha \to +\infty$ , $c=0$ , $\phi \left(x,y\right)=\sin \left(\pi x\right)\sin \left(\pi y\right)$ and $f=-div\left(A\nabla \phi \right)$ , the diffusion operator $A$ is given by

$A=\frac{1}{x^2+y^2}\left[\begin{array}{cc}{10}^{-3}{x}^2+y^2& \left({10}^{-3}-1\right)xy\\ \left({10}^{-3}-1\right)xy& x^2+{10}^{-3}{y}^2\end{array}\right]$

Table 1. Test 1: Error and convergence orders.

level of refinement	Er-inf error	L2 error	H1 error
1	5.2e−02	5.200e−02	1.9e−01
2	1.4e−02	1.200e−02	6.2e−02
3	3.7e−03	3.1e−02	2.08e−02
4	9.5e−04	7.9e−04	7.1e−02
5	2.52e−04	2.1e−04	2.6e−03
ocv1er	1.92	1.98	1.56
ocv2er	1.88	1.87	1.51

Test 2. We consider a case: where $\alpha \to 0$ , $c=0$ $\phi \left(x,y\right)=\sin \left(2\pi y\right){\text{e}}^{-2\times {10}^{\frac{-5}{2}}\pi x}$ and $f=-div\left(A\nabla \phi \right)$ operator $A$ is given by

$A=\left[\begin{array}{cc}1& 0\\ 0& {10}^5\end{array}\right]$

and

$g=\{\begin{array}{ll}2\pi \times {10}^5\cos \left(2\pi y\right)\exp \left(-2\pi x\times {10}^{\frac{-5}{2}}\right)\hfill & \text{if}0\le x\le 1\text{and}y=1\hfill \\ 2\pi \times {10}^5\sin \left(2\pi y\right)\exp \left(-2\pi x\times {10}^{\frac{-5}{2}}\right)\hfill & \text{if}0\le y\le 1\text{and}x=0\hfill \\ -2\pi \times {10}^5\sin \left(2\pi y\right)\exp \left(-2\pi x\times {10}^{\frac{-5}{2}}\right)\hfill & \text{if}0\le y\le 1\text{and}x=1\hfill \\ -2\pi \times {10}^5\cos \left(2\pi y\right)\exp \left(-2\pi x\times {10}^{\frac{-5}{2}}\right)\hfill & \text{if}0\le x\le 1\text{and}y=0\hfill \end{array}$

It is important to note that in this special case the exact solution must satisfies the following null average condition

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

and the following compatibility condition to ensure the uniqueness of the solution

${\displaystyle {\int }_{\Omega }f\left(x,y\right)\text{d}x\text{d}y}-{\displaystyle {\int }_{\Gamma }g\left(x\right)\text{d}\gamma \left(x\right)}=0$

Table 2. Test problem 2.

level of refinement	Er-inf error	L2 error	H1 error
1	2.72e−02	2.49e−02	2.67e−01
2	1.44e−02	1.21e−02	9.77e−02
3	7.5e−03	6.00e−03	4.41e−02
4	3.82e−03	2.99e−03	1.53e−02
5	1.9e−03	1.49e−03	4.36e−03
6	9.7e−04	7.45e−04	1.18e−03
ocv1er	1.92	1.98	1.56
ocv2er	1.88	1.87	1.51

It follows by observing Table 1 and Table 2 of error estimates and convergences that, numerical tests above show a convergence of order two in L²-norm and Er-inf norm, a convergence of order 1.5 in relative norm H¹.

5.3. Conclusion and Perspectives

In summary, we proposed a numerical scheme for anisotropic diffusion problems with prescribed Robin boundary conditions. We showed an existence and uniqueness result of the solution of the discrete problem. Numerical tests confirmed the theoretical cases. Unfortunately, the numerical results for general cases (case where $\alpha \ne 0$ ) were not performed due to the robustness of calculations are designed for further work. Our challenge is to present similar tests on an unstructured mesh.

Funding

This work was supported by ongoing institutional funding. No additional grants to carry out or direct this particular research were obtained.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Njifenjou, A. and Moukouop-Nguena, I. (2001) Traitement des anisotropies de per-meabilite en simulation d’ecoulement en milieu poreux par les volumes finis. Sys-temes 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 Omnes, 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.[CrossRef]
[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.[CrossRef]
[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.[CrossRef]
[7]	Andreianov, B., Boyer, F. and Hubert, F. (2007) Discrete Duality Finite Volume Schemes for Leray-Lions-Type Elliptic Problems on General 2D Meshes. Numerical Methods for Partial Differential Equations, 23, 145-195.[CrossRef]
[8]	Boyer, F. and Hubert, F. (2008) Finite Volume Method for 2D Linear and Nonlinear Elliptic Problems with Discontinuities. SIAM Journal on Numerical Analysis, 46, 3032-3070.[CrossRef]
[9]	Donfack, H. and Jeutsa, A.K. (2022) Numerical Solution of 2-D Diffusion Problems Using Discrete Duality Finite Volume Method on General Boundary Conditions. Journal of Applied Mathematics and Physics, 10, 1968-1997.[CrossRef]
[10]	Gander, M.J., Halpern, L., Hubert, F. and Krell, S. (2021) Discrete Optimization of Robin Transmission Conditions for Anisotropic Diffusion with Discrete Duality Finite Volume Methods. Vietnam Journal of Mathematics, 49, 1349-1378.[CrossRef]
[11]	Le Dret, H. and Lucquin, B. (2016) Partial Differential Equations: Modeling, Analysis and Numerical Approximation. International Series of Numerical Mathematics, 168, Article 40340.