point-value formulation

Accepted Manuscript A multi-moment finite volume method for incompressible Navier–Stokes equations on unstructured grids: volume-average/point-value formulation Bin Xie, Satoshi Ii, Akio Ikebata, Feng Xiao

PII: DOI: Reference:

S0021-9991(14)00557-9 10.1016/j.jcp.2014.08.011 YJCPH 5421

To appear in:

Journal of Computational Physics

Received date: 11 October 2013 Revised date: 30 July 2014 Accepted date: 4 August 2014

Please cite this article in press as: B. Xie et al., A multi-moment finite volume method for incompressible Navier–Stokes equations on unstructured grids: volume-average/point-value formulation, J. Comput. Phys. (2014), http://dx.doi.org/10.1016/j.jcp.2014.08.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A multi-moment finite volume method for incompressible Navier-Stokes equations on unstructured grids: volume-average / point-value formulation Bin Xiea , Satoshi Iib , Akio Ikebatac , Feng Xiaoa,∗ a Department

of Energy Sciences, Tokyo Institute of Technology, 4259 Nagatsuta Midori-ku, Yokohama, 226-8502, Japan. b Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka, 560-8531, Japan. c Production Technology Division, TOTO LTD., 1-1 Nakashima 2-chome, Kokurakita-ku, Kitakyushu, 802-8601, Japan.

Abstract A robust and accurate finite volume method (FVM) is proposed for incompressible viscous fluid dynamics on triangular and tetrahedral unstructured grids. Different from conventional FVM where the volume integrated average (VIA) value is the only computational variable, the present formulation treats both VIA and the point value (PV) as the computational variables which are updated separately at each time step. The VIA is computed from a finite volume scheme of flux form, and is thus numerically conservative. The PV is updated from the differential form of the governing equation that does not have to be conservative but can be solved in a very efficient way. Including PV as the additional variable enables us to make higher-order reconstructions over compact mesh stencil to improve the accuracy, and moreover, the resulting numerical model is more robust for unstructured grids. We present the numerical formulations in both two and three dimensions on triangular and tetrahedral mesh elements. Numerical results of several benchmark tests are also presented to verify the proposed numerical method as an accurate and robust solver for incompressible flows on unstructured grids. Keywords: Finite volume method, unstructured grid, robustness, accuracy, triangular/tetrahedral mesh, incompressible flow, multi-moment, compact-stencil.

1. 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. Nearly all commercial CFD codes are based on the “FVM+unstructured grid” prevailing paradigm. However, developing high-order finite volume model on unstructured grids is not a trivial task. Conventional FVM requires wide stencil to generate highorder reconstructions as in the k-exact finite volume method [3] and the weighted essentially non-oscillatory method [13, 19]. However, the choice of the cell stencil is not straightforward. As mentioned in [13], 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, including the major commercial codes, are limited to the second order on unstructured grids, where a linear reconstruction such as the MUSCL( Monotone Upstreamcentered Schemes for Conservation Laws) scheme [38, 39] is used. The numerical schemes with compact stencil, on the other hand, add the degrees of freedoms (DOFs) locally on each cell for high order reconstructions, and are more flexible and easier to implement on unstructured meshes. Among many others, the representative methods of this sort are the discontinuous Galerkin (DG) method [7, 8, 9, 10] and the spectral finite volume (SV) method [40, 41]. These high order schemes were originally devised for convectiondominant flow problems, and make use of the Riemann solvers for the full Euler flux functions across cell boundaries ∗ Corresponding

author: Dr. F. Xiao (Email: [email protected])

Preprint submitted to Journal of Computational Physics

August 8, 2014

where the reconstructed physical fields are allowed to be discontinuous. These schemes turned out to be a great success in solving convection-dominated flows. However, not as many as implementations of these methods are found for the incompressible Navier-Stokes equations where one has substantial difficulty because of the singularity in the hyperbolic system as the sound speed can be infinite. 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[6], MAC(Marker and Cell)[17] method, SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) method[29], SIMPLEC(SIMPLE Corrected) method[37] and PISO(Pressure Implicit with Splitting Operators) method[22]. 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 [32, 18, 12, 34]. 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 pointwisely updated by a differential form 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. These 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. As justified in this paper, without large increase in algorithmic complexity and computational cost, significant improvements are achieved in (1) the convergence rate in advection computation, (2) the accuracy of the whole fluid solver in terms of numerical errors and (3) the robustness to the computational mesh. The present formulation can be seen as an extension of the CIP multi-moment finite volume methods [47, 48, 46, 43, 45, 20, 44, 21, 1, 4] 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. This paper is organized as follows. The numerical formulation is presented with details in section 2. Numerical tests are given in section 3 to verify the accuracy and the robustness of the present method in comparison with some standard FVM schemes. We end this paper by some conclusion remarks in section 4. 2. Numerical formulation on unstructured mesh We consider the incompressible unsteady Navier-Stokes equations, ∇ · u = 0, 1 ∂u + ∇ · (u ⊗ u) = − ∇p + ν∇2 u, ∂t ρ

(1) (2)

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 density and ν the kinematic viscosity. 2

2.1. The projection solution procedure The projection method of Chorin [6] provides a simple and robust solution procedure for incompressible flow and is adopted in this paper. For the completeness, we briefly summaries the numerical steps to update the velocity from time level m (t = tm ) to m + 1 (t = tm + Δt) as follows. 1. Given the velocity field um at step m, compute the intermediate velocity u∗ from the momentum equation (2) without the pressure gradient term, ∂u = −∇ · (um ⊗ um ) + ν∇2 um . ∂t

(3)

2. The intermediate velocity u∗ usually does not satisfy the mass conservation equation (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,   1 m+1 1 ∇ · u∗ . (4) ∇ · ∇p = ρ Δt 3. Correct the velocity by the projection step, um+1 − u∗ 1 = − ∇pm+1 . Δt ρ

(5)

It is straightforward to show that the updated velocity satisfies ∇ · um+1 = 0. Our major interest in this paper is to develop a novel approach for the spatial discretization for the solution procedure shown above. 2.2. The 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, ..., Ne ) with three vertices θi j ( j = 1, 2, 3) located at (xi1 , yi1 ), (xi2 , yi2 ) and (xi3 , yi3 ) respectively as shown in Fig.1(left). We also denote the mass center of Ωi by θic . The three boundary line segments of element Ωi are denoted by Γi1 = θi2 θi3 , Γi2 = θi3 θi1 and Γi3 = θi1 θi2 . For the later use, we denote the middle points and the outward normal unit vector of segment Γi j by θ˜i j and ni j = (n xi j , nyi j ). Shown in Fig.1 (right), the volume-integrated average (VIA) and point-value (PV) at the vertices of a physical variable φ(x, y, t) are defined as:  1 φi (t) ≡ φ(x, y, t)dΩ, |Ωi | Ωi (6) φi j (t) ≡ φ(xi j , yi j , t), j = 1, 2, 3; 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 ω ≡ [0 ≤ ξ, η ≤ 1, ξ + η ≤ 1] with the local coordinate (ξ, η) defined by ⎧ ⎪ ⎪ ⎨ x = xi1 + ξ(xi2 − xi1 ) + η(xi3 − xi1 ), (7) ⎪ ⎪ ⎩y = yi1 + ξ(yi2 − yi1 ) + η(yi3 − yi1 ). We give some numerical formula to transform the first-order derivatives between the global coordinate system (x, y) and the local coordinate (ξ, η), ∂φi (x, y) ∂φi (x, y) + yξi , ∂x ∂y ∂φi (x, y) ∂φi (x, y) φηi (x, y) = xηi + yηi ∂x ∂y

φξi (x, y) = xξi

3

(8)

θi 3

Γi 2

φi 3

Γi1

Ωi θi1

Γi 3

φi θi 2

φi1

φi 2

Figure 1: The triangular mesh element (left) and the definition of discrete moments (right).

and

  1 ∂φi (ξ, η) ∂φi (ξ, η) − yξi (ξ, η) φ xi (ξ, η) = yηi (ξ, η) , |Ji (ξ, η)| ∂ξ ∂η   1 ∂φi (ξ, η) ∂φi (ξ, η) + xξi (ξ, η) φyi (ξ, η) = −xηi (ξ, η) , |Ji (ξ, η)| ∂ξ ∂η

