‘Discrete Maximum Principle Honored by Conventional Finite Volume Schemes for Diffusion-Convection-Reaction Problems: Proof with Geometrical Arguments

I. INTRODUCTION

Let $\mathrm{\Omega }$ be a bounded connected open subset of ${\mathbb{R}}^2$ whose boundary denoted by $\mathrm{\Gamma }$ is the union of polygonal lines ${\mathrm{\Gamma }}_k{]}_{k\in K}$ where $K$ is a finite subset of $\mathbb{N}$ which denotes the set of positive integers (see Figure 1 below). Note that if $K$ is a singleton then $\mathrm{\Omega }$ is a polygon (and so simply connected). Given the scalar functions $D(\cdot )$, $\mu (\cdot )$ and $f(\cdot )$ together with a vector field $(\cdot )$, all being defined in $\mathrm{\Omega }$, we consider the elliptic problem that consists in finding a scalar function $u(.)$ in an adequate function space such that

with the following homogeneous Dirichlet boundary conditions:

Under reasonable assumptions on the previous data i.e.

with

it is easy to prove that (see [1] for instance): The second order elliptic problem (1.1)-(1.2) gets a unique weak solution in the sense that

where we have set:

and

Following [2] one can prove that if the given function $f$ is positive almost everywhere in $\mathrm{\Omega }$ then the weak solution of the system (1.1)-(1.2) is also positive almost everywhere in $\mathrm{\Omega }$. That is the weak form of the Maximum Principle. Several works on construction of positivity-preserving numerical methods for diffusion, diffusion-convection, diffusion-reaction and diffusion-convection-reaction problems are available in the literature (see for instance [5, 6, 11, 13]). Such numerical methods are sometimes called monotone schemes.

The main objective of this work is to expose geometrical arguments for proving the well-known discrete version of the Maximum Principle satisfied by the conventional finite volume solution to the system $(1.1)-(1.2)$.

II. PRELIMINARY TOOLS

Definition 2.1 (Partition of $\mathrm{\Omega }$ ). Let $\bar{\mathrm{\Omega }}$ be the closure of $\mathrm{\Omega }$ in the sense of the standard topology of ${\mathbb{R}}^2$ and let $J$ be a finite subset of $\mathbb{N}$ which is the set of positive integers. A family $\{{\mathrm{\Omega }}_j{\}}_{j\in J}$ made up of subsets of $\bar{\mathrm{\Omega }}$ defines a partition of $\bar{\mathrm{\Omega }}$ if the following conditions are satisfied:

where $Int(\diamond )$ denotes the interior of $\diamond$ in the sense of standard topology of ${\mathbb{R}}^2$.

Let us consider a partition P over $\mathrm{\Omega }$ consisting in a finite family of closed convex polygons (named also polygonal elements) generically denoted by T. These polygonal elements are the so-called control volumes in the language of Finite Volume theory. The control volumes from the partition P defines a conforming Finite Volume mesh over $\bar{\mathrm{\Omega }}$ if (in addition to conditions (i)-(iii) from Definition 2.1) the following conditions are satisfied:

where $\mathrm{\varnothing }$ denotes the empty set. Let us denote by $\partial P$ the set of boundary edges (viewed as degenerate control volumes) and we briefly define the conventional finite volume mesh T as it follows: $T=\{P,\partial P\}$.

We should use intensively in what follows a notion of characteristic function slightly different from the usual one and defined as follows.

Definition 2.2 Let T be a control volume either from P or from $\partial P$. We call in this work the characteristic function of T denoted by $1_T$ the function defined almost everywhere either in $\mathrm{\Omega }$ (with respect to Lebesgue measure in 2-D) or on $\mathrm{\Gamma }=\bigcup _{k\in K}{\mathrm{\Gamma }}_k$ (with respect to Lebesgue measure in 1-D) by:

where $Ext(\diamond )$ denotes the exterior of a subset $\diamond$ from $R^2$ (with respect to the natural topology of $R^2$ ). Recall that $Int(\diamond )$ stands for the interior of $\diamond$ from $R^2$. ☐

Let us introduce the following discrete function spaces that play a key-role in the sequel.

Definition 2.3 We set:

where $0_{S^{\partial P}}$ is the zero-function (denoted simply 0 if there is no risk of confusion) from the discrete function space $S^{\partial P}$. ☐

The Finite Volume method is based on the fundamental idea that the exact solution u could be approximated inside any control-volume T with a constant $U_T$ corresponding to either the mean-value of u or its approximation at a given point located inside T, with Cartesian coordinates $x_T$. In the context of conventional Finite Volumes the choice of that point is not arbitrary as we will be seeing in assumption $({\mathcal{A}}_3)$ below. Let us denote by ${\mathrm{\Gamma }}_T$ the boundary of any control-volume T. We need to specify the following assumptions that make the conventional Finite Volumes very attractive and realistic for certain engineering problems as subsurface flow problems (notice that [3, 4, 15] are among distinguished references on fluid flow in porous media):

$({\mathcal{A}}_1)$ The diffusion coefficient $D(\cdot )$ is a piecewise constant function i.e.

where $\{{\mathrm{\Omega }}_s{\}}_{s\in S}$ defines a partition $\mathcal{P}$ of the domain $\bar{\mathrm{\Omega }}$ in the sense of Definition 2.1.

Denote by T the Finite Volume mesh corresponding to the partition P. Let us make the following assumption on T.

