A generalized formulation for gradient schemes in unstructured finite volume method

The computation of gradients of flow properties plays a pivotal role in obtaining accurate numerical solutions of Navier-Stokes equations on unstructured meshes. This is because these gradients are necessary to evaluate both inviscid and viscous fluxes in finite volume flow solvers. The solution gradients are employed for higher-order convective schemes for inviscid flux calculations while viscous fluxes are directly dependent on the gradients of velocities and temperature. The construction of robust and accurate flow solvers for incompressible and compressible flows therefore necessitate an accurate and robust methodology for gradient computation among other requirements.

Cell-centered finite volume codes which typically employ unstructured grids reconstruct the centroidal gradients using reconstruction schemes that can be generally classified as belonging to either least squares (LSQ) or Green-Gauss (GG) family of approaches [1, 2]. While the former category of techniques are based on the principle of least squares minimization, the latter class of methods rely on the Gauss divergence theorem to calculate the gradients. These reconstruction strategies have been in widespread use, with LSQ approaches largely favoured for Euler flows while GG-based schemes are more common in viscous computations. The schemes of the two families have their own advantages and drawbacks. While the linear least squares approach gives consistent gradients and accurate estimates on isotropic irregular meshes, it is quite inaccurate on stretched meshes over curved geometries [3]. The use of weights mitigates the problem but the weighted least squares reconstruction requires limiting to stabilise the numerical solution on these meshes even for subsonic inviscid flows [4] and also results in incorrect gradients on triangulated meshes with skewed elements [3]. The choice of the explicit neighbourhood necessary for least squares family of schemes affects both the stability and accuracy and requires special attention. The Green-Gauss family of schemes do not always lead to consistent gradients but can produce accurate estimates of gradients on curved meshes with high aspect ratio [3]. Coupled with their ease of implementation, they are a popular choice in viscous flow solvers where they estimate aerodynamic coefficients quite accurately [4]. The detailed investigations of Diskin and Thomas [5] for high Reynolds number applications also demonstrate that the gradient accuracy on high aspect ratio meshes depends both on the grid and the solution.

It is evident that accuracy, robustness and computational cost are conflicting requirements for any gradient scheme. While there exists no known “all-in-one” scheme with all these attributes for a wide range of flow conditions on arbitrary polyhedral meshes, there have been advancements in devising new improved schemes as well as constructing hybrid reconstruction approaches. The latter class of methods attempt to blend the GG and LSQ strategies with a view to jointly harness their advantages and the GLSQ scheme of Shima and co-workers [4] is a notable effort in this direction. Among the recently proposed schemes for gradient computations are the variational reconstruction (VR) approach of Wang et al. [6], the Implicit Green-Gauss (IGG) [7] scheme, the Taylor-Gauss gradients (TG) [8, 9] and the Modified Green-Gauss (MGG) [10] reconstruction. The first two schemes are inherently implicit and may be regarded effectively as an extension of compact schemes for unstructured meshes, with the VR scheme categorised as a least squares based method and IGG reconstruction belonging to Green-Gauss family. The TG gradients are constructed similar to least squares but have dependence on the face normal vector similar to GG, thereby having a unique flavor of hybridization. The MGG reconstruction, derived from a different variant of the Gauss-divergence theorem, reconstructs centroidal gradients from face normal derivatives of the scalar which makes it an implicit scheme on non-orthogonal meshes. All the schemes however result in consistent gradients on arbitrary mesh topologies and exhibit advantages compared to their conventional counterparts for fluid flow problems. A good overview of various gradient techniques in the finite volume framework including recent developments may be found in [11].