(9)

where |Ji (ξ, η)| = xξi (ξ, η)yηi (ξ, η) − xηi (ξ, η)yξi (ξ, η) is the determinant of the Jacobian matrix. It is straightforward from (7) that xξi = xi2 − xi1 , xηi = xi3 − xi1 , yξi = yi2 − yi1 , yηi = yi3 − yi1 ,

(10)

|Ji (ξ, η)| = 2|Ωi | = (xi2 − xi1 )(yi3 − yi1 ) − (xi3 − xi1 )(yi2 − yi1 ).

(11)

and

2.3. The multi-moment reconstruction The piecewise reconstruction polynomial for physical field φ(x, y) on cell Ωi is written in the local coordinate system as, (12) Φi (ξ, η) = c00 + c10 ξ + c01 η + c11 ξη + c20 ξ2 + c02 η2 . The unknown coefficients are uniquely determined by ⎧ ⎪ c00 = φi1 ⎪ ⎪ ⎪ ⎪ ⎪ 10 4 2 ⎪ ⎪ c 10 = − 3 φi1 − 3 φi2 + 3 φi3 + 4φi + φξi − 2φηi ⎪ ⎪ ⎪ ⎪ 10 2 ⎪ ⎪ ⎨c01 = − 3 φi1 + 3 φi2 − 43 φi3 + 4φi − 2φξi + φηi ⎪ ⎪ 2 2 ⎪ c11 = 16 ⎪ ⎪ 3 φi1 − 3 φi2 − 3 φi3 − 4φi + 2φξi − 2φηi ⎪ ⎪ ⎪ 7 7 2 ⎪ ⎪ c20 = 3 φi1 + 3 φi2 − 3 φi3 − 4φi − φξi + 2φηi ⎪ ⎪ ⎪ ⎪ ⎩c02 = 7 φi1 − 2 φi2 + 7 φi3 − 4φ + 2φξi − 2φηi , i 3 3 3 which are obtained from the following constraint conditions in terms of the multiple moments, ⎧ ⎪ Φi (0, 0) = φi1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Φ ⎪ i (1, 0) = φi2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Φ ⎪ ⎪  i(0, 1) = φi3 ⎪ ⎪ ⎪ ⎨ Φi (ξ, η)dξdη = φi /2 ⎪ ⎪ ⎪ ω ⎪ ⎪ ⎪ ∂Φi (ξ, η) ⎪ ⎪ ⎪ ( 1 , 1 ) = φξi ⎪ ⎪ ⎪ 3 3 ∂ξ ⎪ ⎪ ⎪ ⎪ ∂Φi (ξ, η) ⎪ ⎪ ⎪ = φηi ⎩ ( 13 , 13 ) ∂η 4

(13)

(14)

where the PVs at the vertices, φi j , 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. Interpolation function (12) can be equivalently cast into a basis function form as, Φi (ξ, η) = ψ1 φi1 + ψ2 φi2 + ψ3 φi3 + ψ φi + ψξ φξi + ψη φηi ,

(15)

where the basis functions read ⎧ 10 16 7 2 7 2 ⎪ ⎪ ψ1 = 1 − 10 ⎪ 3 ξ − 3 η + 3 ξη + 3 ξ + 3 η ⎪ ⎪ ⎪ 2 4 2 7 2 2 2 ⎪ ⎪ ψ2 = 3 η − 3 ξ − 3 ξη + 3 ξ − 3 η ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ψ3 = 23 ξ − 43 η − 23 ξη − 23 ξ2 + 73 η2 ⎪ ⎪ ⎪ ψ = 4ξ + 4η − 4ξη − 4ξ2 − 4η2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ψξ = ξ − 2η + 2ξη − ξ 2 + 2η2 ⎪ ⎪ ⎪ ⎪ ⎩ψη = η − 2ξ + 2ξη + 2ξ 2 − η2

(16)

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 a more superior trade-off regarding the numerical accuracy and computational cost. Next, we discuss the spatial discrete schemes at each sub-step in details. 2.4. The scheme for advection-diffusion equation The advection-diffusion equation in 2D is written in conservative form as ∂φ = −∇ · (φu) + ν∇2 φ, ∂t

(17)

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. dφi 1 =− dt |Ωi |

Γi

(φu · n) dΓ +

for VIA moment and

ν |Ωi |

Γi

(n · ∇φ) dΓ

(18)

dφi j = −ui j · (∇φ)i j + ν(∇2 φ)i j , (19) dt 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).

5

2.4.1. The advection flux in (18) The total advection flux for cell Ωi is approximated by the summation of those across each boundary segment,

Γi

(φu · n) dΓ =

3 

3



|Γi j |φuΓi j · ni j ≈ |Γi j |φΓi j uΓi j · ni j ,

j=1

(20)

j=1

where |Γi j | represents the length of edge Γi j , φΓi j the line-integrated average value of φ along Γi j , and uΓi j = (uΓi j , vΓi j ) the velocity averaged over edge Γi j . The line-integrated average of φ along edge Γi j is not continuous and has different values for the two neighboring cells sharing Γi j . We choose the upwinding side value according to the hyperbolic feature of the advection equation, ⎞ ⎛ 3  ⎟⎟⎟ 1 1 ⎜⎜⎜⎜ φΓi j = Φiup (x, y)dΓ = ⎜⎝⎜ φik + 4Φiup (θ˜i j )⎟⎟⎠⎟ , (21) |Γi j | Γi j 6 k j where Φiup denotes the reconstruction function (12) or (15) over the upwinding cell, i.e. ⎧ ⎪ ⎪ ⎨= i, for ni j · uΓi j > 0; iup = ⎪ ⎪ ⎩= i j, otherwise.

(22)

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 Γi j as

⎧ ⎪ ⎪ φΓi1 = 19 φi1 + 19 φi2 + 19 φi3 + 23 φi + 16 φξi + 16 φηi = 16 φi2 + φi3 + 4Φi (θ˜i1 ) , ⎪ ⎪

⎪ ⎨ (23) φΓi2 = 19 φi1 + 19 φi2 + 19 φi3 + 23 φi − 13 φξi + 16 φηi = 16 φi1 + φi3 + 4Φi (θ˜i2 ) , ⎪ ⎪

⎪ ⎪ ⎪ ⎩φΓ = 1 φi2 + 1 φi1 + 1 φi3 + 2 φi + 1 φξi − 1 φηi = 1 φi1 + φi2 + 4Φi (θ˜i3 ) . 9 9 9 3 6 3 6 i3 2.4.2. The diffusion flux in (18) In a similar way, the total diffusion flux for cell Ωi is approximated by the summation of those across each boundary edge,

3 

(n · ∇φ) dΓ = |Γi j |(φ xΓi j n xi j + φyΓi j nyi j ) , (24) Γi

j=1

where φ xΓi j and φyΓi j represent the line-integrated average values of φ x and φy along edge Γi j , which are the average of the values computed from the reconstructions, Φi and Φi j , over two neighboring cells sharing edge Γi j , i.e. ⎛ 3 ⎞



⎟⎟⎟ 1 ⎜⎜⎜⎜ ⎜⎜ Φ xi (θik ) + Φ xi j (θik ) + 4 Φ xi (θ˜i j ) + Φ xi j (θ˜i j ) ⎟⎟⎟⎠ φ xΓi j = (25) 12 ⎝ k j and φyΓi j

⎛ 3 ⎞



⎟⎟⎟ 1 ⎜⎜⎜⎜ ⎜⎜ Φyi (θik ) + Φyi j (θik ) + 4 Φyi (θ˜i j ) + Φyi j (θ˜i j ) ⎟⎟⎟⎠ . = 12 ⎝ k j

(26)

