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


 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, we explored a novel global-direction stencil and combine it with face-area-weighted centroid on unstructured finite volume methods from differential form to realize the skewness reduction and a better reflection of flow anisotropy. Note, however, that the differential form is hard to achieve higher-order accuracy, and in order to set stage for the method promotion on higher-order numerical simulation, in this research, we demonstrate that it is also feasible to extend this novel method to the unstructured finite volume discretization in integral form. 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 a better numerical performance.

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 [Chinese Physics B. 2020, In press], 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 7 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.

Unstructured finite volume method
In this section, we give unstructured finite volume discretizations from both differential and integral forms at first. In addition, the k-exact reconstruction method based on any local origins is also derived in this section.

Finite volume discretization from differential form
Usually, the differential inviscid governing equation could be formulated as, where, j u is a solution vector. c j F is the convective flux and j s is a source (or forcing) term vector. In differential unstructured finite volume method, both solution and source term vectors are evaluated as point values at the geometric centroid, rather than cell-averaged values. Therefore, the source term integration is totally avoided.
Approximating the flux divergence by the integral average, we obtain Therefore, the governing equation has already committed the second-order errors because of the linear exactness, and it could not be eliminated directly.
Note, however, that higher-order accurate discretization from the differential form is also possible to realize. In Ref. [29], where Masatsuka provided an idea by reconstructing the primitive function of flux to guarantee the governing equation is analytically accurate. The main process is shown as follows.
The original inviscid governing equation from differential form is, 0, c tj  +   = uf (4) where, flux c f could be formulated by the primitive function c F as, (5) Taking this equation to Eq. (4), we obtain 1 0.
This equation is accurate, and the feasibility of extending the differential form to higher-order accuracy is demonstrated.
As mentioned in introduction, the 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 9 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. The governing equation and basic discretization process of integral form is given in the following section.

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, According to the divergence theorem, Eq. (7) is transformed as, On this basis, Eq. (8) could be discretized as follows, {} 11 , (9) where,  As shown in Fig. 1, where u L and u R are two state vectors at the gauss point obtained by owner and neighbor cells of the common face. For the second-order unstructured finite volume solver, the gauss point is just located at the face centroid and the weight is equaling to 1, while for higher-order discretization, there should be at least two points to guarantee the flux quadrature is 10 higher-order accurate. The gauss point location as well as the corresponding weights will be specifically given in Section 5.2. Besides, r j and r k 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 j u and k u stored at the reference point 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-value solution at any local origins is derived, and the main process will be given in the following section. As a result, j u and k u these two point values could be written in terms of ori j u and ori k u , x y y xy (10) where point value solution of ori j u and ori k u , as well as the solution gradients are all derived and obtained by the k-exact reconstruction algorithm, which will be specifically described in Section 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.
As demonstrated in Eq. (12), 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 stored at cells: jj → uu , 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 orig j u is unknown at any local origins in every time step. Besides, least-square 12 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 stored 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 k-exact reconstruction [30][31][32][33][34][35][36], the geometric centroid is always used to formulate the reconstructed function where, gc j x and gc j y are coordinates of the geometric centroid in Cartesian-coordinate system.
Taking the integral average of this expansion over the control volume gives (14) where nm j xy is the integration of geometric quantities and is called Moments [34,35]. Besides, conservation of the mean requires that the average of the reconstructed function .
This is the precondition of the k-exact reconstruction algorithm, and if the expansion is based on the geometric centroid, the mean constraint could be written in terms of 13 . gc j j j j j j xy xy  = + +  uu uu (16) To meet this requirement, the point value stored at the geometric centroid gc j u need to be computed according to this constraint, while in order to obtain the point value stored at any local origins, the expansion should be formulated based on any local origins as well, where, j ori u is the point value stored at any local origins. Integrating Eq. (17) over the control volume, we obtain  (19) which could be further transformed as where the first row is the mean constraint, and geometric terms could be expressed as,  (25) and the weights are set to emphasize the geometrically adjacent data, where j x and k x are local origins of central and stencil cells. When solving the Eq. (24), the unconstraint reconstructed equations are solved at first and after obtaining the gradients, the mean constraint must be used to compute the constant term j ori u in the expansion as, 15 . (27) Thus, the mean constraint is always exactly satisfied and the obtained In short, by solving Eq. (24), solution vectors at any local origins as well as derivatives could be obtained. Therefore, the k-exact reconstruction method based on any local origins is derived, and it sets stage for the employment of novel reference points on unstructured finite volume discretization from integral form. Higher-order, i.e., the third-order derivation of k-exact reconstruction will be given in Section 5.

