0 Introduction

The finite volume method (FVM) has gained a great popularity in computational fluid dynamics (CFD) because of its better conservativeness and flexibility to adapt to both structured and unstructured grids. However,conventional FVM requires wide stencil to generate high-order reconstructions,and the choice of the cell stencil is not straightforward. As mentioned in ,particular attention must be paid to choose the admissible stencils and there is not a simple rule to guide this procedure for reconstructions higher than first order. As a matter of fact,most operational codes are limited to the second order on unstructured grids,where a linear reconstruction such as the MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) scheme is used.

On the other hand,the majority of numerical methods for incompressible flows are essentially based on the pressure-correction (or pressure-projection) approach,even they appear in the literature as different variants,such as the projection method[4],MAC(Marker and Cell)[5] method,SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) method[6],SIMPLEC(SIMPLE Corrected) method [7] and PISO(Pressure Implicit with Splitting Operators) method[8]. All these methods work very well under the second-order FVM framework where the velocity can be directly coupled with the pressure through a pressure-projection step when staggering or co-volume arrangement is adopted.

Although the second-order FVM on unstructured meshes proves itself a good trade-off between computational complexity and numerical accuracy,and has been widely accepted as a practical CFD framework for real-case applications,some remaining problems,for example the dependency of solution quality on computational mesh and the poor accuracy in advection computation,still deserve further efforts for more accurate and robust formulations. Some high-order DG methods have been devised for incompressible unsteady Navier-Stokes equations [9, 10, 11, 12]. A common practice is to use the fractional-step pressure projection for the temporal updating and use DG for the spatial discretization,which is usually more complex and computationally expensive compared to the conventional FVM.

In this paper,we propose a novel numerical formulation to solve incompressible unsteady Navier-Stokes equations by including the point value (PV) of the velocity field at the cell vertices as a new DOF in addition to the volume integrated average (VIA) that is the only computational variable in the conventional FVM. The present method is called Volume integrated average and Point value based Multi-moment (VPM) method,in which the VIA moment is computed by a finite volume formulation of flux form,and thus exactly conservative,while the PV moment is updated by point-wise Riemann solver that can be computed very efficiently in unstructured grids. The reconstruction that is needed to evaluate the numerical fluxes and the derivatives in the spatial discretization is determined from the constraint conditions in terms of both VIA and PV as the computational variables. More precisely,we construct piecewise incomplete quadratic polynomials of 6 DOFs in 2 and 8 DOFs in 3 dimensions. In addition to one VIA and three (2D) or four (3D) PVs of a single grid cell,the derivatives at the cell center are also used to determine the polynomials. The derivatives are approximated from the VIA and PVs of the cells sharing the boundary edges with the target cell. The stencil is compact,and more importantly there is no arbitrariness in choosing the stencil for reconstruction. Compared to the conventional finite volume method,significant improvements in accuracy and robustness can be achieved by the present method without large increase in algorithmic complexity and computational cost.

The present formulation can be seen as an extension of the CIP multi-moment finite volume methods [13, 14, 15, 16, 17, 18, 19, 20, 21, 22] to incompressible Navier-Stokes equations on unstructured grids with triangular and tetrahedral elements. The multi-moment finite volume methods have been developed on structured grids as accurate and robust fluid solvers and applied so far to various applications. The underlying idea to make use of multiple types of moments facilitates novel numerical models of greater efficiency and flexibility and is consequently much beneficial for implementations in unstructured grids.

1 Numerical Formulation on Unstructured Mesh

We consider the incompressible unsteady Navier-Stokes equations,

where u =(u,v,w) is the velocity vector with components u,v and w in x,y and z directions respectively. p is the pressure,ρ the kinematic viscosity.

The projection method of Chorin [4] provides a simple and robust solution procedure for incompressible flow and is adopted in this paper. We briefly summarize the numerical steps to update the velocity from time level m(t=t^m) to m+1(t=t^m+Δt) as follows.

