# The Finite Volume Formulation for 2D Second-Order Elliptic

Authors: Stefano Ferraris; Ivan Bevilacqua; Davide Canone; Davide Pognant; Maurizio Previati
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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ℝ𝕋0-ℙ0 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 Section2, 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 fieldu given by(2.1)-div⁡(Λ∇u)=finΩ,(2.2)u=gon⁡Γ,
where Λ is the diffusion tensor describing the material properties, f is the forcing term and g the boundary function that defines nonhomogeneous Dirichlet conditions on the boundary Γ. Under suitable regularity assumptions on f, g, Λ, 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 Λ(x) is normally represented by a strictly positive definite matrix for every x 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 {Ωh}h of the domain Ω. For any mesh Ωh, the subscripted label h is the mesh size and is defined by(2.3)h=max⁡K∈ΩhhK,hK=sup⁡x′,x''∈K|x′-x′′|,wherehK is the diameter of the polygonal cell  K∈Ωh. We also label the numerical solution uh calculated on the mesh Ωh by h. Also,(i)
a genericmesh vertex is denoted by  V  and its coordinate vector by xV. We will also find it convenient to introduce a local numbering of the vertices, for example, Vi, and to ease notation to denote the vertex position of the i-th vertex by xi (instead of xVi);(ii)
a genericcell interface or a boundary edge is denoted by σ, its center (i.e., its edge midpoint) by xσ, and its measure (edge length) by |σ|;(iii)
a genericpolygonal cell is denoted by K its measure (area) by |K|, its center of gravity by xK, and its boundary by ∂K.The orientation of each mesh interfaceσ is reflected by its unit normal vector nσ, which is fixed once and for all. For any mesh cell K and any face σ of the polygonal boundary ∂K, we define the unit normal vector nKσ that points out of K and we also use the notation NKσ=|σ|nKσ.Let the algebraic vectoruh=(uK)K∈Ωh be the numerical solution, where each uK approximates the cell average of the scalar solution u over the cell K. The finite volume scheme for uh on the computational mesh Ωh reads(2.4)-∑σ∈∂KΛσ∇σuh⋅NKσ=fK|K|∀K∈Ωh,
where  Λσ is an evaluation of Λ at face σ taken from the side of K;  ∇σ is the discrete face gradient built inside the diamond cell centered at face σ; NK,σ is the geometric vector perpendicular to σ, pointing from K to L, and with lenght |NKσ|=|σ|;  fK is the cell average of the right-hand side term f over the cell K:(2.5)fK:=1|K|∫Kf(x)dV.LetΛK and ΛL be first-order approximations of the diffusion tensor Λ within the control volumes K and L. For example, we can take ΛK:=Λ(xK) and ΛL:=Λ(xL), where xK and xL are internal points (not necessarily the barycenters) of  K and L, respectively. Now, we can consider either Λσ=ΛK inside K and Λσ=ΛL inside L or a suitable mean of these two tensors. In particular, we can deal with the arithmetic mean:(2.6)Λσ=12(ΛK+ΛL),or with theharmonic mean:(2.7)Λσ-1=12(ΛK-1+ΛL-1).Both approaches deserve a special care whenΛ is discontinuous across σ and a proper definition of the discrete gradient ∇σuh 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 D centered at the mesh face σ with vertices xi and xi+1 and shared by the cells with centers xK and xL. The diamond cell D can also be seen as the union of the subdiamonds DK and DL, which are the triangular cells sharing σ as common base and having, respectively, vertices  xK and xL, the centers of gravity of cells K and L. The four vectors NKi, NLi, NLi+1, and NKi+1 shown in Figure 1 are respectively orthogonal to the four boundary faces σKi, σLi, σLi+1, and σKi+1, and their lenght is equal to the lenght of the corresponding face. Instead, when σ is a boundary face, thus defined by σ=∂K∩Γ (Γ being the boundary of the computational domain Ω), the diamond cell associated to σ coincides with DK, that is, it is the triangle defined by σ and the vertex xK.Figure 1
Geometry of the2D 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 K and L. All these formulas can be easily adapted to the case of a boundary face by taking xL equal to the center of gravity of face σ.Using the Green-Gauss theorem yields:(2.8)∇σuh≈1|D|∫D∇u(x)dV=1|D|∫∂Du(x)n(x)dS=1|D|∑σ∈∂D∫σu(x)n(x)dS,
where n is the unit vector orthogonal to σ∈∂D and pointing out of D. If the restriction of u to the face σ of ∂D is an affine function, the boundary integral on ∂D in (2.8) only depends on the values of u at the vertexes of D and the constant vector provided by the formulas that we derive below must coincide with the gradient of u.The Gauss-Green theorem applied to the solution gradient∇u integrated on the 2D diamond cell D of Figure 1 yields:(2.9)∑σ∈∂D∫σu(x)n(x)dS=12((u(xK)+u(xi))NKi(u(xK)+u(xi))NKi+(u(xi)+u(xL))NLi+(u(xL)+u(xi+1))NLi+1+(u(xi+1)+u(xK))NKi+1(u(xK)+u(xi))NKi).
Let us introduce the vectors NKL and Nii+1, which are, respectively, orthogonal to the edge connecting xK to xL and xi to xi+1, and whose lenght are equal to the length of these edges. Using the geometric relations NKL=-(NKi+NLi) and NKL=-(NKi+NLi) into (2.9) and rearranging the terms yields: (2.10)∫∂D∇u(x)⋅n(x)dS=12((u(xi+1)-u(xi))(u(xi+1)-u(xi))NKL+(u(xL)-u(xK))Nii+1(u(xi+1)-u(xi))).To derive the gradient formula, we replace the function valuesu(xK) and u(xL) with the cell unknowns uK and uL, respectively, and the function values u(xi) and u(xi+1) with the corresponding nodal unknowns ui and ui+1. We obtain(2.11)(uh):=12|D|(∇KL(uL-uK)NKL+(ui+1-ui)Nii+1).
which is a constant approximation of the solution gradient ∇u within the quadrilateral cell D.With any vertexxi of the mesh Ωh, we associate the reconstructed vertex value ui=∑K∈ΩiωKiuK,  where  Ωi={K:xi∈K} denotes the subset of the mesh cells which share the vertex xi. The interpolation weights ωKi are assumed to verify the consistency relations [3]:(2.12)∑K∈ΩiωKi=1,∑K∈ΩiωKi(xi-xK)=0.
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ωKi are obtained by solving the reconstruction problem that approximates the cell-averaged data set {(xK,uK)K∈Ωi} by the affine function:(2.13)ũi(x)=α+β⋅(xK-xi)forx∈Vi
on the covolume 𝒱i=⋃K∈ΩiK and in a least square sense, compared to [3, 41]. The reconstructed value at vertex xi is now given by taking ui=ũi(xi)=α. The coefficients (α,β) are the minimizers of the least squares functional:(2.14)Ji(α,β)=∑K∈Ωi(α+β⋅(xK-xi)-uK)2.
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 β=(βx,βy); a straightforward calculation yields:(2.15)∂Ji∂α=2∑K∈Ωi(α+β⋅(xK-xi)-uK)=0,∂Ji∂βx=2∑K∈Ωi(α+β⋅(xK-xi)-uK)(xK-xi)=0,∂Ji∂βy=2∑K∈Ωi(α+β⋅(xK-xi)-uK)(yK-yi)=0.Letmi be the number of cells of Ωi; we reformulate the linear system written above as(2.16)∑K∈Ωimiα+∑K∈Ωi(xK-xi)βx+∑K∈Ωi(yK-yi)βy=∑K∈ΩiuK,(2.17)∑K∈Ωi(xK-xi)α+∑K∈Ωi(xK-xi)2βx+∑K∈Ωi(xK-xi)(yK-yi)βy=∑K∈Ωi(xK-xi)uK,(2.18)∑K∈Ωi(yK-yi)α+∑K∈Ωi(xK-xi)(yK-yi)βx+∑K∈Ωi(yK-yi)2βy=∑K∈Ωi(yK-yi)uK.
For convenience, we consider a suitable local numbering of the cells of 𝒱i, for example, 𝒱i=K1∪K2∪…Kmi, and we introduce the two matrices:(2.19)A=(mi∑K∈Ωi(xK-xi)∑K∈Ωi(yK-yi)∑K∈Ωi(xK-xi)∑K∈Ωi(xK-xi)2∑K∈Ωi(xK-xi)(yK-yi)∑K∈Ωi(yK-yi)∑K∈Ωi(xK-xi)(yK-yi)∑K∈Ωi(yK-yi)2),B=(11…1xK1-xixK2-xi…xKmi-xiyK1-yiyK2-yi…yKmi-yi),and the column vectoru=(u1,u2,…,umi)T, which collects the solution averages to be imposed. Using such definitions, the linear system (2.16)–(2.18) becomes(2.20)A(αβxβy)=Bu,whichimpliesthat(αβxβy)=A-1Bu,since matrix  A is, obviously, symmetric and positive definite and, thus, nonsingular. The coefficient α that provides the vertex value is given by(2.21)α=(1,0,0)(αβxβy)=(1,0,0)A-1Bu.
We require that such a coefficient be an average of the data in u with coefficients ω=(ωK1i,ωK2i,…,ωKmii)T, that is, α=ωTu. By comparison with (2.21), we have the final weight formula:(2.22)ω=BTA-T(100).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 ui at the vertexes xi∈Γ on the Dirichlet boundary are constrained to the boundary data, for instance ui=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 K and L, that is, σ=K∣L. 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)∫σΛ∇(u)⋅ndS≈ΛKσ∇σ(uh)⋅NKσcannot be consistent wheneverΛKσ is some kind of average between ΛK and ΛL and the constant vector ∇σ(uh)  approximates ∇(u)  over all the diamond cellhull(xK,xL,σ). To tackle this problem, we consider two distinct approximations of the solution gradient within the subcells DK and DL, denoted by  ∇Kσ(uh)  and  ∇Lσ(uh), respectively, and impose directly the flux conservation. To derive an expression for such one-sided gradients, we introduce an additional unknown uσ along face σ and we apply the Gauss-Green Theorem as for the derivation of ∇KLuh in the previous section. Also, we recall that NKσ=NKL and |NKσ|=|σ| and we use a similar notation for the normal vector from the other side of σ, so that NLσ=-NKL and |NLσ|=|σ|. The two one-sided gradient formulas are(3.2)∇Kσ(uh):=12|DK|((uσ-uK)NKσ(uσ-uK)NKσ+(ui+1-ui)Nii+1(uσ-uK)NKσ),∇Lσ(uh):=12|DL|((uσ-uK)NKσ(uσ-uL)NLσ-(ui+1-ui)Nii+1(uσ-uK)NKσ).Now, we search for a tensorΛKσ, which is an average of the diffusion tensors Λk and ΛL, and a gradient vector ∇KL, which is an average of the gradient vectors ∇Kσ(uh)  and  ∇Lσ(uh), that makes it possible to express directly the flux conservation as(3.3)ΛKσ∇KL(uh)⋅nσ:=ΛK∇Kσ(uh)⋅nσ=ΛL∇Lσ(uh)⋅nσ.
Since we expect that the normal component of the vectors  ∇Kσ(uh) and ∇Lσ(uh) be the same, we impose directly this condition by requiring that a scalar coefficient φ exists such that(3.4)∇Lσ(uh)=∇Kσ(uh)+φnσ,
(recall that nσ=nKL). We also assume that ∇KL(uh) be a weighted average of ∇Kσ(uh) and ∇Lσ(uh). Let us be given two nonnegative coefficients μK and  μL such that(3.5)μK+μL=1,
and such that(3.6)∇KL(uh)=μK∇Kσ(uh)+μL∇Lσ(uh).To this purpose, a number of possible choices have been proposed and widely investigated in the literature forμK and μL. For example, a popular option is to take the volume of the two subdiamonds DK and DL normalized by the total volume of the diamond, for example, μK=|DK|/|D| and μL=|DL|/|D|. The diffusion tensor can be chosen as the simple average of ΛK and ΛL, for example, ΛKL=μKΛK+μLΛL, or as the harmonic average, for example, ΛKL-1=μKΛK-1+μLΛL-1. 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 K and L. Our aim in the next developments is at determining a value for φ and of the average diffusion tensor ΛKL in terms of ΛK and ΛL that ensures flux conservation (3.3) for any given pair of coefficients μK and μL. 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 ΛLnσ and use flux conservation to obtain(3.7)nσ⋅ΛL∇Lσ(uh)=nσ⋅(ΛL∇Kσ(uh)+φΛLnσ)=nσ⋅ΛK∇Kσ(uh),from which we have(3.8)nσ⋅(ΛK-ΛL)∇Kσ(uh)=φnσ⋅ΛLnσ.
Likewise, we multiply (3.4) by ΛKnσ and use flux conservation to obtain(3.9)nσ⋅ΛK∇Lσ(uh)=nσ⋅(ΛK∇Kσ(uh)+φΛKnσ)=nσ⋅(ΛL∇Lσ(uh)+φΛKnσ),from which we have(3.10)nσ⋅(ΛK-ΛL)∇Lσ(uh)=φnσ⋅ΛKnσ.
We multiply (3.8) by μK and (3.10) by  μL, we sum the resulting relations and use (3.5) to have(3.11)nσ⋅(ΛK-ΛL)(μK∇Kσ(uh)+μL∇Lσ(uh))=φnσ⋅(μKΛL+μLΛK)nσ.Solving this last equation forφ gives us the formula:(3.12)φ=nσ⋅(ΛK-ΛL)∇KL(uh)nσ⋅(μKΛL+μLΛK)nσ.
Now, we derive an expression for ΛKL. In view of (3.5) and (3.3), it holds that(3.13)nσ⋅ΛKL∇KL(uh)=nσ⋅(μKΛK∇Kσ(uh)+μLΛL∇Lσ(uh)).
First, in (3.13), we substitute the expression for ∇Lσ(uh) provided by (3.4), we collect the factor ∇Kσ(uh) and we obtain:(3.14)nσ⋅ΛKL∇KL(uh)=nσ⋅(μKΛK+μLΛL)∇Kσ(uh)+φμKnσ⋅ΛLnσ.
Second, in (3.13), we substitute the expression for ∇Kσ(uh) provided by (3.4), we collect the factor ∇Lσ(uh) and we obtain(3.15)nσ⋅ΛKL∇KL(uh)=nσ⋅(μKΛK+μLΛL)∇Lσ(uh)-φμLnσ⋅ΛKnσ.
Third, we multiply (3.14) by μL, (3.15) by μK, we add the two resulting equations, we use again (3.5) and (3.6), and we obtain(3.16)nσ⋅ΛKL∇KL(uh)=nσ⋅(μKΛK+μLΛL)∇KL(uh)+φμKμLnσ⋅(ΛL-ΛK)nσ.
Finally, in (3.16) we substitute the scalar factor φ given by (3.12), we collect the vectors nF and ∇KL(uh) and we obtain:(3.17)nσ⋅ΛKL∇KL(uh)=nσ⋅(μKΛK+μLΛL+μKμL+nσ⋅(ΛL-ΛK)nσnσ⋅(μKΛK+μLΛL)nσ(ΛK-ΛL))∇KL(uh).
Equation (3.17) suggests us to set(3.18)ΛKL=μKΛK+μLΛL+βKL(ΛK-ΛL),
where(3.19)βKL=nσ⋅(ΛL-ΛK)nσnσ⋅(ΛK/μK+ΛL/μL)nσ.