$({\mathcal{A}}_2)$ $\mathcal{T}$ is compatible with the discontinuities of $D(\cdot )$ in the sense that the discontinuity points of $D(\cdot )$ belong to the mesh interfaces ${\mathrm{\Gamma }}^T=\bigcup _{T\in \bar{P}}{\mathrm{\Gamma }}_T$, where we have set $\bar{P}=P\cup \partial P$. In other words any discontinuity point of the function $D(\cdot )$ is located in a control volume boundary.

(A3) For all $(T^{\prime },T^{″})\in \bar{\mathcal{P}}\times \mathcal{P}$ such that $T^{\prime }$ and $T^{″}$ are adjacent (that is ${\mathrm{\Gamma }}_{T^{\prime }}\cap {\mathrm{\Gamma }}_{T^{″}}$ is a common edge for control volumes $T^{\prime }$ and $T^{″}$ ), the vector $x_{T^{\prime }}-x_{T^{″}}$ is orthogonal to the common edge. This is the so-called orthogonality condition required for conventional Finite Volume meshes (see [5, 6]).

An immediate consequence of the assumption (1.3) is that $D(\cdot )$ is a nonnegative constant function in each control volume T. We denote by $D^T$ the constant value of $D(\cdot )$ in the control volume T.

Let us give a brief description of the different steps for getting a conventional finite volume scheme. We start with introducing some useful notations: E is the set of all mesh edges, $E^{int}$ is the subset of E made of interior mesh edges and $E^{ext}$ is the subset of E made of exterior mesh edges i.e. mesh edges lying on the domain boundary.

Step 1: Integrate the two sides of the balance equation (1.1) in each control volume $T$ from the family $\mathcal{P}$. So we get what follows (thanks to Ostrogradski's theorem):

where ${\nu }_T$ stands for outward unit vector normal to the control-volume boundary ${\mathrm{\Gamma }}_T$.

Step 2: Re-write the first two integral terms from the left-hand side of (3.2) as follows for all $T\in P$:

where $E_T$ is the set of mesh edges $\sigma$ lying in ${\mathrm{\Gamma }}_T$ and where ${\nu }_{\sigma ,T}$ stands for outward unit vector normal to the portion $\sigma$ of the control-volume boundary ${\mathrm{\Gamma }}_T$, called again mesh edge associated with ${\mathrm{\Gamma }}_T$. Integrals from the first summation are diffusion fluxes while integrals from the second summation are convection fluxes (called sometimes advection fluxes).

Step 3: Perform the approximation of the unknown function u in the control-volume T with the unknown real constant $u(x_T)$. So one could set what follows concerning approximation of the reaction term:

where $I_T(\diamond )$ is the integral of a function $\diamond$ defined in the control volume T.

Step 4: Look for reasonable approximations of flux integral terms from the left-hand side of (3.3). What should one understand by reasonable approximations? We mean that the flux approximations should take account of the following constraints:

• Perform the upwind approximation of the convection flux in view to ensure the stability of the global finite volume scheme. For that purpose, let us start with setting:

Definition 3.1

Definition 3.2 (Upwind approximation of the convection flux over $\sigma \in E^{int}$ ) Let $\sigma$ in $E_T\cap E_L$, with T and L from the set P. We set:

In other words the upwind approximation of the convective flux across the interior edge $\sigma$ in $E_T\cap E_L$ could be defined as follows:

Since (according to the flux continuity principle over grid-block interfaces)

the preceding approximation of the convective flux is equivalent to the following one

• The flux continuity across interior edges $\sigma$, i.e. $\sigma \in {\mathcal{E}}^{int}$, is a fundamental physical principle to be met. So we have necessarily for all $\sigma \in {\mathcal{E}}^{int}$

Integrating the right-hand and the left-hand sides of (3.9) over $\sigma \in E^{int}$ leads to the following "weak formulation" of flux continuity:

Since the weak solution $u$ of the system (1.1)-(1.2) lies in $H_0^1(\mathrm{\Omega })$, the trace $u_{|\sigma }$ exists (in $H^{\frac{1}{2}}(\sigma )$ for instance) in a unique manner. In consequence we naturally get what follows:

Thus, the previous "weak formulation" of flux continuity (3.10) is reduced to

In the context of conventional Finite Volumes the family P satisfies the so-called orthogonality condition (see assumption $({\mathcal{A}}_3)$ above at the beginning of the current Section). So there exists a family of points $\{x_T;T\in P\}$, such that for any pair $(T,L)\in \mathcal{P}\times \bar{\mathcal{P}}$, with T and L adjacent, the orthogonal projections of $x_T$ and $x_L$ on their common edge $\sigma$ coincides and let call it $x_{\sigma }$. We make the following convention:

"If $T$ is adjacent to the domain boundary we set: $L\stackrel{def}{=}\sigma$, where $\sigma$ is the boundary edge associated with $T$, and $x_L$ coincides with $x_{\sigma }$ ".

This being said, from the following diffusion flux approximation (assuming the exact solution restriction $u_{{|}_T}$ in $C^0(\bar{T})$ for any $T\in \mathcal{P}$; it is the case if $u_{{|}_T}\in H^2(T)$ ):

where $mes(.)$ stands for Lebesgue measure in one-space dimension, $dist(..)$ represents the Euclidean distance and where $x_{\sigma ,T}$ is in fact the point $x_{\sigma }$ seen as from the boundary of T by an observer standing inside T. The principle of continuity of u on grid-block interfaces is expressed at the discrete level by the relation:

where $(\mathcal{P}\times \mathcal{P}{)}_{adj}$ is the subset of $P\times P$ made of $(T,L)$ such that T and L are adjacent. So it is reasonable to set:

With the above notation the diffusion flux approximation could read as follows