1) Given the velocity field u^m at step m,compute the intermediate velocity u^* from the momentum Eq.(2) without the pressure gradient term,

2) The intermediate velocity u^* usually does not satisfy the mass conservation Eq.(1),and must be corrected by the projection shown in the next step. For that purpose,we solve the pressure field at step (m+1) from the following Poisson equation,

3) Correct the velocity by the projection step,

It is straight forward to show that the updated velocity satisfies ▽· u^m+1=0.

Our major interest in this paper is to develop a novel approach for the spatial discretization for the solution procedure shown above.

1.1 Computational Mesh

We firstly describe the numerical formulation in 2D by discarding the component in z direction. The computational domain is divided into non-overlapping triangular cells Ω_i(i=1,…,N_e) with three vertices θ_ij(j=1,2,3) located at (x_i1,y_i1),(x_i2,y_i2) and (x_i3,y_i3) respectively as shown in Fig. 1(a). We also denote the mass center of Ω_i by θ_ic.

The three boundary line segments of element Ω_i are denoted by For the later use,we denote the middle points and the outward normal unit vector of segment Γ_ij by and n_ij=(n_xij,n_yij). Shown in Fig. 1(b),the volume-integrated average (VIA) and point-value (PV) at the vertices of a physical variable ϕ(x,y,t) are defined as:

where |Ω_i| denotes the volume of cell Ω_i.

Given both VIA and PV,a polynomial can be constructed in terms of VIA and PV. For simplicity,mesh cell Ω_i is mapped onto a standard element

We give some numerical formula to transform the first-order derivatives between the global coordinate system (x,y) and the local coordinate (ξ,η) and where |J_i(ξ,η)|=x_ξi(ξ,η)y_ηi(ξ,η)-x_ηi(ξ,η)y_ξi(ξ,η) is the determinant of the Jacobian matrix.

It is straightforward from (7) that

and 1.2 Multi-moment Reconstruction

The piecewise reconstruction polynomial for physical field ϕ(x,y)on cell Ω_i is written in the local coordinate system as,

The unknown coefficients are uniquely determined by which are obtained from the following constraint conditions in terms of the multiple moments, where the PVs at the vertices,ϕ_ij,and the VIA,ϕ_i,as defined in (6),are the computational variables to be predicted at every step. The derivatives at the cell center,ϕ_ξi and ϕ_ηi,are computed from the PVs and VIAs of the neighboring cells as described in appendix A in [23].

Interpolation function (12) can be equivalently cast into a basis function form as,

where the basis functions read which are identical for all mesh cells.

In the present model,both VIA and PV are treated as the computational variables and updated at every time step through the pressure-projection solution procedure shown above. It should be noted that inclusion of the PVs as new computational variables causes increase in computational cost. For a single physical variable,we have to memorize both the PV for each vertex and VIA for each grid cell. An increase in CPU is also observed. All these depend on the topology of the grid elements. A detailed comparison with the conventional finite volume method will be shown later. It reveals that the present method is more superior regarding the compromise between numerical accuracy and computational cost.

Next,we discuss the spatial discrete schemes at each sub-step in details.

1.3 Scheme for Advection-diffusion Equation

The advection-diffusion equation in 2D is written in conservative form as

where ϕ is the physical quantity,which stands for either u or v in 2D Navier-Stokes equations. Remembering that all terms on the right hand side of (17) are known at each time level,we omit the superscript m hereafter unless it is needed.

In the present method,two types of moments,i.e. the volume integrated average (VIA) and the point value (PV) of the physical field ϕ,are treated as the prognostic variables and updated separately. Their governing equations can be immediately obtained from definition (6) and (17),i.e.

for VIA moment and for PV moment. Note that we have used Gauss divergence theorem and assumed a constant kinematic viscosity to yield (18).

Next,we describe the spatial discretization for the terms on the right hand sides of (18) and (19).

1.3.1 Advection Flux

The total advection flux for cell Ω_i is approximated by the summation of those across each boundary segment,