The least squares and Green-Gauss reconstruction schemes have been traditionally considered as two different families for gradient computation. This is arguably because they adopt vastly different mathematical philosophies to compute the same quantity of interest. However, efforts to hybridize the constructional philosophies of both classes have yielded some interesting insights into the methods. For instance, it is possible to derive the Green-Gauss formula from a face-based least squares reconstruction using a specific choice of weights [12]. In another work by Syrakos et al. [9], a unified framework for gradient reconstruction schemes is proposed which encompasses both GG and LSQ classes of schemes. The basic idea is to use Taylor expansion to express the scalar at a neighbour point (which may or may not be a neighbour cell center) in terms of the value and gradients of the scalar at the cell center where the gradients are sought. The difference in values is then multiplied by a vector and the summation across a suitable stencil yields a system of equations that can be solved to obtain the centroidal gradients. The choice of the neighbour point and the vector used to left multiply the system, distinguishes the different gradients schemes, with LSQ class of methods generated by an orthogonal projection and the GG class of methods by an oblique projection.

In this work, we propose a generalized formulation for nominally first order accurate gradient schemes in unstructured finite volume method, in the similar spirit as Syrakos et al. [9]. Our proposed framework reconstructs centroidal gradients from the directional derivatives of the scalar in a specified neighbourhood, which brings in more generality in the formulation. The mathematical framework of the proposed formulation is presented in section 2. Following that in section 3, we demonstrate how existing gradient schemes of both GG and LSQ families can be constructed using the generalized formulation. The framework also enables us to interpret and possible derive hybrid schemes that could retain some benefits of both GG and LSQ family of schemes. This is discussed in section 4. Finally, we extend the proposed framework to construct a generalized version of flexible gradient scheme [13] in section 5. The flexible scheme has a free-parameter that can be tuned to increase robustness of finite volume computations without compromising the gradient reconstruction accuracy. Some general discussions on the novelties and drawbacks of the proposed formalism are presented in section 6.

We present a generic mathematical framework for nominally first order accurate gradient reconstruction schemes in unstructured cell-centered finite volume method. Traditional approaches of deriving gradient schemes are generally through some variants of the Gauss-divergence theorem or through a least squares minimization problem. We, however, adopt a different approach wherein the development of a gradient reconstruction formula is through a tensor identity, that requires the definition of two key entities - the ‘neighbourhood’ of a cell, and, a geometric vector pair ($\mathbf{a},\mathbf{b}$) defined at each element forming the neighbourhood of a cell.

To illustrate, let us consider a typical unstructured mesh with a characteristic dimension ‘$h$’¹¹1For a general unstructured mesh in a $d$-dimensional domain, $h=(V/N)^{1/d}$, where $V$ is the total volume of the computational domain and $N$ is the number of cells discretizing the domain., where the gradients of a scalar are required to be computed in a representative cell, with the cell center “$i$” (see Fig. 1). We define a ‘neighbourhood’ of the cell comprising of a set of points that maintain a compact support about the cell centroid. For the sake of generality, the location of these points are kept arbitrary but typical choices could be (but not restricted to) the cell centres of face-sharing or node-sharing cells, face centers of the faces forming the cell, mid-point of the line joining the face-sharing cell centers, etc. We use the subscript “$f$” to denote a representative point in the defined neighbourhood of the cell “$i$”, as shown in Fig. 1. With the neighbourhood set, we define a pair of “geometric” vectors $\mathbf{a}_{f}$ and $\mathbf{b}_{f}$ at each point forming the cell neighbourhood. The term “geometric” here implies that the vectors shall be function of the mesh metrics alone with a constraint on the magnitudes of these vectors,

where, $p$ and $q$ are positive integers²²2It is presumed that $h\ll 1$ as a requirement for sufficient resolution of the calculations..

For a continuous vector field $\bf{u}$, we can write the following tensor identity,

where ${\bf u}_{f}$ denotes the value of the vector at the point “$f$”. Since $\bf u$ is a continuous vector field, we use a Taylor expansion to express ${\bf u}_{f}$ on the left hand side of Eq. 2 in terms of the value at the cell centroid ${\bf u}_{i}$, which after rearranging leads to the expression,

