The Finite Volume Formulation for 2D Second-Order Elliptic Problems with Discontinuous Diffusion/Dispersion Coefficients

Abstract

We propose a finite volume method for the numerical resolution of two-dimensional steady diffusion problems with possibly discontinuous coefficients on unstructured polygonal meshes. Our numerical method is cellcentered, secondorder accurate on smooth solutions and based on a special numerical treatment of the diffusion/dispersion coefficients that makes its application possible also when such coefficients are discontinuous. Numerical experiments confirm the convergence of the numerical approximation and show a good behavior on a set of benchmark problems in two space dimensions.

1. Introduction

In the last decades, finite volume methods have been greatly successful in solving engineering models of flows in porous media on complex geometrical domain because the finite volume formulation works on general polygonal and polyhedral meshes. This great mesh flexibility is now combined with a strong theoretical foundation, that is, convergence analysis and error estimates are available on the simplest mathematical models [1]. In this work, we address a version of the finite volume method that is popularly known as “the diamond scheme” and was originally presented for the advection-diffusion equation in two dimensions in [2–4] and successively extended to convection-dominated problems in [5, 6] and nonlinear flow problems in partially saturated porous media [7].

The numerical treatment of both diffusion/dispersion flux and of the partially saturated flow is based on averaging one-sided gradients and diffusion (conductivity) tensors. Many papers addressed the issue of conductivity averaging both in finite differences and in finite elements. A seminal work [8] started with a finite difference application to the unsaturated flow equation. Later on they, dealt with the problem in a finite volume scheme [9]. Recent reviews on this topic can be found in [10–12]. They also made performance comparisons using different averages. The latter paper shows in conclusion that their scheme, namely, an averaging schemes based on Darcian mean principle used in the framework of either vertex-centered or cell-centered approach compare favorably to other methods for a range of test cases. They compared their method with the one of [11]. Aside of the numerical comparisons, in [13] it is remarked that the usage of the harmonic average can be motivated on the basis of physical arguments.

Another reason for dealing with the tensor averaging issue is related to the scale problems. Typically, the Darcy equation is written at a scale much smaller than the scale of practical interest [14, 15]. This gives boost to new numerical techniques to able to deal with the heterogeneity of hydrodynamic parameters. One possibility is to obtain either analytically or empirically explicit equations for the scale of interest, eliminating other scales in the problem [16].

Several methods have been recently developed, such as, the heterogeneous multiscale method (HMM) [16] and the multiscale finite volume element method [17]. Babuška in the 70s motivated the multiscale finite element method [18–20]. Multiscale methods have been proposed also in the case of saturated flows in heterogeneous porous media and also applied the multiscale finite volume element method to the Richards equation [21, 22]. In [14], a multiscale method based on the framework of HMM was recently extended to solve the Richards equation with random coefficients.

Another feature of their numerical scheme is the formula estimating the macroscopic flux in which the unsaturated hydraulic conductivity can be calculated as a diagonal tensor. In particular, it is worth noting that a finite volume high-order accurate approximation of the pressure head field also allows one to achieve a better resolution in the approximation of other important fields like the components of the hydraulic conductivity tensor at the mesh vertices [7]. For the sake of simplicity, in this work it was considered the steady flow in a layered porous medium in the presence of a source term. It is the Poisson equation which is suitable to test new schemes which will be applied in further work to the more general unsaturated transient cases.

We propose a novel technique that automatically adjusts when the diffusion tensor is discontinuous across a mesh interface shared by two adjacent cells. This technique is general and can be easily implemented in any finite volume scheme that has an explicit numerical flux and may result in a particularly efficient strategy for the numerical resolution of problems where both the diffusive/dispersive and convective phenomena are simultaneously significant. In fact, in such problems high-resolution finite volume methods are normally coupled with mixed finite elements following the criterion of choosing the best available technique in accordance with the nature of the equation to be discretized. In fact, the ℝ𝕋₀−ℙ₀ method is more suitable to the numerical discretization of the diffusive part of the model, and the finite volume method gives an accurate and stable discretization of the convective part of the model, even in presence of strong convective fields. This approach was proved successful in the numerical modeling of oil reservoir problems [23] and of groundwater flow and transport of contaminants in porous media, combared with [24–26]. Moreover, other different engineering areas may benefit from this new technique, such as [27–38].

The outline of the paper is as follows. In Section 2, we introduce the mathematical model and discuss the formulation of the finite volume methods based on vertex reconstruction. In Section 3, we present the numerical treatment of the diffusion tensor. In Section 4, we confirm the theoretical results with numerical experiments. In Section 5 we offer final remarks and conclusions.

2. The Mathematical Model and the Finite Volume Formulation

Let Ω be a polygonal domain with boundary Γ.

We consider the steady diffusion problem for the scalar solution field 𝑢 given by−div(Λ∇𝑢)=𝑓inΩ,𝑢=𝑔onΓ, (2.1)(2.2) where Λ is the diffusion tensor describing the material properties, 𝑓 is the forcing term and 𝑔 the boundary function that defines nonhomogeneous Dirichlet conditions on the boundary Γ. Under suitable regularity assumptions on 𝑓, 𝑔, Λ, and Ω, it turns out that the diffusion problem is mathematically wellposed, a unique solution exists [39] and such a solution is continuously dependent on the model data. In particular, due to the elliptic nature of the model, the diffusion tensor Λ(𝐱) is normally represented by a strictly positive definite matrix for every 𝐱 in Ω. The components of the diffusion tensor may be discontinuous in the computational domain; in such a case and without loss of generality, we assume that the mesh is conforming with each discontinuity of Λ.

The numerical approximation to (2.1)–(2.2) is performed on a sequence of polygonal partitions {Ω_ℎ}_ℎ of the domain Ω. For any mesh Ω_ℎ, the subscripted label ℎ is the mesh size and is defined byℎ=max_𝖪∈Ωℎℎ_𝖪,ℎ_𝖪=sup_{𝐱′,𝐱′′∈𝖪}𝐱′−𝐱^′′, (2.3)