where |Γ_ij| represents the length of edge Γ_ij,ϕ_Γij the line-integrated average value of ϕ along Γ_ij,and u_Γij=(u_Γij,ν_Γij) the velocity averaged over edge Γ_ij.

The line-integrated average of ϕ along edge Γ_ij is not continuous and has different values for the two neighboring cells sharing Γ_ij. We choose the upwinding side value according to the hyperbolic feature of the advection equation,

where Φ_{i up} denotes the reconstruction function (12) or (15) over the upwinding cell,i.e.

It should be noted that the Simpson quadrature formula in (21) is exact for reconstruction (12) or (15). To show this,we give explicitly the line-integrated average of Φ_i along edge Γ_ij as

1.3.2 Diffusion Flux

In a similar way,the total diffusion flux for cell Ω_i is approximated by the summation of those across each boundary edge,

where ϕ_xΓij and ϕ_yΓij represent the line-integrated average values of ϕ_x and ϕ_y along edge Γ_ij,which are the average of the values computed from the reconstructions,Φ_i and Φ_ij,over two neighboring cells sharing edge Γ_ij,i.e. 1.3.3 Gradient Operator in the Advection Term The PV moment at the vertices (θ_ij) of the triangular mesh element Ω_i,ϕ_ij(j=1,2,3),is predicted from the differential form equation (19),where the gradients of the advection terms are computed from a weighted upwinding scheme, for derivative with respect to x and for derivative with respect to y. Here,Ω_k(k=1,2,…,K) represent the union of cells that share vertex θ_ij,and Φ_k the reconstruction functions (12) or (15) on cell Ω_k. The weight Ω_k is computed by where denotes the vector from the center of cell k to vertex θ_lj. It is clear from (29) that the gradient computed in the upwinding cell is used,which is effectively a point-wise Riemann solution. 1.3.4 Diffusion Operator

We denote again by Ω_k(k=1,2,…,K) the union of cells that share vertex (θ_ij). Given the net diffusion fluxes for these cells calculated by (24),the point wise value of the diffusion term at θ_ij can be simply obtained from the Time-evolution converting (TEC) formula described in [16, 19]. To be more precise,the time tendency of VIA in this case is computed by

which is already obtained by (24).

The time tendency of the PV due to the Laplace operator is then computed from the TEC formula based on the least square minimization described in appendix B in [23],which yields

1.3.5 Time Integration Scheme

The semi discrete equations (18) and (19) are used to update the VIA and PV moments. We summarize these equations into a form of ordinary differential equation,

where L(ϕ) stands for the spatial approximations discussed above.

Given the solution ϕ^m at time level t^m,we use the Runge-Kutta time integration schemes to predict the solution ϕ^* for the advection-diffusion equation. Both second and third order Runge-Kutta schemes can be used. The third order TVD Runge-Kutta time integration method [24] is applied for time integration in this paper. We give the scheme as follows.

where Δt=t^m+1-t^m is the time integration interval.

In the context of incompressible Navier-Stokes equations,the above procedure applies component-wisely to get the intermediate velocity field u^*.

1.4 Scheme for Pressure Poisson Equation

In the present formulation,we only use the VIA moment for pressure. The Poisson equation (4) is recast in an integral form.

Discretizing (34) over cell Ω_i yields where u_Γij^* and ν_Γij^* on the right hand side are obtained from reconstruction function in the same manner as described before. The pressure gradient is appro ximated by a linear interpolation spanned over two neighboring cells sharing edge Γ_ij. We denote (∂p/∂r)_ije_ij as the pressure gradient along line segment r_ij= r(θ_cij)-r(θ_ci) with r(θ_cij)=(x_cij,y_cij) and r(θ_ci)=(x_ci,y_ci) being the position vectors of the centroids of cells Ω_ij and Ω_i respectively,i.e.

The unit vector of r_ij is denoted by e_ij=(e_xij,e_yij)= r_ij/|r_ij|.

The linear interpolation function over two neighboring cells Ω_i and Ω_ij is then written in the following form,

