Convergence analysis of a cell centered finite volume diffusion operator on non-orthogonal polyhedral meshes

by   L. Bonaventura, et al.
Politecnico di Milano

A simple but successful strategy for building a discrete diffusion operator in finite volume schemes of industrial use is to correct the standard two-point flux approximation with a term accounting for the local mesh non-orthogonality. Practical experience with a variety of different mesh typologies, including non-orthogonal tetrahedral, hexahedral and polyhedral meshes, has shown that this discrete diffusion operator is accurate and robust whenever the mesh is not too distorted and sufficiently regular. In this work, we show that this approach can be interpreted as equivalent to introducing an anisotropic operator that accounts for the preferential directions induced by the local mesh non-orthogonality. This allows to derive a convergence analysis of the corrected method under a quite weak global assumption on mesh distortion. This convergence proof, which is obtained for the first time for this finite volume method widely employed in industrial software packages such as OpenFOAM, provides a reference framework on how to interpret some of its variants commonly implemented in commercial finite volume codes. Numerical experiments are presented that confirm the accuracy and robustness of the results. Furthermore, we also show empirically that a least square approach to the gradient computation can provide second order convergence even when the mild non-orthogonality condition on the mesh is violated.



There are no comments yet.


page 32


A second-order face-centred finite volume method on general meshes with automatic mesh adaptation

A second-order face-centred finite volume strategy on general meshes is ...

Weak Consistency of Finite Volume Schemes for Systems of Non Linear Conservation Laws: Extension to Staggered Schemes

We prove in this paper the weak consistency of a general finite volume c...

Boundedness-Preserving Implicit Correction of Mesh-Induced Errors for VoF Based Heat and Mass Transfer

Spatial discretisation of geometrically complex computational domains of...

Interplay between diffusion anisotropy and mesh skewness in Hybrid High-Order schemes

We explore the effects of mesh skewness on the accuracy of standard Hybr...

A fully local hybridised second-order accurate scheme for advection-diffusion equations

In this paper, we present a fully local second-order upwind scheme, appl...

Construction and analysis of the quadratic finite volume methods on tetrahedral meshes

A family of quadratic finite volume method (FVM) schemes are constructed...

A new Hodge operator in Discrete Exterior Calculus. Application to fluid mechanics

This article introduces a new and general construction of discrete Hodge...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Finite volume methods have been extremely popular in computational fluid dynamics (CFD) in the past and they still are an area of active research in numerical mathematics. Among the many different developments in this field, we recall the finite volume element scheme [9], the multi-point flux approximation schemes (MPFA) [1], [2], [3], or more recent variants, such as the mixed finite volume scheme (MFV) [16], the hybrid finite volume scheme (HFV) [23], [25] and the discrete duality finite volume schemes (DDFV) [4],[5], [11], [12], [32], [33], [34]. All these methods share indeed many common features, as discussed in [18].

In this work, however, we will focus on cell centered schemes, in which a single unknown is associated to each mesh cell. Cell centered finite volume methods are widely employed in industrial codes [7], [40] for a number of practical reasons. Indeed, they rely on relatively simple data structures, even for general unstructured meshes, and they allow for easy treatment of boundary conditions at singular boundary points, such as inner corners, while effectively handling general shapes of the computational domain. Cell centered methods can be naturally parallelized by domain decomposition techniques, guaranteeing minimal interprocessor communications, especially in their low order variants, due to the use of discrete operators built from local stencils. They allow an easy implementation of locally adaptive multilevel refinement strategies and they can be easily equipped with very efficient geometric multigrid procedures [43]. Finally, cell centered finite volume methods also allow an immediate extension to nonlinear coupled problems [24].

Known drawbacks of cell centered schemes are the reduced accuracy in strongly heterogeneous diffusion problems [25] with respect to MPFA, MFV or HFV schemes, as well as the only asymptotical recovery of the discrete Stokes formula, in contrast with the exact discrete property provided for example by DDFV schemes [11]. On the other hand, MPFA, MFV, HFV and DDFV schemes achieve such properties by introducing additional unknowns at selected mesh locations, thus implying an additional cost with respect to cell centered discretizations. It is still an open question if similar accuracy improvements can be obtained from cell centered schemes by introducing additional unknowns through local mesh refinement.