## 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 λKL is equal to the weighted average μKΛK+μLΛL 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:(4.1)Eu=(∑K∈Ωh|K||u(xK)-uK|2∑K∈Ωh|K||u(xK)|2)1/2,error on the gradient:(4.2)E∇u=(∑σ∈Ωh|σ||∇u(xσ)-∇KLuh|2∑K∈Ωh|σ||∇u(xσ)|2)1/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
(4.3)u(x,y)=sin⁡⁡(2πx)sin⁡⁡(2πy)+x3+xy2.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:
(4.4)Λiso=(10.250.251),Λani=(10.10.10.25).
We solve this test case on mesh family ℳ1. Each mesh of such family is built as follows. First, we remap the position (x̂,ŷ) of the nodes of an n×n uniform partition by the smooth coordinate transformation:
(4.5)x=x̂+(110)sin⁡(2πx̂)sin⁡(2πŷ),y=ŷ+(110)sin⁡(2πx̂)sin⁡(2πŷ).The corresponding mesh of ℳ1 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 ℳ1 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 h and h2. 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ℳ1; l is the mesh label, 𝒩K is the number of cells, 𝒩σ is the number of mesh edges, 𝒩V is the number of mesh nodes, and h is the mesh size parameter.
l𝒩K𝒩σ𝒩Vh11214002802.00810-1244114009601.07110-131681520035205.42210-24656120000134402.71910-252592178400524801.36110-2First and second mesh of mesh familyℳ1. This mesh family is used in test case 1.
(a)(b)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Ω=∪i=12Ωi where
(4.6)Ω1={(x,y)∈Ω:0≤x≤1,0≤y≤12},Ω2={(x,y)∈Ω:0≤x≤1,12<y≤1}.The diffusion tensor is discontinuous across the horizontal line y=1/2 and is given by
(4.7)Λ(x,y)={λ1Λ̃foranyx∈[0,1],y∈[0,12]λ2Λ̃foranyx∈[0,1],y∈[12,1],
where
(4.8)Λ̃=(1+x2+y2-(x-12)(y-12)-(x-12)(y-12)1+x2+y2),and λ1=1, λ2=10. The exact solution is continuous across the line y=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
(4.9)u(x,y)={1+x+3y3for anyx∈[0,1],y∈[0,12]3λ2-λ123λ2+x+3λ1λ2y3for anyx∈[0,1],y∈[12,1].
This test case is solved on the two mesh families described below.Mesh Familyℳ2
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 Figure3; mesh data are reported in Table 2.Mesh Familyℳ3
Belowy=1/2, we consider a regular mesh formed by 2n×n squares, while above y=1/2 we consider a regular mesh formed by 4n×2n squares. All the mesh nodes except those located at the boundaries and those located at the internal discontinuity line y=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/n (note that 1/n is the distance between two consecutive nodes in each direction). Instead, the nodes lying at y=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 y=1/2: these polygons are treated as degenerate pentagons with two parallel consecutive edges. The first mesh is built by taking n=10 and mesh refinement is implemented by doubling the value of n 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ℳ2 are shown in the two plots of Figure 8 and on mesh family ℳ3 in the two plots of Figure 9. In the bottom-left corner of both plots, we show the exact slopes proportional to h1/2 (gradient errors) and to h and h2 (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ℳ2; l is the mesh label, 𝒩K is the number of cells, 𝒩σ is the number of mesh edges, 𝒩V is the number of mesh nodes, and h  is the mesh size parameter.
l𝒩K𝒩σ𝒩Vh12505352861.86510-121000207010719.05110-234000814041414.69310-241600032280162812.40710-2564000128560645611.21510-2Table 3
Data of mesh familyℳ3; l is the mesh label, 𝒩K is the number of cells, 𝒩σ is the number of mesh edges, 𝒩V is the number of mesh nodes, and h is the mesh size parameter.
l𝒩K𝒩σ𝒩Vh15692372.50010-122243521291.25010-1389613764816.25010-243584544018573.12510-25143362163272971.56310-265734486272289297.81310-372293763445761152013.90610-3First and second mesh of mesh familyℳ2. This mesh family is used in test case 2.
(a)(b)First and second mesh of mesh familyℳ3. This mesh family is used in test case 2.
(a)(b)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Ω=∪i=14Ωi where
(4.10)Ω1={(x,y)∈Ω:0≤x≤12,0≤y≤12},Ω2={(x,y)∈Ω:12≤x≤1,0≤y≤12},Ω3={(x,y)∈Ω:12≤x≤1,12<y≤1},Ω4={(x,y)∈Ω:0≤x≤12,12<y≤1}.The diffusion tensor is discontinuous across the horizontal line y=1/2 and the vertical line x=1/2 and is given by
(4.11)Λi=1αi(1+(x-12)(y-12)-(x-1/2)(y-12)-(x-12)(y-12)1+(x-12)(y-12)for(x,y)∈Ωi,i=1,…,4,),
where α1=α3=1 and α2=α4=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
(4.12)u(x,y)=2αi(1-x)(12-x)cos⁡(2πy)for(x,y)∈Ωi,i=1,…,4.This test case is solved on the mesh family ℳ4 5 that is built by displacing randomly all the nodes except those at the subdomain interfaces x=1/2 and y=1/2 of a n×n regular partition of Ω. We consider five meshes with n=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ℳ4 are shown in the two plots of Figure 10. In the bottom-left corner of these plots, we show the exact slopes proportional to h1/2 (gradient errors) and to h and h2 (solution errors). As in test case 2, the normal component of the solution gradient is discontinuous across the internal lines x=1/2 and y=1/2. Therefore, the convergence rate that is numerically measured in this test case is proportional to h1/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ℳ4; l is the mesh label, 𝒩K is the number of cells, 𝒩σ is the number of mesh edges, 𝒩V is the number of mesh nodes, and h is the mesh size parameter.
l𝒩K𝒩σ𝒩Vh11002201211.79610-124008404419.15010-231600328016814.78810-2464001296065612.42310-252560051520259211.21610-2First and second mesh of mesh familyℳ4. This mesh family is used in test case 3.
(a)(b)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 h and h2.
(a)(b)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 h and h2.
(a)(b)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 h1/2 (gradient errors) and to h and h2 (solution errors).
(a)(b)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 h1/2 (gradient errors) and to h and h2 (solution errors).
(a)(b)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 h1/2 (gradient errors) and to h and h2 (solution errors).
(a)(b)


## 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.



---

*Source: 187634-2012-04-03.xml*