where ℎ_𝖪 is the diameter of the polygonal cell 𝖪∈Ω_ℎ. We also label the numerical solution 𝑢_ℎ calculated on the mesh Ω_ℎ by ℎ. Also,(i)a generic mesh vertex is denoted by 𝖵 and its coordinate vector by 𝐱_𝖵. We will also find it convenient to introduce a local numbering of the vertices, for example, 𝖵_𝑖, and to ease notation to denote the vertex position of the 𝑖-th vertex by 𝐱_𝑖 (instead of 𝐱_𝖵𝑖);(ii) a generic cell interface or a boundary edge is denoted by 𝜎, its center (i.e., its edge midpoint) by 𝐱_𝜎, and its measure (edge length) by |𝜎|;(iii) a generic polygonal cell is denoted by 𝖪 its measure (area) by |𝖪|, its center of gravity by 𝐱_𝖪, and its boundary by 𝜕𝖪.

The orientation of each mesh interface 𝜎 is reflected by its unit normal vector 𝐧_𝜎, which is fixed once and for all. For any mesh cell 𝖪 and any face 𝜎 of the polygonal boundary 𝜕𝖪, we define the unit normal vector 𝐧_𝖪𝜎 that points out of 𝖪 and we also use the notation 𝐍_𝖪𝜎=|𝜎|𝐧_𝖪𝜎.

Let the algebraic vector 𝑢_ℎ=(𝑢_𝖪)_𝖪∈Ωℎ be the numerical solution, where each 𝑢_𝖪 approximates the cell average of the scalar solution 𝑢 over the cell 𝖪. The finite volume scheme for 𝑢_ℎ on the computational mesh Ω_ℎ reads−_{𝜎∈𝜕𝖪}Λ_𝜎∇_𝜎𝑢_ℎ⋅𝐍_𝖪𝜎=𝑓_𝖪|𝖪|∀𝖪∈Ω_ℎ, (2.4) where Λ_𝜎 is an evaluation of Λ at face 𝜎 taken from the side of 𝖪; ∇_𝜎 is the discrete face gradient built inside the diamond cell centered at face 𝜎; 𝐍_𝖪,𝜎 is the geometric vector perpendicular to 𝜎, pointing from 𝖪 to 𝖫, and with lenght |𝐍_𝖪𝜎|=|𝜎|; 𝑓_𝖪 is the cell average of the right-hand side term 𝑓 over the cell 𝖪:𝑓_𝖪1∶=|𝖪|_𝖪𝑓(𝐱)𝑑𝑉. (2.5)

Let Λ_𝖪 and Λ_𝖫 be first-order approximations of the diffusion tensor Λ within the control volumes 𝖪 and 𝖫. For example, we can take Λ_𝖪∶=Λ(𝐱_𝖪) and Λ_𝖫∶=Λ(𝐱_𝖫), where 𝐱_𝖪 and 𝐱_𝖫 are internal points (not necessarily the barycenters) of 𝖪 and 𝖫, respectively. Now, we can consider either Λ_𝜎=Λ_𝖪 inside 𝖪 and Λ_𝜎=Λ_𝖫 inside 𝖫 or a suitable mean of these two tensors. In particular, we can deal with the arithmetic mean:Λ_𝜎=12Λ_𝖪+Λ_𝖫, (2.6)

or with the harmonic mean:Λ_𝜎⁻¹=12Λ_𝖪⁻¹+Λ_𝖫⁻¹. (2.7)

Both approaches deserve a special care when Λ is discontinuous across 𝜎 and a proper definition of the discrete gradient ∇_𝜎𝑢_ℎ is also needed in order to preserve the property of flux conservation.

To define the numerical gradient, a special control volume is built around this interface, which has a quadrilateral shape in two dimensions and is named “diamond cell”. The geometry of a diamond cell is shown in Figure 1, which plots the diamond cell 𝖣 centered at the mesh face 𝜎 with vertices 𝐱_𝑖 and 𝐱_𝑖+1 and shared by the cells with centers 𝐱_𝖪 and 𝐱_𝖫. The diamond cell 𝖣 can also be seen as the union of the subdiamonds 𝐷_𝖪 and 𝐷_𝖫, which are the triangular cells sharing 𝜎 as common base and having, respectively, vertices 𝐱_𝖪 and 𝐱_𝖫, the centers of gravity of cells 𝖪 and 𝖫. The four vectors 𝐍_𝖪𝑖, 𝐍_𝖫𝑖, 𝐍_𝖫𝑖+1, and 𝐍_𝖪𝑖+1 shown in Figure 1 are respectively orthogonal to the four boundary faces 𝜎_𝖪𝑖, 𝜎_𝖫𝑖, 𝜎_𝖫𝑖+1, and 𝜎_𝖪𝑖+1, and their lenght is equal to the lenght of the corresponding face. Instead, when 𝜎 is a boundary face, thus defined by 𝜎=𝜕𝖪∩Γ (Γ being the boundary of the computational domain Ω), the diamond cell associated to 𝜎 coincides with 𝐷_𝖪, that is, it is the triangle defined by 𝜎 and the vertex 𝐱_𝖪.

Figure 1 

Geometry of the 2D diamond cell.

The numerical diffusive flux is built by using a constant approximation of the solution gradient and an evaluation of the diffusion tensor Λ within the diamond cell, as discussed before. Let us give the formulas for an internal mesh face 𝜎, that we suppose to be shared by two cells 𝖪 and 𝖫. All these formulas can be easily adapted to the case of a boundary face by taking 𝐱_𝖫 equal to the center of gravity of face 𝜎.

Using the Green-Gauss theorem yields:∇_𝜎𝑢_ℎ≈1|𝖣|_𝖣(𝐱)1∇𝑢𝑑𝑉=|𝖣|_𝜕𝖣𝑢(𝐱)𝐧(𝐱)1𝑑𝑆=|𝖣|_{𝜎∈𝜕𝖣𝜎}𝑢(𝐱)𝐧(𝐱)𝑑𝑆, (2.8) where 𝐧 is the unit vector orthogonal to 𝜎∈𝜕𝖣 and pointing out of 𝖣. If the restriction of 𝑢 to the face 𝜎 of 𝜕𝖣 is an affine function, the boundary integral on 𝜕𝖣 in (2.8) only depends on the values of 𝑢 at the vertexes of 𝖣 and the constant vector provided by the formulas that we derive below must coincide with the gradient of 𝑢.