For these reasons, it is important to understand the analytical behaviour of cell centered finite volume discretizations on the typical non-orthogonal meshes practicaly required for industrial applications [7], [40]. For these applications, the so-called Gauss corrected scheme, widely adopted by finite volume practitioners [30], [36], [37], [39], [42], appears to be a simple, robust and sufficiently accurate option. Notice that this scheme can also be interpreted as a specific realization of the recently introduced asymmetric gradient discretization method [17].

To the best of the authors’ knowledge, the convergence properties of this finite volume method have never been analyzed in the case of non-orthogonal meshes. Indeed, convergence analyses of finite volume schemes for diffusion operators on unstructured mesh types are usually limited to polyhedral meshes satisfying an orthogonality condition [21], [22]. This is quite restrictive in practice, since none of the robust mesh generators usually adopted for pre-processing of industrial configurations are able to guarantee this condition.

In this work, we show that it is possible to prove the convergence of the Gauss corrected scheme on unstructured meshes satisfying a global and rather weak mesh regularity condition. This goal is achieved adapting the approach used in [22] for the convergence analysis of a cell-centered finite volume scheme for anisotropic diffusion problems on orthogonal meshes. A preliminary version of these results has been presented in [14]. Furthermore, we also show empirically that a least square approach to the gradient computation can provide second order convergence even when the mild non-orthogonality condition on the mesh is violated. It is to be remarked that existing convergence proofs for finite volume methods on non-orthogonal meshes either involve discretization schemes not guaranteeing local flux conservativity [24], [25]

, or DDFV schemes employing additional degrees of freedom

[4], [5], or two-dimensional diamond schemes on meshes satisfying more restrictive regularity conditions [13]. We will focus here on the isotropic steady state diffusion equation


We will assume that is a measurable function, , such that for a.e. , with , and for . The classical weak formulation of problem (1) consists in finding such that


Rather than proving convergence directly for the finite volume scheme associated to the strong problem formulation (1), we will identify a discrete weak formulation underlying the finite volume scheme and then to prove convergence of its solution to that of the continuous weak problem (26).

The paper is organized as follows. In section 2, several fundamental definitions of mesh related quantities and discrete functional spaces are introduced. In section 3, the cell centered finite volume method that is the focus of our analysis is presented. In section 4, the discrete weak formulation is recovered and in section 5, the convergence analysis of the Gauss corrected scheme is presented. In section 6, the results of some numerical experiments are reported. A proposal to overcome the constraints on the mesh for some specific three-dimensional mesh types is introduced in section 7. Finally, in section 8 some conclusions are drawn and some future developments are outlined.

2 Meshes and discrete spaces

The finite volume method is a mesh-based discretization technique suitable for any number of space dimensions, but in this work we only consider the case. Since the computational domains of practical interest are usually of complex geometry, the focus here is on meshes composed of arbitrarily shaped polyhedral cells, in the sense of the formal definition below, see also [19], [25].

Definition 2.1.

(Polyhedral mesh): Let be a bounded, open polyhedral subset of A polyhedral mesh for is denoted by , where the quadruple includes:

  1. is a finite family of non-empty, connected, polyhedral, open, disjoint subsets of called cells (or control volumes), such that . For any , is the boundary of , denotes the measure of , and is the diameter of , that is the maximum distance between two points in .

  2. is a finite family of disjoint subsets of representing the faces. Let be the set of interior faces such that, for all ,

    is a non-empty open subset of a hyperplane in

    with , and let be the set of boundary faces such that, for all , is a non-empty open subset of . It is assumed that, for any , there exists a subset such that . The set of cells sharing one face is . It is assumed that, for all , either has exactly two elements and then , or has exactly one element and then . For all , denotes the -dimensional measure of , and is the barycenter of .

  3. is a family of points of indexed by , such that for all , and it is called the center of , possibly corresponding to its barycenter. It is assumed that all cells are -star-shaped, in the sense that if , then the line segment .

  4. is the finite set of vertices of the mesh. For , collects all the vertices belonging to , while for , collects all the vertices belonging to .

