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 multipoint 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 nonorthogonal meshes practicaly required for industrial applications [7], [40]. For these applications, the socalled 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 nonorthogonal 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 preprocessing 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 cellcentered 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 nonorthogonality condition on the mesh is violated. It is to be remarked that existing convergence proofs for finite volume methods on nonorthogonal meshes either involve discretization schemes not guaranteeing local flux conservativity [24], [25]
, or DDFV schemes employing additional degrees of freedom
[4], [5], or twodimensional diamond schemes on meshes satisfying more restrictive regularity conditions [13]. We will focus here on the isotropic steady state diffusion equation(1a)  
(1b) 
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
(2) 
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 threedimensional 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 meshbased 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:

is a finite family of nonempty, 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 .

is a finite family of disjoint subsets of representing the faces. Let be the set of interior faces such that, for all ,
is a nonempty open subset of a hyperplane in
with , and let be the set of boundary faces such that, for all , is a nonempty 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 . 
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 starshaped, in the sense that if , then the line segment .

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(3) 
Additionally denotes the orthogonal distance between and
(4) 
which is constant for all From the assumption that is starshaped, it follows that and that it also holds:
(5) 
For all and , denotes the cone with vertex and basis , also called halfdiamond, that is the volume defined by
(6) 
For all , denotes the diamond associated to face , as in Figure 1.
Definition 2.1 covers a wide range of meshes, including meshes with nonconvex cells, with nonplanar 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.
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(7) 
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
(8) 
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 nonorthogonal 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
(9) 
where the numerical flux through face is computed as
(10) 
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 nonhomogeneous [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 twopoint flux approximation [21]
(11) 
Even though unconditionally monotone and coercive [15], it is of limited accuracy on unstructured meshes, where mesh nonorthogonality may lead to severe errors in the approximation of the diffusion fluxes [20]. In order to compensate for the unavoidable nonorthogonality of realistic unstructured meshes, a simple but effective solution is provided by the Gauss corrected scheme, which introduces a nonorthogonal correction term [36] in the twopoint flux scheme, thus obtaining the approximation
(12)  
Here, the first term corresponds to the twopoint 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 nonorthgonality 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 twopoint 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 twopoint 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
(13) 
for Since for any closed control volume the geometrical relations
(14) 
hold, with , Eq.(13) is also equal to
(15) 
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
(16) 
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 nonorthogonal 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
(17) 
where the diffusive fluxes take the forms
(18a)  
(18b) 
with the shorthand notation and using the unit vectors
(19a)  
(19b) 
The diffusivity in Eq.(18) is defined as
(20a)  
(20b) 
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
(21) 
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
(22) 
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
(23) 
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 nonorthogonal polyhedral meshes, it is useful to consider the associated isotropic diffusion problem
(24) 
where
is isotropic diffusivity tensor associated to the scalar diffusivity
. This allows to reformulate problem (1) as(25a)  
(25b) 
with naturally verifying the usual assumptions [22]. Similarly, the associated weak formulation is given by
(26)  
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
(27) 
In order to allow for the treatment of nonorthogonal meshes, the following diffusivity tensor decomposition can be applied
(28) 
with anisotropic (directional) diffusivity tensors
(29a)  
(29b) 
by following the natural directions locally identified from the nonorthogonal polyhedral mesh. Notice also that both diffusivity tensors are symmetric, since and , and that
(30) 
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
(31)  
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 twopoint 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
(32)  
that generates a directional derivative which can be easily approximated via a twopoint 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 twopoint flux portion increases its dominance for increasing mesh nonorthogonality, 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
(33a)  
(33b) 
for the internal and external faces, respectively. In these formulae, the transmissivities
(34) 
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
(35a)  
(35b) 
where the vector quantities
(36a)  
(36b) 
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
(37)  
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
(38) 
where
(39a)  
(39b)  
(39c) 
Notice that term defines a symmetric bilinear form
(40)  
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
(41)  
expressing the fact that the test function gradient cannot be separated from the diffusivity tensor , since the latter is a facebased quantity defined from the local mesh nonorthogonality (i.e., from the angle between and unit vectors). In this case, the term can be rewritten as a nonsymmetric discrete bilinear form
(42) 
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
(43)  
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 nonorthogonality, 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
(44) 
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 nonorthogonality symmetric tensor has been defined. As a consequence, the diffusive flux (33) can be recast into the form
from which, after introducing the vectors
(45a)  
(45b) 
one obtains that
Comments
There are no comments yet.