Global-direction stencil based on the face-area-weighted centroid
In Section 2.3, 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.

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 16 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.
(a) Vertex-neighbor stencil (b) Face-neighbor stencil 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.
(a) Minor-aspect-ratio (b) High-aspect-ratio 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.

Fig. 5
Global-direction stencil cells on minor and high-aspect-ratio triangular grids.
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 is 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. 18 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 stored 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.

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.

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 ˆj k e is a unit vector pointing from centroid of cell j to that of cell k, ˆj k n is the unit outward normal vector of common face. The grid skewness is measured by ˆĵ k jk  en, 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 19 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.
For high-aspect-ratio triangular grids, a novel face-area-weighted centroid was proposed by Nishiakwa [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, where, ( ) , mk mk xy is coordinate of the k-th face centroid. Besides, we generalize Eq. (28) to the face-area-weighted centroid formula [28] : where, jk A 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. 20 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. For these two common faces, if the geometric centroid is utilized, the grid skewness measure is [28]. 12 23 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 equaling to 2, the skewness measure becomes [28], 22 12 23 2 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 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.

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 stored at cells reference points, 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.3, 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 22 cells and local origins. Global-direction stencils with geometric and face-area-weighted centroids on grids with straight and curved boundaries are displayed in Fig. 9 and Fig. 10.
(a) Geometric centroid (b) Face-area-weighted centroid 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 23 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.

Numerical examples
To

Manufactured boundary layer (Governed by Linear convective equation)
In this section, we first use the Method of Manufactured Solutions (MMS) [43][44][45][46] 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 .
(a) μ = 10 -6 (b) μ = 10 -8 On this basis, the manufactured solution represents the analytical solution of this modified governing equation, and we can calculate both L 2 and L ∞ errors of different stencils.
( ) where ˆy h is the first layer vertical spacing, and  is the stretching factor that could be computed by the known condition 3 10 Besides, distribution of five sets of background quadrilateral grids from the coarsest to finest is shown in Table 2. On this basis, computational accuracy and discretization errors of four different stencils on integral unstructured finite volume solver are counted and given in Fig. 13 and Fig. 14.  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.

Computational results on regular grids
Besides, combining the specific data listed in Table 3, we find both L 2 and L ∞ errors of G-Stencil (F-a-w centroid) are the lowest among all stencils we tested, and when AR = 10 4 , 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.

Computational results on randomly perturbed grids
On randomly perturbed triangular grids, we also test the case of AR = 10 3 and 10 4 . For simplicity, results of AR = 10 4 are given here.

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.  (37) and   Combining the results shown in 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-averaged errors are almost identical and both of them are close to errors counted at the geometric centroid.
In addition, both cell-averaged errors and errors counted at the geometric centroid are always higher than that counted 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 ori j u , which could be obtained by the k-exact reconstruction. 31 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 Therefore, combining the analysis and numerical verification, we notice that compared with the traditional finite volume scheme from integral form, where the result is counted by the cell-averaged value, 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 stored at the face-area-weighted centroid is more 32 appropriate to reflect the computational result. Hence, in the following test, errors are all counted at the local origin.

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 r i = 1 and r 0 = 1.384 located in the first quadrant. These two circular arcs represent the inviscid wall boundary, and the flow at both inlet and outlet are supersonic.
Analytical solution [47] of this numerical example could be derived by isentropic relation and is given as follows, (38) where the value of Mach number at the inner radius is 2.25 i M = and the density 1 Besides, the sound speed is calculated as, Flow structure of this numerical case is shown as follows. 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.  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.   As shown in Fig. 20 and Fig. 21, for three stencils with the geometric centroid, both L 2 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 L 2 and L ∞ errors is achieved by G-Stencil (F-a-w centroid). Particularly, L 2 errors can reach 2.418 and 2.472 order accurate. As a 35 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.

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

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     Conclusions related to lift an drag coefficients mentioned above can be better illustured by residuals of different stencils. As exhibited in Fig. 28, we can eaisly find that for three stencil with geometric centroid, the convergence speed of G-Stencil is colose to the commonly used V-Stencil, 39 and in 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

Dissipative term evaluation (Governed by Laplacian equation)
From But, the numerical performance on dissipative term is also need to be further tested to set stage for the extension to viscous problems. In this section, the MMS method is also used on 41 Laplacian equation to examine discretization errors and convergence speed of different stencils.
Laplacian model equation could be formulated as, where  is a constant equaling to 1. Integrating this model equation over the control volume, we obtain ( ) According to the divergence theorem, Eq. (41) could be transformed as From Eq. (42), we can clearly find the difference between convective term and dissipative term.
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 42 where,  is a damping coefficient, and here, a special value is 43  = , which has been known to provide superior accuracy and robustness [28,48,51]. ˆj k n is the unit face normal vector, and ˆj k e is the unit vector between two cell centroids. L jk is the length from one centroid to another.
Besides, L   and R   are solutions linearly reconstructed at integral point, On the basis of determining the discretization method of dissipative term, the manufactured  From Eq. (44), 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. 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
Geometric and face-area-weighted centroids on grids with regular and 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 Fig 34,   Fig. 35 and Table 8.   From Fig. 34, Fig. 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.
45 On this basis, residuals on the coarsest and finest regular grids of different stencils are shown in Fig. 36.  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. 46

Accuracy analysis of the method on the third-order discretization
In Section 2.3, the k-exact reconstruction based on any local origins is derived, and from the Section 4.1 to Section 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.

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,  (49) where ( )   For k-exact reconstruction, no matter the reconstructed function is formulated based on which reference point, the mean constraint requires [34,35], Besides, integrating the reconstructed polynomial over the control volume, we have This mean constraint could be simplified as    (53) which is utilized for calculating the point value ori j u , and after that, the mean constraint is always satisfied. If we take j T to represent the sum of some terms 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 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.

Higher-order finite volume spatial discretization
In section 2.2, the second-order finite volume discretization has been introduced, and as show 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.
Accordingly, the numerical flux could also be obtained by the Roe scheme [37].

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, Hp     = − + v 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] are utilized again for the source term integration. 51 Besides, regular and randomly perturbed triangular grids with 5×10 2 and 10 3 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.  Computational results on both regular and randomly perturbed triangular grids are shown as follows.   As shown in Fig. 38 and Fig. 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 are 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×10 2 or 10 3 . 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. 53   From Fig. 40 and Fig. 41, we find that on randomly perturb 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 54 this novel local origin to higher-order unstructured finite volume method from integral form are verified.

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-averaged errors are also counted to further contrast with the mentioned point value errors counted at the local origin.   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.
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.

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. 56 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 discretization errors are counted at the local origin rather than the integral average, since we realize when the global-direction stencil is utilized, point-value solution ori j u 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.
Besides, in subsonic flow over NACA0012 airfoil, 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 57 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. The next work will be carried out from two aspects. Firstly, we will study the higher-order finite volume discretization from differential form to further extend the face-area-weighted centroid and its combination with the global-direction stencil. In addition, applications on viscous and high-mach-number flows with complex surface are also necessary to examine its numerical performance. Face-tangent scheme.