The size of the polyhedral mesh is defined as .

Furthermore, for any and for any ,

is the constant unit vector normal to

and outward to . For any , the set of neighbors of is denoted by


Additionally denotes the orthogonal distance between and


which is constant for all From the assumption that is -star-shaped, it follows that and that it also holds:


For all and , denotes the cone with vertex and basis , also called half-diamond, that is the volume defined by


For all , denotes the diamond associated to face , as in Figure 1.

Definition 2.1 covers a wide range of meshes, including meshes with non-convex cells, with non-planar faces requiring triangulation, or with hanging nodes. Furthermore, Definition 2.1 also includes tetrahedral and hexahedral meshes as particular cases, as well as meshes with wedge and pyramidal cells.

Figure 1: Non-orthogonal generally polyhedral mesh: the line segment is not aligned with the face normal unit vector . Furthermore, the intersection point between the face and the vector is not necessarily coincident with the face centroid . The half diamonds and are represented with dashed and dotted lines, respectively.

Finite volume methods are traditionally introduced in discrete functional spaces of piecewise constant functions [21]. In recent analyses [27], [28]

, associated inner products, norms and seminorms are exploited to recast the discrete flux balance equations into an equivalent variational form, which naturally allows to derive stability estimates and to investigate the numerical convergence of specific schemes

[25]. In the classical finite volume framework, the discrete flux balance equation corresponding to problem (1) takes the form


where the face flux is such that

and denotes the infinitesimal face area element. A relevant feature of the scheme is the flux conservativity property


which is assumed to hold for all interior faces , where and are the cells sharing the face .

The convergence analysis of cell centered finite volume schemes on arbitrary polyhedral meshes [25] may also require to introduce the space , which consists of the functions that are piecewise constant on each cell . For all and for all , the constant value of in is denoted by . Consequently, discrete functional analysis results for the convergence of finite volume schemes [21], [25] can be exploited.

In addition, in order to introduce proper test functions to check the convergence of the discrete solution to the continuous solution of the weak formulation, for all functions a projection operator is defined, such that .

3 A cell centered diffusion scheme for non-orthogonal meshes

The vast majority of finite volume schemes for diffusive problems are based on the application of the discrete Gauss theorem. The numerical approximation is derived as


where the numerical flux through face is computed as


and depends on the definition of the the face normal gradient. Usually,

is approximated by the surface value interpolation

, obtained from standard interpolation schemes. Linear interpolation is often chosen to preserve second order accuracy, while harmonic interpolation is sometimes selected, especially when the scalar diffusivity field is strongly non-homogeneous [21]. A variety of alternative schemes can be constructed to approximate the face normal gradient , each with its own specific features. Most of them are traditionally studied empirically, by directly testing them on specific meshes and representative flow problems [29], [41].

The simplest scheme for the face normal gradient is represented by the two-point flux approximation [21]


Even though unconditionally monotone and coercive [15], it is of limited accuracy on unstructured meshes, where mesh non-orthogonality may lead to severe errors in the approximation of the diffusion fluxes [20]. In order to compensate for the unavoidable non-orthogonality of realistic unstructured meshes, a simple but effective solution is provided by the Gauss corrected scheme, which introduces a non-orthogonal correction term [36] in the two-point flux scheme, thus obtaining the approximation