2.4.3. The gradient operator in the advection term in (19) The PV moment at the vertices (θi j ) of the triangular mesh element Ωi , φi j ( 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,   K  ∂φ = ωk Φ xk (θi j ) (27) ∂x i j k=1 6

for derivative with respect to x and



∂φ ∂y

 = ij

K  ωk Φyk (θi j )

(28)

k=1

for derivative with respect to y. Here, Ωk (k = 1, 2, · · · , K) represent the union of cells that share vertex θi j , and Φk the reconstruction functions (12) or (15) on cell Ωk . The weight ωk is computed by   −−−−→ max 0, ui j · θkc θi j  , (29) ωk =  −−−−→ K k=1 max 0, ui j · θkc θi j −−−−→ where θkc θi j denotes the vector from the center of cell k to vertex θi j . It is clear from (29) that the gradient computed in the upwinding cell is used, which is effectively a point-wise Riemann solution. 2.4.4. The diffusion operator in the diffusion term in (19) We denote again by Ωk (k = 1, 2, · · · , K) the union of cells that share vertex (θi j ). Given the net diffusion fluxes for these cells calculated by (24), the point wise value of the diffusion term at θi j can be simply obtained from the Time-evolution converting (TEC) formula described in appendix B. To be more precise, the time tendency of VIA in this case is computed by

ν (n · ∇φ) dΓ, (30) δt φk = |Ωk | Γk 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 (B.1) described in appendix B, which yields ν(∇2 φ)i j = δt φi j = T EC(δt φk ).

(31)

2.4.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, dφ = L (φ) dt

(32)

where L(φ) stands for the spatial approximations discussed above. Given the solution φm at time level tm , 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 [33] is applied for time integration in this paper. We give the scheme as follows. φ(1) = φn + ΔtL(φn ) 3 1 φ(2) = φn + φ(1) + 4 4 1 2 φ∗ = φn + φ(2) + 3 3

1 ΔtL(φ(1) ) 4 2 ΔtL(φ(2) ), 3

(33)

where Δt = tm+1 − tm 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∗ .

7

2.5. The 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. 



1 1 (n · u∗ ) dΓ. (34) n · ∇p dΓ = Δt Γi Γi ρ Discretizing (34) over cell Ωi yields ⎞   3 ⎛ 3 

⎜⎜⎜ |Γi j | ∂p ⎟⎟ 1  |Γi j |(u∗Γi j n xi j + v∗Γi j nyi j ) , ei j · ni j ⎟⎟⎠ = ⎝⎜ ρi j ∂r i j Δt j=1 j=1

(35)

where u∗Γi j and v∗Γi j on the right hand side are obtained from reconstruction function in the same manner as described before. The pressure gradient is approximated by a linear interpolation spanned over two neighboring cells sharing edge Γi j . We denote ( ∂p ∂r )i j ei j as the pressure gradient along line segment ri j = r(θci j ) − r(θci ) with r(θci j ) = (xci j , yci j ) and r(θci ) = (xci , yci ) being the position vectors of the centroids of cells Ωi j and Ωi respectively, i.e. 1 1 1 1 xik , xci j = xi jk , yci = yik , yci j = yi jk . 3 k=1 3 k=1 3 k=1 3 k=1 3

xci =

3

3

3

The unit vector of ri j is denoted by ei j = (e xi j , eyi j ) = ri j /|ri j |. The linear interpolation function over two neighboring cells Ωi and Ωi j is then written in the following form,   ∂p P(x, y) = p0 + (e xi j x + eyi j y). (36) ∂r i j The unknowns p0 and ( ∂p ∂r )i j in (36) are uniquely determined by requiring  1 P(x, y)dΩ = pi , |Ωi | Ωi  1 P(x, y)dΩ = pi j . |Ωi j | Ωi j From (37) and (38), with some algebraic manipulation, we have   (pi j − pi ) ∂p . = ∂r i j e xi j (xci j − xci ) + eyi j (yci j − yci )

(37) (38)

(39)

Following [31, 25], we consider a correction to the non-orthogonality of the mesh. We decompose the normal vector of edge Γi j as follows, ei j ni j = + n i j , (40) e i j · ni j where n i j represents a vector co-linear with edge Γi j . The pressure gradient term in (35) is then split into two parts, i.e. the orthogonal part and the non-orthogonal correction part,       ∂p 1 ∂p ∂p ei j · n i j = + ei j · n i j (41) ∂r i j ei j · ni j ∂r i j ∂r i j As adopted in [23], we calculate the orthogonal component implicitly using the values at level m + 1 and the nonorthogonal correction explicitly with the values at level m,    m+1  m 1 ∂p ∂p ∂p ei j · n i j ≈ + ei j · n i j (42) ∂r i j ei j · ni j ∂r i j ∂r i j =

m (pm (pm+1 − pm+1 ) 1 ij i i j − pi ) + n i j · ei j . ei j · ni j e xi j (xci j − xci ) + eyi j (yci j − yci ) e xi j (xci j − xci ) + eyi j (yci j − yci j )

8

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

αi j (pm+1 − pm+1 ) = βi ij i

(43)

j=1

where αi j = and βi =

|Γi j | 1

ei j · ni j ρi j e xi j (xci j − xci ) + eyi j (yci j − yci j )

⎛ ⎞⎞ m 3 ⎛  (pm ⎜1 ⎜⎜⎜ ⎟⎟⎟⎟ 1 i j − pi ) ⎜⎝|Γi j | ⎜⎜⎜⎝ (u∗Γi j n xi j + v∗Γi j nyi j ) − n i j · ei j ⎟⎟⎠⎟⎟⎠ . ( Δt ρi j e xi j (xci j − xci ) + eyi j (yci j − yci j ) j=1

and can be solved by iterative Equation (43) is a simultaneously linked linear algebraic equation for pressure pm+1 i solvers. 2.6. 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



p xΓ

pyΓ

 ij

ij

= |Γi j |e xi j n xi j (

∂p m+1 ) , ∂r i j

(44)

= |Γi j |eyi j nyi j (

∂p m+1 ) . ∂r i j

(45)

The VIA of velocity components are then updated by, ⎛ ⎞ 3 

⎟⎟⎟ Δt ⎜⎜⎜⎜ 1 m+1 ∗ ⎜⎜  |n xi j |p xΓ ⎟⎟⎟⎠ , ui = ui − ij ρ|Ωi | ⎝ 3j=1 |n xi j | j=1 and vm+1 i

=

v∗i

⎞ ⎛ 3 

⎟⎟⎟ Δt ⎜⎜⎜⎜ 1 ⎜⎜  |nyi j |pyΓ ⎟⎟⎟⎠ . − ij ρ|Ωi | ⎝ 3j=1 |nyi j | j=1

(46)

(47)

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

(48) um+1 = u∗i j + T EC um+1 − u∗k k ij and

= v∗i j + T EC vm+1 − v∗k . vm+1 k ij

(49)

2.7. 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. m 1. Given the velocity field in terms of VIA (um i ) and PV (ui j ) at step m, compute the intermediate velocity for both ∗ ∗ VIA, ui , and PV, ui j , by (18) and (19) using the numerical algorithms detailed in 2.4.

at step (m + 1). 2. Solve Poisson equation (43) to get pressure pm+1 i 3. Update the VIA and PV of velocity to the new time level, um+1 and um+1 i i j , by projection corrections (46)-(49). 9

2.8. Extension to three dimensions 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 components in 2D. We use the tetrahedral element for three dimensions. Shown in Fig. 2, cell Ωi (i = 1, ..., Ne ) has four vertices θi j ( j = 1, 2, 3, 4) located at (xi1 , yi1 ), (xi2 , yi2 ), (xi3 , yi3 ) and (xi4 , yi4 ). The four boundary surfaces of element Ωi are denoted by Γi1 = (θi3 , θi2 , θi4 ), Γi2 = (θi4 , θi1 , θi3 ), Γi3 = (θi1 , θi4 , θi2 ), and Γi4 = (θi2 , θi3 , θi1 ).

θi 3

Γi 2

Γi 4

Ωi

Γi1

θi 4 Γi 3

θi1

θi 2 Figure 2: The tetrahedral mesh element for 3D.

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:  1 φi (t) ≡ φ(x, y, z, t)dΩ, |Ωi | Ωi (50) φi j (t) ≡ φ(xi j , yi j , zi j , t), j = 1, 2, 3, 4. 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, Φi (ξ, η, ζ) = ψ1 φi1 + ψ2 φi2 + ψ3 φi3 + ψ4 φi4 + ψ φi + ψξ φξi + ψη φηi + ψζ φζi ,

(51)

where φξi , φηi and φζi are the partial derivatives at the cell center, and the basis functions read ⎧ 23 11 2 2 2 ⎪ ψ1 = 1 − 17 ⎪ ⎪ 6 (ξ + η + ζ) + 6 (ξη + ξζ + ηζ) + 6 (ξ + η + ζ ) ⎪ ⎪ ⎪ 1 1 1 2 2 2 ⎪ ψ2 = 6 (η + ζ − 5ξ) − 6 (ξη + ξζ + ηζ) + 6 (11ξ − η − ζ ) ⎪ ⎪ ⎪ ⎪ ⎪ 1 1 1 2 2 2 ⎪ ⎪ ψ 3 = 6 (ξ + ζ − 5η) − 6 (ξη + ξζ + ηζ) + 6 (11η − ξ − ζ ) ⎪ ⎪ ⎪ ⎪ 1 1 1 2 2 ⎪ ⎪ ⎨ψ4 = 6 (ξ + η − 5ζ) − 6 (ξη + ξζ + ηζ) + 6 (11ζ − ξ − η2 ) ⎪ ⎪ ⎪ ψ = 5(ξ + η + ζ) − 5(ξη + ξζ + ηζ) − 5(ξ 2 + η2 + ζ 2 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ψξ = ξ − η − ζ + ξη + ξζ + ηζ − ξ 2 + η2 + ζ 2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ψη = η − ξ − ζ + ξη + ξζ + ηζ + ξ 2 − η2 + ζ 2 ⎪ ⎪ ⎪ ⎪ ⎩ψ = ζ − ξ − η + ξη + ξζ + ηζ + ξ 2 + η2 − ζ 2 . ζ

(52)

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

10

moments, ⎧ ⎪ Φi (0, 0, 0) = φi1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Φi (1, 0, 0) = φi2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Φi (0, 1, 0) = φi3 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ Φ ⎪ ⎪  i(0,  0, 1) = φi4 ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ Φi (ξ, η, ζ)dξdηdζ = φi /4 ⎪ ⎪ ω ⎪ ⎪ ⎪ (ξ, η, ζ) ∂Φ i ⎪ ⎪ ( 1 , 1 , 1 ) = φξi ⎪ ⎪ ⎪ 4 4 4 ⎪ ∂ξ ⎪ ⎪ ⎪ ∂Φ (ξ, η, ζ) ⎪ i ⎪ ⎪ ( 1 , 1 , 1 ) = φηi ⎪ ⎪ ⎪ 4 4 4 ∂η ⎪ ⎪ ⎪ ⎪ ∂Φi (ξ, η, ζ) ⎪ ⎪ ⎪ = φζi , ⎩ ( 14 , 14 , 14 ) ∂η

(53)

where the PVs at the vertices, φi j , 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. 3. Numerical results In this section, extensive numerical tests are presented to evaluate the accuracy, robustness and computational cost of the numerical method described above. In order to quantify the numerical errors, we define L1 and L2 errors as follows,   N Ne e 2 i=1 (|φni − φei ||Ωi |) i=1 ((φni − φei ) |Ωi |) E(L1 ) = and E(L2 ) = , (54) Ne Ne 2 i=1 (|φei ||Ωi |) i=1 (φei |Ωi |) where φni and φei denote numerical and exact solutions respectively. 3.1. Advection of a sine wave To evaluate the convergence rate of the advection scheme presented in section 2.4, we computed the advection transport of a sine function as in [19] with gradually refined grids. The advected profile is initially defined as φ0 (x, y) = sin 2π(x + y)

(55)

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 as u = v = 1.0. Numerical experiments were conducted on gradually refined grids generated by Delaunay algorithm with doubly increased resolution. The time step of Δt = 0.001 is used throughout all computations. The L1 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. Table 1: Comparison of L1 errors and convergence rates of different advection schemes.

Cell number 226 894 3588 14412 57518

QUICK 1.916 × 10−1 4.821 × 10−2 1.254 × 10−2 3.605 × 10−3 1.103 × 10−3

Rate – 2.01 1.94 1.79 1.71

Superbee 1.027 × 10−1 7.228 × 10−2 3.327 × 10−2 1.162 × 10−2 3.458 × 10−3

Rate – 0.51 1.12 1.51 1.75 11

MUSCL 2.409 × 10−1 6.548 × 10−2 1.887 × 10−2 5.618 × 10−3 1.682 × 10−3

Rate – 1.89 1.79 1.74 1.74

VPM 1.297 × 10−1 1.772 × 10−2 2.188 × 10−3 2.735 × 10−4 3.522 × 10−5

Rate – 2.90 3.01 2.99 2.95

It is observed that the presented scheme is of 3th order accuracy. Because the multi-moment reconstruction is of second order, the VIA is theoretically expected to be of third order accuracy, which is justified by the numerical output. For comparison, we also include the numerical results from other standard advection schemes used in conventional FVM codes, such as the QUICK scheme [27], superbee-TVD scheme [30] and MUSCL scheme [38, 39], which are available in the OpenFOAM [42, 23] platform. These standard schemes are widely used for unstructured grids and 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 shows a large improvement in numerical accuracy by VPM scheme. On the grid of 57518 elements, for example, the L1 error of VPM is only 3% of the QUICK scheme.

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

X

0.4

0.6

0.8

1

X

Figure 3: The computational meshes for advection test. Left (grid A): a mesh generated by Delaunay algorithm. Right (grid B): a mesh where a region has heavily distorted triangular elements of bad conditions in orthogonality, skewness and uniformity.

3.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   (x − 0.5)2 + (y − 0.75)2 φ0 (x, y) = exp (56) 0.015 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. Figs. 5, 6 and 7 display the corresponding outputs of other standard schemes on the two grids. Compared against Fig.4, it is found 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 L2 errors on both A and B grids for VPM method and other standard finite volume schemes with the “degradation rate” defined by Degradation rate =

L2 error on grid B − 1. L2 error on grid A 12

(57)

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

exact

0.2

exact

0.2

numrical

0

0

0.2

0.4

0.6

0.8

numrical

0

1

0

0.2

0.4

X

0.6

0.8

1

X

Figure 4: The numerical results of VPM scheme for the solid rotation advection tests on grid A (left) and grid B (right).

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. Table 2: Numerical errors on different grids and degradation rate.

Scheme VPM QUICK Superbee MUSCL

L2 error on grid A 0.07220 0.11919 0.13150 0.13588

L2 error on grid B 0.09014 0.24073 0.24995 0.32268

Degradation rate 0.25 1.02 0.90 1.37

We can see that having the smallest degradation rate, the VPM scheme is the most robust one in 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. 3.3. Taylor vortex problem To evaluate the accuracy of the whole incompressible fluid solver, we computed the Taylor vortex problem [32, 12]. An analytical solution for this problem is available as   2π2 t , (58) u(x, y, t) = − cos(πx) sin(πy) exp − Re   2π2 t v(x, y, t) = sin(πx) cos(πy) exp − , (59) Re   4π2 t 1 p(x, y, t) = − (cos(2πx) + sin(2πy)) exp − , (60) 4 Re where Re is the Reynolds number. 13

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

exact

0.2

exact

0.2

numrical

0

0

0.2

0.4

0.6

0.8

numrical

0

1

0

0.2

X

0.4

0.6

0.8

1

X Figure 5: Same as Fig.4 but by the QUICK scheme.

We computed this test problem over [−1, 1] × [−1, 1] domain with the periodical boundary condition, using grids of different resolutions generated by Delaunay algorithm. We use Re = 100 in this example. 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 2.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 (order) higher than 2.5. In the results on the refined mesh (57670 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. This benchmark problem also provides a test bed to evaluate kinetic energy change of numerical solution[16]. The total kinetic energy in the computational domain depends on Reynolds number and varies in time as  1 1 4π2 t 1 1 ). (61) (u2 + v2 )dxdy = exp(− Kexact (t) = 2|Ω| −1 −1 4 Re The total kinetic energy of the numerical solution is obtained by   1  u¯ i (t)2 + v¯ i (t)2 Knumer (t) =  |Ωi | . |Ωi | i 2

(62)

In order to see the numerical dissipation of the model, we computed a “nearly inviscid” case with a Reynolds number of 2 × 106 on different grid resolutions and show the total kinetic energy with respect to time in Fig.9 (left). The numerical results indicate that the VPM incompressible fluid model is stable and the numerical dissipation can be effectively reduced by using refined grid resolution. We also experimented with different Reynolds numbers on a 57670-cell grid. It is shown in Fig.9 (right) that the present method is able to accurately reproduce the kinetic energy for a wide range of Reynolds numbers. The proposed method thus represents a promising numerical framework for accurate and robust RANS simulations. Although a thorough implementation of RANS models which makes full use of both VIA and PV is a new topic 14

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

exact

0.2

exact

0.2

numrical

0

0

0.2

0.4

0.6

0.8

numrical

0

1

0

0.2

X

0.4

0.6

0.8

1

X Figure 6: Same as Fig.4 but by the superbee TVD scheme.

for us to explore in future, a practical implementation can be conducted by using the existing RANS models to VIA variables first under the conventional FVM framework, and update the PV moments through the TEC formula which are proved to be consistent and stable for any forcing term in the governing equations. 3.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 used in the conventional FVM. It also leads to an increase in CPU consumption. As described before, the spatial discretization can be computed efficiently. Being the computational variables, the PVs at cell vertices are always ready for use in the calculation of the numerical fluxes, which largely simplifies the operations and saves CPU time. 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.10GHZ). 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%. Shown in Table 4, the elapse time increase may reach 1.5 fold in comparison with the QUICK scheme. As shown in §3.1, large improvement on numerical accuracy can be obtained by VPM scheme. From Table 1, we see that the numerical error of VPM on a 5463-element grid is smaller then that of the QUICK scheme on a 14412-element grid. It shows 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. The elapse time is measured by integrating the incompressible fluid model for 1000 steps with a small time step interval. In order to compare the VPM model with the QUICK conventional FVM model in terms of both numerical accuracy and computational cost, we also include the numerical errors obtained in the previous test. Shown in Table 5, about 34% increase is seen in the memory requirement due to 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 104 . 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 common understanding that the pressure Poisson equation consumes a large portion of the 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 15

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

exact

0.2

exact

0.2

numrical

0

0

0.2

0.4

0.6

0.8

numrical

0

1

0

0.2

X

0.4

0.6

0.8

1

X Figure 7: Same as Fig.4 but by the MUSCL scheme.

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. 3.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. [14], 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 shown in Fig.10. The Reynolds number used in the following four test cases is 1000. 3.5.1. Lid driven square cavity problem This is the most widely used benchmark test for numerical codes designed to solve viscous incompressible flows[14]. 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 = 1 is moving at a constant speed (u = 1, v = 0). Numerical tests were carried out with triangular unstructured grids (upper-left panel of Fig.10) of different resolutions. It is known that a −9 − pm steady solution is reached with this Reynolds number, which is identified by the criterion of max |pm+1 i i | < 10 . We plot the numerical results of different grid resolutions in Fig.11. 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.12 for comparison. It is obvious that the present model is more accurate than the conventional FVM. 3.5.2. Lid driven forward-step cavity problem A cavity with a forward-step shown in the upper-right panel of Fig.10 is used in [28] to evaluate numerical solvers developed in general coordinates. Numerical results are given in Fig.13 with the reference velocity profiles provided in [28] which were computed on a 512 × 512 curvilinear grids. Shown in Fig.13 (middle and right panels), the present model is able to reproduce adequate solutions even with much coarser mesh resolutions. The solutions on meshes of more than 2704 cells agree well with the reference solution.

16

Table 3: Numerical errors and convergence rates of VPM and conventional FVM models for Taylor vortex problem.

Method VPM

FVM (QUICK)

FVM (Superbee)

FVM (MUSCL)

Cell number 902 3604 14362 57670 902 3604 14362 57670 902 3604 14362 57670 902 3604 14362 57670

L2 error of |u| 7.842 × 10−3 1.409 × 10−3 2.183 × 10−4 3.598 × 10−5 3.744 × 10−2 1.443 × 10−2 6.987 × 10−3 3.029 × 10−3 4.322 × 10−2 1.704 × 10−2 7.188 × 10−3 3.027 × 10−3 3.778 × 10−2 1.569 × 10−2 7.283 × 10−3 3.087 × 10−3

Order in |u| – 2.48 2.70 2.59 – 1.38 1.05 1.20 – 1.34 1.25 1.24 – 1.27 1.11 1.23

L2 error of p 2.361 × 10−2 5.523 × 10−3 1.173 × 10−3 3.015 × 10−4 3.284 × 10−1 1.423 × 10−1 6.710 × 10−2 4.147 × 10−2 3.581 × 10−1 1.469 × 10−1 6.756 × 10−2 4.113 × 10−2 3.619 × 10−1 1.469 × 10−1 6.609 × 10−2 4.083 × 10−2

Order in p – 2.10 2.24 1.95 – 1.21 1.09 0.69 – 1.29 1.12 0.71 – 1.30 1.16 0.69

Table 4: Computational cost estimation for advection calculation of VPM in comparison with the QUICK scheme.

Elements 226 894 3588 14412 57518

DOFs 226 894 3588 14412 57518

QUICK Time(s) 0.21 0.35 0.84 2.66 9.95

L1 error 1.916 × 10−1 4.821 × 10−2 1.254 × 10−2 3.605 × 10−3 1.103 × 10−3

DOFs 360 1382 5463 21779 86598

VPM Time(s) 0.22 0.55 1.67 6.27 25.95

L1 error 1.297 × 10−1 1.772 × 10−2 2.188 × 10−3 2.735 × 10−4 3.522 × 10−5

3.5.3. Lid driven skew cavity problem Numerical examples of viscous flow in skew cavities (shown in the bottom-left panel of Fig.10) with inclined angles of 30o and 45o are found in [11] and [28]. Our numerical results are shown in Fig 14 and Fig 15. The reference solution in [11] is computed on 320 × 320 structured grids. In the case of 30o , the numerical solution converges to the reference solution with a much coarser grid compared to the case of 45o , which is common to other existing works in the literature. Even though the 45o cavity flow requires a relatively higher resolution mesh for converged solution, it is still much less in comparison with other conventional FVM schemes. 3.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.10) which has been extensively investigated in [15]. We use the results in [15] for comparison which were computed on a triangular mesh with a representative size of 1/128. Shown in Fig.16, our results converge to the reference results even with a relatively low mesh resolution. 3.6. 2D viscous flow around a circular cylinder In order to validate the capability of present model to solve time-dependent problem, we simulate the viscous flow with Re = 100 passing a circular cylinder as proposed in [35, 36]. For this Reynolds number, the flow is time dependent and characterized by periodical vortex shedding in the downstream region of the cylinder. It provides a good test bed of transient flow to evaluate numerical models. 17

-3

-1

1.30

1.34 1.38

1.25

ln(L2 Error)

2.48

1.20

2.70

-8

0.69

-4

2.10

-5 2.24

-6 2.59

-9

-7

VPM QUICK superBee vanLeer

-10 -11

0.69

-3

1.24

-6 -7

1.16 1.09

1.05

-5

1.21

-2

ln(L2 Error)

-4

3.5

4

1.95

VPM QUICK superBee vanLeer

-8

4.5

5

-9

5.5

3.5

ln(√Elements)

4

4.5

5

5.5

ln(√Elements)

Figure 8: The numerical errors and convergence rates for velocity (left) and pressure (right).

Table 5: Computation cost and numerical accuracy for incompressible fluid solvers.

Elements 902 3604 14362 57670

DOFs 2706 10812 43086 173010

Time(s) 4.24 12.43 47.78 226.32

QUICK L2 error of |u| 3.74 × 10−2 1.44 × 10−2 6.99 × 10−3 3.03 × 10−3

L2 error of p 3.28 × 10−1 1.42 × 10−1 6.71 × 10−2 4.15 × 10−2

DOFs 3690 14578 57770 231322

Time(s) 5.03 16.27 68.25 317.60

VPM L2 error of |u| 7.84 × 10−3 1.40 × 10−3 2.18 × 10−4 3.60 × 10−5

L2 error of p 2.36 × 10−2 5.52 × 10−3 1.17 × 10−3 3.01 × 10−4

The computational condition was the same as in [36] (see case III therein for details). We use an unstructured grid with 81394 triangle elements locally refined near the cylinder. Fig.17 shows the contours of velocity norm, pressure and vorticity fields at the instant when the lift force reaches the maximum. We can clearly observe the Karman vortex street behind the cylinder. All the fields recover the major features observed in [36]. The time histories of drag and lift coefficients are plotted in Fig.18 which look in agreement with other standard solutions. An enlarged view of drag and lift coefficients are plotted in Fig.19. The periodical oscillations are accurately reproduced, and the amplitudes remain perfectly constant. To quantitatively verify the numerical solutions, we computed the maximum drag coefficient (Cd )max , the maximum lift coefficient (C L )max and the Strouhal number which are summarized in Table 6 with the reference data from [35] and [36]. It is seen that the present results are among the most accurate solutions of the benchmark test.

Present model Tu and Aliabadi[36] Turek and Shahrouz[35]

Lower bound Upper bound

(Cd )max 3.23 3.215 3.22 3.24

(C L )max 1.01 1.0374 0.99 1.01

St 0.298 0.2959 0.295 0.305

Table 6: Comparison between the present solution and the reference solutions.

3.7. 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 [2], we imposed a constant velocity (v(−0.5, y, z) = 1) on surface x = −0.5 and computed 18

0.26 exact result with Re=2000000 numerical result with 3604 cells numerical result with 14362 cells numerical result with 57670 cells

0.24

K

0.25

K

0.22

0.248 exact result with Re=200 exact result with Re=2000 exact result with Re=20000 numerical result for Re=200 numerical result for Re=2000 numerical result for Re=20000

0.2

0.246

0.18 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

t

t

Figure 9: The total kinetic energy of Taylor vortex benchmark test with a Reynolds number of 2 × 106 on different grid resolutions (left), and the total kinetic energy changes with different Reynolds numbers on a 57670-cell mesh (right).

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.20. For comparison, we also include the numerical results in [2] 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.8. Viscous flows passing a sphere with different Reynolds numbers As the final test, we further verified the present model by simulating viscous flows passing a sphere with different Reynolds numbers, which is widely used to evaluate numerical codes for incompressible viscous flows [24, 26]. The computational domain is set as 15 × 10 × 10 in streamwise (x), spanwise (y) and vertical(z) directions respectively. A unit sphere is located at (5, 5, 5). The computational grid is composed of 661780 tetrahedral elements. We set the inlet flow velocity u∞ to be unit, The Reynolds number is adjusted by altering the viscosity. The time step is variable to maintain a constant CFL number of 0.5. For small Reynolds numbers (Re= 50,100,150 and 200), we reach a steady axisymmetric flow which is realized when the variation in pressure is less than 10−9 in our numerical tests. Streamlines on (x, z)-plane cutting through the sphere center are shown in Fig.21 for Reynolds numbers of 50,100,150 and 200, which visually agree well with those in [24]. The quantitative comparison has been also conducted with the available data in [24]. Shown in Fig.22, we plot the polar separation angle, separation length, vortex position and drag coefficient against the Reynolds numbers up to 200 and compare our results with the existing numerical and experimental results summarized in [24]. It reveals that the present method is in good agreement with the well-accepted results for these steady viscous flows.

Present model Kim et al. [26] Johnson and Patel [24] Constantinescu and Squires[5]

Re 250 300 250 300 250 300 250 300

CD 0.71 0.659 0.701 0.657 0.70 0.656 0.70 0.655

CL 0.062 0.069 0.059 0.067 0.062 0.069 0.062 0.065

Table 7: Numerical results of viscous flows over a sphere.

19

St 0.135 0.134 0.137 0.136

0.8

0.8

0.6

0.6

Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.2

0.4

X

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

X 0

0.6

-0.2

Y

Y

0.4

0.2

-0.4

0

0

0.5

X

1

0

1.5

0.2

0.4

X

0.6

0.8

1

Figure 10: Geometrical configurations and computational meshes for 2D incompressible flow benchmark tests.

It is known that with the Reynolds number increasing, flow turns into a steady nonaxisymmetric flow regime (210 ≤ Re ≤ 270), and then an unsteady flow regime (Re ≥ 280). We computed two cases with the Reynolds number of 250 and 300 to reproduce these two flow regimes. The numerical results are given in Table 7. The drag and lift coefficients for both cases are in excellent agreement with the existing researches [24, 5, 26]. As expected a steady flow is reproduced in the case of Re = 250. A periodical flow is observed in the case of Re = 300. The Strouhal number at Re = 300 agrees well with the reported results aforementioned. Four snapshots of the instant verticalcomponent vorticity field on the x − y cross section cutting through the sphere center are plotted in Fig.23. We set the start point φ = 0 to the instant when the lift force reaches the minimum value. It is found that the vorticity structure evolves in time and looks quite similar to those reported in the existing researches [24]. 4. Conclusion remarks We have presented and evaluated a novel numerical solver for incompressible 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. 20

1

0.4

0.8

0.2

0.6

V

Y

0

0.4

-0.2

0.2

0 -0.4

226 triangles 894 triangles 5612 triangles Ghia

-0.2

0

0.2

0.4

0.6

U

0.8

226 triangles 894 triangles 5612 triangles Ghia

-0.4

1

0

0.2

0.4

0.6

0.8

1

X

Figure 11: Numerical results of square cavity flow problem. Displayed are stream line on a 5612-cell mesh (left), u profiles along x = 0.5 line (middle) and v profiles along y = 0.5 line (right) on grids of different resolutions.

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 practice. 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. Acknowledgment This work was supported in part by JSPS KAKENHI(24560187).

21

1

0.4

0.8

0.2

0.6

V

Y

0

0.4

-0.2 present scheme QUICK superBee vanLeer Ghia

0.2

0 -0.4

-0.2

0

0.2

0.4

0.6

0.8

present scheme QUICK superBee vanLeer Ghia

-0.4

1

0

0.2

0.4

U

0.6

0.8

1

X

Figure 12: Numerical results of square cavity flow problem on a 5612-cell mesh. Displayed are the u profile along x = 0.5 line (left) and v profile along y = 0.5 line (right).

0.8

0.2

0.6

0

V

0.4

Y

1

-0.2

0.4

-0.4 0.2

0

682 triangles 2704 triangles 7538 triangles Oosterlee

-0.4

-0.2

0

0.2

U

0.4

0.6

0.8

682 triangles 2704 triangles 7538 triangles Oosterlee

-0.6 1

0

0.2

0.4

0.6

0.8

1

X

Figure 13: Numerical results of forward-step cavity flow problem. Displayed are stream line on a 7538-cell mesh (left), u profile along x = 0.75 line (middle) and v profile along y = 0.75 line (right).

22

1 0.9 0.05 0.8 0.7 0

0.5

V

Y

0.6

-0.05 0.4 0.3 -0.1 0.2

132 triangles 460 triangles 1872 triangles Demirdzic

0.1

132 triangles 460 triangles 1872 triangles Demirdzic

-0.15 0

-0.2

0

0.2

0.4

0.6

0.8

1

0.2

0.4

U

0.6

0.8

1

X

Figure 14: Numerical results of a 30o inclined skew cavity flow problem. Displayed are stream line on a 1872-cell mesh (left), u profile along y = tan(30o )(x − 0.5) line (middle) and v profile along y = tan(30o )/2 line (right).

1 0.05 654 triangles 2634 triangles 7352 triangles Demirdzic

0.9 0.8 0.7 0

0.5

V

Y

0.6

0.4 0.3 0.2

0

-0.05

654 triangles 2654 triangles 7352 triangles Demirdzic

0.1 -0.2

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

U

0.8

1

X

Figure 15: Numerical results of a 45o inclined skew cavity flow problem. Displayed are stream line on a 7352-cell mesh (left), u profile along y = tan(45o )(x − 0.5) line (middle) and v profile along y = tan(45o )/2 line (right).

0

0.2

-0.1

-0.2

V

Y

0

-0.3 -0.2

-0.4

-0.5

92 triangles 372 triangles 2262 triangles Glowinski

-0.4

-0.2

0

0.2

0.4

U

0.6

0.8

92 triangles 372 triangles 2262 triangles Glowinski

-0.4

1

0.2

0.4

0.6

0.8

X

Figure 16: Numerical results of semi-circular cavity flow problem. Displayed are stream line on a 2266-cell mesh (left), u profile along x = 0.5 line (middle) and v profile along y = −0.25 line (right).

23

Figure 17: The physical fields at the instant when the lift force becomes largest. Displayed are velocity norm (top), pressure (middle) and vorticity (bottom) fields.

4

Cd 3

2

1

CL 0

-1 0

2

4

6

8

Time(sec) Figure 18: Time history of drag and lift coefficient.

24

10

3.24 1

3.22

0.5

CL 0

Cd 3.2

-0.5

3.18

-1 3.16 6

6.2

6.4

5.5

6

Time(sec)

6.5

Time(sec)

Figure 19: Zoomed view of time history of drag (left) and lift (right) coefficient.

1 0.4

7660 tetrahedrals 58783 tetrahedrals 411202 tetrahedrals Albensoeder

0.8

7660 tetrahedrals 58783 tetrahedrals 411202 tetrahedrals Albensoeder

0.3

0.6 0.2

Y

V

0.4

0.1

0.2

0

0

-0.1

-0.2

-0.2 -0.4

-0.2

0

0.2

0.4

-0.4

X

-0.2

0

0.2

0.4

U

Figure 20: Numerical results of 3D lid driven cubic cavity. Displayed are u profile along centerline (x, 0, 0) (left) and v profile along centerline (0, y, 0) (right). The reference results of [2] was computed by a Chebyshev-collocation method with 96 × 96 × 64 Gauss-Lobatto points.

25

Figure 21: Streamlines on the central cross section of flows passing a sphere. Displayed are the steady cases with Reynolds numbers of 50 (top-left), 100 (top-right), 150 (bottom-left) and 200 (bottom-right).

26

180 1.5 Johnson et al. Pruppacher et al. Taneda Present results

160

xs

Θs

1

140 Johnson et al. Taneda Tomboulides Magnaudet et al. Present results

0.5

120

0

50

100

150

0

200

0

50

Re

100

150

200

Re

1

3

2.5

0.8

Johnson et al. Roos & Willmarth Present results

xc 2

Cd

Johnson et al. Taneda Present results

0.6

1.5

0.4 1

yc 0.2

0

0.5

0

50

100

150

200

Re

0

0

50

100

150

200

Re

Figure 22: Flow geometry as a function of the Reynolds number. Displayed are the polar separation angle Θ s (top-left), the separation length x s (top-right), the vortex position (xc , yc ) (bottom-left) and the drag coefficient Cd (bottom-right). See Figs. 4 and 5 in [24] as well as the references therein for the reproduced data.

Figure 23: The contours of ωz at Re=300 for a quarter period on xy-plane at φ = 0 (top-left), φ = π/2 (top-right), φ = π (bottom-left) and φ = 3π/2 (bottom-right). Contours are drawn at 0.5 increment with the solid and dot lines for positive and negative values respectively.

27

Appendix A. Computation of the average slope at cell center The first-order derivatives (slope) at the centroid of cell Ωi , φξi and φηi , are computed by the following steps. ˆ i j (ξ, η), j = 1, 2, 3, over cell Ωi and its three ˆ i (ξ, η) and Φ • Construct the piecewise interpolation functions, Φ neighboring cells Ωi j , j = 1, 2, 3 , by

and

ˆ i (ξ, η) = ψˆ 1 φi1 + ψˆ 2 φi2 + ψˆ 1 φi3 + ψˆ φi Φ

(A.1)

ˆ i j (ξ, η) = ψˆ 1 φi j1 + ψˆ 2 φi j2 + ψˆ 1 φi j3 + ψˆ φi j , j = 1, 2, 3; Φ

(A.2)

where the basis functions are

⎧ ⎪ ⎪ ψˆ 1 = 1 − ξ − η − 4ξη ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ψ ⎨ ˆ 2 = ξ − 4ξη ⎪ ⎪ ⎪ ψˆ 3 = η − 4ξη ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ψˆ = 12ξη

(A.3)

Note that we only use the VIA and PV in reconstructions (A.1) and (A.2) . • The derivatives at the middle point θ˜i j of edge Γi j are obtained by averaging the values computed separately from (A.1) and (A.2) for the two neighbor cells. ⎛ ⎞ ⎛   ˆ i j ∂yi j ∂Φ ˆ i j ⎞⎟⎟ ⎟⎟ ˆ i ∂yi ∂Φ ˆi 1 ⎜⎜⎜⎜ 1 ∂yi ∂Φ 1 ⎜⎜⎜ ∂yi j ∂Φ ⎟ ⎟ ˜ ⎜ ⎟ φ x (θi j ) = ⎝⎜ − − + ⎝ ⎠ ⎟⎟ , 2 |Ji | ∂η ∂ξ ∂ξ ∂η θ˜i j |Ji j | ∂η ∂ξ ∂ξ ∂η θ˜i j ⎠ ⎛ ⎞ (A.4) ⎛   ˆ i j ∂xi j ∂Φ ˆ i j ⎞⎟⎟ ⎟⎟ ˆ i ∂xi ∂Φ ˆi 1 ⎜⎜ 1 ∂xi ∂Φ 1 ⎜⎜⎜ ∂xi j ∂Φ ⎟⎟⎠ ⎟⎟⎟ . ⎜⎝ − − + φy (θ˜i j ) = ⎜⎜⎜⎝ 2 |Ji | ∂ξ ∂η ∂η ∂ξ θ˜i j |Ji j | ∂ξ ∂η ∂η ∂ξ θ˜i j ⎠ • Given the derivatives at the middle points of the three edges, φ x (θ˜i j ) and φy (θ˜i j ) ( j = 1, 2, 3), the derivatives at the cell center, φ xi and φyi are computed by a simple averaging, 1 φ x (θ˜i j ), 3 j=1 3

φ xi =

1 ˜ φyi = φy (θi j ). 3 j=1 3

(A.5)

The derivatives at the cell center in the local coordinate can be immediately obtained by (8). The cell averaged derivatives computed from (A.5) can be further modified to suppress the numerical oscillations as in [46]. Appendix B. Time-evolution converting (TEC) formula Given the time variation in VIA moment, the time variation in PV moment can easily obtained via the Timeevolution converting (TEC) formula, which are in fact the spatial interpolation of time variation of a physical variable [43, 44]. In unstructured grids, we make use of the least square minimization technique described below. Provided the time variation of VIA δt φk , k = 1, 2, · · · , K, at the centers of cell union which share vertex θi j , the increment of PV at θi j , δt φi j , is computed by least square minimization [36] as ⎤⎡ ⎤ ⎡  ⎤ ⎡   δφ ⎥ ⎢⎢⎢ K x y ⎥⎢ δ φ ⎥ ⎢ k 2k  k k ⎥⎥⎥⎥ ⎢⎢⎢⎢ t i j ⎥⎥⎥⎥ ⎢⎢⎢⎢  k t k ⎥⎥⎥⎥ ⎢⎢⎢  ⎥ x y ⎥⎢ δ φ ⎥ ⎢ (B.1) ⎢⎢⎣ k xk  k xk k k 2 k ⎥⎥⎦ ⎢⎢⎣ t xi j ⎥⎥⎦ = ⎢⎢⎣ k δt φk xk ⎥⎥⎦ δt φyi j k yk k xk y k k yk k δ t φ k yk 28

where (xk , yk ) is the centroid of cell k. Not only δt φi j , δt φ xi j and δt φyi j can be also obtained from (B.1). The TEC formula is in fact an interpolation for an operator rather than a physical variable. It applies to any “forcing” term in the equations, thus it largely simplifies the computations. The TEC formula also allows the VPM model to easily makes use of the existing modules in the finite volume framework. We denote a TEC formula that computes δt φi j from δt φk by δt φi j = T EC(δt φk ).

29

(B.2)

[1] R. Akoh, S. Ii and F. Xiao, A CIP/multi-moment finite volume method for shallow water equations with source terms, Int. J. Numer. Method in Fluid, 56 (2008) 2245-2270. [2] S. Albensoeder and H.C. Kuhlmann Accurate three-dimensional lid-driven cavity flow, J. Comp. Phys., 206 (2005) 536-558. [3] T. J. Barth and P. O. Frederickson, High-Order Solution of the Euler Equations on Unstructured Grids using Quadratic Reconstruction, AIAA Paper No. 90-0013 (1990). [4] C.G. Chen and F. Xiao, Shallow water model on spherical-cubic grid by Multi-moment finite volume method, J. Comput. Phys. 227 (2008) 5019-5044. [5] G.S.Constantinescu and K.D.Squires, LES and DES Investigations of Turbulent Flow over a Sphere, AIAA Paper 2000-0540 (2000). [6] A.J. Chorin, ”Numerical Solution of the Navier-Stokes Equations”, Math. Comp., 22 (1968), 745-762. [7] B. Cockburn and C.W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Math. Comput. 52 (1989) 411-435. [8] B. Cockburn, S.Y. Lin and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems, J. Comput. Phys. 84 (1989) 90-113. [9] B. Cockburn, S. Hou and C.W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case, Math. Comput. 54 (1990) 545-581. [10] B. Cockburn and C.W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V: The multidimensional systems, J. Comput. Phys. 141 (1998) 199-224. ˘ Lilek and M. Peri´c, Fluid flow and heat transfer test problems for non-orthogonal grids: bench-mark solutions, Int. J. Numer. [11] I. Demird˘zi´c, Z. Method in Fluid, 15 (1992) 329-354. [12] E. Ferrer and R. Willden, A high order discontinuous Galerkin finite element solver for the incompressible Navier-Stokes equations, Comput. and Fluids 46 (2011) 224-230. [13] O. Friedrich, Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids, J. Comput. Phys. 144 (1998) 194-212. [14] U. Ghia, K.N. Ghia and C.T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comp. Phys., 48 (1982), 387-411. [15] R. Glowinski, G. Guidoboni and T.W. Pan, Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity, J. Comput. Phys., 216 (2006) 76-91. [16] F. Ham and G. Iaccarino, G. 2004 A hybrid vertex-centered finite volume/element method for viscous incompressible flows on non-staggered unstructured meshes, Annual Research Briefs, Center for Turbulence Research, Stanford University and NASA-Ames, 2004, 3-14. [17] F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface, Physics of Fluids, 8 (1965) 2182-2189. [18] J. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods, Springer, 2008. [19] C. Hu and C.W. Shu, Weighted essentially non-oscillatory schemes on triangular meshes, J. Comput. Phys. 150 (1999) 97-127. [20] S. Ii, M. Shimuta and F. Xiao, A 4th-order and single-cell-based advection scheme on unstructured grids using multi-moments, Comput. Phys. Comm. 173 (2005) 17-33. [21] S. Ii and F. Xiao, CIP/multi-moment finite volume method for Euler equations: A semi-Lagrangian characteristic formulation, J. Comput. Phys. 222 (2007) 849-871. [22] R.I. Issa, Solution of Implicitly Discretized Fluid Flow Equations by Operator Splitting, J. Comput. Phys., 62 (1986) 40-65. [23] H. Jasak, OpenFOAM: Open source CFD in research and industry, Int. J. Naval Archit. and Ocean engin. 1 (2009) 89 -94. [24] T.Johnson and V.C. Patel, Flow past a sphere up to Reynolds number of 300, J.Fluid Mech., 378 (1999) 19-70. [25] F. Jureti´c and A.D. Gosman, Error analysis of the finite-volume method with respect to mesh type, Numerical Heat Transfer (Part B), 57 (2010) 414-439. [26] J. Kim, D. Kim and H. Choi, An Immersed-Boundary Finite-Volume Method for Simulations of Flow in Complex Geometries, J. Comput. Phys., 171 (2001) 132-150. [27] B.P. Leonard, A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Methods Appl. Mech. Eng., 19 (1979) 59-98. [28] C.W. Oosterlee, P. Wesseling and E. Brakkee, Benchmark solutions for the incompressible Navier-Stokes equations in general co-ordinates on staggered grids, Int. J. Numer. Method in Fluid, 17 (1993) 301-321. [29] S.V. Patankar and D.B. Spalding, A calculation procedure for heat and mass transfer in three-dimensional parabolic flows. Int. J. Heat and Mass Transfer, 15 (1972) 1787-1806. [30] P.L. Roe, Some contributions to the modeling of discontinuous flows, Lectures in Applied Mathematics, 22 (1985), Springer Verlag, 163-193. [31] R. Rossi, Direct numerical simulation of scalar transport using unstructured finite-volume schemes, J. Comput. Phys 229 (2009) 1639-1657. [32] K. Shahbazi, P. Fischer and C. Ethier, A high-order discontinuous Galerkin method for the unsteady incompressible Navier-Stokes equations, J. Comput. Phys 222 (2007) 391-407. [33] C.W. Shu, Total-variation-diminishing time discretizations, SIAM J. Sci. Stat. Comput. 9 (1988) 1073-1084. [34] D.T. Steinmoeller, M. Stastna and K.G. Lamb, A short note on the discontinuous Galerkin discretization of the pressure projection operator in incompressible flow, J. Comput. Phys. 251 (2013) 480-486. [35] S. Turek and M. Schafer, Benchmark computations of laminar flows around a cylinder. In Flow Simulation with High-performance Computers II, Hirschel EH (ed.), Notes on Numerical Fluid Mechanics, 52 (1996) 547-566. [36] S. Tu and S. Aliabadi, Development of a hybrid finite volume/element solver for incompressible flows. Int. J. Numer. Meth. Fluids, 55 (2007), 177-203. [37] J.P. Van Doormaal and G.D. Raithby, Enhancements of the SIMPLE method for predicting incompressible fluid flows, Num. Heat Transfer, 7 (1984) 147-163. [38] B. Van Leer, Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys., 23 (1977) 276-299.

30

[39] B. Van Leer, Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov’s Method, J. Comp. Phys., 32(1979) 101-136. [40] Z.J. Wang, Spectral (finite) volume method for conservation laws on unstructured grids: basic formulation, J. Comput. Phys. 178 (2002) 210-251. [41] Z.J. Wang, Y. Liu, Spectral (finite) volume method for conservation laws on unstructured grids II: extension to two-dimensional scalar equation, J. Comput. Phys. 179 (2002) 665-697. [42] H.G. Weller, G.Tabor, H.Jasak and C.Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Computers in Physics 12 (1998) 620-631. [43] F. Xiao, Unified formulation for compressible and incompressible flows by using multi integrated moment method: one-dimensional inviscid compressible flow, J. Comput. Phys. 195 (2004) 629-654. [44] F. Xiao, R. Akoh and S. Ii, Unified formulation for compressible and incompressible flows by using multi-integrated moments II: Multidimensional version for compressible and incompressible flows, J. Comput. Phys. 213 (2006) 31-56. [45] F. Xiao, A. Ikebata and T. Hasegawa, Numerical simulations of free-interface fluids by a multi integrated moment method, Comput. Struct. 83 (2005) 409-423. [46] F. Xiao and T. Yabe, Completely conservative and oscillation-less semi-Lagrangian schemes for advection transportation, J. Comput. Phys. 170 (2001) 498-522. [47] T. Yabe and T. Aoki, A universal solver for hyperbolic-equations by cubic-polynomial interpolation. 1. One-dimensional solver, Comput. Phys. Commun. 66 (1991) 219-232. [48] T. Yabe, F. Xiao and T. Utsumi, The constrained interpolation profile method for multiphase analysis, J. Comput. Phys. 169 (2001) 556-593.

31