The unknowns p₀ and (∂p/∂r)_ij in (36) are uniquely determined by requiring From (37) and (38),with some algebraic manipulation,we have

Following [25, 26],we consider a correction to the non-orthogonality of the mesh. We decompose the normal vector of edge Γ_ij as follows,

where n_‖ij represents a vector co-linear with edge Γ_ij. The pressure gradient term in (35) is then split into two parts,i.e. the orthogonal part and the non-orthogonal correction part, As adopted in (17),we calculate the orthogonal component implicitly using the values at level m+1 and the non-orthogonal correction explicitly with the values at level m,

Thus,Poisson equation (35) for pressure is finally written as

where

and

Equation (43) is a simultaneously linked linear algebraic equation for pressure p_i^m+1 and can be solved by iterative solvers.

1.5 Projection of the Velocity Field

As the final step of the solution procedure,the velocity field must be projected onto a divergence-free subset. From the pressure field computed above,we firstly calculate the pressure gradients over each boundary edge by

and

The VIA of velocity components are then updated by,

and

Given the increments of VIA from (46) and (47) on the union of cells (Ω_k,k=1,2,…,K) that share vertex (θ_ij),we update the PVs of velocity at the vertices again by the time-evolution converting (TEC) formula in appendix B as

and 1.6 Summary of the Numerical Procedure

In order to allow the users to follow the algorithm easily,we summarize again the solution procedure as the following steps.

1) Given the velocity field in terms of VIA (u_i^m) and PV (u_ij^m) at step m,compute the intermediate velocity for both VIA,u_i^*,and PV,u_ij^*,by (18) and (19) and the numerical algorithms detailed in section 1.3.

2) Solve Poisson equation (43) to get pressure ^m+1_i at step (m+1).

3) Update the VIA and PV of velocity to the new time level,u_i^m+1 and u_ij^m+1,by projection correction (46)-(49).

1.7 Extension of the Numerical Procedure

All the numerical procedures in 2D described above can be extended to 3D in a straightforward manner. Here,we only present the 3D multi-moment reconstruction for tetrahedral mesh element,presuming that the rest of the 3D formulation can be derived by analogizing the corresponding component in 2D.

We use the tetrahedral element for three dimensions. Shown in Fig. 2,cell Ω_i(i=1,…,N_e) has four vertices θ_ij(j=1,2,3,4) located at (x_i1,y_i1),(x_i2,y_i2),(x_i3,y_i3) and (x_i4,y_i4).

The four boundary surfaces of element Ω_i are denoted by .

In a similar way,the volume-integrated average (VIA) and point-value (PV) at the vertices of a physical variable ϕ(x,y,z,t) are defined as:

The 3D piecewise reconstruction polynomial for physical field ϕ(x,y,z) in cell Ω_i in the local coordinate (ξ,η,ζ) system is cast into a basis function form in terms of multiple moments as,

where ϕ_ξi,ϕ_ηi and ϕ_ξi are the partial derivatives at the cell center,and the basis functions read

It is noted that the basis function (52) are determined from the following constraint conditions in terms of the multiple moments,

where the PVs at the vertices,ϕ_ij,and the VIA,ϕ_i,as defined in (50),are the computational variables to be predicted at every step. The derivatives at the cell center,are computed from the PVs and VIAs of the four neighboring cells in a similar way of the 2D case. 2 Numerical Results

We present some numerical tests to show the accuracy and robustness of the numerical method presented in this paper.

In order to quantify the numerical errors,we define L₁ and L₂ errors as follows,

and where ϕ_ni and ϕ_ei denote numerical and exact solutions respectively. 2.1 Advection of a Sine Wave

To evaluate the convergence rate of the advection scheme presented in section 1.3,we computed the advection transport of a sine function as in [27] with gradually refined grids. The advected profile is initially defined as

over a [0, 1]×[0, 1] computational domain,the periodic boundary condition is imposed in both x and y directions.

A constant and uniform advection velocity is specified