Here, the first term corresponds to the two-point flux contribution in Eq.(11), expressing the diffusion flux component in the direction of the line segment , while the second term accounts for the local mesh non-orthgonality across the face , expressed as the difference between the correct face normal diffusive flux estimated from a proper face gradient and the diffusive flux along the direction of . In order to avoid oscillatory solutions [41], it is important that the gradient at face is evaluated using a different scheme from the one employed in the first term of Eq.(12). Thus, the gradient is usually estimated at face by interpolation of the neighbouring cells gradients and for , which can be either the standard linear interpolation or, for increased simplicity, the midpoint rule. Indeed, if the cell derivatives are linear approximations, the diffusion flux will be more accurate than first order on very regular meshes [36]. The Gauss corrected scheme allows more accurate approximations than the two-point flux approximation (11), but it is not, in general, unconditionally coercive on arbitrary unstructured meshes. As a consequence, on irregular meshes it may become a source of numerical instability. On orthogonal grids, this scheme reduces to the classical two-point flux scheme, since the correction term vanishes.

Following [10], on a general unstructured polyhedral mesh like that of Definition 2.1, the centered discrete gradient operator is defined as the piecewise constant function


for Since for any closed control volume the geometrical relations


hold, with , Eq.(13) is also equal to


which is easily recognized as the finite volume discretization of the gradient based on the Gauss theorem [30]. For this reason, the gradient approximation in Eq.(13) is often identified as the Gauss gradient scheme.

The consistency of the discrete gradient in Eq.(13) has been analyzed in [25]. It stems directly from the geometrical identity


where is the transpose of the vector , see Figure 1, and

is the identity matrix. For any affine function

defined by , with and , assuming that and , it results that . Hence, expression (13) leads to , which amounts to linear exactness for any affine function on , provided that , which is verified whenever , , see Figure 1.

Finally, if the face gradient in Eq.(12) is computed using a linear interpolation operator applied to the cell gradients reconstructed via the Gauss scheme (15) from both cells sharing the face , the non-orthogonal correction term in Eq.(12) is associated to a large stencil which includes, besides cells and sharing face , all their neighbouring cells .

By applying the Gauss corrected scheme from Eq.(12) to the diffusion problem (1), one obtains the finite volume scheme


where the diffusive fluxes take the forms


with the shorthand notation and using the unit vectors


The diffusivity in Eq.(18) is defined as


which define piecewise constant functions over the diamond cells and dual to internal and external mesh faces respectively. Finally, it is important to notice that the fluxes (18) are locally conservative, since


In the definition of the fluxes, a reconstruction of the face gradient must be employed. For this purpose, a linear interpolation operator is selected at internal faces


while at boundary faces the simplest choice is , for all . Here, we will use the approximation

in order to recover at boundaries. The linear interpolation makes use of (15) with linear interpolation of the face values


while the boundary face values follow directly from the homogeneous Dirichlet conditions in problem (1). Additionally, in practical implementations it is customary to compute the scalar diffusivity at internal faces by a linear interpolation operator as .

To allow for the treatment of non-orthogonal polyhedral meshes, it is useful to consider the associated isotropic diffusion problem



is isotropic diffusivity tensor associated to the scalar diffusivity

. This allows to reformulate problem (1) as


with naturally verifying the usual assumptions [22]. Similarly, the associated weak formulation is given by


It is possible to derive a finite volume scheme for diffusion problems with tensorial diffusivity by constructing a local discrete gradient [22], in order to obtain at cell face a consistent approximation of the diffusive flux , with usual notation for finite volume schemes. To this purpose, it is beneficial to rewrite the diffusive flux for an internal face using the diffusivity tensor from Eq.(24). Since is symmetric, it follows that


In order to allow for the treatment of non-orthogonal meshes, the following diffusivity tensor decomposition can be applied


with anisotropic (directional) diffusivity tensors


by following the natural directions locally identified from the non-orthogonal polyhedral mesh. Notice also that both diffusivity tensors are symmetric, since and , and that


A similar flux decomposition can be carried out at boundary faces by substituting the unit vector with .

By substituting the tensor decomposition from Eqs.(28)-(29) into the finite volume fluxes (27), one obtains that


which directly corresponds to the terms of the Gauss corrected scheme appearing in Eq.(18). In particular, the first term in Eq.(31), corresponding to the anisotropic diffusivity tensor , is amenable to approximation by a two-point flux scheme, in a manner similar to what is done in the perpendicular bisection method in [31]. This term, when inserted into the finite volume diffusive flux, yields


