Skip to main content

Extending the global-direction stencil with “face-area-weighted centroid” to unstructured finite volume discretization from integral form

Abstract

Accuracy of unstructured finite volume discretization is greatly influenced by the gradient reconstruction. For the commonly used k-exact reconstruction method, the cell centroid is always chosen as the reference point to formulate the reconstructed function. But in some practical problems, such as the boundary layer, cells in this area are always set with high aspect ratio to improve the local field resolution, and if geometric centroid is still utilized for the spatial discretization, the severe grid skewness cannot be avoided, which is adverse to the numerical performance of unstructured finite volume solver. In previous work [Kong, et al. Chin Phys B 29(10):100203, 2020], we explored a novel global-direction stencil and combined it with the face-area-weighted centroid on unstructured finite volume methods from differential form to realize the skewness reduction and a better reflection of flow anisotropy. Greatly inspired by the differential form, in this research, we demonstrate that it is also feasible to extend this novel method to the unstructured finite volume discretization from integral form on both second and third-order finite volume solver. Numerical examples governed by linear convective, Euler and Laplacian equations are utilized to examine the correctness as well as effectiveness of this extension. Compared with traditional vertex-neighbor and face-neighbor stencils based on the geometric centroid, the grid skewness is almost eliminated and computational accuracy as well as convergence rate is greatly improved by the global-direction stencil with face-area-weighted centroid. As a result, on unstructured finite volume discretization from integral form, the method also has superiorities on both computational accuracy and convergence rate.

1 Introduction

Compared with the block-structured grid, unstructured grid is highly automated in grid generation [1,2,3] and convenient for the grid adaptation [4, 5], while the computational accuracy and stabilities are hard to be guaranteed [6], especially on high-aspect-ratio triangular grids [7,8,9]. Regarding the current difficulties, scholars are continued trying to improve the discretization algorithms to break through the bottleneck of accuracy loss and stability deterioration on highly anisotropic grids, and realize the unification of automated grid generation and accurate numerical simulation [10, 11].

As for the research about unstructured finite volume discretization on high-aspect-ratio triangular grids, Diskin and Thomas [12, 13] tested the accuracy of gradient reconstruction, and results show that the accuracy of reconstructed gradient on such grid is determined by both solution and grid, and the existing discretization schemes cannot meet the requirements of grid quality and computational accuracy simultaneously. On this basis, Diskin and Thomas et al. [14, 15], systematically compared the numerical performance of inviscid and viscous fluxes on different node-centered and cell-centered unstructured finite volume methods. Research suggests that for high-aspect-ratio triangular grids, there are few discretization schemes able to preserve the simulation accuracy, and after adding disturbance to grid points, the magnitude of computational errors obtained by all schemes is proportional to cell aspect ratios. Similar conclusions are obtained in Ref. [16].

On this basis, we notice that although there are numerous discretization approaches for high-aspect-ratio triangular grids, studies related to the stencil selection are quite few, and within the framework of unstructured finite volume method, both inviscid and viscous fluxes are greatly influenced by the gradient reconstruction [17], where different stencils play a crucial role. Commonly used stencils are face-neighbor and vertex-neighbor stencils that consist of cells sharing faces or vertices with the central cell. Apart from these two commonly used stencils, an ingenious stencil augmentation method is proposed in Ref. [18], which is named as the Smart Augmentation (SA) stencil. For each vertex, an extra cell closest to geometric centroid is selected from vertex-adjacent cells and appended to the face-neighbor stencil. This method improves the stabilities of face-neighbor stencil to a certain extent. In order to further improve stabilities of the face-neighbor stencil, efficient gradient stencils are proposed by Nishikawa [19], where the F-decreasing augmentation in combination with the symmetric augmentation stencil augments a face-neighbor stencil with extra cells to attempt to increase the symmetry of the stencil. On this basis, extra cells are added to decrease the reciprocal of the Frobenius norm of a scaled least-squares matrix to minimize the lower bound of the magnitude of the gradient. This stencil augmentation approach is much more robust as well as efficient, and the instability of face-neighbor stencil is effectively avoided even on extremely irregular grids.

Apart from the mentioned stencil augmentation methods, Xiong et al. [20], proposed a local-direction stencil selection method, where stencil cells are augmented along two local directions. In this method, characteristics of the flowfield are taken into account during the process of determining local directions. On isotropic triangles, two local directions are close to the normal and tangential directions of the wall, while on high-aspect-ratio triangular grids, stencil cells selected by this method will deviate a lot from the boundary normal, and the numerical performance on such grids will get deteriorated [21, 22].

Compared with the local-direction stencil selection method, global direction stencil selection method [22, 23] effectively overcomes the deviation on grids with high aspect ratios, and cells are always selected along two global directions, that is, normal and tangential directions of the wall. Hence, the flow anisotropy could be well reflected, and after verification, on high-aspect-ratio triangular grids, global-direction stencil has a better numerical performance on both computational accuracy and efficiency than commonly used face-neighbor and vertex-neighbor stencils as well as the local-direction stencil.

Although global-direction stencil preliminary exhibits a better numerical performance, the distribution of reference points is still much more skewed, especially on high-aspect-ratio triangular grids. You and Mittal et al. [24], first proposed the grid skewness on such grids, and concluded that grid skewness is adverse to the computational accuracy and stabilities of CFD solvers. Regarding the grid skewness, there are various definitions of that, such as the angle between face normal and the vector pointing form face centroid to the cell centroid, the minimal internal angle of grid cell [25], ratio of the max diagonal to the minimum [26], etc. Nevertheless, from different definitions, the same conclusion could be drawn that on high-aspect-ratio triangular grids, the grid skewness is always quite evident [27].

Different from above skewness measures, Nishikawa [28] proposed a novel definition of that, and it is defined at a face shared by two neighbor cells, say, A and B, as the dot product of the unit vector pointing from the centroid of cell A to that of cell B. Therefore, a non-skewed grid has the measure one and highly-skewed grid nearly zero. Specifically, on highly-skewed triangular grids, the reference point distribution is irregular and exhibits deflective phenomenon. Thus, although global-direction stencil cells mentioned above are along the normal and tangential directions of the wall, the irregular distribution of reference points is not changed, and it is hard to say whether flowfield characteristics are well reflected or not. In order to reduce the grid skewness and optimize the reference point distribution, Nishikawa proposed a novel face-area-weighted centroid [28] for the second-order unstructured finite volume discretization from differential form. After verification, the second-order accuracy is also achieved, and the convergence rate is greatly improved. Besides, compared with the traditional geometric centroid, the distribution of this novel reference point is more regular, and the connection of that is almost parallel to the boundary normal vector.

Based on this phenomenon, in previous work, we combined this novel centroid and global-direction stencil for the second-order unstructured finite volume method [23], and it is verified that the global-direction stencil with face-area-weighted centroid has a lower discretization error than stencils with the geometric centroid including face-neighbor and vertex-neighbor stencils as well as the global-direction stencil. What’s more, this novel method has a more stable numerical performance on high-mach-number flow. As a result, the current situation related to accuracy loss and stability deterioration on high-aspect-ratio triangular grids is greatly ameliorated by the employment of this improved global-direction stencil. This method is designed for the second-order differential finite volume solver [28], where both solution and source term vectors are evaluated as point values, and therefore, the source term integration is totally avoided. Besides, higher-order accuracy of the differential finite volume solver is also able to be achieved, and professor Nishikawa gave an idea in Ref. [28, 29] about the primitive function reconstruction of the flux, by which the governing equation is higher-order accurate without the second-order errors.

Greatly inspired by the face-area-weighted centroid and efficient differential finite volume solver, we explore the novel local origin in the finite-volume discretization based on the integral form and investigate its impact not only on a second-order method but also on a third-order method, since on the integral form, the grid skewness is also existing and needs to be eliminated to guarantee the accuracy and convergence rate.

In this research, the second-order accurate integral finite volume discretization is considered at first, where the k-exact reconstruction [30,31,32,33,34,35,36] based on any local origins is derived, and we combine the global-direction with the face-area-weighted centroid to further reduce the grid skewness and capture the variation of flowfield more accurately. In addition, four representative numerical examples are designed to examine the effectiveness of this extension, and on this basis, the third-order k-exact reconstruction and finite volume discretization based on any local origins are also derived to further extend the global-direction stencil with face-area-weighted centroid to higher-order accurate unstructured finite volume method, and the accuracy as well as the discretization errors is examined by a numerical example with the method of manufactured solutions.

The paper is organized as follows. In Section 2, governing equations and spatial discretization from both differential and integral forms are given at first. And then, we theoretically derive the second-order k-exact reconstruction process based on any local origins. Different stencil selection methods and global-direction stencil with face-area-weighted centroid are presented in Section 3. In Section 4, numerical examples governed by linear convective, Euler and Laplacian equations are employed to verify the effectiveness, feasibilities as well as superiorities of the method on the second-order integral unstructured finite volume solver. The third-order k-exact reconstruction and finite volume discretization based on any local origins, as well as the accuracy analysis are given in Section 5. Concluding remarks and future work will be summarized in Section 6.

2 Unstructured finite volume method

As mentioned in introduction, the finite volume method from differential form and the discretization base on the face-area-weighted centroid [28] give us great inspirations, and in this paper, to further extend the applicable scope of this novel local origin, we promote it on unstructured finite volume method from integral form, since when unstructured finite volume discretization from integral form is employed, the grid skewness is existing and needs to be eliminated as well.

As a result, in this section, we give unstructured finite volume discretization from integral form at first. In addition, the k-exact reconstruction method based on any local origins is also derived in this section.

2.1 Finite volume discretization from integral form

The finite volume discretization from integral form is utilized in this paper, where the governing equation in integral form could be written as

$$ \underset{V_j}{\int }{\partial}_t{\boldsymbol{u}}_j\kern0.5em dV+\underset{V_j}{\int}\mathit{\operatorname{div}}{\boldsymbol{F}}_j^c\kern0.5em dV=\underset{V_j}{\int }{\boldsymbol{s}}_j\kern0.5em dV, $$
(1)

where, uj is a solution vector. \( {\boldsymbol{F}}_j^c \) is the convective flux and sj is a source (or forcing) term vector. Besides, Vj and ∂Vj are the volume and boundaries of cell j. According to the divergence theorem, Eq. (1) is transformed as

$$ \underset{V_j}{\int }{\partial}_t{\boldsymbol{u}}_j\kern0.5em dV+\underset{\partial {V}_j}{\oint }{\boldsymbol{f}}_n^c\kern0.5em dS=\underset{V_j}{\int }{\boldsymbol{s}}_j\kern0.5em dV, $$
(2)

where, \( {\boldsymbol{f}}_n^c \) is the convective flux along the face normal vector. On this basis, Eq. (2) could be discretized as follows.

$$ {\partial}_t{\overline{\boldsymbol{u}}}_j+\frac{1}{V_j}\sum \limits_{k\in \left\{{k}_j\right\}}{\boldsymbol{\varPhi}}_{jk}{A}_{jk}=\frac{1}{V_j}\underset{V_j}{\int }{\boldsymbol{s}}_j\kern0.5em dV, $$
(3)

where, \( {\overline{\boldsymbol{u}}}_j \) is the cell-averaged solution vector and {kj} are faces of cell j. Besides, Φjk is the numerical flux along the normal vector at the k-th face, and Ajk is the k-th face area or length (in 2D). The diagram of unstructured finite volume discretization is shown as follows.

As shown in Fig. 1, where uL and uR are two state vectors at the gauss point obtained by owner and neighbor cells of the common face, and Cj, Ck are two reference points of cell j and cell k respectively. In fact, for traditional unstructured finite volume discretization, these two reference points are located at the geometric centroid of the control volume, while in this paper, the discretization based on any local origins are derived and it will be specifically described in the following section.

Fig. 1
figure 1

Diagram of unstructured finite volume discretization

For the second-order unstructured finite volume solver, the gauss point is just located at the face centroid and the weight is equal to 1, while for higher-order discretization, there should be at least two points to guarantee the flux quadrature is higher-order accurate. The gauss point location as well as the corresponding weights will be specifically given in Section 5.2. Besides, rj and rk are two vectors pointing from the cell centroid of owner and neighbor stencil cells to gauss point respectively. Two state vectors could be interpolated by the point value \( {\boldsymbol{u}}_j^{ori} \) and \( {\boldsymbol{u}}_k^{ori} \) at the local origin obtained by the k-exact reconstruction algorithm [30,31,32,33,34,35,36].

Usually, the reference point is chosen at the geometric centroid, while in this paper, the point-valued solution at any local origins is derived, and the main process will be given in the following section. As a result, uL and uR these two state vectors could be formulated by