where, $\epsilon_{1}\sim O(h^{2})$ represents the higher order terms. Summing up the above equation over the entire neighbourhood of the cell, we get the relation,

The sum of the dyadic product of the vectors $\mathbf{a}_{f}$ and $\mathbf{b}_{f}$ appearing in the left hand side of the above equation, is a second order tensor, referred to from here on as the ‘dyadic sum’,

To construct the gradient reconstruction scheme, we place another constraint on the choice of the neighbourhood and the geometric vectors, which is that the dyadic sum, $\mathbb{P}$, is invertible. Now, we simply substitute the vector field as the gradient of a scalar field, i.e. $\bf u=\nabla\phi$ in Eq. 4 and obtain the following expression for centroidal gradients,

In obtaining the above expression, the contribution from the second term on the right hand side of Eq. 4 is ignored. It can be shown that the contribution from that term is $O(h)$ (see appendix A.1) and hence can be ignored to obtain an expression for nominally first order accurate gradients.

The generalized gradient scheme of Eq. 5 reconstructs centroidal gradients from the directional derivatives of $\phi$ along the $\mathbf{b}_{f}$ vectors. It is possible to generalize existing gradient reconstruction schemes using this framework, as shall be demonstrated in the next few sections. With the neighbourhood and geometric vectors defined, the next crucial step in the formalism is the evaluation of quantity $\mathbf{b}_{f}\cdot\nabla\phi_{f}$ on the right hand side of Eq. 5, from which the centroidal gradients are to be reconstructed. For a cell-centered finite volume discretization, we have the values of $\phi$ only at the cell centers. Therefore, we can express the quantity $(\mathbf{b}_{f}\cdot{\nabla\phi}_{f})$, in terms of centroidal values of $\phi$. To do that, we write a general expression for the directional derivative using Taylor series expansion,

where, the $\mathbf{x}_{f,1}=\mathbf{x}_{f}+\delta_{f,1}\mathbf{b}_{f}$ and $\mathbf{x}_{f,2}=\mathbf{x}_{f}+\delta_{f,2}\mathbf{b}_{f}$ are points along $\mathbf{b}_{f}$, at a distance of $\delta_{f,1}/|\mathbf{b}_{f}|$ and $\delta_{f,2}/|\mathbf{b}_{f}|$ from $\mathbf{x}_{f}$, respectively (see Fig. 1). The first term on the right hand side of Eq. 6 essentially expresses $(\mathbf{b}_{f}\cdot{\nabla\phi}_{f})$ in terms of $\phi$ at specific locations along the $\mathbf{b}_{f}$ vector. While the choice of the locations $\mathbf{x}_{f,1}$ and $\mathbf{x}_{f,2}$ are generic, to ensure compactness, we impose a restriction that $|\delta_{f,1}\mathbf{b}_{f}|\sim O(h)$ and $|\delta_{f,2}\mathbf{b}_{f}|\sim O(h)$, i.e., the points $\mathbf{x}_{f,1}$ and $\mathbf{x}_{f,2}$ are at an $O(h)$ distance from $\mathbf{x}_{f}$. Therefore, substituting Eq. 6 into Eq. 5, we obtain a nominally first order accurate gradient reconstruction formula,

where, $\phi_{f,1}=\phi(\mathbf{x}_{f,1})$ and $\phi_{f,2}=\phi(\mathbf{x}_{f,2})$. The contribution from the second term on the right hand side of Eq. 6 is ignored as it can be shown to be $O(h)$ in magnitude (see appendix A.2).