Numerical experiments were conducted on gradually refined grids generated by Delaunay algorithm with doubly increased resolution. The time step of Δt is used throughout all computations.

The L₁ errors and the convergence rates of VPM and other standard schemes at t=1.0 are given in Table 1 for grids of different resolutions.

It is observed that the presented scheme is of 3th order accuracy. Because the multi-moment reconstructions is of second order,the VIA is theoretically expected to be of third order accuracy,which is justified by the numerical output. Other standard schemes widely used for unstructured grids show convergence rates less than 2nd order,among which the QUICK scheme is more accurate for this smooth-profile test problem. Moreover,a comparison between the errors of VPM and standard schemes reveal a large improvement in numerical accuracy. On the grid of 57 518 elements,for example,the L₁ error of VPM is only 3% of the QUICK scheme.

2.2 Solid Rotation of a 2D Cone

It is well known that numerical solutions on unstructured grid depend on the quality of computational grid. In this test,we examine the robustness of the advection scheme to the quality of computational grid. A 2D cone of the initial distribution specified by

over [0, 1]×[0, 1] is transported with a rotational velocity field defined by u=-2π(y-0.5) and v=2π(x-0.5).

We use two unstructured grids with triangular elements. As shown in Fig. 3,grid A is generated by the Delaunay algorithm and has good quality,while grid B has a region where the mesh elements are distorted with large skewness and high aspect ratio. The numerical results of the VPM scheme after one revolution of rotation are shown in Fig. 4. It is found that the degradation in mesh quality doesn’t significantly affect the numerical results.

For comparison,we also include the numerical results from other standard advection schemes used in conventional FVM codes,such as the QUICK scheme [28],Superbee-TVD scheme [29] and MUSCL scheme [2, 3],which are available in the OpenFoam [30] platform. Figs. 5,6 and 7 display the corresponding outputs of these schemes on two grids.

Compared against Fig. 4,we observe that the VPM scheme is the most accurate on both grids. More important,the conventional schemes are more sensitive to the mesh quality,whose numerical results look remarkably worse on grid B.

To quantify this observation,we give in Table 2 the L₂ errors on both A and B grids for VPM method and other standard finite volume schemes with the “degradation rate” defined by

As an indicator of the robustness to the computational grid,a smaller value of the degradation rate shows a better robustness in respect to the quality of mesh elements.

We can see that having the smallest degradation rate,the VPM scheme is the most robust one with regard to the computational grid,and thus retains the numerical accuracy even on a grid of worse quality. It is very sensible if we consider that in the present scheme both the cell average and vertex value are updated,which compensates the loss of the information due to the worse nonorthogonality and skewness of the mesh.

2.3 Taylor Vortex Problem

To evaluate the accuracy of the whole incompressible fluid solver,we computed the Taylor vortex problem [9, 11]. An analytical solution for this problem is available as

where Re is Reynolds number. We use Re=100 in this example.

We computed this test problem over [0, 1]×[0, 1] domain with the periodical boundary condition,using grids of different resolutions. Our interest here is to evaluate the numerical accuracy in spatial discretization,and a relatively small time stepping increment (Δt=10^-4) is used. We show the numerical errors and convergence rates for both velocity and pressure fields at t=0.1 in Table 3 and Fig. 8. We also include the results of the conventional finite volume method with the QUICK,Superbee-TVD and MUSCL schemes as the advection solvers. It should be noted that the conventional FVM model was generated by merging the templates provided by the OpenFOAM source code into the projection solution procedure shown in section 1.

It is observed that the VPM method is more accurate than the conventional FVM models regarding both numerical errors and convergence rate. Numerical solution in velocity of VPM shows a convergence rate higher 2.5. In the results on the refined mesh (57 670 cells),due to the high convergence rate the numerical errors of VPM model in both velocity and pressure are reduced down to only 1% of those in the conventional FVM models. It is also found that the errors in velocity are approximately one-tenth of those in pressure in all results.

2.4 Remarks on Computational Cost