The Gauss-Green theorem applied to the solution gradient ∇𝑢 integrated on the 2D diamond cell 𝖣 of Figure 1 yields:_{𝜎∈𝜕𝖣𝜎}𝑢(𝐱)𝐧(𝐱)1𝑑𝑆=2𝑢𝐱_𝖪+𝑢𝐱_𝑖𝐍_𝖪𝑖+𝑢𝐱_𝑖+𝑢𝐱_𝖫𝐍_𝖫𝑖+𝑢𝐱_𝖫+𝑢𝐱_𝑖+1𝐍_𝖫𝑖+1+𝑢𝐱_𝑖+1+𝑢𝐱_𝖪𝐍_𝖪𝑖+1. (2.9) Let us introduce the vectors 𝐍_𝖪𝖫 and 𝐍_𝑖𝑖+1, which are, respectively, orthogonal to the edge connecting 𝐱_𝖪 to 𝐱_𝖫 and 𝐱_𝑖 to 𝐱_𝑖+1, and whose lenght are equal to the length of these edges. Using the geometric relations 𝐍_𝖪𝖫=−(𝐍_𝖪𝑖+𝐍_𝖫𝑖) and 𝐍_𝖪𝖫=−(𝐍_𝖪𝑖+𝐍_𝖫𝑖) into (2.9) and rearranging the terms yields: _𝜕D(𝐱)(𝐱)1∇𝑢⋅𝐧𝑑𝑆=2𝑢𝐱_𝑖+1−𝑢𝐱_𝑖𝐍_𝖪𝖫+𝑢𝐱_𝖫−𝑢𝐱_𝖪𝐍_𝑖𝑖+1. (2.10)

To derive the gradient formula, we replace the function values 𝑢(𝐱_𝖪) and 𝑢(𝐱_𝖫) with the cell unknowns 𝑢_𝖪 and 𝑢_𝖫, respectively, and the function values 𝑢(𝐱_𝑖) and 𝑢(𝐱_𝑖+1) with the corresponding nodal unknowns 𝑢_𝑖 and 𝑢_𝑖+1. We obtain𝑢_ℎ1∶=2|𝖣|∇_𝖪𝖫𝑢_𝖫−𝑢_𝖪𝐍_𝖪𝖫+𝑢_𝑖+1−𝑢_𝑖𝐍_𝑖𝑖+1. (2.11) which is a constant approximation of the solution gradient ∇𝑢 within the quadrilateral cell 𝖣.

With any vertex 𝐱_𝑖 of the mesh Ω_ℎ, we associate the reconstructed vertex value 𝑢_𝑖=∑_{𝖪∈Ω𝑖}𝜔_𝖪𝑖𝑢_𝖪, where Ω_𝑖={𝖪∶𝐱_𝑖∈𝖪} denotes the subset of the mesh cells which share the vertex 𝐱_𝑖. The interpolation weights 𝜔_𝖪𝑖 are assumed to verify the consistency relations [3]:_{𝖪∈Ω𝑖}𝜔_𝖪𝑖=1,_{𝖪∈Ω𝑖}𝜔_𝖪𝑖𝐱_𝑖−𝐱_𝖪=0. (2.12) Alternative constructions of the interpolation weights are available in the literature, see, for instance, [40]. The main advantage offered by these alternative choices concerns the control of positivity of the weights, which turns out to be significant when we are aimed at developing a numerical method with a maximum or minimum principle. However, as this topic is out of the scope of the current work, we will not pursue it anymore.

The interpolation weights 𝜔_𝖪𝑖 are obtained by solving the reconstruction problem that approximates the cell-averaged data set {(𝐱_𝖪,𝑢_𝖪)_{𝖪∈Ω𝑖}} by the affine function:𝑢_𝑖(𝐱)=𝛼+𝜷⋅𝐱_𝖪−𝐱_𝑖for𝐱∈𝒱_𝑖 (2.13) on the covolume 𝒱_𝑖=⋃_{𝖪∈Ω𝑖}𝖪 and in a least square sense, compared to [3, 41]. The reconstructed value at vertex 𝐱_𝑖 is now given by taking 𝑢_𝑖=𝑢_𝑖(𝐱_𝑖)=𝛼. The coefficients (𝛼,𝜷) are the minimizers of the least squares functional:𝒥_𝑖(𝛼,𝜷)=_{𝖪∈Ω𝑖}𝛼+𝜷⋅𝐱_𝖪−𝐱_𝑖−𝑢_𝖪². (2.14) Imposing the zero gradient condition, that is, ∇_𝛼,𝜷𝒥(𝛼,𝜷)=0, yields a linear system for the coefficients (𝛼,𝜷), whose solution returns the interpolation weights. For completeness of exposition, we briefly review the derivation of the weight formula given in [41]. Let 𝜷=(𝛽_𝑥,𝛽_𝑦); a straightforward calculation yields:𝜕𝒥_𝑖𝜕𝛼=2_{𝖪∈Ω𝑖}𝛼+𝜷⋅𝐱_𝖪−𝐱_𝑖−𝑢_𝖪=0,𝜕𝒥_𝑖𝜕𝛽_𝑥=2_{𝖪∈Ω𝑖}𝛼+𝜷⋅𝐱_𝖪−𝐱_𝑖−𝑢_𝖪𝑥_𝖪−𝑥_𝑖=0,𝜕𝒥_𝑖𝜕𝛽_𝑦=2_{𝖪∈Ω𝑖}𝛼+𝜷⋅𝐱_𝖪−𝐱_𝑖−𝑢_𝖪𝑦_𝖪−𝑦_𝑖=0. (2.15)