that generates a directional derivative which can be easily approximated via a two-point flux scheme. On the other hand, the second term in the diffusive flux corresponding to the anisotropic diffusivity tensor must be treated via a reconstruction of the cell gradient. it is important also to notice that, from the tensor decomposition in Eqs.(28)-(29), the two-point flux portion increases its dominance for increasing mesh non-orthogonality, due to the increasing angle between the unit vectors and . This property is beneficial in guaranteeing diagonal dominance of the linear system matrix and thus numerical stability, as will be clear from the rest of the discussion.

4 Discrete weak formulation

Returning to the diffusive fluxes from Gauss corrected scheme (18), using the diffusivity tensor decomposition in Eqs.(28)-(29), the finite volume fluxes can be rewritten in the form


for the internal and external faces, respectively. In these formulae, the transmissivities


are introduced to simplify the notation and denotes a generic discrete gradient operator, still to be defined, that is piecewise constant on the diamond cells for all and . If one defines the diamond cell gradient from the linear interpolation of cell gradients as in Eq.(22), then the fluxes become


where the vector quantities


have been introduced, which are such that generally. Notice also the approximation introduced in the boundary term . It is now possible to derive the weak formulation underlying the finite volume scheme (17). By multiplying Eq.(17) by the test function and summing the result for all , one obtains

that, after discrete integration by parts, produces

from which, due to flux conservativity (21), one obtains that


By substituting into Eq.(37) the fluxes (35) with the face gradient from the linear interpolation of the Gauss scheme (15), it is possible to identify two terms and in the expression




Notice that term defines a symmetric bilinear form


which is a discretization of the term , directly corresponding to the portion of diffusive fluxes that can be ascribed to the anisotropic diffusion tensor . This expresses the flux component that is parallel to the local vector associated to internal faces , or to the vector associated to boundary faces. The term contains instead the vectors and , related to the diffusivity tensor , but it is not yet in a form readily corresponding to a discrete weak formulation. To this purpose, it is convenient to rewrite as

The two summation terms between brackets contained in the last expression correspond to the internal faces and the boundary faces contributions, respectively. They can be interpreted as a discretization of the term for all , which allows to introduce the piecewise constant function that is defined on each cell as


expressing the fact that the test function gradient cannot be separated from the diffusivity tensor , since the latter is a face-based quantity defined from the local mesh non-orthogonality (i.e., from the angle between and unit vectors). In this case, the term can be rewritten as a non-symmetric discrete bilinear form


which is thus associated to the diffusivity tensor , expressing the contribution of diffusion from a local direction not aligned with at internal faces, or with at boundary faces.

Thus, the discrete weak formulation implied when using the Gauss corrected scheme for the heterogeneous isotropic diffusion problem (1) takes the form


Several remarks are in order on the basis of the previously introduced formulation. Similarly to [22], cell gradients can lead to a discrete inner product whenever the mesh geometry allows for a direct estimation of face normal fluxes, e.g., in the case of an orthogonal polyhedral mesh. When instead anisotropic effects (directional bias) emerge locally on cell faces due to mesh non-orthogonality, the construction of face gradients becomes inevitable, as done in [25]. In this latter case, the diffusivity tensor is necessarily defined on diamond cell support.

Secondly, the linear interpolation operator used to obtain the face gradient in the Gauss corrected scheme from a linear combination of cell gradients calculated via the Gauss gradient scheme for all implies that


which defines the face gradient as the diamond cell gradient obtained via an inverse volume weighting procedure. The same conclusion is also valid for the scalar diffusivity defined from Eq.(20).

Finally, when the scalar diffusivity appearing inside the fluxes (33) is computed by a linear interpolation procedure, the anisotropic diffusivity tensor in the Gauss corrected scheme becomes

where the face non-orthogonality symmetric tensor has been defined. As a consequence, the diffusive flux (33) can be recast into the form

from which, after introducing the vectors


one obtains that