Using both VIA and PV as the computational variables causes an increase in memory requirement. For a single physical field,we need to keep the PVs at the element vertices as the new degrees of freedom (DOFs) in addition to the VIAs as in the conventional FVM. It also leads to an increase in CPU consumption.

In order to quantify the computational cost in comparison with the conventional FVM,we run the 2D advection test with different grid resolutions for 1000 steps on a PC with a single CPU (Intel(R) Xeon E5-2687W,3.10 GHz). We give in Table 4 the DOFs for the memory requirement and elapse time for CPU consumption. It is observed that for the triangular mesh the number of DOFs increases about 51%. As described before,the spatial discretization can be computed by firstly finding and saving the numerical fluxes and the derivatives from the multi-moment reconstructions for all cells,and then updating VIA cell-wisely and PV point-wisely. It is noted that the PVs at cell vertices are readily for use in the calculation of the numerical fluxes,which simplifies the operations and saves CPU time. Shown in Table 4,the elapse time increase may reach 1.5 fold in comparison with the QUICK scheme. As shown in section 2.1,large improvement can be obtained by VPM scheme. From Table 1,we see that the numerical error of VPM on a 5 463-element grid is smaller than that of the QUICK scheme on a 14 412-element grid. It leads to an observation that for a given error level the VPM scheme is much more efficient in terms of both memory requirement and CPU consumption compared to the standard schemes due to its high accuracy and faster convergence.

The computational cost of the full fluid solver is also evaluated using Taylor vortex test shown above. About 34% increase is seen in the memory requirement due the fact that we only use the VIA for pressure field. The increase in elapse time varies with grid resolution,which is smaller for a coarse grid and becomes larger toward a saturation value of 42% when the number of grid cells is larger than 10⁴. The increase in CPU consumption of the whole fluid solver is less significant than that of the pure advection computation. This result is consistent with the well-known observation that the pressure Poisson equation consumes the large part of CPU time.

It is observed that the numerical errors of VPM can be smaller than the conventional FVM by two orders of magnitude. Similar to the advection case,we can conclude that the VPM model is much superior in computational efficiency and requires much less hardware resource to reach a given error level in comparison with conventional FVM model. This justifies the practical utility of the VPM method.

2.5 Lid Driven Cavity Flows with Different Shapes

Incompressible viscous flows in closed cavities under the forcing of the upper lid moving at a constant speed provide a benchmark test set to verify numerical solvers. Following the standard test of Ghia et al. [31],other tests of different geometric configurations were suggested as well to evaluate the capability of numerical codes to deal with complex geometries. We show next the results of the VPM model for four cavities of different shapes,i.e. square,forward-step duct,skewed rectangle and semi-circle. The Reynolds number used in the following four test cases is 1000.

2.5.1 Lid Driven Square Cavity Problem

This is the most widely used benchmark test for numerical codes designed to solve viscous incompressible flows [31]. Incompressible flow is enclosed in a square [0, 1]×[0, 1] domain. Nonslip boundary conditions are imposed on the 4 walls. The upper wall y =1is moving at a constant speed (u=1,v=1). Numerical tests were carried out with triangular unstructured grids (upper-left panel of Fig. 9) of different resolutions. It is known that a steady solution is reached with this Reynolds number,which is identified by the criterion of |p_i^m+1-p_i^m|＜1×10^-9. We plot the numerical results of different grid resolutions in Fig. 10. It is observed that the numerical solution of present model converges to the reference solution with a modest mesh resolution. We also show the results of the conventional FVM with different advection schemes in Fig. 11 for comparison. It is obvious that the present model is more accurate than the conventional FVM.

2.5.2 Lid Driven Square Cavity Problem

A cavity with a forward-step shown in the upper-right panel of Fig. 9 is used in [32] to evaluate numerical solvers developed in general coordinates. Numerical results are given in Fig. 12 with the reference velocity profiles provided in [32] which were computed on a 512×512 curvilinear grids. Shown in Fig. 12(b)-(c),the present model is able to reproduce adequate solutions even with much coarser mesh resolutions. The solutions on meshes of more than 2 704 cells agree well with the reference solution.