Let 𝑚_𝑖 be the number of cells of Ω_𝑖; we reformulate the linear system written above as_{𝖪∈Ω𝑖}𝑚_𝑖𝛼+_{𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖𝛽_𝑥+_{𝖪∈Ω𝑖}𝑦_𝖪−𝑦_𝑖𝛽_𝑦=_{𝖪∈Ω𝑖}𝑢_𝖪,_𝖪∑_∈Ω𝑖𝑥_𝖪−𝑥_𝑖𝛼+∑_{𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖²𝛽_𝑥+∑_{𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖𝑦_𝖪−𝑦_𝑖𝛽_𝑦=∑_{𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖𝑢_{𝖪,𝖪∈Ω𝑖}𝑦_𝖪−𝑦_𝑖𝛼+_{𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖𝑦_𝖪−𝑦_𝑖𝛽_𝑥+_{𝖪∈Ω𝑖}𝑦_𝖪−𝑦_𝑖²𝛽_𝑦=_{𝖪∈Ω𝑖}𝑦_𝖪−𝑦_𝑖𝑢_𝖪. (2.16)(2.17)(2.18) For convenience, we consider a suitable local numbering of the cells of 𝒱_𝑖, for example, 𝒱_𝑖=𝖪₁∪𝖪₂∪…𝖪_𝑚𝑖, and we introduce the two matrices:𝑚𝖠=^{𝑖𝖪∈Ω𝑖}𝑥_𝖪−𝑥_{𝑖𝖪∈Ω𝑖}𝑦_𝖪−𝑦_{𝑖𝖪∈Ω𝑖}𝑥_𝖪−𝑥_{𝑖𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖²_{𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖𝑦_𝖪−𝑦_{𝑖𝖪∈Ω𝑖}𝑦_𝖪−𝑦_{𝑖𝖪∈Ω𝑖}𝑥_𝖪−𝑥_𝑖𝑦_𝖪−𝑦_{𝑖𝖪∈Ω𝑖}𝑦_𝖪−𝑦_𝑖²,𝖡=𝑥11…1_𝖪1−𝑥_𝑖𝑥_𝖪2−𝑥_𝑖…𝑥_𝖪𝑚𝑖−𝑥_𝑖𝑦_𝖪1−𝑦_𝑖𝑦_𝖪2−𝑦_𝑖…𝑦_𝖪𝑚𝑖−𝑦_𝑖, (2.19)

and the column vector 𝐮=(𝑢₁,𝑢₂,…,𝑢_𝑚𝑖)^𝑇, which collects the solution averages to be imposed. Using such definitions, the linear system (2.16)–(2.18) becomes𝖠𝛼𝛽_𝑥𝛽_𝑦𝛼=𝖡𝐮,whichimpliesthat𝛽_𝑥𝛽_𝑦=𝖠⁻¹𝖡𝐮, (2.20)since matrix 𝖠 is, obviously, symmetric and positive definite and, thus, nonsingular. The coefficient 𝛼 that provides the vertex value is given by()𝛼𝛼=1,0,0𝛽_𝑥𝛽_𝑦=(1,0,0)𝖠⁻¹𝖡𝐮. (2.21) We require that such a coefficient be an average of the data in 𝐮 with coefficients 𝝎=(𝜔_K1𝑖,𝜔_K2𝑖,…,𝜔_{K𝑚𝑖𝑖})^𝑇, that is, 𝛼=𝝎^𝑇𝐮. By comparison with (2.21), we have the final weight formula:𝝎=𝖡^𝑇𝖠^−𝑇100. (2.22)

Conditions (2.12) can be checked by a direct calculations.

The functional in (2.14) can be suitably modified to take into account different kind of boundary conditions, for example, Neumann or Robin, compared to[41]. The values 𝑢_𝑖 at the vertexes 𝐱_𝑖∈Γ on the Dirichlet boundary are constrained to the boundary data, for instance 𝑢_𝑖=0 for a homogeneous condition. We also mention that this choice of weights provides a very robust technique, compared to the study of locking effect due to the accuracy of the reconstruction [42].

3. Treatment of Discontinuous Diffusion Tensors

We discuss, here, how to treat the case of a diffusion tensor that is discontinuous across an interface 𝜎 shared by the cells 𝖪 and 𝖫, that is, 𝜎=𝖪∣𝖫. Since we suppose that Λ is discontinuous across 𝜎 and the normal flux of the exact solution is continuous across 𝜎, the normal component of the solution gradient must be discontinuous. Consequently, a numerical approximation of the diffusive flux across 𝜎 likeΛ∇(_𝜎𝑢)⋅𝐧𝑑𝑆≈Λ_𝖪𝜎∇_𝜎𝑢_ℎ⋅𝐍_𝖪𝜎 (3.1)

cannot be consistent whenever Λ_𝖪𝜎 is some kind of average between Λ_𝖪 and Λ_𝖫 and the constant vector ∇_𝜎(𝑢_ℎ) approximates ∇(𝑢)  over all the diamond cell hull(𝐱_𝖪,𝐱_𝖫,𝜎). To tackle this problem, we consider two distinct approximations of the solution gradient within the subcells 𝐷_𝖪 and 𝐷_𝖫, denoted by ∇_𝖪𝜎(𝑢_ℎ) and ∇_𝖫𝜎(𝑢_ℎ), respectively, and impose directly the flux conservation. To derive an expression for such one-sided gradients, we introduce an additional unknown 𝑢_𝜎 along face 𝜎 and we apply the Gauss-Green Theorem as for the derivation of ∇_𝖪𝖫𝑢_ℎ in the previous section. Also, we recall that 𝐍_𝖪𝜎=𝐍_𝖪𝖫 and |𝐍_𝖪𝜎|=|𝜎| and we use a similar notation for the normal vector from the other side of 𝜎, so that 𝐍_𝖫𝜎=−𝐍_𝖪𝖫 and |𝐍_𝖫𝜎|=|𝜎|. The two one-sided gradient formulas are∇_𝖪𝜎𝑢_ℎ1∶=2𝐷_𝖪𝑢_𝜎−𝑢_𝖪𝐍_𝖪𝜎+𝑢_𝑖+1−𝑢_𝑖𝐍_𝑖𝑖+1,∇_𝖫𝜎𝑢_ℎ1∶=2𝐷_𝖫𝑢_𝜎−𝑢_L𝐍_𝖫𝜎−𝑢_𝑖+1−𝑢_𝑖𝐍_𝑖𝑖+1. (3.2)

Now, we search for a tensor Λ_𝖪𝜎, which is an average of the diffusion tensors Λ_k and Λ_𝖫, and a gradient vector ∇_𝖪𝖫, which is an average of the gradient vectors ∇_𝖪𝜎(𝑢_ℎ) and ∇_𝖫𝜎(𝑢_ℎ), that makes it possible to express directly the flux conservation asΛ_𝖪𝜎∇_𝖪𝖫𝑢_ℎ⋅𝐧_𝜎∶=Λ_𝖪∇_𝖪𝜎𝑢_ℎ⋅𝐧_𝜎=Λ_𝖫∇_𝖫𝜎𝑢_ℎ⋅𝐧_𝜎. (3.3) Since we expect that the normal component of the vectors ∇_𝖪𝜎(𝑢_ℎ) and ∇_𝖫𝜎(𝑢_ℎ) be the same, we impose directly this condition by requiring that a scalar coefficient 𝜑 exists such that∇_𝖫𝜎𝑢_ℎ=∇_𝖪𝜎𝑢_ℎ+𝜑𝐧_𝜎, (3.4) (recall that 𝐧_𝜎=𝐧_𝖪𝖫). We also assume that ∇_𝖪𝖫(𝑢_ℎ) be a weighted average of ∇_𝖪𝜎(𝑢_ℎ) and ∇_𝖫𝜎(𝑢_ℎ). Let us be given two nonnegative coefficients 𝜇_𝖪 and 𝜇_𝖫 such that𝜇_𝖪+𝜇_𝖫=1, (3.5) and such that∇_𝖪𝖫𝑢_ℎ=𝜇_𝖪∇_𝖪𝜎𝑢_ℎ+𝜇_𝖫∇_𝖫𝜎𝑢_ℎ. (3.6)

To this purpose, a number of possible choices have been proposed and widely investigated in the literature for 𝜇_𝖪 and 𝜇_𝖫. For example, a popular option is to take the volume of the two subdiamonds 𝐷_𝖪 and 𝐷_𝖫 normalized by the total volume of the diamond, for example, 𝜇_𝖪=|𝐷_𝖪|/|𝖣| and 𝜇_𝖫=|𝐷_𝖫|/|𝖣|. The diffusion tensor can be chosen as the simple average of Λ_𝖪 and Λ_𝖫, for example, Λ_𝖪𝖫=𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫, or as the harmonic average, for example, Λ⁻¹_𝖪𝖫=𝜇_𝖪Λ_𝖪⁻¹+𝜇_𝖫Λ_𝖫⁻¹. Much more attention must be paid when Λ is discontinuous across 𝜎 as none of the previous choices may provide the correct flux information across the common interface of cells 𝖪 and 𝖫. Our aim in the next developments is at determining a value for 𝜑 and of the average diffusion tensor Λ_𝖪𝖫 in terms of Λ_𝖪 and Λ_𝖫 that ensures flux conservation (3.3) for any given pair of coefficients 𝜇_𝖪 and 𝜇_𝖫. Our derivation is similar to that considered in [43], and as pointed out therein, the resulting matrix is the same obtained for other purposes in the field of upscaling of conductivity, either by means of the large average techniques [44] or using the homogenization theory in the case of layered materials [45].

To this purpose, we first derive an expression for 𝜑. We multiply (3.4) by Λ_𝖫𝐧_𝜎 and use flux conservation to obtain𝐧_𝜎⋅Λ_𝖫∇_𝖫𝜎𝑢_ℎ=𝐧_𝜎⋅Λ_𝖫∇_𝖪𝜎𝑢_ℎ+𝜑Λ_𝖫𝐧_𝜎=𝐧_𝜎⋅Λ_𝖪∇_𝖪𝜎𝑢_ℎ, (3.7)

from which we have𝐧_𝜎⋅Λ_𝖪−Λ_𝖫∇_𝖪𝜎𝑢_ℎ=𝜑𝐧_𝜎⋅Λ_𝖫𝐧_𝜎. (3.8) Likewise, we multiply (3.4) by Λ_𝖪𝐧_𝜎 and use flux conservation to obtain𝐧_𝜎⋅Λ_𝖪∇_𝖫𝜎𝑢_ℎ=𝐧_𝜎⋅Λ_𝖪∇_𝖪𝜎𝑢_ℎ+𝜑Λ_𝖪𝐧_𝜎=𝐧_𝜎⋅Λ_𝖫∇_𝖫𝜎𝑢_ℎ+𝜑Λ_𝖪𝐧_𝜎, (3.9)

from which we have𝐧_𝜎⋅Λ_𝖪−Λ_𝖫∇_𝖫𝜎𝑢_ℎ=𝜑𝐧_𝜎⋅Λ_K𝐧_𝜎. (3.10) We multiply (3.8) by 𝜇_𝖪 and (3.10) by 𝜇_𝖫, we sum the resulting relations and use (3.5) to have𝐧_𝜎⋅Λ_𝖪−Λ_𝖫𝜇_𝖪∇_𝖪𝜎𝑢_ℎ+𝜇_𝖫∇_𝖫𝜎𝑢_ℎ=𝜑𝐧_𝜎⋅𝜇_𝖪Λ_𝖫+𝜇_𝖫Λ_𝖪𝐧_𝜎. (3.11)

Solving this last equation for 𝜑 gives us the formula:𝐧𝜑=^𝜎⋅Λ_𝖪−Λ_𝖫∇_𝖪𝖫𝑢_ℎ𝐧_𝜎⋅𝜇_𝖪Λ_𝖫+𝜇_𝖫Λ_𝖪𝐧_𝜎. (3.12) Now, we derive an expression for Λ_𝖪𝖫. In view of (3.5) and (3.3), it holds that𝐧_𝜎⋅Λ_𝖪𝖫∇_𝖪𝖫𝑢_ℎ=𝐧_𝜎⋅𝜇_𝖪Λ_𝖪∇_𝖪𝜎𝑢_ℎ+𝜇_𝖫Λ_𝖫∇_𝖫𝜎𝑢_ℎ. (3.13) First, in (3.13), we substitute the expression for ∇_𝖫𝜎(𝑢_ℎ) provided by (3.4), we collect the factor ∇_𝖪𝜎(𝑢_ℎ) and we obtain:𝐧_𝜎⋅Λ_𝖪𝖫∇_𝖪𝖫𝑢_ℎ=𝐧_𝜎⋅𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫∇_𝖪𝜎𝑢_ℎ+𝜑𝜇_𝖪𝐧_𝜎⋅Λ_𝖫𝐧_𝜎. (3.14) Second, in (3.13), we substitute the expression for ∇_𝖪𝜎(𝑢_ℎ) provided by (3.4), we collect the factor ∇_𝖫𝜎(𝑢_ℎ) and we obtain𝐧_𝜎⋅Λ_𝖪𝖫∇_𝖪𝖫𝑢_ℎ=𝐧_𝜎⋅𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫∇_𝖫𝜎𝑢_ℎ−𝜑𝜇_𝖫𝐧_𝜎⋅Λ_𝖪𝐧_𝜎. (3.15) Third, we multiply (3.14) by 𝜇_𝖫, (3.15) by 𝜇_𝖪, we add the two resulting equations, we use again (3.5) and (3.6), and we obtain𝐧_𝜎⋅Λ_𝖪𝖫∇_𝖪𝖫𝑢_ℎ=𝐧_𝜎⋅𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫∇_𝖪𝖫𝑢_ℎ+𝜑𝜇_𝖪𝜇_𝖫𝐧_𝜎⋅Λ_𝖫−Λ_𝖪𝐧_𝜎. (3.16) Finally, in (3.16) we substitute the scalar factor 𝜑 given by (3.12), we collect the vectors 𝐧_𝐹 and ∇_𝖪𝖫(𝑢_ℎ) and we obtain:𝐧_𝜎⋅Λ_𝖪𝖫∇_𝖪𝖫𝑢_ℎ=𝐧_𝜎⋅𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫+𝜇_𝖪𝜇_𝖫+𝐧_𝜎⋅Λ_𝖫−Λ_𝖪𝐧_𝜎𝐧_𝜎⋅𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫𝐧_𝜎Λ_𝖪−Λ_𝖫∇_𝖪𝖫𝑢_ℎ. (3.17) Equation (3.17) suggests us to setΛ_𝖪𝖫=𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫+𝛽_𝖪𝖫Λ_𝖪−Λ_𝖫, (3.18) where𝛽_KL=𝐧_𝜎⋅Λ_𝖫−Λ_𝖪𝐧_𝜎𝐧_𝜎⋅Λ_𝖪/𝜇𝖪+Λ_𝖫/𝜇𝖫𝐧_𝜎. (3.19)

4. Numerical Experiments

In this section, we present and discuss a set of numerical results to show the convergence behavior of the finite volume method when we compute the diffusive flux using the technique described in the previous sections. In all the test case that we present in this section, we compare the behavior of the new numerical treatment that we consider in this paper with the weighted average, which is the standard approach in the method [5, 41, 42]. In fact, from the formulation given in equation (3.18), it follows that the average diffusion tensor 𝜆_𝖪𝖫 is equal to the weighted average 𝜇_𝖪Λ_𝖪+𝜇_𝖫Λ_𝖫 plus a correction term. This correction term is specifically designed to take care of possible discontinuities in the diffusion coefficients.

The finite volume formulation based on the least squares reconstruction of vertex values leads to a linear system for the cell-average unknowns whose coefficients matrix is generally unsymmetric although displaying a symmetric nonzero pattern. The positive definiteness of such system is still an open issue and this fact is also the major difficulty for the development of a full theory of convergence of such method. Theoretical results that prove coercivity are available only for meshes of slightly deformed quadrilaterals [2–4].

We solve such linear system using the MA41 routine of the HSL software collection, which implements an unsymmetric multifrontal sparse LU factorization technique especially designed for matrices with a symmetric nonzero pattern and unsymmetric values. Different software packages like UMFPACK and standard Krylov methods like BiCG-Stab and GMRES can be used in alternative. Efficiency issue is beyond the scope of our investigation but more details and comparison of performance when implementing different linear algebra solvers are found in the benchmark of the FVCA-6 Conference held in Prague, Czech Republic, in June 2011, [46].

We solve (2.1)-(2.2) on the domain Ω=]0,1[×]0,1[ for the data specified in the three test cases reported below. For all calculations, we measure the following relative errors:

error on the solution:ℰ_𝑢=∑_𝖪∈Ωℎ|𝖪|𝑢𝐱_𝖪−𝑢_𝖪²∑_𝖪∈Ωℎ|𝖪|𝑢𝐱_𝖪^21/2, (4.1)

error on the gradient:ℰ_∇𝑢=∑_𝜎∈Ωℎ|𝜎|∇𝑢𝐱_𝜎−∇_𝖪𝖫𝑢_ℎ²∑_𝖪∈Ωℎ|𝜎|∇𝑢𝐱_𝜎^21/2. (4.2)

Test Case 1
This test case is taken from [47] and is devoted to confirm that the new treatment of tensor coefficients and the standard weighted arithmetic average show the same behavior when the diffusion coefficients are regular, for example, constant or smoothly varying functions of position, even if small anisotropies along the principal direction of diffusion are present.The exact solution that we want to approximate is given by 𝑢𝑥,𝑦=sin(2𝜋𝑥)sin2𝜋𝑦+𝑥³+𝑥𝑦². (4.3)We consider two different constant diffusion tensor. The first one, called Λ_iso, is isotropic, while the second one, called Λ_ani, is anisotropic.The two diffusion tensors are given by: Λ_iso=10.250.251,Λ_ani=10.10.10.25. (4.4) We solve this test case on mesh family ℳ₁. Each mesh of such family is built as follows. First, we remap the position (𝑥,𝑦) of the nodes of an 𝑛×𝑛 uniform partition by the smooth coordinate transformation: 1𝑥=𝑥+(),110sin2𝜋𝑥sin2𝜋𝑦𝑦=𝑦+().10sin2𝜋𝑥sin2𝜋𝑦 (4.5)The corresponding mesh of ℳ₁ is built from this “primal” mesh by splitting each quadrilateral cell into two triangles and connecting the barycenters of adjacent triangular cells by a straight segment. The mesh construction is completed at the boundary Γ by connecting the barycenters of the triangular cells close to Γ to the midpoints of the boundary edges and these latters to the boundary vertices of the “primal” mesh. The first and the second mesh of mesh family ℳ₁ are displayed in Figure 2; mesh data are given in Table 1. Approximation errors for the isotropic diffusion tensor Λ_iso are shown in Figure 6; approximation errors for the anisotropic tensor Λ_ani are shown in Figure 7. In both plots, we also show the exact slopes proportional to ℎ and ℎ². From Figures 6 and 7, we deduce that in the case of constant diffusion tensors the performance of the new diffusion average proposed in this work and the weighted arithmetic average are almost the same. This behavior is confirmed both in the case of an isotropic diffusion tensor and in the case of an anisotropic diffusion tensor.

Table 1 

Data of mesh family ℳ₁; l is the mesh label, 𝒩_𝖪 is the number of cells, 𝒩_𝜎 is the number of mesh edges, 𝒩_𝖵 is the number of mesh nodes, and ℎ is the mesh size parameter.

(a)

(b)

Figure 2 

First and second mesh of mesh family ℳ1. This mesh family is used in test case 1.

Test Case 2
In this test case, we show the behavior of the new technique that is proposed in this paper when the diffusion coefficients are discontinuous along an internal line of the computational domain. To such purpose, we split Ω=∪²_𝑖=1Ω_𝑖 where Ω₁=1𝑥,𝑦∈Ω∶0≤𝑥≤1,0≤𝑦≤2,Ω₂=1𝑥,𝑦∈Ω∶0≤𝑥≤1,2<𝑦≤1. (4.6)The diffusion tensor is discontinuous across the horizontal line 𝑦=1/2 and is given by Λ=𝜆𝑥,𝑦¹[]1Λforany𝑥∈0,1,𝑦∈0,2𝜆₂[]1Λforany𝑥∈0,1,𝑦∈2,1, (4.7) where Λ=1+𝑥²+𝑦²−1𝑥−21𝑦−2−1𝑥−21𝑦−21+𝑥²+𝑦², (4.8)and 𝜆₁=1, 𝜆₂=10. The exact solution is continuous across the line 𝑦=1/2 and is designed to ensure flux conservation, that is, continuity of the normal component of the flux field. The exact solution is given by 𝑢𝑥,𝑦=1+𝑥+3𝑦³[]1forany𝑥∈0,1,𝑦∈0,23𝜆₂−𝜆₁2³𝜆₂𝜆+𝑥+3¹𝜆₂𝑦³[]1forany𝑥∈0,1,𝑦∈2,1. (4.9) This test case is solved on the two mesh families described below.Mesh Family ℳ₂
This mesh sequence is the first of the mesh collection of the two-dimensional benchmark of the conference “Finite Volumes for Complex Applications—V’’ held in Aussois (France) in 2008. The first mesh and the first refined mesh of this mesh suite are shown in Figure 3; mesh data are reported in Table 2.Mesh Family ℳ₃
Below 𝑦=1/2, we consider a regular mesh formed by 2𝑛×𝑛 squares, while above 𝑦=1/2 we consider a regular mesh formed by 4𝑛×2𝑛 squares. All the mesh nodes except those located at the boundaries and those located at the internal discontinuity line 𝑦=1/2 are then displaced by perturbing their position to a random position inside a square box centered at that original node position. The sides of this square box are aligned with the coordinate axis and their length is equal to 0.8/𝑛 (note that 1/𝑛 is the distance between two consecutive nodes in each direction). Instead, the nodes lying at 𝑦=1/2 are allowed to change only the abscissa in order not to modify the shape of the interface line while the nodes at the boundary are not displaced. We treat a nonmatching mesh as a conformal polygonal mesh modifying the shape of the polygons that are immediately below the line 𝑦=1/2: these polygons are treated as degenerate pentagons with two parallel consecutive edges. The first mesh is built by taking 𝑛=10 and mesh refinement is implemented by doubling the value of 𝑛 at each refinement step, thus implying that the mesh size parameter is approximately halved when passing from one mesh to the next one in the refinement process.The first mesh and the first refined mesh of this mesh suite are shown in Figure 4; mesh data are reported in Table 3.
The numerical results for the gradient and the solution approximation on mesh family ℳ₂ are shown in the two plots of Figure 8 and on mesh family ℳ₃ in the two plots of Figure 9. In the bottom-left corner of both plots, we show the exact slopes proportional to ℎ^1/2 (gradient errors) and to ℎ and ℎ² (solution errors). The solution gradient has a discontinuous normal component across such lines, and, due to this lack of regularity, the convergence rate cannot be expected to be better than that of a first-order scheme. This behavior is reflected in both left plots of Figures 8 and 9. Even if the convergence rate of the gradient is the same in the two cases, the approximation errors of the formulation using the new average are smaller than those obtained when using the standard weighted arithmetic average. Concerning the solution approximation, the results are even more spectacular because adopting the new average allows us to recover the second-order convergence rate, thus confirming the superior behavior of the new method.

Table 2 

Data of mesh family ℳ₂; l is the mesh label, 𝒩_𝖪 is the number of cells, 𝒩_𝜎 is the number of mesh edges, 𝒩_𝖵 is the number of mesh nodes, and ℎ is the mesh size parameter.

Table 3 

Data of mesh family ℳ₃; l is the mesh label, 𝒩_𝖪 is the number of cells, 𝒩_𝜎 is the number of mesh edges, 𝒩_𝖵 is the number of mesh nodes, and ℎ is the mesh size parameter.

(a)

(b)

Figure 3 

First and second mesh of mesh family ℳ2. This mesh family is used in test case 2.

(a)

(b)

Figure 4 

First and second mesh of mesh family ℳ3. This mesh family is used in test case 2.

Test Case 3
In this test case we aim at confirming the behavior of test case 2 for a more difficult case in which the diffusion tensor has intersecting discontinuities with an internal cross point.To such purpose, we split the computational domain as Ω=∪⁴_𝑖=1Ω_𝑖 where Ω₁=1𝑥,𝑦∈Ω∶0≤𝑥≤21,0≤𝑦≤2,Ω₂=1𝑥,𝑦∈Ω∶21≤𝑥≤1,0≤𝑦≤2,Ω₃=1𝑥,𝑦∈Ω∶21≤𝑥≤1,2,Ω<𝑦≤1₄=1𝑥,𝑦∈Ω∶0≤𝑥≤2,12<𝑦≤1. (4.10)The diffusion tensor is discontinuous across the horizontal line 𝑦=1/2 and the vertical line 𝑥=1/2 and is given by Λ_𝑖=1𝛼_𝑖11+𝑥−21𝑦−2−()1𝑥−1/2𝑦−2−1𝑥−21𝑦−211+𝑥−21𝑦−2for𝑥,𝑦∈Ω_𝑖,𝑖=1,…,4,, (4.11) where 𝛼₁=𝛼₃=1 and 𝛼₂=𝛼₄=100. The exact solution is continuous across such line and is designed to ensure flux conservation, that is, continuity of the normal component of the flux field. The exact solution is given by 𝑢𝑥,𝑦=2𝛼_𝑖(11−𝑥)2−𝑥cos2𝜋𝑦for𝑥,𝑦∈Ω_𝑖,𝑖=1,…,4. (4.12)This test case is solved on the mesh family ℳ₄ 5 that is built by displacing randomly all the nodes except those at the subdomain interfaces 𝑥=1/2 and 𝑦=1/2 of a 𝑛×𝑛 regular partition of Ω. We consider five meshes with 𝑛=10,20,40,80,160; mesh data are shown in Table 4, and the first mesh and the first refined mesh of this mesh suite are shown in Figure 5.
The numerical results for the gradient and the solution approximation on mesh family ℳ₄ are shown in the two plots of Figure 10. In the bottom-left corner of these plots, we show the exact slopes proportional to ℎ^1/2 (gradient errors) and to ℎ and ℎ² (solution errors). As in test case 2, the normal component of the solution gradient is discontinuous across the internal lines 𝑥=1/2 and 𝑦=1/2. Therefore, the convergence rate that is numerically measured in this test case is proportional to ℎ^1/2, cf. the left plot of Figure 10. This fact is independent of the average technique that is used for the diffusion coefficients. Nonetheless, the approximation errors given by the new average are smaller than those obtained when using the standard weighted arithmetic average. Concerning the solution approximation, such results confirms the behavior seen in test case 2. Consequently, the finite volume approximation of the scalar solution recovers the optimal second-order convergence rate when we apply the new average technique.

Table 4 

Data of mesh family ℳ₄; l is the mesh label, 𝒩_𝖪 is the number of cells, 𝒩_𝜎 is the number of mesh edges, 𝒩_𝖵 is the number of mesh nodes, and ℎ is the mesh size parameter.

(a)

(b)

Figure 5 

First and second mesh of mesh family ℳ4. This mesh family is used in test case 3.

(a)

(b)

Figure 6 

Test case 1: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)-(2.2) is solved using mesh family ℳ1 and the isotropic diffusion tensor Λiso in (4.4).The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. Note that the two error curves are superimposed. In the bottom-left corners of both plots, we show the exact slopes proportional to ℎ and ℎ2.

(a)

(b)

Figure 7 

Test case 1: error curves for the gradient approximation (left plot) and the solution approximation (right plot). We use the anisotropic diffusion tensor Λani of (4.4). The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. Note that the two error curves are superimposed. In the bottom-left corners of both plots, we show the exact slopes proportional to ℎ and ℎ2.

(a)

(b)

Figure 8 

Test case 2: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)-(2.2) is solved using mesh family ℳ2 and the discontinuous diffusion tensor (4.7).The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. In the bottom-left corner of both plots, we show the exact slopes proportional to ℎ1/2 (gradient errors) and to ℎ and ℎ2 (solution errors).

(a)

(b)

Figure 9 

Test case 2: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)–(2.2) is solved using mesh family ℳ3 and the discontinuous diffusion tensor (4.7). The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. In the bottom-left corner of both plots we show the exact slopes proportional to ℎ1/2 (gradient errors) and to ℎ and ℎ2 (solution errors).

(a)

(b)

Figure 10 

Test case 3: error curves for the gradient approximation (left plot) and the solution approximation (right plot). Problem (2.1)-(2.2) is solved using mesh family ℳ4 and the discontinuous diffusion tensor (4.11). The errors obtained by average (3.18)-(3.19) are marked by circles, the errors obtained by using the weighted arithmetic average are marked by squares. In the bottom-left corner of both plots we show the exact slopes proportional to ℎ1/2 (gradient errors) and to ℎ and ℎ2 (solution errors).

5. Conclusions

The numerical treatment of diffusion coefficients is an open problem and major problem, as, for example, the conductivity of different layers of soil is abruptly spatially variable. Moreover, it is particularly important in view of the problems of different scales involved and of random fields of the coefficients themselves. To this purpose, we proposed and tested a new technique for the numerical treatment of the conductivity tensor that interpolates the information coming from the different sides of a mesh interface shared by two adjacent cells. We point it out that this design is particularly suitable to the case of discontinuous conductivity tensors since the interpolation automatically adjusts its value by introducing an appropriate correction to the standard weighted average. Numerical experiments confirm this behavior. Furthermore, it may result in a particularly efficient strategy for the numerical resolution of problems where both the diffusive/dispersive and convective phenomena are simultaneously significant.

Acknowledgments

This work was funded by the PRIN 2007 project “Experimental measurement of the soil-vegetation-atmosphere interaction processes and numerical modelling of their response to climate change”. The authors also thank Dr. G. Manzini for the meshes used in test case 1 and test case 3.