The expression in Eq. 7 represents an explicit gradient reconstruction formula if the values of $\phi$ are available at $\mathbf{x}_{f,1}$ and $\mathbf{x}_{f,2}$. In general, $\mathbf{b}_{f}$ need not pass through cell centers, implying that there need not exist a point along the vector $\mathbf{b}_{f}$ where the values of $\phi$ are available. Therefore, in order to obtain these values, one may need to interpolate. If that is the case, then it can be shown (see appendix A.3) that the interpolation needs to be at least second order accurate in order to obtain first order accurate centroidal gradients. Instead of obtaining $\phi_{f,1}$ and $\phi_{f,2}$, using the values from neighbour cells, one may also use the values and the gradients from neighbour cells in the interpolation scheme, thereby creating an implicit gradient scheme. Implicit schemes can also be constructed by directly expressing the directional derivative $(\mathbf{b}_{f}\cdot\nabla\phi_{f})$ in terms of the values and gradients in the neighbourhood of the cell.

In this section, we demonstrate that existing gradient schemes including the popularly employed Green-Gauss reconstruction [14, 15] and the linear least squares reconstruction [3, 16] can be generalized under the proposed formalism. In all the following discussions, we assume a two-dimensional domain discretized by non-overlapping arbitrary polygonal volumes where the gradient of the scalar field needs to be reconstructed at the cell-centers. The extension of ideas presented herein to three dimensions is relatively straightforward.

Besides the GG and LSQ schemes discussed in the previous section, there are some schemes which cannot be clearly categorised into either of these families but have essence of both. One such scheme is the recently proposed Taylor-Gauss (TG) gradients [8, 9], which is obtained through an oblique projection and hence is not a minimization problem but leads to a least squares like system of equations. In our proposed framework, the same scheme can be obtained by first choosing the neighbourhood as the points of intersection of each face with the lines joining the cell center “$i$” to the cell center of the respective face-sharing neighbour cell. The geometric vectors are then chosen as, $\mathbf{a}_{f}=\mathbf{S}_{f}$, and $\mathbf{b}_{f}=\mathbf{x}_{j}-\mathbf{x}_{i}$, where $\mathbf{S}_{f}$ is the face area vector. The resulting dyadic sum is,

With the directional derivative evaluated the same way as in the linear least squares scheme, we get,

and the corresponding system of equations to obtain the TG gradients as,

The weighted variants of Taylor-Gauss scheme [9] can be easily obtained by scaling the geometric vectors similar to the weighted variants of LSQ.

It is remarked upon by the authors of [8] that the TG gradient schemes have a flavor of GG due to the use of face normals in reconstructing the centroidal gradients but with the additional ability to create weighted variants like an LSQ scheme. If we notice closely, to construct the TG scheme through the generalized formulation, we effectively borrowed $\mathbf{a}_{f}$ and $\mathbf{b}_{f}$ from GG and LSQ schemes, respectively. Therefore, the scheme is, in some sense, a hybrid of GG and LSQ class of schemes, despite not being derivable from either routes clearly. This has been explained in [9] through their unification philosophy based on orthogonal and oblique projections. Within the context of our proposed formalism, the construction of TG gradients show a logical way of devising hybrid schemes out of traditional schemes, which is through different permutations of their respective geometric vectors. The efficacy of such a hybrid scheme can only be determined through numerical studies which is beyond the scope of this paper. Nevertheless, the proposed generalization provides an alternate route to creating nominally first order accurate hybrid gradient schemes where the advantages of the constituent schemes can possible be gained.

Gradients of scalar fields are required during finite volume computation of fluid flow problems, where inconsistent gradients could lead to incorrect solutions. However, certain gradients reconstruction schemes like the linear least squares, despite being consistent, give rise to numerical instabilities on highly distorted meshes [3, 4]. Therefore, the efficacy of a gradient reconstruction scheme cannot be judged only in terms of its accuracy. For least squares reconstruction, the choice of the weights provide an indirect means to introduce robustness into the scheme. But a better idea is to directly introduce a damping parameter into the scheme, which was proposed in a flexible gradient scheme by Nishikawa [13]. The key idea in [13] is first to obtain gradients using conventional least squares approach and then use these gradients to compute the “$\alpha$-damped” gradients, which are then averaged to obtain the centroidal gradients. The parameter ‘$\alpha$’ can be controlled to directly introduce damping without compromising the accuracy of the computed gradients. Additionally, it is also shown that for a specific value of ‘$\alpha$’, the gradients become fourth order accurate on uniform Cartesian meshes. In this section, we propose a generalization of the flexible gradient idea using the generalized gradient reconstruction formalism of section 2. To motivate this, we first demonstrate the construction of a flexible variant of the Modified Green-Gauss (MGG) scheme [10].