2.5.3 Lid Driven Skew Cavity Problem

Numerical example of viscous flow in skew cavities with inclined angles of 30° is found in [33] and [32]. Our numerical results are shown in Fig. 13. The reference solution in [33] is computed on 320×320 structured grids. Our numerical solution converges to the reference solution with much coarser grids compared to other existing conventional FVM schemes.

2.5.4 Lid Driven Semi-circular Cavity Problem

We computed the lid-driven viscous incompressible flow in a semi-circular cavity (bottom-right panel of Fig. 9) which has been extensively investigated in [34]. We use the results in [34] for comparison which were computed on a triangular mesh with a representative size of 1/128. Shown in Fig. 14,our results converge to the reference results even with a relatively low mesh resolution.

2.6 3D Lid Driven Cubic Cavity Problem

As a verification of the 3D code,we simulated a 3D lid driven flow in a unit cube[-0.5,0.5]×[-0.5,0.5]×[-0.5,0.5]. Following [35],we imposed a constant velocity v(-0.5,y,z)=1 on surface x=-0.5 and computed this test using tetrahedral grids of different resolutions. The velocity component in y direction (v) along centerline (x,0,0) and that in x direction (u) along centerline (0,y,0) are plotted in Fig. 15. For comparison,we also include the numerical results in [35] as a reference solution which are solved by a Chebyshev-collocation method with 96×96×64 Gauss-Lobatto points. It is seen that the present solver converges quickly to the reference solution.

3 Conclusion Remarks

We have presented and evaluated a novel numerical solver for i ncompressible Navier-Stokes equations on triangular and tetrahedral unstructured meshes by using the multi-moment finite volume formulation. The numerical method,so-called VPM,treats the point value (PV) at the cell vertex as a new computational variable in addition to the volume integrated average (VIA) value that is the only computational variable in the conventional finite volume method. VIA is solved via a finite volume formulation of flux form,and is thus numerically conserved. PV is computed point-wisely based on the differential form of the governing equations which need not to be numerically conservative but can be solved efficiently. The numerical formulation can be very easily implemented on unstructured grids. With multiple moments,more DOFs can be used for each grid cell and high order reconstruction can be built on a compact stencil to calculate the numerical flux and derivatives for spatial discretization,which appears to be particularly important to unstructured grids. Another advantage of the present method is that there is not any arbitrariness in choosing the neighboring cells for the reconstructions.

As shown in this paper,when applied to incompressible fluid dynamics,VPM model causes about 30%-40% increase in the memory requirement and 40%-50% increase in CPU time. Nevertheless,a large leap in improving numerical accuracy and robustness is achieved. Our numerical results show that the numerical errors of VPM can be smaller than the conventional FVM by two orders of magnitude. The numerical experiments reveal that the numerical results of the present method converge to the reference solutions even with the relatively coarse grid resolutions. Consequently,the VPM model is much superior in computational efficiency and requires much less computer resource to reach a given error level in comparison with conventional FVM models. Moreover,the numerical solutions of the VPM model are less dependent on the quality of computational grids. All these justify the practicability of the VPM method in real-case applications.

The present formulation uses only cell average and point values at cell vertices as the unknown DOFs for triangular (2D) and tetrahedral (3D) elements. We have also applied the same idea to unstructured grids with other elements,such as quadrilateral (2D) and hexahedral (3D) elements,and obtained very promising results,which will be reported soon.

Higher order formulations for unstructured grids with other types of moments need further exploration. Given the DG formulation as a standard track,how to balance the algorithmic cost and the solution quality,which is always problem-dependent,should be a key in future practices.

To sum up,with a modest increase in computational complexity and cost,the present VPM method is able to achieve significant improvements in accuracy and robustness in comparison with conventional finite volume method. Thus,it provides a practical solver for incompressible viscous flows which well balances the numerical accuracy and computational complexity.