$$ \left\{\begin{array}{l}{\boldsymbol{u}}_L={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_G\hbox{-} {x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_G\hbox{-} {y}_j^{ori}\right),\\ {}{\boldsymbol{u}}_R={\boldsymbol{u}}_k^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_k\left({x}_G\hbox{-} {x}_k^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_k\left({y}_G\hbox{-} {y}_k^{ori}\right),\end{array}\right. $$
(4)

where point value solution of \( {\boldsymbol{u}}_j^{ori} \) and \( {\boldsymbol{u}}_k^{ori} \), as well as the solution gradients are all derived and obtained by the k-exact reconstruction algorithm at the local origins Cj and Ck, which will be specifically described in Section 2.2, and \( \left({x}_j^{ori},\kern0.5em {y}_j^{ori}\right) \), \( \left({x}_k^{ori},\kern0.5em {y}_k^{ori}\right) \) are local origins coordinates of cell j and face-adjacent cell k in Cartesian-coordinate system. Besides, (xG,  yG) is the gauss point coordinate at per face. On this basis, the numerical flux along the face normal vector could be formulated as

$$ {\boldsymbol{\varPhi}}_{jk}\left({\boldsymbol{u}}_L,\kern0.5em {\boldsymbol{u}}_R\right)=\left(\begin{array}{c}{\rho u}_n\\ {}{\rho u}_n\boldsymbol{v}+P{\boldsymbol{n}}_{jk}\\ {}{\rho u}_nH\end{array}\right), $$
(5)

where un = vnjk is computed by the Roe flux [37], and H = (γp/(γ – 1) + ρv2/2)/ρ is the specific total enthalpy.

The above content is the second-order accurate unstructured finite volume discretization from integral form, which is employed in this paper and in Section 5, we will specifically exhibit the higher-order, i.e., the third-order accurate discretization process.

2.2 k-exact reconstruction based on any local origins

Least Square (LSQR) gradient reconstruction [38,39,40] is commonly used in the second-order unstructured finite volume methods [13, 14, 18, 41].

As shown in Fig. 2, when calculating the gradient of cell Ci, three stencil cells C1, C2 and C3 are employed to compose the reconstructed equations.

$$ \left(\begin{array}{cc}{\omega}_{i1}\left({x}_1-{x}_i\right)& {\omega}_{i1}\left({y}_1-{y}_i\right)\\ {}{\omega}_{i2}\left({x}_2-{x}_i\right)& {\omega}_{i2}\left({y}_2\hbox{-} {y}_i\right)\\ {}{\omega}_{i3}\left({x}_3-{x}_i\right)& {\omega}_{i3}\left({y}_3\hbox{-} {y}_i\right)\end{array}\right)\left(\begin{array}{l}{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_i\\ {}{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_i\end{array}\right)=\left(\begin{array}{l}{\boldsymbol{u}}_1-{u}_i\\ {}{\boldsymbol{u}}_2-{u}_i\\ {}{\boldsymbol{u}}_3-{u}_i\end{array}\right), $$
(6)

where ωi1, ωi2 and ωi3 are weights set to emphasize geometrically nearby data.

$$ {\omega}_{ij}=\frac{1}{\left|{\boldsymbol{x}}_i-{\boldsymbol{x}}_j\right|}, $$
(7)

and, xi and xj are two position vectors from the origin to geometric centroid of cell i and its stencil cell j respectively.

Fig. 2
figure 2

Diagram of reconstructed triangular stencil cells

As demonstrated in Eq. (6), solution utilized for the least-square reconstruction is the point value at the geometric centroid, and for second-order finite volume methods, this point value is always written in terms of the cell average: \( {\boldsymbol{u}}_j\to {\overline{\boldsymbol{u}}}_j \), but if the reference point is moved anywhere, the point value used for the least-square reconstruction could not be given by the cell average again, and therefore, the reconstruction process could not be implemented, since the point value \( {\boldsymbol{u}}_j^{ori} \) is unknown at any local origins in every time step. Besides, least-square reconstruction is second-order accurate and considering the method extension to higher-order accuracy, the k-exact reconstruction method is employed in this manuscript.

Besides, compared with LSQR, the point value at the reference point could be obtained by the k-exact reconstruction, and this point value is utilized for the flux construction. As a result, although this research is focused on the second-order unstructured finite volume discretization, we also employ the k-exact reconstruction method to solve the solution gradient and obtain the point value at the corresponding reference point. The specific process is shown as follows.

In order to obtain the point value at any local origins, the expansion should be formulated based on any local origins as well,

$$ {R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left(x-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left(y-{y}_j^{ori}\right), $$
(8)

where \( {\boldsymbol{u}}_j^{ori} \) is the point valued solution recovered at a local origin from cell averaged solutions as we will describe later. Integrating Eq. (8) over the control volume, we obtain

$$ \frac{1}{V_j}\underset{V_j}{\int }{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV={\boldsymbol{u}}_j^{ori}+\operatorname{}\frac{\partial \boldsymbol{u}}{\partial x}{\left|{}_j\underset{V_j}{\int}\left(x-{x}_j^{ori}\right)\kern0.5em dV+\operatorname{}\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\underset{V_j}{\int}\left(y-{y}_j^{ori}\right)\kern0.5em dV, $$
(9)

where, \( {x}_j^{ori} \) and \( {y}_j^{ori} \) are the coordinate of any local origins. Replacing \( \left(x-{x}_j^{ori}\right) \) and \( \left(y-{y}_j^{ori}\right) \) with \( \left(x-{x}_j^{gc}\right)+\left({x}_j^{gc}-{x}_j^{ori}\right) \) and \( \left(y-{y}_j^{gc}\right)+\left({y}_j^{gc}-{y}_j^{ori}\right) \) respectively, Eq. (9) could be written as

$$ {\displaystyle \begin{array}{l}\frac{1}{V_j}\underset{V_j}{\int }{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\underset{V_j}{\int}\left(x-{x}_j^{gc}\right)\kern0.5em dV+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\underset{V_j}{\int}\left(y-{y}_j^{gc}\right)\kern0.5em dV\\ {}{\left.+\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\underset{V_j}{\int}\left({x}_j^{gc}-{x}_j^{ori}\right)\kern0.5em dV+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\underset{V_j}{\int}\left({y}_j^{gc}-{y}_j^{ori}\right)\kern0.5em dV,\end{array}} $$
(10)

which could be further transformed as

$$ \frac{1}{V_j}\underset{V_j}{\int }{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV\equiv \kern0.5em {\overline{\boldsymbol{u}}}_j={\boldsymbol{u}}_j^{ori}+\operatorname{}\frac{\partial \boldsymbol{u}}{\partial x}{\left|{}_j{\overline{x}}_j+\operatorname{}\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\overline{y}}_j+\operatorname{}\frac{\partial \boldsymbol{u}}{\partial x}{\left|{}_j\varDelta {x}_j+\operatorname{}\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\varDelta {y}_j, $$
(11)

where, Δxj and Δyj are abbreviations of \( \left({x}_j^{gc}-{x}_j^{ori}\right) \) and \( \left({y}_j^{gc}-{y}_j^{ori}\right) \), and if the reconstructed function is expanded based on any local origins, there are another two terms in its integral average, and the mean constraint is rewritten as

$$ {\overline{\boldsymbol{u}}}_j={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j{\overline{x}}_j+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\overline{y}}_j+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\varDelta {x}_j+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\varDelta {y}_j. $$
(12)

On this basis, if the reconstructed function is integrated over a stencil cell k, we obtain

$$ {\overline{\boldsymbol{u}}}_k={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\underset{V_k}{\int}\left(x-{x}_j^{ori}\right)\kern0.5em dV+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\underset{V_k}{\int}\left(y-{y}_j^{ori}\right)\kern0.5em dV. $$
(13)

Similarly replace \( \left(x-{x}_j^{ori}\right) \), \( \left(y-{y}_j^{ori}\right) \) by \( \left(x-{x}_k^{gc}\right)+\left({x}_k^{gc}-{x}_j^{ori}\right) \) and \( \left(y-{y}_k^{gc}\right)+\left({y}_k^{gc}-{y}_j^{ori}\right) \) respectively, and reconstructed equation on stencil cell k could be written as

$$ {\overline{\boldsymbol{u}}}_k={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({\overline{x}}_k+\left({x}_k^{gc}-{x}_j^{ori}\right)\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({\overline{y}}_k+\left({y}_k^{gc}-{y}_j^{ori}\right)\right). $$
(14)

This equation is written for every stencil cell, of which there should be more than the number of derivatives to be solved to create an overconstrained system. If we write the novel mean constraint together with the Eq. (12) for each stencil cell, we have

$$ \left(\begin{array}{ccc}1& \overline{x}& \overline{y}\\ {}{\omega}_{j1}& {\omega}_{j1}{\hat{x}}_{j1}& {\omega}_{j1}{\hat{y}}_{j1}\\ {}{\omega}_{j2}& {\omega}_{j2}{\hat{x}}_{j2}& {\omega}_{j1}{\hat{y}}_{j2}\\ {}\vdots & \vdots & \vdots \\ {}{\omega}_{jN}& {\omega}_{jN}{\hat{x}}_{jN}& {\omega}_{j1}{\hat{y}}_{jN}\end{array}\right)\times \left(\begin{array}{c}{\boldsymbol{u}}_j^{ori}\\ {}{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\\ {}{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\end{array}\right)=\left(\begin{array}{c}{\overline{\boldsymbol{u}}}_j-\left({\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\varDelta {x}_j+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\varDelta {y}_j\right)\\ {}{\omega}_{j1}{\overline{\boldsymbol{u}}}_1\\ {}{\omega}_{j2}{\overline{\boldsymbol{u}}}_2\\ {}\vdots \\ {}{\omega}_{jN}{\overline{\boldsymbol{u}}}_N\end{array}\right), $$
(15)

where the first row is the mean constraint, and geometric terms could be expressed as

$$ {\displaystyle \begin{array}{l}{\hat{x^n{y}^m}}_{jk}\kern0.5em \equiv \frac{1}{A_k}{\int}_{V_k}{\left(\left(x-{x}_k^{gc}\right)+\left({x}_k^{gc}-{x}_j^{ori}\right)\right)}^n{\left(\left(y-{y}_k^{gc}\right)\operatorname{}+\operatorname{}\left({y}_k^{gc}-{y}_j^{ori}\right)\right)}^m dA\\ {}\kern1em =\sum \limits_{p=0}^m\sum \limits_{q=0}^n\frac{m!}{p!\left(m-p\right)!}\frac{n!}{q!\left(n-q\right)!}{\left({x}_k^{gc}-{x}_j^{ori}\right)}^q{\left({y}_k^{gc}-{y}_j^{ori}\right)}^p{\overline{x^{n-q}{y}^{m-p}}}_k,\end{array}} $$
(16)

and the weights are set to emphasize the geometrically adjacent data,

$$ {\omega}_{jk}=\frac{1}{\left|{\boldsymbol{x}}_j-{\boldsymbol{x}}_k\right|} $$
(17)

where xj and xk are local origins of central and stencil cells. When solving the Eq. (15), the unconstraint reconstructed equations are solved at first,

$$ \left(\begin{array}{ccc}{\omega}_{j1}& {\omega}_{j1}{\hat{x}}_{j1}& {\omega}_{j1}{\hat{y}}_{j1}\\ {}{\omega}_{j2}& {\omega}_{j2}{\hat{x}}_{j2}& {\omega}_{j1}{\hat{y}}_{j2}\\ {}\vdots & \vdots & \vdots \\ {}{\omega}_{jN}& {\omega}_{jN}{\hat{x}}_{jN}& {\omega}_{j1}{\hat{y}}_{jN}\end{array}\right)\times \left(\begin{array}{c}{\boldsymbol{u}}_j^{ori}\\ {}{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\\ {}{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\end{array}\right)=\left(\begin{array}{c}{\omega}_{j1}{\overline{\boldsymbol{u}}}_1\\ {}{\omega}_{j2}{\overline{\boldsymbol{u}}}_2\\ {}\vdots \\ {}{\omega}_{jN}{\overline{\boldsymbol{u}}}_N\end{array}\right). $$
(18)

Assuming that the current reconstructed problem is Ax = b, transposition matrix AT is multiplied at both sides of the equation as ATAx = ATb at first, and then, since the matrix ATA is the square matrix, we obtain x = (ATA)−1ATb.

On this basis, solution gradients as well as point-valued solution \( {\boldsymbol{u}}_j^{ori} \) could be obtained. Note, however, that the mean constraint is not satisfied currently, and therefore, the point-valued solution \( {\boldsymbol{u}}_j^{ori} \) should be recalculated as

$$ {\boldsymbol{u}}_j^{ori}={\overline{\boldsymbol{u}}}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j{\overline{x}}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\overline{y}}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\varDelta {x}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\varDelta {y}_j. $$
(19)

Thus, the mean constraint is always exactly satisfied and the obtained \( {\boldsymbol{u}}_j^{ori} \), as well as \( {\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j \) and \( {\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j \) is utilized for calculating the pointwise solution at per gauss point as shown in Eq. (4).

In short, we perform the solution reconstruction at a Gauss point over a face (i.e., obtain uL and uR) by first computing the point valued solutions and gradients at the local origins with the LSQR method and then use Eq. (4) or a higher-order variant to be presented later. It sets a stage for the employment of novel reference points on unstructured finite volume discretization from integral form without degrading the design accuracy.

3 Global-direction stencil based on the face-area-weighted centroid

In Section 2.2, we theoretically derive the feasibility of k-exact reconstruction based on any local origins. In this section, some commonly used stencil selection methods are introduced at first. And then, we briefly discuss the grid skewness on high-aspect-ratio triangular grids as well as the effect of skewness reduction by the employment of face-area-weighted centroid. Finally, to reduce the grid skewness and achieve a better reflection of flow anisotropy, the global-direction stencil with this novel centroid is introduced. The main content in this section has been analyzed in Ref. [23], while considering the completeness of the article, it is also discussed here.

3.1 Stencil selection methods

Commonly used stencils are vertex-neighbor and face-neighbor stencils. As Fig. 3 exhibits, face-neighbor stencil includes entire neighbor cells that share faces with the central cell, and vertex-neighbor stencil is similarly constructed by cells that share vertices with the central cell. Besides, both of them are topological-dependent, and therefore are limited by the original mesh topology. Furtherly, the stencil size of them is hard to accurately control, especially for vertex-neighbor stencil, and characteristics of flowfield cannot be well reflected.

Fig. 3
figure 3

(a) Vertex-neighbor and (b) face-neighbor stencils, where different numbers represent stencil layers (e.g., for face-neighbor stencil, the first layer stencil is composed of all cells that share faces with the central cell, and the second layer stencil consists of cells that share faces with the first layer stencil)

Apart from two commonly used stencils, in 2018, Xiong et al. [20], put forward the local-direction stencil selection method, by which selected stencil cells are along two local directions. As shown in Fig. 4a, on isotropic grid, two local directions are close to the normal and tangential directions of the wall, while on high-aspect-ratio triangular grids, as Fig. 4b displays, one of the local directions has severely deviated from the normal direction of the wall, and flow anisotropy is not well reflected. Besides, it is verified that on this grid type, accuracy loss and stability deterioration cannot be avoided [22], and the implementation process of this stencil selection method is quite complicated.

Fig. 4
figure 4

Local-direction stencil cells on (a) minor and (b) high-aspect-ratio triangular grids

In previous work, based on the existing problems on local-direction stencil, a novel global-direction stencil selection method [22, 23] was proposed. Compared with the local-direction stencil, problems mentioned above are well solved by this novel stencil. Specifically, for this method, two global directions, that is normal and tangential directions of the wall, are determined at first. And then, for each central cell, two lines which are parallel to global directions respectively and pass the cell centroid are generated. Finally, cells in a given set, such as the vertex-adjacent cells that intersect with these two lines are selected to construct the new stencil, and the stencil size is governed by layer of vertex-adjacent cells.

As Fig. 5 demonstrates, stencil cells selected by this method are always along two global directions no matter the grid with high aspect ratio or not. Therefore, flow anisotropy can be well reflected, and in corresponding numerical examples [23], global-direction stencil has a better numerical performance than commonly used vertex-neighbor and face-neighbor stencils as well as local-direction stencil. Besides, on cases with simple shapes, two global directions could be easily determined, while for complex surface, there are no analytical expressions to obtain the normal and tangential directions of the wall. Hence, in this situation, we can refer to the method of computing wall distance to get the corresponding normal vector [42], and construct the global-direction stencil.

Fig. 5
figure 5

Global-direction stencil cells on (a) minor and (b) high-aspect-ratio triangular grids

But after analysis, although global-direction stencil cells are always along the normal and tangential directions of the wall, no matter on grid with high aspect ratio or not, the only data obtained by k-exact reconstruction and required for the flux evaluation are solution vectors at the reference point rather than stencil cells themselves, and it is hard to guarantee whether flow anisotropy is well reflected or not. In the next section, we will focus on different locations of reference points.

3.2 Face-area-weighted centroid and global-direction stencil

In this section, we give a brief analysis about grid skewness and face-area-weighted centroid on triangular mesh, and explain the reason as well as the idea of the combination of global-direction stencil and this novel reference point in detail.

3.2.1 Grid skewness measure and face-area-weighted centroid

In introduction, we have analyzed that although there are various definitions of grid skewness, the same conclusion could be drawn that on high-aspect-ratio triangles, the grid skewness is always evident. Here, a typical skewness measure is taken into account.

As shown in Fig. 6, the grid skewness is defined at common face shared by two neighbor cells, where \( {\hat{\boldsymbol{e}}}_{jk} \) is a unit vector pointing from centroid of cell j to that of cell k, \( {\hat{\boldsymbol{n}}}_{jk} \) is the unit outward normal vector of common face. The grid skewness is measured by \( \left|{\hat{\boldsymbol{e}}}_{jk}\cdot {\hat{\boldsymbol{n}}}_{jk}\right| \), and from Fig. 6, it could be easily found that the non-skewed grid has the measure one, and with the increase of cell aspect ratio, the skewness measure of highly-skewed grid is close to zero. On this basis, the most direct intuition is using isotropic or minor-aspect-ratio triangles to ensure the computing process is carried out on non-skewed grids, while for some typical flows, such as boundary-layer-type flow, solutions in boundary layer are changed dramatically, particularly along the normal direction of the wall. Therefore, to enhance the resolution in this local field, highly-anisotropic grids cannot be avoided, and we can only rely on a novel local origin to reduce the grid skewness and improve the numerical performance.

Fig. 6
figure 6

Grid skewness measure on grids with different aspect ratios

For high-aspect-ratio triangular grids, a novel face-area-weighted centroid was proposed by Nishikawa [28] on the second-order differential unstructured finite volume solver. By the employment of this reference point, the grid skewness is almost eliminated. In the following analysis, we will briefly summarize the face-area-weighted formula and the effect about skewness reduction. The specific derivation has been given in Ref. [28], and the main results are listed in this section.

A typical choice for the local origin is the geometric centroid that could be written for a triangle by the arithmetic average of face midpoints,

$$ \left({x}_j,\kern0.5em {y}_j\right)=\frac{1}{3}\sum \limits_{k=1}^3\left({x}_{mk},\kern0.5em {y}_{mk}\right), $$
(20)

where, (xmk,  ymk) is coordinate of the k-th face centroid. Besides, we generalize Eq. (20) to the face-area-weighted centroid formula [28]:

$$ \left({x}_j,\kern0.5em {y}_j\right)=\frac{\sum \limits_{k=1}^3{\hat{A}}_{jk}^p\left({x}_{mk},\kern0.5em {y}_{mk}\right)}{\sum \limits_{k=1}^3{\hat{A}}_{jk}},{\hat{A}}_{jk}=\frac{A_{jk}}{\underset{k\in \left\{1,2,3\right\}}{\max }{A}_{jk}}, $$
(21)

where, Ajk is the area or length (in 2D) of the face shared by two neighbor cells, and p (> 0) is a real value that controls the skewness degree. Note when p = 0, face-area-weighted centroid is just consistent with the geometric centroid.

As shown in Fig. 7, where h is the grid spacing in y-direction and R is cell aspect ratio. For a typical triangular grid in Cartesian-coordinate system, skewness measure along the x-direction approaches one, while in y-direction is nearly zero. Thus, in the following analysis, we will focus on common faces shared by cell 1 and cell 2 as well as cell 2 and cell 3.

Fig. 7
figure 7

A typical high-aspect-ratio triangular grid in Cartesian-coordinate system

For these two common faces, if the geometric centroid is utilized, the grid skewness measure is [28]

$$ {d}_{12}=\frac{2R}{R^2+1},\kern4.5em {d}_{23}=\frac{2}{\sqrt{R^2+4}}. $$
(22)

With the increase of cell aspect ratio, both measures nearly equal to zero. Hence the grid is highly skewed. But if we employ the face-area-weighted centroid and when parameter p is equal to 2, the skewness measure becomes [28]

$$ {d}_{12}^{p=2}=1,\kern4.5em {d}_{23}^{p=2}=1-\frac{1}{2{R}^2}. $$
(23)

As a result, for high-aspect-ratio triangular grid, the grid skewness could be eliminated, and note, moreover, that with the increase of cell aspect ratio, the grid skewness, at the face like between cell 2 and cell 3, could be further reduced. Besides, it is demonstrated in Fig. 8 when the face-area-weighted centroid is employed, line connecting the novel local origins is almost parallel to the normal direction of the wall, and the serrated phenomenon exhibited on geometric centroid is effectively avoided.

Fig. 8
figure 8

(a) Geometric and (b) face-area-weighted centroids on high-aspect-ratio triangular grids

This special distribution just coincides with our original motivation of designing the global-direction stencil, and in previous work, we combined the global-direction stencil and face-area-weighted centroid on the second-order unstructured finite volume solver in differential form to capture the flow anisotropy more accurately. After verification, better computational accuracy and stabilities are obtained by this novel method. In this work, we further investigate feasibilities of extension to the integral form.

3.2.2 Combination of global-direction stencil and face-area-weighted centroid

In Section 3.1, we have demonstrated that compared with commonly used vertex-neighbor and face-neighbor stencils, although global-direction stencil cells are along the normal and tangential directions of the wall, the only data obtained by k-exact reconstruction and required for the flux evaluation are solution vectors at cells local origins, rather than stencil cells themselves, and if we still use the geometric centroid, grid skewness is unable to be eliminated. Therefore, it is hard to guarantee whether the flowfield characteristics are well captured or not. But when face-area-weighted centroid is employed, the distribution of reference points is more regular, and the line connecting them is almost along the normal direction of the wall.

Besides, in Section 2.2, we have derived the k-exact reconstruction process based on any local origins, and therefore, in this section, we give the method of combining the global-direction stencil and face-area-weighted centroid in detail to realize the unified direction of both stencil cells and local origins. Global-direction stencils with geometric and face-area-weighted centroids on grids with straight and curved boundaries are displayed in Figs. 9 and 10.

Fig. 9
figure 9

Global-direction stencil cells combined with (a) geometric centroids and (b) face-area-weighted centroids on grid with the straight boundary

Fig. 10
figure 10

Global-direction stencil cells combined with (a) geometric centroids and (b) face-area-weighted centroids on grid with the curved boundary

From these two figures, we can easily find although same stencil cells are selected, geometric centroids are deflective and exhibit serrated phenomenon. Particularly on high-aspect-ratio triangular grids, the mentioned phenomenon is much more evident. By comparison, when face-area-weighted centroids are employed, line connecting them is close to the normal direction of the wall no matter on grid with straight or curved boundaries, and it is consistent with one of the global directions. As a result, the flow anisotropy could be well reflected, and the grid skewness is reduced.

On unstructured finite volume method from differential form, both computational accuracy and stabilities are greatly improved by this novel method, and in the next section, four representative numerical examples are designed to verify the effectiveness and superiorities of this novel method on integral unstructured finite volume solver.

4 Numerical examples

To examine the effectiveness of global-direction stencil with face-area-weighted centroid on unstructured finite volume method from integral form, in this section, numerical examples governed by linear convective, Euler and Laplacian equations are utilized. For comparison, these numerical examples are simulated with four different stencils, including vertex-neighbor and face-neighbor stencils, as well as the global direction stencil with geometric centroid and face-area-weighted centroid. Note that from Section 3.2.1, we can find the face-area-weighted centroid could be further distinguished by different p values, and when p = 2, skewness has been almost eliminated. Hence, the face-area-weighted centroid in this section refers to the result of p = 2. In addition, to simplify the presentation in the following analysis, different stencils are abbreviated as follows (Table 1).

Table 1 Abbreviation of different stencils

More importantly, in this section, we make a discussion about the result exhibition of unstructured finite volume method from integral form, and notice that the point-valued solution \( {\boldsymbol{u}}_j^{ori} \) at the face-area-weighted centroid obtained by the mean constraint of k-exact reconstruction algorithm with global-direction stencil is more accurate than the cell-averaged solution \( {\overline{\boldsymbol{u}}}_j \). The specific discussion will be given in Section 4.1.3.

4.1 Manufactured boundary layer (governed by linear convective equation)

In this section, we first use the Method of Manufactured Solutions (MMS) [43,44,45,46] on linear convective equation to verify feasibilities of the employment of face-area-weighted centroid on unstructured finite volume discretization from integral form, and examine the numerical performance of global-direction stencil with this novel reference point. The model equation can be written as

$$ \frac{\partial u}{\partial t}+\boldsymbol{\alpha} \cdot \nabla u=0, $$
(24)

where, \( \boldsymbol{\alpha} =\left(\cos \frac{\pi }{16},\kern0.5em \sin \frac{\pi }{16}\right) \) is the constant convective velocity, and to simulate characteristics of boundary layer, the manufactured solution is

$$ u\left(x,y\right)=1-{e}^{\frac{-\left(y-{y}_0\right)}{\sqrt{c\mu \left(x-{x}_0\right)}}}, $$
(25)

where c = 0.59 is a constant, and parameter μ is utilized to control the thickness of boundary layer. Flowfields corresponding to different μ values are shown in Fig. 11, and in the following test, μ is set as 10− 6.

Fig. 11
figure 11

(a) μ = 10-6 (b) μ = 10-8. Flowfields with different μ values

By bringing the manufactured solution to Eq. (24), the modified equation with source term could be written as

$$ \frac{\partial u}{\partial t}+\boldsymbol{\alpha} \cdot \nabla u=\frac{a_y\left(y-{y}_0\right){e}^{\frac{-\left(y-{y}_0\right)}{\sqrt{c\mu \left(x-{x}_0\right)}}}}{\sqrt{c\mu \left(x-{x}_0\right)}}-\frac{a_x\left(y-{y}_0\right){e}^{\frac{-\left(y-{y}_0\right)}{\sqrt{c\mu \left(x-{x}_0\right)}}}}{2\sqrt{c\mu {\left(x-{x}_0\right)}^3}}. $$
(26)

On this basis, the manufactured solution represents the analytical solution of this modified governing equation, and we can calculate both L2 and L errors of different stencils.

$$ \left\{\begin{array}{l}{L}_2=\sqrt{\frac{\sum_{j=1}^N\left[{\left({\tilde{u}}_j-{u}_{analy}\right)}^2\cdot {A}_j\right]}{\sum_{j=1}^N{A}_j}},\\ {}{L}_{\infty }=\max {\left|{\tilde{u}}_j-{u}_{analy}\right|}_{j=1,N},\end{array}\right. $$
(27)

where \( {\tilde{u}}_j \), uanaly and Aj are numerical and analytical solutions and area of cell j respectively. Note, however, that errors are computed for the point valued solutions recovered at the local origin from the cell averaged solutions stored at cells rather than for the cell-averaged solutions for the reason explained in Section 4.1.3. Therefore, the exact solution used to compute the error is a point value evaluated at the local origin. Besides, as shown in Fig. 12, both regular and randomly perturbed triangular grids are used in this numerical example, and two levels of grid stretching, including 103 and 104 these two wall cell aspect ratios (AR), are tested. In each level, five sets of triangular grids from the coarsest to finest are generated within x, y [0.05,  1.05] × [0,  0.001].

Fig. 12
figure 12

(a) Regular and (b) randomly perturbed triangular grids

During the grid generation process, nodes in x-direction are equidistantly distributed, while the y-coordinates of different nodes are determined by

$$ {y}_{j+1}={\left.{y}_j+{\hat{h}}_y{\beta}^{j+1}\right|}_{j=1,2,\dots, {N}_y}, $$
(28)

where \( {\hat{h}}_y \) is the first layer vertical spacing, and β is the stretching factor that could be computed by the known condition yN = 10−3. Besides, distribution of five sets of background quadrilateral grids from the coarsest to finest is shown in Table 2.

Table 2 Distribution of background quadrilateral grids in x and y directions

On this basis, computational accuracy and discretization errors of four different stencils on integral unstructured finite volume solver are counted and given in Figs. 13 and 14.

Fig. 13
figure 13

(a) L2 and (b) L errors. Errors of different stencils on regular grids with AR = 103

Fig. 14
figure 14

(a) L2 and (b) L errors. Errors of different stencils on regular grids with AR = 104

4.1.1 Computational results on regular grids

From Figs. 13 and 14, we can easily find the second-order accuracy is achieved by all stencils we tested for both L2 and L errors. Therefore, the effectiveness of employing face-area-weighted centroid on unstructured finite volume method from integral form is verified. In addition, it is also proved that the derivation of k-exact reconstruction based on any local origins is feasible, and its correctness has been well demonstrated by computational results.

Besides, combining the specific data listed in Table 3, we find both L2 and L errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested, and when AR = 104, computational accuracy of L errors could be obviously improved. What’s more, from Table 3, we notice that G-Stencil requires the least number of stencil cells, and the efficiency can also be improved.

Table 3 Computational errors and average stencil size of different stencils on the finest grid

4.1.2 Computational results on randomly perturbed grids

On randomly perturbed triangular grids, we also test the case of AR = 103 and 104. For simplicity, results of AR = 104 are given here.

From Fig. 15, we find errors on randomly perturbed grids demonstrate the similar trends to results on regular grids, and when face-area-weighted centroid is employed on the global-direction stencil, both L2 and L errors are greatly reduced. As a result, on linear convective governing equation, the effectiveness as well as superiorities of this novel method is well verified.

Fig. 15
figure 15

(a) L2 and (b) L errors. Errors of different stencils on randomly perturbed grids with AR = 104

4.1.3 Discussion about results exhibition of finite volume discretization from integral form

In this work, we also realize that no matter what numerical scheme is adopted, the main target is obtaining a more accurate result or flowfield, and as mentioned above, for unstructured finite volume method from integral form, the obtained solution vector is the cell-averaged value. Consequently, the result is also exhibited by the cell average, and there is no relevant research to further find the point-valued solution and compare it with the cell average, since the point-valued solution could not be obtained by solving the integral governing equation directly. But we notice that by k-exact reconstruction, the solution at the local origin could be obtained by the mean constraint shown as follows,

$$ {\boldsymbol{u}}_j^{ori}={\overline{\boldsymbol{u}}}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j{\overline{x}}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\overline{y}}_j-{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)-{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right), $$
(29)

and computational results shown in Sections 4.1.1 and 4.1.2 are just counted at the current local origin. In order to further contrast with the result counted by the cell average, we additionally test the discretization errors counted by the cell-averaged value for the second-order finite volume solver, where the spatial discretizations are based on both geometric centroid and face-area-weighted centroid, and compare it with errors counted at the geometric centroid and face-area-weighted centroid respectively.

Besides, when counting the cell-averaged errors, the source term and analytical solution need to be integrated over the control volume, and in this section, three-point quadrature rule is employed [34, 35]. For a triangular cell, the integral point is just per face centroid and the corresponding weights are all equal to \( \frac{1}{3} \).

For the convenience of analysis, the cell average and point value at the geometric centroid and face-area-weighted centroid are all based on the identical global-direction stencil, i.e., G-Stencil. Besides, “Avg” and “Ref” these two labels in Figs. 16 and 17 represent the cell-averaged and point-valued errors respectively.

Fig. 16
figure 16

(a) L2 and (b) L errors. Cell-averaged and point-valued discretization errors of G-Stencil with the geometric centroid and face-area-weighted centroid on grids with AR = 103

Fig. 17
figure 17

(a) L2 and (b) L errors. Cell-averaged and point-valued discretization errors of G-Stencil with the geometric centroid and face-area-weighted centroid on grids with AR = 104

Combining the results shown in the two figures and Table 4, we notice that no matter the spatial discretization is based on the geometric centroid or face-area-weighted centroid, cell averages are almost identical and both of them are close to point value at the geometric centroid.

Table 4 Computational errors on the finest grid and L2, L accuracy between the last two grids of different stencils on the second- and third-order unstructured finite volume solver

In addition, both cell-averaged errors and errors at the geometric centroid are always higher than that at the face-area-weighted centroid. On this basis, we wonder whether the computational result could be exhibited by the value at a certain local origin, that is \( {\boldsymbol{u}}_j^{ori} \), which could be obtained by the k-exact reconstruction.

By further consideration, we realize there may be three requirements that should be satisfied if the computational result exhibited by the value at a certain local origin. The first one is that, in every time step, the point value at the local origin could be obtained. Secondly, the designed order of accuracy is achieved by the point-valued evaluation. Finally, solution at the local origin should be more accurate than the cell-averaged value, and it means errors counted at this specific local origin are lower than the cell-averaged errors. If these three requirements are all satisfied, we believe the result could be evaluated by the point value at this special local origin.

Firstly, we derive the k-exact reconstruction algorithm based on any local origins, and the corresponding point-valued solution \( {\boldsymbol{u}}_j^{ori} \) in reconstructed polynomial at every time step is obtained by the mean constraint. Hence, the first requirement is satisfied. Besides, as results shown in Sections 4.1.1 and 4.1.2, where point values are computed at the geometric centroid and face-area-weighted centroid, the second-order accuracy is always achieved. More importantly, when the global-direction stencil is utilized, cell averages are close to point value at the geometric centroid, and both of them are not as accurate as the point value at the face-area-weighted centroid. Consequently, both the second and third requirements are satisfied.

Therefore, combining the analysis and numerical verification, we notice that compared with the traditional finite volume scheme from integral form, result exhibited at the specific face-area-weighted centroid with the global-direction stencil is more accurate, and the designed order of accuracy could also be achieved.

In short, although the spatial discretization is from integral form, it is not necessary to exhibit the result by the cell average, and the point value at the face-area-weighted centroid is more appropriate to reflect the computational result. Hence, in the following test, errors are all computed for the point-valued solutions recovered at the local origin.

4.2 Supersonic vortex flow (governed by Euler equations)

To further test the effectiveness and feasibilities of this novel method on Euler equations, in this section, the supersonic vortex flow is introduced. Computational domain is two concentric circular arcs with radius ri = 1 and r0 = 1.384 located in the first quadrant. These two circular arcs represent the inviscid wall boundary, and the flow at both inlet and outlet is supersonic. Analytical solution [47] of this numerical example could be derived by isentropic relation and is given as follows,

$$ \left\{\begin{array}{l}\rho ={\rho}_i{\left(1+\frac{\gamma -1}{2}{M}_i^2\left(1-{\left(\frac{r_i}{r}\right)}^2\right)\right)}^{\left(1/\left(\gamma -1\right)\right)},\\ {}P=\frac{\rho^{\gamma }}{\gamma },\kern1em \left\Vert \boldsymbol{v}\right\Vert =\frac{c_i{M}_i}{r},\end{array}\right. $$
(30)

where the value of Mach number at the inner radius is Mi = 2.25 and the density ρi = 1. Besides, the sound speed is calculated as

$$ {c}_i=\sqrt{\gamma {P}_i/{\rho}_i}=1, $$
(31)

Flow structure of this numerical case is shown in Fig. 18.

Fig. 18
figure 18

(a) Mach number (b) density. Flowfields of supersonic vortex flow

For this numerical example, both regular and randomly perturbed triangular grids are utilized. As shown in Fig. 19, where the regular grid is generated by splitting the background quadrilateral grid with right diagonals, and randomly perturbed grid is generated by introducing the random node perturbation to the regular grid with topology and the number of cells unchanged. Besides, two levels of grid aspect ratios are employed, and in each level, five sets of grids from the coarsest to finest are generated.

Fig. 19
figure 19

(a) Regular and (b) randomly perturbed triangular grids

Considering that the aspect ratio of grid with the curved boundary is not a fixed value, similarly, the wall cell aspect ratio is utilized, and two aspect ratios are approximately equal to 4 and 8 respectively. The distribution of background quadrilateral grids in the radical and circumferential is shown in Table 5.

Table 5 Distribution of background quadrilateral grids in radical and circumferential directions

4.2.1 Computational errors on regular grids

As shown in Figs. 20 and 21, for three stencils with the geometric centroid, both L2 and L errors of global-direction stencil are the lowest. On this basis, we employ the face-area-weighted centroid, and discretization errors are further reduced. Besides, according to the specific data listed in Table 6, higher-order accuracy on both L2 and L errors is achieved by G-Stencil (F-a-w centroid). Particularly, L2 errors can reach 2.418 and 2.472 order accuracy. As a result, the employment of face-area-weighted centroid on unstructured finite volume method from integral form will not deteriorate the designed order of accuracy but greatly improve the computational accuracy of the finite volume solver after being introduced in global-direction stencil.

Fig. 20
figure 20

(a) L2 and (b) L errors. Errors of different stencils on regular grids with AR ≈ 4

Fig. 21
figure 21

(a) L2 and (b) L errors. Errors of different stencils on regular grids with AR ≈ 8

Table 6 Computational errors on the finest grid and L2, L accuracy between the last two grids of different stencils

4.2.2 Computational errors on randomly perturbed grids

In order to adequately illustrate the numerical performance of different methods, errors on randomly perturbed grids are also counted. For simplicity, results of AR ≈ 8 are given in Fig. 22.

From Fig. 22, errors exhibit the similar trends to that of regular grids, and both average and max errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested. Therefore, correctness of k-exact reconstruction based on any local origins as well as the effectiveness and superiorities of the global-direction stencil with face-area-weighted centroid is verified on Euler equations as well.

Fig. 22
figure 22

(a) L2 and (b) L errors. Errors of different stencils on randomly perturbed grids with AR ≈ 8

4.3 Subsonic flow over a NACA0012 airfoil (governed by Euler equations)

In Section 3.1, we give a brief introduction about the determination of global directions based on wall distance [42], and in this section, the subsonic flow over a NACA0012 airfoil governed by Euler equations is simulated to examine the effectiveness of the mentioned method. The angle of attack is α = 0, and the initial condition is Ma = 0.5.

Regular and randomly perturbed triangular grids with both O-type and C-type topologies are utilized in this numerical example. For O-type regular grid, the first layer spacing in normal direction is 10−3, and there are 201 and 71 grid points distributed in circumferential and normal directions respectively. Besides, for C-type regular grid, the background quadrilateral grid is split with random diagonals, and the distribution is 226 × 66 (129 points on airfoil surface). Regular and randomly perturbed grids with both O-type and C-type grids near the airfoil are shown in Figs. 23 and 24.

Fig. 23
figure 23

(a) Regular and (b) randomly perturbed triangular grids with O-type over a NACA0012 airfoil

Fig. 24
figure 24

(a) Regular and (b) randomly perturbed triangular grids with C-type over a NACA0012 airfoil

Besides, to intuitively contrast the effect of geometric centroid and face-area-weighted centroid on cells adjacent to the airfoil, we display these two reference points on randomly perturbed grid with O-type respectively.

From Fig. 25, we can easily find that compared with the geometric centroid, the distribution of face-area-weighted centroid is more regular, and the reference points are almost spread along the normal direction of the wall, especially on grid adjacent to the airfoil surface. On this basis, we count lift and drag coefficients as well as the residual of different stencils with O-type and C-type grids respectively. Before calculating these two coefficients, pressures at the per gauss point of airfoil boundary face are needed at first, and they are interpolated by the point value at the local origin, and the corresponding formula is given as follows,

$$ {p}_{iG}={p}_i^{ori}+{\left.\frac{\partial p}{\partial x}\right|}_i\left({x}_{iG}-{x}_i^{ori}\right)+{\left.\frac{\partial p}{\partial y}\right|}_i\left({y}_{iG}-{y}_i^{ori}\right), $$
(32)

where, piG is the pressure at the gauss point of wall-adjacent face. Similarly, for the second-order accurate unstructured finite volume discretization, the gauss point is just the face centroid. Besides, (xiG,  yiG) is the coordinate of the current gauss point, and \( \left({x}_i^{ori},\kern0.5em {y}_i^{ori}\right) \) is the coordinate of the local origin Ci.

Fig. 25
figure 25

(a) Geometric centroid and (b) face-area-weighted centroid of perturbed triangular grid over a NACA0012 airfoil

4.3.1 Results on regular and randomly perturbed O-type grid

From results of lift and drag coefficients shown in Figs. 26 and 27, we can easily find from local amplifications of both regular and randomly perturbed grids that the lift coefficient of G-Stencil (F-a-w centroid) is the lowest among four different stencils, and in the same time steps, there is almost no oscillations on V-Stencil and G-Stencil (Cell centroid) as well as G-Stencil (F-a-w centroid). However, compared with these three stencils, oscillations of the F-Stencil are quite evident, and have not yet converged.

Fig. 26
figure 26

Cl on (a) regular and (b) perturbed triangular grids

Fig. 27
figure 27

Cd on (a) regular and (b) perturbed triangular grids

Similar phenomenon could be concluded on drag coefficient, where the result of G-Stencil (F-a-w centroid) is lower than V-Stencil and G-Stencil (Cell centroid) but slightly higher than F-Stencil. Nevertheless, oscillations of F-Stencil are quite obvious, and the convergence property of this stencil is much poorer than another three stencils.

Conclusions related to lift and drag coefficients mentioned above can be better illustrated by residuals of different stencils. As exhibited in Fig. 28, we can easily find that for three stencils with geometric centroid, the convergence speed of G-Stencil is close to the commonly used V-Stencil, and in the same time steps, residuals of these two stencils are decreased by 7 orders of magnitude, while only 4 orders of magnitude is decreased on F-Stencil. On this basis, with the introduction of face-area-weighted centroid on global-direction stencil, the residual is further decreased, and therefore, the convergence property of integral unstructured finite volume solver is greatly improved by global-direction stencil with face-area-weighed centroid. In the following section, we will further examine the numerical performance on C-type grids.

Fig. 28
figure 28

(a) Residuals on regular and (b) randomly perturbed grids

4.3.2 Results on regular and randomly perturbed C-type grid

From the results shown in Figs. 29 and 30, we can easily find on C-type grid, more accurate lift and drag coefficients are obtained by the global-direction stencil with face-area-weighted centroid as well. However, oscillations of the face-neighbor stencil are extremly evident compared with another three stencils we tested, and the result does not converge at all. It could be well illustrated by the residual shown in Fig. 31.

Fig. 29
figure 29

Cl on (a) regular and (b) perturbed triangular grids

Fig. 30
figure 30

Cd on (a) regular and (b) perturbed triangular grids

Fig. 31
figure 31

(a) Residuals on regular and (b) randomly perturbed grids

As shown in Fig. 31, vertex-neighbor and global-direction stencils have a similar convergence rate, and it is improved by the global-direction stencil with face-area-weighted centroid. However, the result of face-neighbor stencil is almost diverged within the same time steps.

In short, feasibilities as well as superiorities of global-direction stencil with face-area-weighted centroid on unstructured finite volume method from integral form are demonstrated again. And according to results displayed in Sections 4.1, 4.2 and 4.3, we can easily find computational accuracy, efficiency and convergence speed are all improved by this novel method.

4.4 Dissipative term evaluation (governed by Laplacian equation)

From Section 4.1 to 4.3, three numerical examples govern by linear convective and Euler equations are simulated, and the correctness of k-exact reconstruction based on any local origins is verified. Moreover, on unstructured finite volume method from integral form, the global-direction stencil with face-area-weighted centroid also has a better computational accuracy, efficiency as well as convergence rate for convective flux evaluation.

But, the numerical performance on dissipative term also needs to be further tested to set stage for the extension to viscous problems. In this section, the MMS method is also used on Laplacian equation to examine discretization errors and convergence speed of different stencils. Laplacian model equation could be formulated as

$$ \frac{\partial \phi }{\partial t}-\nabla \cdot \left(\nu \nabla \phi \right)=0, $$
(33)

where ν is a constant equaling to 1. Integrating this model equation over the control volume, we obtain

$$ \frac{\partial {\overline{\phi}}_j}{\partial t}-\frac{1}{V_j}{\int}_{V_j}\nabla \cdot \left(\nu \nabla {\phi}_j\right)\kern0.5em dV=0. $$
(34)

According to the divergence theorem, Eq. (34) could be transformed as

$$ \frac{\partial {\overline{\phi}}_j}{\partial t}+\frac{1}{V_j}{\oint}_{\partial {V}_j}{\left(\nu \nabla {\phi}_j\right)}_n dS=0. $$
(35)

From Eq. (35), we can clearly find the difference between convective term and dissipative term. For dissipative term computation, what we need to compute at the face integral point is not two state vectors obtained by owner and neighbor cells, but solution gradient, which could be evaluated by the arithmetic average of left and right cells for the second-order unstructured finite volume method,

$$ \nabla {\phi}_{int}^f=\frac{1}{2}\left(\nabla {\phi}_{int}^L+\nabla {\phi}_{int}^R\right). $$
(36)

On this basis, in order to improve the stability and reduce the truncation errors, the average face gradient is always added a solution jump term [28, 48, 49]. Here are many different schemes for adding the jump term, such as the edge-normal (EN) and face-tangent (FT) schemes [50]. Based on the FT scheme, a novel α-damping scheme is proposed by Nishikawa [51,52,53,54]. The accuracy as well as robustness of second-order unstructured finite volume solver is greatly improved by this scheme, and it could be formulated as

$$ \nabla {\phi}_{int}^f=\frac{1}{2}\left(\nabla {\phi}_{int}^L+\nabla {\phi}_{int}^R\right)+\frac{\alpha }{\left|{\hat{\boldsymbol{e}}}_{jk}\cdot {\hat{\boldsymbol{n}}}_{jk}\right|{L}_{jk}}\left({\phi}_R^{\ast }-{\phi}_L^{\ast}\right){\hat{\boldsymbol{n}}}_{jk}, $$
(37)

where, α is a damping coefficient, and here, a special value is α = 4/3, which has been known to provide superior accuracy and robustness [28, 48, 51]. \( {\hat{\boldsymbol{n}}}_{jk} \) is the unit face normal vector, and \( {\hat{\boldsymbol{e}}}_{jk} \) is the unit vector between two cell centroids. Ljk is the length from one centroid to another.

$$ {\hat{\boldsymbol{e}}}_{jk}=\frac{{\boldsymbol{x}}_k-{\boldsymbol{x}}_j}{L_{jk}},\kern3em {L}_{jk}=\left|{\boldsymbol{x}}_k-{\boldsymbol{x}}_j\right|. $$
(38)

Besides, \( {\phi}_L^{\ast } \) and \( {\phi}_R^{\ast } \) are solutions linearly reconstructed at integral point,

$$ {\phi}_L^{\ast }={\phi}_j+\frac{1}{2}\nabla {\phi}_j\cdot \left({\boldsymbol{x}}_k-{\boldsymbol{x}}_j\right),\kern2.5em {\phi}_R^{\ast }={\phi}_k+\frac{1}{2}\nabla {\phi}_k\cdot \left({\boldsymbol{x}}_j-{\boldsymbol{x}}_k\right) $$
(39)

On the basis of determining the discretization method of dissipative term, the manufactured solution of Eq. (33) is

$$ \phi =\cos \left(\omega x\right)\sin \left(10\omega y\right), $$
(40)

where ω = 2π is a constant, and the source term could be derived as

$$ s=\left({10}^2+1\right){\omega}^2\cos \left(\omega x\right)\sin \left(10\omega y\right). $$
(41)

As shown in Fig. 32, the simulation is carried out on five sets of triangular grids generated by splitting the quadrilateral grids with regular or random diagonals, and the cell aspect ratio is AR = 10. The distribution of these five sets of grids is from 10 × 10 to 80 × 80, and the computational domain is x, y [−0.5,  0.5] × [−5 × 10−2,  5 × 10−2].

Fig. 32
figure 32

Quadrilateral grids split by (a) regular and (b) random diagonals

From Eq. (37), we notice that the damping term is greatly influenced by the grid skewness measure. On this basis, we first count it on both two grid types, and list the results in Table 7.

Table 7 Grid skewness measure of geometric centroid and face-area-weighted centroid on two grid types

As a result, superiorities of face-area-weighted centroid are illustrated again, and with the employment of this novel centroid, the grid skewness is almost eliminated. In order to reflect the relative position of two different centroids on grid with regular and random diagonals, the distribution of these two local origins is shown in Fig. 33.

Fig. 33
figure 33

Geometric and face-area-weighted centroids on grids with (a) regular and (b) randomly diagonals

From Fig. 33, we can easily find that face-area-weighted centroids are almost distributed along y-axis no matter on grid with regular or random diagonals, but the distribution of geometric centroid is irregular, and according to the data listed in Table 7, it has been demonstrated that the evident skewness will be introduced by this traditional centroid. Based on above analysis, we count discretization errors of four different stencils, and corresponding results are shown in Figs. 34 and 35 and Table 8.

Fig. 34
figure 34

(a) L2 and (b) L errors. Errors of different stencils on grids with regular diagonals

Fig. 35
figure 35

(a) L2 and (b) L errors. Errors of different stencils on grids with random diagonals

Table 8 Computational errors on the finest regular grid and L2, L accuracy between the last two grids of different stencils

From Figs. 34 and 35 and combining the specific data listed in Table 8, we find that discretization errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested no matter on grid with regular or random diagonals, and the computational accuracy is greatly improved by this novel method.

On this basis, residuals on the coarsest and finest regular grids of different stencils are shown in Fig. 36.

Fig. 36
figure 36

Residuals on the (a) coarsest and (b) finest regular grids

From Fig. 36, we can easily find that among three stencils with the geometric centroid, the G-Stencil (Cell centroid) has a faster convergence speed than two commonly used V-Stencil and F-Stencil. What’s more, when face-area-weighted centroid is employed, the convergence rate is further promoted. Hence, except for convective term, on dissipative term, a better computational accuracy and faster convergence rate are also realized by this novel method. In short, feasibilities and superiorities of employing face-area-weighted centroid on unstructured finite volume discretization from integral form are verified again on Laplacian model equation.

5 Accuracy analysis of the method on the third-order discretization

In Section 2.2, the k-exact reconstruction based on any local origins is derived, and from Section 4.1 to 4.4, we notice that on the second-order finite volume solver, the result of global-direction stencil with face-area-weighted centroid is much more accurate than two commonly used vertex-neighbor and face-neighbor stencils, as well as the global-direction stencil with the geometric centroid. In this section, the k-exact reconstruction and spatial discretization based on any local origins are derived at first. On this basis, a typical numerical example based on the Method of Manufactured Solutions (MMS) [43,44,45,46] is employed for the accuracy analysis.

5.1 Derivation of higher-order k-exact reconstruction based on any local origins

For higher-order k-exact reconstruction, e.g., the third-order, the reconstructed polynomial based on any local origins could be formulated as

$$ {\displaystyle \begin{array}{l}{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left(x-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left(y-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left(x-{x}_j^{ori}\right)}^2\\ {}\kern1em +\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\left(x-{x}_j^{ori}\right)\left(y-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left(y-{y}_j^{ori}\right)}^2,\end{array}} $$
(42)

where \( \left({x}_j^{ori},\kern0.5em {y}_j^{ori}\right) \) is the coordinate of the current local origin, which is not a necessary geometric centroid, and if the \( \left(x-{x}_j^{ori}\right) \) and \( \left(y-{y}_j^{ori}\right) \) are replaced by \( \left(x-{x}_j^{gc}\right)+\left({x}_j^{gc}-{x}_j^{ori}\right) \) and \( \left(y-{y}_j^{gc}\right)+\left({y}_j^{gc}-{y}_j^{ori}\right) \) respectively, we obtain

$$ {\displaystyle \begin{array}{l}{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left(x-{x}_j^{gc}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left(y-{y}_j^{gc}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left(x-{x}_j^{gc}\right)}^2\\ {}\kern1em +\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\left(x-{x}_j^{gc}\right)\left(y-{y}_j^{gc}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left(y-{y}_j^{gc}\right)}^2+{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)\left(x-{x}_j^{gc}\right)\\ {}\kern1em +\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\left({y}_j^{gc}-{y}_j^{ori}\right)\left(x-{x}_j^{gc}\right)\kern1em +\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\left({x}_j^{gc}-{x}_j^{ori}\right)\left(y-{y}_j^{gc}\right)+{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)\left(y-{y}_j^{gc}\right)\\ {}\kern1em +{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left({x}_j^{gc}-{x}_j^{ori}\right)}^2\\ {}+\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\left({x}_j^{gc}-{x}_j^{ori}\right)\left({x}_j^{gc}-{x}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left({y}_j^{gc}-{y}_j^{ori}\right)}^2.\end{array}} $$
(43)

For k-exact reconstruction, no matter the reconstructed function is formulated based on which reference point, the mean constraint requires [34, 35]

$$ \frac{1}{V_j}\underset{V_j}{\int }{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV=\underset{V_j}{\int }u\left(\boldsymbol{x}\right)\kern0.5em dV={\overline{\boldsymbol{u}}}_j, $$
(44)

Besides, integrating the reconstructed polynomial over the control volume, we have

$$ {\displaystyle \begin{array}{l}{\int}_{V_j}{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV\equiv {\overline{\boldsymbol{u}}}_j={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\underset{V_j}{\int}\left(x-{x}_j^{gc}\right)\kern0.5em dV+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\underset{V_j}{\int}\left(y-{y}_j^{gc}\right)\kern0.5em dV\\ {}+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\underset{V_j}{\int }{\left(x-{x}_j^{gc}\right)}^2\kern0.5em dV+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\underset{V_j}{\int}\left(x-{x}_j^{gc}\right)\left(y-{y}_j^{gc}\right)\kern0.5em dV+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\underset{V_j}{\int }{\left(y-{y}_j^{gc}\right)}^2\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)\underset{V_j}{\int}\left(x-{x}_j^{gc}\right)\kern0.5em dV+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)\underset{V_j}{\int}\left(x-{x}_j^{gc}\right)\kern0.5em dV\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)\underset{V_j}{\int}\left(y-{y}_j^{gc}\right)\kern0.5em dV\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)\underset{V_j}{\int}\left(y-{y}_j^{gc}\right)\kern0.5em dV\\ {}\kern1em +{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)\kern1em +\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left({x}_j^{gc}-{x}_j^{ori}\right)}^2\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)\left({x}_j^{gc}-{x}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left({y}_j^{gc}-{y}_j^{ori}\right)}^2.\end{array}} $$
(45)

This mean constraint could be simplified as

$$ {\displaystyle \begin{array}{l}{\overline{\boldsymbol{u}}}_j={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j{\overline{x}}_j+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\overline{y}}_j+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\overline{x^2}}_j+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j{\overline{xy^2}}_j+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\overline{y^2}}_j\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right){\overline{x}}_j+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right){\overline{x}}_j+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right){\overline{y}}_j\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right){\overline{y}}_j\kern1em +{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left({x}_j^{gc}-{x}_j^{ori}\right)}^2\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)\left({x}_j^{gc}-{x}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left({y}_j^{gc}-{y}_j^{ori}\right)}^2,\end{array}} $$
(46)

which is utilized for calculating the point value \( {\boldsymbol{u}}_j^{ori} \), and after that, the mean constraint is always satisfied. If we take Tj to represent the sum of some terms

$$ {\displaystyle \begin{array}{l}{\boldsymbol{T}}_j={\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right){\overline{x}}_j+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right){\overline{x}}_j+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right){\overline{y}}_j\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right){\overline{y}}_j\kern1em +{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_j^{gc}-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left({x}_j^{gc}-{x}_j^{ori}\right)}^2\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_j^{gc}-{x}_j^{ori}\right)\left({x}_j^{gc}-{x}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left({y}_j^{gc}-{y}_j^{ori}\right)}^2,\end{array}} $$
(47)

we have

$$ {\overline{\boldsymbol{u}}}_j={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j{\overline{x}}_j+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\overline{y}}_j+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\overline{x^2}}_j+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j{\overline{xy^2}}_j+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\overline{y^2}}_j+{\boldsymbol{T}}_j. $$
(48)

To achieve kth-order accuracy requires that we compute the kth derivatives by minimizing the error in predicting the mean value of the reconstructed function for the stencil k, that is, by minimizing the difference between the actual cell average \( {\overline{\boldsymbol{u}}}_k \) and the average of \( {R}_j\left(x-{\boldsymbol{x}}_j^{ori}\right) \) over stencil cell k. The mean value for a single cell Vk of the reconstructed function Rj is

$$ {\displaystyle \begin{array}{l}\underset{V_k}{\int }{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\underset{V_k}{\int}\left(x-{x}_j^{ori}\right)\kern0.5em dV+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\underset{V_k}{\int}\left(y-{y}_j^{ori}\right)\kern0.5em dV\\ {}+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\underset{V_k}{\int }{\left(x-{x}_j^{ori}\right)}^2\kern0.5em dV+{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\underset{V_k}{\int}\left(x-{x}_j^{ori}\right)\left(y-{y}_j^{ori}\right)\kern0.5em dV+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\underset{V_k}{\int }{\left(y-{y}_j^{ori}\right)}^2\kern0.5em dV.\end{array}} $$
(49)

Likewise, \( \left(x-{x}_j^{ori}\right) \) and \( \left(y-{y}_j^{ori}\right) \) could be written as \( \left(x-{x}_k^{gc}\right)+\left({x}_k^{gc}-{x}_j^{ori}\right) \) and \( \left(y-{y}_k^{gc}\right)+\left({y}_k^{gc}-{y}_j^{ori}\right) \) respectively,

$$ {\displaystyle \begin{array}{l}\underset{V_k}{\int }{R}_j\left(\boldsymbol{x}-{\boldsymbol{x}}_j^{ori}\right)\kern0.5em dV={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({\overline{x}}_k+\left({x}_k^{gc}-{x}_j^{ori}\right)\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j{\int}_{V_k}\left({\overline{y}}_k+\left({y}_k^{gc}-{y}_j^{ori}\right)\right)\kern0.5em dV\\ {}\kern1em +\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j\left({\overline{x^2}}_k+2{\overline{x}}_k\left({x}_k^{gc}-{x}_j^{ori}\right)+{\left({x}_k^{gc}-{x}_j^{ori}\right)}^2\right)\\ {}\kern1em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({\overline{xy}}_k+{\overline{x}}_k\left({y}_k^{gc}-{y}_j^{ori}\right)+{\overline{y}}_k\left({x}_k^{gc}-{x}_j^{ori}\right)+\left({x}_k^{gc}-{x}_j^{ori}\right)\left({y}_k^{gc}-{y}_j^{ori}\right)\right)\\ {}\kern1em +\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j\left({\overline{y^2}}_k+2{\overline{y}}_k\left({y}_k^{gc}-{y}_j^{ori}\right)+{\left({y}_k^{gc}-{y}_j^{ori}\right)}^2\right).\end{array}} $$
(50)

Combining the mean constraint given in Eq. (48) and writing Eq. (50) for every stencil cell, we obtain

$$ \left(\begin{array}{ccccccc}1& {\overline{x}}_j& {\overline{y}}_j& {\overline{x^2}}_j& {\overline{xy}}_j& {\overline{y^2}}_j& \cdots \\ {}{\omega}_{j1}& {\omega}_{j1}{\hat{x}}_{j1}& {\omega}_{j1}{\hat{y}}_{j1}& {\omega}_{j1}{\hat{x^2}}_{j1}& {\omega}_{j1}{\hat{xy}}_{j1}& {\omega}_{j1}{\hat{y^2}}_{j1}& \cdots \\ {}{\omega}_{j2}& {\omega}_{j2}{\hat{x}}_{j2}& {\omega}_{j2}{\hat{y}}_{j2}& {\omega}_{j2}{\hat{x^2}}_{j2}& {\omega}_{j2}{\hat{xy}}_{j2}& {\omega}_{j2}{\hat{y^2}}_{j2}& \cdots \\ {}{\omega}_{j3}& {\omega}_{j3}{\hat{x}}_{j3}& {\omega}_{j3}{\hat{y}}_{j3}& {\omega}_{j3}{\hat{x^2}}_{j3}& {\omega}_{j3}{\hat{xy}}_{j3}& {\omega}_{j3}{\hat{y^2}}_{j3}& \cdots \\ {}\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ {}{\omega}_{jN}& {\omega}_{jN}{\hat{x}}_{jN}& {\omega}_{jN}{\hat{y}}_{jN}& {\omega}_{jN}{\hat{x^2}}_{jN}& {\omega}_{jN}{\hat{xy}}_{jN}& {\omega}_{jN}{\hat{y^2}}_{jN}& \cdots \end{array}\right)\left(\begin{array}{c}{\left.{\boldsymbol{u}}^{ori}\right|}_j\\ {}{\left.\partial \boldsymbol{u}/\partial x\right|}_j\\ {}{\left.\partial \boldsymbol{u}/\partial y\right|}_j\\ {}{\left.\frac{1}{2}{\partial}^2\boldsymbol{u}/\partial {x}^2\right|}_j\\ {}{\left.{\partial}^2\boldsymbol{u}/\partial x\partial y\right|}_j\\ {}{\left.\frac{1}{2}{\partial}^2\boldsymbol{u}/\partial {y}^2\right|}_j\\ {}\vdots \end{array}\right)=\left(\begin{array}{c}\overline{\boldsymbol{u}}-{\boldsymbol{T}}_j\\ {}{\omega}_{j1}{\overline{\boldsymbol{u}}}_1\\ {}{\omega}_{j2}{\overline{\boldsymbol{u}}}_2\\ {}{\omega}_{j3}{\overline{\boldsymbol{u}}}_3\\ {}\vdots \\ {}{\omega}_{jN}{\overline{\boldsymbol{u}}}_N\end{array}\right), $$
(51)

where, the first row is the mean constraint, and the unconstraint equations are solved at first. After that, the point value at any local origins should be computed by the mean constraint to ensure Eq. (48) is exactly satisfied,

$$ {\left.{\boldsymbol{u}}_j^{ori}={\overline{\boldsymbol{u}}}_j-\operatorname{}\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j{\overline{x}}_j-\operatorname{}\frac{\partial \boldsymbol{u}}{\partial y}{\left|{}_j{\overline{y}}_j-\frac{1}{2}\operatorname{}\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\overline{x^2}}_j-\operatorname{}\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}{\left|{}_j{\overline{x{y}^2}}_j-\frac{1}{2}\operatorname{}\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\overline{y^2}}_j-{\boldsymbol{T}}_j, $$
(52)

which is utilized to calculate the solution at per gauss point of cell face. Eventually, the k-exact reconstruction based on any local origins is derived, and in the next section, higher-order spatial discretization will be specifically described.

5.2 Higher-order finite volume spatial discretization

In Section 2.1, the second-order finite volume discretization has been introduced, and as shown in Fig. 1, we mentioned that only one gauss point, i.e., the face centroid, is required for the second-order accurate discretization, while there should be at least two gauss points to guarantee the flux integration is higher-order accurate. The required order of accuracy, gauss point location as well as corresponding weights are given in Table 9, where A and B represent two edge points of a certain cell face.

Table 9 Required order of accuracy and the gauss point location as well as the corresponding weights

As a result, for the third-order accurate discretization, two gauss points are employed, and the left and right state vectors at per gauss point could be interpolated by the point value \( {\boldsymbol{u}}_j^{ori} \) and \( {\boldsymbol{u}}_k^{ori} \), and both of them are obtained by the mean constraint of k-exact reconstruction.

$$ \left\{\begin{array}{l}{\boldsymbol{u}}_L={\boldsymbol{u}}_j^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_j\left({x}_G-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_j\left({y}_G-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_j{\left({x}_G-{x}_j^{ori}\right)}^2\\ {}\kern2em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_j\left({x}_G-{x}_j^{ori}\right)\left({y}_G-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_j{\left({y}_G-{y}_j^{ori}\right)}^2,\\ {}{\boldsymbol{u}}_R={\boldsymbol{u}}_k^{ori}+{\left.\frac{\partial \boldsymbol{u}}{\partial x}\right|}_k\left({x}_G-{x}_j^{ori}\right)+{\left.\frac{\partial \boldsymbol{u}}{\partial y}\right|}_k\left({y}_G-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {x}^2}\right|}_k{\left({x}_G-{x}_j^{ori}\right)}^2\\ {}\kern2em +{\left.\frac{\partial^2\boldsymbol{u}}{\partial x\partial y}\right|}_k\left({x}_G-{x}_j^{ori}\right)\left({y}_G-{y}_j^{ori}\right)+\frac{1}{2}{\left.\frac{\partial^2\boldsymbol{u}}{\partial {y}^2}\right|}_k{\left({y}_G-{y}_j^{ori}\right)}^2.\end{array}\right. $$
(53)

Accordingly, the numerical flux could also be obtained by the Roe flux [37].

5.3 Accuracy analysis

In the first two sections, the third-order k-exact reconstruction as well as the spatial discretization is presented in detail. In this section, a numerical example, which is governed by the Euler equations, based on the Method of Manufactured Solutions (MMS) [43,44,45,46] is employed for the accuracy analysis. The manufactured solution [28] is shown as follows,

$$ \left\{\begin{array}{l}\rho =1.12+0.15\sin \left[\pi \left(3.12x+1895.92y\right)\right],\\ {}u=1.32+0.06\sin \left[\pi \left(2.09x+2099.21y\right)\right],\\ {}v=1.18+0.03\sin \left[\pi \left(2.15x+2001.32y\right)\right],\\ {}p=1.62+0.31\sin \left[\pi \left(3.79x+1973.98y\right)\right],\end{array}\right. $$
(54)

Taking the analytic solution to Euler equations, the source term is obtained

$$ \boldsymbol{s}=\nabla \cdot \boldsymbol{F}=\left(\begin{array}{c}u\frac{\partial \rho }{\partial x}+\rho \frac{\partial u}{\partial x}+v\frac{\partial \rho }{\partial y}+\rho \frac{\partial v}{\partial y}\\ {}2\rho u\frac{\partial u}{\partial x}+{u}^2\frac{\partial \rho }{\partial x}+\rho \left(v\frac{\partial u}{\partial y}+u\frac{\partial v}{\partial y}\right)+ uv\frac{\partial \rho }{\partial y}+\frac{\partial p}{\partial x}\\ {}2\rho v\frac{\partial v}{\partial y}+{v}^2\frac{\partial \rho }{\partial y}+\rho \left(v\frac{\partial u}{\partial x}+u\frac{\partial v}{\partial x}\right)+ uv\frac{\partial \rho }{\partial x}+\frac{\partial p}{\partial y}\\ {}u\frac{\partial \left(\rho H\right)}{\partial x}+\rho H\frac{\partial u}{\partial x}+v\frac{\partial \left(\rho H\right)}{\partial y}+\rho H\frac{\partial v}{\partial y}\end{array}\right), $$
(55)

where H = (γp/(γ − 1) + ρv2/2)/ρ is the specific total enthalpy. For higher-order unstructured finite volume discretization, the source term should be integrated over the control volume rather than being evaluated at the reference point to avoid the second-order error. Three-point quadrature [34, 35] is utilized again for the source term integration.

Besides, regular and randomly perturbed triangular grids with 5 × 102 and 103 these two aspect ratios are employed, and in each cell aspect ratio, five sets of triangular grids from the coarsest to finest are utilized to count the discretization errors and computational accuracy. The coarsest regular and randomly perturbed triangular grids and the distribution of background quadrilateral grids in x and y directions are shown in Fig. 37 and Table 10.

Fig. 37
figure 37

(a) Regular and (b) randomly perturbed triangular grids

Table 10 Distribution of background quadrilateral grids in x and y directions

Computational results on both regular and randomly perturbed triangular grids are shown as follows.

5.3.1 Computational results on regular triangular grids

As shown in Figs. 38 and 39, the third-order accuracy is achieved by employing the global-direction stencil with face-area-weighted centroid. Therefore, correctness as well as feasibilities of k-exact reconstruction and spatial discretization based on any local origins is demonstrated. Besides, combining the result given in Table 11, we can easily find that discretization errors of global-direction stencil with face-area-weighted centroid are always the lowest among all methods we tested, no matter on grids with AR = 5 × 102 or 103.

Fig. 38
figure 38

(a) L2 and (b) L errors. Errors of different stencils on regular triangular grids with AR = 5 × 102

Fig. 39
figure 39

(a) L2 and (b) L errors. Errors of different stencils on regular triangular grids with AR = 103

Table 11 Computational errors on the finest grid and L2, L accuracy between the last two grids of different stencils

But note that errors are also counted at the local origin rather than the cell-averaged result, and for higher-order accurate finite volume method, the comparison between the cell-averaged result and that counted at the local origin is listed in Section 5.3.3.

5.3.2 Computational results on randomly perturbed triangular grids

From Figs. 40 and 41, we find that on randomly perturbed triangular grids, there is a similar trend to results on regular grids, that is, the global-direction stencil with face-area-weighted centroid is lower than another three stencils. As a result, computational accuracy is improved, which is beneficial for obtaining a more accurate flowfield. In short, feasibilities and superiorities of extending the face-area-weighted centroid and the combination of global-direction stencil with this novel local origin to higher-order unstructured finite volume method from integral form are verified.

Fig. 40
figure 40

(a) L2 and (b) L errors. Errors of different stencils on randomly perturbed triangular grids with AR = 5 × 102

Fig. 41
figure 41

(a) L2 and (b) L errors. Errors of different stencils on randomly perturbed triangular grids with AR = 103

5.3.3 Discussion about results exhibition of higher-order discretization from integral form

In Section 4.1.3, it is verified by analysis and numerical example that for the second-order unstructured finite volume discretization, the point value at the face-area-weighted centroid is more appropriate for the result exhibition, since errors at this certain local origin are lower than the cell average. In this section, cell averages are also utilized to further contrast with the mentioned point values at the local origin.

According to the three requirements mentioned in Section 4.1.3, firstly, the point value at any local origins \( {\boldsymbol{u}}_j^{ori} \) could be obtained by the mean constraint. In addition, as shown in Sections 5.3.1 and 5.3.2, the third-order accuracy is always achieved by the point-value solution at both geometric centroid and the face-area-weighted centroid, and more importantly, as shown in Fig. 42, errors counted at the face-area-weighted centroid are obviously lower than that counted at the geometric centroid as well as the cell-averaged value.

Fig. 42
figure 42

(a) L2 and (b) L errors. Cell-averaged and point-valued discretization errors of G-Stencil with the geometric centroid and face-area-weighted centroid on grids with AR = 103

In summary, for higher-order accurate unstructured finite volume method from integral form, when the global-direction stencil is utilized, we find the result counted at the face-area-weighted centroid is more accurate than the cell-averaged solution. Hence, feasibilities as well as the superiorities of the face-area-weighted centroid are demonstrated not only on the second-order finite volume method, but on discretization with higher-order accuracy.

6 Concluding remarks

Inspired by the novel face-area-weighted centroid utilized for unstructured finite volume discretization from differential form, in this paper, we have shown that the reference point used in the integral form does not have to be the geometric centroid as well, and the face-area-weighted centroid is a better choice, which effectively reduces the grid skewness. Besides, we combine the global-direction stencil with this novel reference point trying to improve the computational accuracy and convergence rate of unstructured finite volume solver.

Specifically, we first illustrate in detail that the traditional LSQR reconstruction method is unable to obtain the point-value solution located at any local origins, and that is the reason why we focus on k-exact reconstruction algorithm. On this basis, the k-exact reconstruction based on any local origins is analytically derived, and during the reconstruction process, the mean constraint can always be satisfied. Besides, we extend the global-direction stencil with novel face-area-weighted centroid from unstructured finite volume method in differential form to the integral form for the grid skewness reduction and a better reflection of flow anisotropy. More importantly, it sets stage for the promotion on higher-order accuracy.

Four representative numerical examples governed by linear convective, Euler and Laplacian equations are simulated to examine the correctness of k-exact derivation and the effectiveness of global-direction stencil with face-area-weighted centroid on the unstructured finite volume discretization from integral form. After verification, in numerical cases governed by linear convective and Euler equations, when the face-area-weighted centroid is utilized on finite volume discretization, the second-order accuracy also could be achieved, and the global-direction stencil with this novel reference point has the lowest discretization errors among all stencils we tested, where, the result is exhibited by the point value rather than the cell average, since we realize when the global-direction stencil is utilized, the point-value solution \( {\boldsymbol{u}}_j^{ori} \) obtained by the mean constraint of k-exact reconstruction at the face-area-weighted centroid is more accurate than the cell average, and this conclusion is also demonstrated on the third-order finite volume discretization. On this basis, the result exhibition of unstructured finite volume method could be replaced by this point value at the face-area-weighted centroid.

Besides, in subsonic flow over NACA0012 airfoil, pressure values utilized to compute the lift and drag coefficients at the gauss point of wall-adjacent face are also interpolated by the point-valued solution at the local origin, and we find this novel method has a better numerical performance on calculating lift and drag coefficients, and has a faster convergence speed than commonly used stencils based on the geometric centroid. On this basis, a Laplacian model equation is utilized to further test its effectiveness on dissipative term evaluation. From the result, similar conclusions are demonstrated that errors of the global-direction stencil with face-area-weighted centroid are still lower than another three stencils, and the convergence rate is also greatly promoted by this novel method. As a result, this conclusion sets stage for the further extension to viscous flows.

Finally, the accuracy analysis is implemented on the third-order unstructured finite volume discretization from integral form as well, where the third-order k-exact reconstruction as well as the spatial discretization based on any local origins is derived. After verification, it is illustrated that the third-order accuracy is also achieved by the global-direction stencil with face-area-weighted centroid, and discretization errors of this method are lower than another three stencils we tested.

In conclusion, the computational accuracy of the second and third-order unstructured finite volume method from integral form will not get deteriorated by the employment of face-area-weighted centroid, and it is feasible to extend the global-direction stencil with face-area-weighted centroid from differential finite volume solver to the integral form. Future work will be carried out from two aspects. Firstly, we will not only further compare the computational accuracy of point-valued solution at the certain local origin with the abstract cell-averaged solution, but figure out the specific approach to display the flowfield by this point value. In addition, extensions of the global-direction stencil selection method and face-area-weighted formula to three dimensions are also quite necessary for practical applications.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.

Abbreviations

CFD:

Computational fluid dynamics

2D:

Two-dimensional

LSQR:

Least-squares

MMS:

Method of manufactured solutions

V-Stencil:

Vertex-neighbor stencil

F-Stencil:

Face-neighbor stencil

G-Stencil (Cell centroid):

Global-direction stencil with cell centroid

G-Stencil (F-a-w centroid):

Global-direction stencil with face-area-weighted centroid

EN:

Edge-normal scheme

FT:

Face-tangent scheme

References

  1. Mavriplis DJ (1997) Unstructured grid techniques. Annu Rev Fluid Mechc 29(1):473–547

    Article  MathSciNet  Google Scholar 

  2. Wang ZJ, Zhang LP, Liu Y (2004) Spectral (finite) volume method for conservation laws on unstructured grids IV: extension to two-dimensional systems. J Comput Phys 194(2):716–741

    Article  MathSciNet  Google Scholar 

  3. Wang NH, Li M, Ma R, Zhang LP (2019) Accuracy analysis of gradient reconstruction on isotropic unstructured meshes and its effects on inviscid flow simulation. Adv Aerodyn 1(18). https://doi.org/10.1186/s42774-019-0020-9

  4. Venkatakrishnan V (1996) Perspective on unstructured grid flow solvers. AIAA J 34(3):533–547

    Article  Google Scholar 

  5. Fidkowski KJ, Darmofal DL (2011) Review of output-based error estimation and mesh adaptation in computational fluid dynamics. AIAA J 49(4):673–694

    Article  Google Scholar 

  6. Snyder DO, Koutsavdis EK, Anttonen JSR (2003) Transonic store separation using unstructured CFD with dynamic meshing. In: 33rd AIAA fluid dynamics conference and exhibit, AIAA 2003-3919, Orlando, 23–26 June 2003

  7. Murayama M, Nakahashi K, Matsushima K (2002) Unstructured dynamic mesh for large movement and deformation. In: 40th AIAA aerospace science meeting & exhibit, AIAA-2002-0122, Reno, 14–17 Jan 2002

  8. Weller HG, Tabor G, Jasak H, Fureby C (1998) A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput Phys 12(6):620–631

  9. Katz A, Sankaran V (2011) Mesh quality effects on the accuracy of CFD solutions on unstructured meshes. J Comput Phys 230(20):7670–7686

    Article  MathSciNet  Google Scholar 

  10. Katz A, Sankaran V (2011) Discretization methodology for high aspect ratio prismatic grids. In: 20th AIAA computational fluid dynamics conference, AIAA 2011-3378, Honolulu, 27–30 June 2011

  11. Katz A, Sankaran V (2012) High aspect ratio grid effects on the accuracy of Navier–Stokes solutions on unstructured meshes. Comput Fluids 65:66–79

    Article  MathSciNet  Google Scholar 

  12. Diskin B, Thomas JL (2008) Accuracy of gradient reconstruction on grids with high aspect ratio. NIA report, No. 2008-12

  13. Diskin B, Thomas JL (2011) Comparison of node-centered and cell-centered unstructured finite-volume discretizations: inviscid fluxes. AIAA J 49(4):836–854

    Article  Google Scholar 

  14. Diskin B, Thomas JL, Nielsen EJ, Nishikawa H, White JA (2010) Comparison of node-centered and cell-centered unstructured finite-volume discretizations: viscous fluxes. AIAA J 48(7):1326–1338

  15. Diskin B, Thomas JL (2012) Effects of mesh regularity on accuracy of finite-volume schemes. In: 50th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, AIAA 2012-0609, Nashville, 09–12 Jan 2012

  16. Jalali A, Sharbatdar M, Ollivier-Gooch C (2014) Accuracy analysis of unstructured finite volume discretization schemes for diffusive fluxes. Comput Fluids 101:220–232

  17. Sozer E, Brehm C, Kiris CC (2014) Gradient calculation methods on arbitrary polyhedral unstructured meshes for cell-centered CFD solvers. In: 52nd aerospace sciences meeting, AIAA 2014-1440, National Harbor, 13–17 Jan 2014

  18. Schwöppe A, Diskin B (2013) Accuracy of the cell-centered grid metric in the DLR TAU-code. In: New results in numerical and experimental fluid mechanics VIII. Springer-Verlag, Berlin, Heidelberg, pp 429–437

  19. Nishikawa H (2019) Efficient gradient stencils for robust implicit finite-volume solver convergence on distorted grids. J Comput Phys 386:486–501

    Article  MathSciNet  Google Scholar 

  20. Xiong M, Deng XG, Gao X et al (2018) A novel stencil selection method for the gradient reconstruction on unstructured grid based on OpenFOAM. Comput Fluids 172:426–442

  21. Xiong M (2019) Research and application of unstructured finite volume gradient reconstruction algorithms. Dissertation, National University of Defense Technology

  22. Kong LF, Dong YD, Liu W (2020) The influence of global-direction stencil on gradient and high-order derivatives of unstructured finite volume methods. Chin J Theor Appl Mech 52(5):1334–1349

  23. Kong LF, Dong YD, Liu W, Zhang HB (2020) An improved global-direction stencil based on the “face-area-weighted centroid” for the gradient reconstruction of unstructured finite volume methods. Chin Phys B 29(10):100203. https://doi.org/10.1088/1674-1056/aba2da

  24. You D, Mittal R, Wang M, Moin P (2006) Analysis of stability and accuracy of finite-difference schemes on a skewed mesh. J Comput Phys 213(1):184–204

    Article  Google Scholar 

  25. Veluri SP (2010) Code verification and numerical accuracy assessment for finite volume CFD codes. Dissertation, Virginia Polytechnic Institute and State University

  26. Dannenhoffer J (2012) Correlation of grid quality metrics and solution accuracy for supersonic flows. In: 50th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, AIAA 2012-0610, Nashville, 09–12 Jan 2012

  27. Kallinderis Y, Fotia S (2015) A priori mesh quality metrics for three-dimensional hybrid grids. J Comput Phys 280:465–488

    Article  Google Scholar 

  28. Nishikawa H (2020) A face-area-weighted ‘centroid’ formula for finite-volume method that improves skewness and convergence on triangular grids. J Comput Phys 401:109001

    Article  MathSciNet  Google Scholar 

  29. Masatsuka K (2013) I do like CFD, vol 1, 2nd edn, p 28 http://www.cfdbooks.com

    Google Scholar 

  30. Barth TJ, Jespersen DC (1989) The design and application of upwind schemes on unstructured meshes. In: 27th aerospace sciences meeting, AIAA 1989-0366, Reno, 09–12 Jan 1989

  31. Barth TJ, Frederickson PO (1990) Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. In: 28th aerospace sciences meeting, AIAA-90-0013, Reno, 08–11 Jan 1990

  32. Barth T (1991) A 3-D upwind Euler solver for unstructured meshes. In: 10th computational fluid dynamics conference, AIAA-91-1548, Honolulu, 24–26 June 1991

  33. Barth TJ (1993) Recent developments in high order k-exact reconstruction on unstructured meshes. In: 31st aerospace sciences meeting, AIAA-93-0668, Reno, 11–14 Jan 1993

  34. Ollivier-Gooch CF, Nejat A, Michalak K (2007) On obtaining high-order finite-volume solutions to the Euler equations on unstructured meshes. In: 18th AIAA computational fluid dynamics conference, AIAA 2007-4464, Miami, 25–28 June 2007

  35. Ollivier-Gooch C, Nejat A, Michalak K (2009) Obtaining and verifying high-order unstructured finite volume solutions to the Euler equations. AIAA J 47(9):2105–2120

  36. Li WA, Ren YX (2012) High-order k-exact WENO finite volume schemes for solving gas dynamic Euler equations on unstructured grids. Int J Numer Meth Fl 70(6):742–763

    Article  MathSciNet  Google Scholar 

  37. Roe PL (1986) Characteristic-based schemes for the Euler equations. Annu Rev Fluid Mech 18(1):337–365

    Article  MathSciNet  Google Scholar 

  38. Shima E, Kitamura K, Haga T (2013) Green–gauss/weighted-least-squares hybrid gradient reconstruction for arbitrary polyhedra unstructured grids. AIAA J 51(11):2740–2747

    Article  Google Scholar 

  39. Deka M, Brahmachary S, Thirumalaisamy R (2018) A new green–gauss reconstruction on unstructured meshes. Part I: gradient reconstruction. J Comput Phys 1:108325

    MathSciNet  Google Scholar 

  40. Nishikawa H (2018) From hyperbolic diffusion scheme to gradient method: implicit green–gauss gradients for unstructured grids. J Comput Phys 372:126–160

    Article  MathSciNet  Google Scholar 

  41. Jameson A, Mavriplis D (1986) Finite volume solution of the two-dimensional Euler equations on a regular triangular mesh. AIAA J 24(4):611–618

    Article  Google Scholar 

  42. Jalali A, Ollivier-Gooch C (2017) Higher-order unstructured finite volume RANS solution of turbulent compressible flows. Comput Fluids 143:32–47

    Article  MathSciNet  Google Scholar 

  43. Roy CJ (2005) Review of code and solution verification procedures for computational simulation. J Comput Phys 205(1):131–156

    Article  Google Scholar 

  44. Roache PJ (2002) Code verification by the method of manufactured solutions. J Fluid Eng 124(1):4–10

    Article  Google Scholar 

  45. Roy CJ, Smith TM, Ober CC (2002) Verification of a compressible CFD code using the method of manufactured solutions. In: 32nd AIAA fluid dynamics conference and exhibit, AIAA 2002-3110, St. Louis, 24–26 June 2002

  46. Park MA, Barral N, Ibanez D et al (2018) Unstructured grid adaptation and solver technology for turbulent flows. In: 2018 AIAA aerospace sciences meeting, AIAA 2018-1103, Kissimmee, 8–12 Jan 2018

  47. Krivodonova L, Berger M (2006) High-order accurate implementation of solid wall boundary conditions in curved geometries. J Comput Phys 211(2):492–512

    Article  MathSciNet  Google Scholar 

  48. Jalali A, Ollivier-Gooch C (2013) Accuracy assessment of finite volume discretizations of convective fluxes on unstructured meshes. In: 51st AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, AIAA 2013-0705, Grapevine, 07–10 Jan 2013

  49. Wang Q, Ren YX, Pan J, Li WA (2017) Compact high order finite volume method on unstructured grids III: Variational reconstruction. J Comput Phys 337:1–26

    Article  MathSciNet  Google Scholar 

  50. Thomas JL, Diskin B, Nishikawa H (2011) A critical study of agglomerated multigrid methods for diffusion on highly-stretched grids. Comput Fluids 41(1):82–93

    Article  Google Scholar 

  51. Nishikawa H, Nakashima Y, Watanabe N (2017) Effects of high-frequency damping on iterative convergence of implicit viscous solver. J Comput Phys 348:66–81

    Article  MathSciNet  Google Scholar 

  52. Nishikawa H (2010) Beyond interface gradient: a general principle for constructing diffusion schemes. In: 40th fluid dynamics conference and exhibit, AIAA 2010-5093, Chicago, 28 June–01 July 2010

  53. Nishikawa H (2011) Robust and accurate viscous discretization via upwind scheme–I: basic principle. Comput Fluids 49(1):62–86

    Article  MathSciNet  Google Scholar 

  54. Nishikawa H (2011) Two ways to extend diffusion schemes to Navier-Stokes schemes: gradient formula or upwind flux. In: 20th AIAA computational fluid dynamics conference, AIAA 2011-3044, Honolulu, 27–30 June 2011

Download references

Acknowledgments

Not applicable.

Funding

This work is supported by National Key Project (No. GJXM92579).

Author information

Authors and Affiliations

Authors

Contributions

The contribution of the authors to this work is equivalent. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Yidao Dong.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kong, L., Dong, Y., Liu, W. et al. Extending the global-direction stencil with “face-area-weighted centroid” to unstructured finite volume discretization from integral form. Adv. Aerodyn. 2, 25 (2020). https://doi.org/10.1186/s42774-020-00048-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42774-020-00048-5

Keywords