In this work, we have presented a generalization of gradient reconstruction schemes in cell-centered unstructured finite volume method. The key components of the proposed framework are a pair of geometric vectors defined in a suitably chosen neighbourhood of a cell, different choices of which yield different gradient schemes. Commonly employed gradient reconstruction schemes generally belong to either the Green-Gauss family or the least squares family, both of which are built from starkly different philosophies. The proposed formalism shows that they can be mathematically unified under a common framework. But as these schemes can be conveniently derived from their own mathematical principles, it might appear that the proposed generalization merely serves as a post facto mathematical ‘fitting’ tool. The authors, however, believe that it is not just the intended purpose of this framework.

The generalization provides some important insights into the basic philosophy of gradient reconstruction. It shows that every nominally first order accurate gradient scheme effectively reconstructs centroidal gradients from the directional derivative of the scalar along specific directions in a defined neighbourhood. This viewpoint is not trivial and, to the best of the authors’ knowledge, has not been perceived before, as it is not seen from either the divergence theorem or the minimization based route. The only known generalization of gradient schemes is based on the idea of orthogonal and oblique projections proposed by Syrakos et al. [9]. Their framework has a similar flavor to the one proposed here as it also has a dyadic sum that requires inversion. But there is an important distinction between the two formalisms that is worth understanding. The framework proposed by Syrakos et al. [9] begins with obtaining a first order accurate estimate of the change in the scalar along a particular direction, using the value and gradients at the centroid of a cell “$i$”, i.e.,

Here, “$f$” denotes a point in a compact neighbourhood about “$i$”. The above expression is now multiplied by another vector $\mathbf{V}_{f}$, and then summed across the neighbourhood to obtain the generalized gradient formula,

Comparing this to Eq. 5, it can be easily seen that the vector $\mathbf{V}_{f}$ is equivalent to $\mathbf{a}_{f}$. Also, $\Delta\phi_{f}$ on the left hand side of Eq. 20 can be seen as an evaluation of the quantity $\mathbf{R}_{f}\cdot\nabla\phi_{i}$ (from Eq. 19), which is effectively the directional derivative of $\phi$ along the vector $\mathbf{R}_{f}$. Therefore, we can draw an equivalence of $\mathbf{R}_{f}$ to $\mathbf{b}_{f}$, but that is not strictly true. The way in which the formulation proposed by Syrakos et al. [9] is built, $\mathbf{R}_{f}$ always points from the cell center “$i$” to some point “$f$” in the defined neighbourhood. The location of “$f$” is kept arbitrary, which provides a degree of freedom in the choice of $\mathbf{R}_{f}$, but the vector always passes through the cell center “$i$”. This restriction is not imposed on the vector $\mathbf{b}_{f}$ in our framework. Physically it means that at conception, the framework proposed by Syrakos et al. [9] reconstructs centroidal gradients from the variation of the scalar between the cell center and an arbitrary point in the neighbourhood, whereas we propose that the same variation could be taken between any two arbitrary points in a compact neighbourhood. Of course, it is always possible to express an arbitrary $\mathbf{b}_{f}$ vector as a sum of $\mathbf{R}_{f}$ and another vector to obtain some form of a generalization like [9], but at a philosophical level, our proposition has a different interpretation. It says that instead of the $\phi$ difference between the cell center and an arbitrary ‘point’, it is, in fact, the change in $\phi$ along an arbitrary ‘direction’ that is required to reconstruct nominally first order accurate centroidal gradients.

The implication of this insight is apparent from the re-derivation of some of the existing gradient schemes. From section 3, one can notice that for the traditional Green-Gauss and least squares reconstruction, the vector $\mathbf{b}_{f}$ is either $\mathbf{x}_{f}-\mathbf{x}_{i}$ (for GG) or $\mathbf{x}_{j}-\mathbf{x}_{i}$ (for LSQ), both of which pass through “$i$”, and therefore can be comfortably derived using the framework of Syrakos et al. [9] as well as ours. But for schemes like the MGG or FLSQ, where $\mathbf{b}_{f}$ is along the face normal (which for an arbitrary unstructured mesh does not pass through cell centers), it not clear how the framework of Syrakos et al. [9] can be used to derive them. The arbitrariness of $\mathbf{b}_{f}$ also allows for construction of hybrid schemes by just using $\mathbf{a}_{f}$ and $\mathbf{b}_{f}$ vectors directly from any two known schemes as shown in section 4. In addition, the perspective of reconstructing centroidal gradients from directional derivatives allows us to create flexible variants of known schemes as shown in section 5. The idea is to introduce an additional damping component in the evaluation of the directional derivative that is dependent on the centroidal gradients. These centroidal gradients could be evaluated apriori from a base scheme that is possibly unstable and then substituted into the $\alpha$-damped version of the directional derivatives. With a suitably chosen value of the parameter $\alpha$, the reconstructed gradients can introduce stability to finite volume computations without compromising the gradient accuracy. From a user perspective, this provides a simple route to increasing the robustness of their solver by modifying their choice of gradient reconstruction scheme into the flexible variant.

While the novelties of the proposed framework have been explained, it is incumbent upon us to highlight some of its shortcomings. The mathematical generalization proposed here does not outline a clear directive as to how the geometric vectors are to be chosen to construct a new gradient reconstruction scheme. The only constraints on the geometric vectors are,

1. 1.

   The geometric vectors should depend on the mesh metrics,

2. 2.

   Their magnitudes should be $O(h^{k})$, where $k>0$, and,

3. 3.

   The sum of their dyadic product over the chosen neighbourhood should be an invertible tensor.

These guidelines are reasonably arbitrary from the point of view of directing possible choices of $\mathbf{a}_{f}$ and $\mathbf{b}_{f}$ to construct a scheme. The generalization, however, does indicate that one need not necessarily follow the traditional routes, like the divergence theorem or a minimization idea, to derive a first order accurate gradient reconstruction scheme. A good example is the Taylor-Gauss (TG) gradient scheme [8, 9] which is derivable neither from the GG or the LSQ route, but can be easily obtained from the unified framework of Syrakos et al. [9] as well as ours. Despite not being derivable from either routes, the scheme has flavor of both GG and LSQ families [8, 9], simply because the geometric vectors that define the scheme are borrowed from GG or LSQ. The generalization, therefore, allows us to see the mathematical hybridization that exists in the TG scheme, and could possibly pave a way for more such hybrid schemes in the future. Nevertheless, the framework does not specify how such a hybridization could be ‘better’ or ‘worse’ that a conventional scheme. As like in most cases, these can only be determined by numerical tests. Similarly in the extension to flexible gradient schemes, the generalization does not specify how to choose the flexible parameter to increase robustness of a given scheme. It is not generally possible to do so on arbitrary unstructured meshes as the degree of damping is highly dependent on the scheme as well as the solution. This is demonstrated in a future work by the same authors, where the proposed generalization of $\alpha$-damped scheme is used to construct flexible variants of some known schemes.

Acknowledgments

The authors would like to express their gratitude to Thibault Pringuey for his write-up on CFD Online detailing the ideas behind fvc::reconstruct in OpenFOAM which motivated the present study. The last author would like to acknowledge the support from Science and Engineering Research Board (SERB), Government of India through the MATRICS grant MTR/2019/000241 during the course of this work.