High-order multi-moment finite volume method with smoothness adaptive fitting reconstruction for compressible viscous flow

Journal of Computational Physics 394 (2019) 559–593

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

High-order multi-moment finite volume method with smoothness adaptive fitting reconstruction for compressible viscous flow Bin Xie a,∗ , Xi Deng b , ShiJun Liao a,∗ , Feng Xiao b,∗ a b

School of Naval Architecture, Department of Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai, 200240, China School of Engineering, Department of Mechanical Engineering, Tokyo Institute of Technology, Tokyo, 152-8550, Japan

a r t i c l e

i n f o

Article history: Received 10 October 2018 Received in revised form 2 June 2019 Accepted 3 June 2019 Available online 10 June 2019 Keywords: Navier–Stokes equations Unstructured grid Compressible flow Multi-moment reconstruction Numerical oscillation Numerical dissipation

a b s t r a c t We proposed a novel multi-moment finite volume method (MMFVM) to compute compressible viscous flows on arbitrary unstructured grids. Unlike the previous volumeaverage/point-value multi-moment (VPM) method (Xie et al. 2014 [60]), the present method, so-called VPM-Ex (VPM with extended-stencil reconstruction), formulates novel spatial reconstructions with substantially improved robustness and accuracy for both smooth and discontinuous solutions. Our efforts in this work have brought at least the following new features to the numerical formulation: 1) a more robust and accurate spatial reconstruction based on the stencil of all immediately adjacent cells that supports higher-order limiting projection; 2) two different options for both explicit and implicit reconstructions; 3) a simplified and more accurate formulation for point-value. In addition, we have proposed a smoothness-adaptive fitting (SAF) method for limiting projection to improve the solution quality for discontinuities. The resulting reconstruction method that combines VPM-Ex scheme with SAF-based limiting projection schemes is able to attain 3rd-order accuracy for smooth solution, and non-oscillatory and less-dissipative numerical results for discontinuous solution on unstructured grids of different types of elements. The numerical method has been extensively verified by various benchmark tests for scalar conservation equation, Euler equations and Navier–Stokes equations in 2 and 3 dimensions. The numerical results substantiate that the present method provides a novel numerical formulation of great practical significance to accurately resolve both continuous and discontinuous solutions. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Motivated by the needs from applications, developing high-order (beyond 2nd order) schemes on unstructured grids to solve compressible flows in presence of both shock waves and other flow structures remains a very active research field. Since high-order schemes have less truncation errors, they can generally converge faster and be more efficient than lower-order schemes to reach the same level of errors in numerical solutions. In spite of tremendous efforts made in the past decades, high-order numerical schemes on unstructured grids that satisfy the important requirements from application perspective, such as computational efficiency, stability and robustness in presence of shocks and other discontinuities, adap-

*

Corresponding author. E-mail addresses: [email protected] (B. Xie), [email protected] (S.J. Liao), [email protected] (F. Xiao).

https://doi.org/10.1016/j.jcp.2019.06.002 0021-9991/© 2019 Elsevier Inc. All rights reserved.

560

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

tivity to geometrical complexity, as well as amenableness to complex physical processes, are far from well-established. As a matter of fact, nearly all commercial codes of computational fluid dynamics (CFD) still dominantly rely on second-order methods which are widely known to be excessively dissipative for resolving flow structures different scales, like vortices, eddies and turbulence, which are essential for many applications. The major stumbling blocks that prevent high-order schemes from being widely adopted are attributed to algorithmic complexities and lack of numerical robustness particularly for discontinuous solutions. Under the conventional finite volume framework, how to represent the physical field with high-order polynomials is far from trivial on unstructured grids. The k−exact least-square method, pioneered by Barth and Frederickson [2,3] and developed in [11] and others, forms a standard approach to construct polynomials of degree k. A 2D k−exact polynomial, for example, has (k + 1)(k + 2)/2 coefficients (degrees of freedom), and thus needs such a number or more surrounding cells which make up of the stencil of support to determine the polynomial. Given k, the topological connection among the cells in the support stencil depends on the shape of mesh elements in an unstructured grid. For a high-order (k ≥ 3) polynomial, a stencil that only includes the immediate neighbors of the target cell might be not enough, and the stencil needs to be extended to involve the farther distanced cells. Another important issue for spatial reconstruction is the limiting projection, which is necessary to suppress spurious oscillations in presence of discontinuous solutions, such as shock waves in compressible flow. The most representative and popularly used high-order limiting projection methods for both smooth and non-smooth solutions are the essentially non-oscillatory (ENO) [21] method and the weighted-ENO (WENO) [29,25] scheme, where polynomials over sub-stencils are combined to realize high-order polynomial approximation for smooth profile and to suppress spurious oscillation in presence of large jump of discontinuity. Although a WENO-type method requires a stencil, the implementation on structured grids is straightforward. However, it is not a trivial task to apply such limiting projections for unstructured grids. Along the path of extending the (W)ENO schemes to multi-dimensions and unstructured grids, tremendous effort has been made to improve the computational efficiency and numerical robustness of the schemes. For example, Abgrall [1] and Sonar [48] devised ENO schemes on unstructured grids. Hu and Shu [22] presented higher order WENO schemes on triangular meshes using a combination of linear or quadratic polynomials. To get rid of the negative linear weights, a variant [45] was proposed to enhance numerical stability, at expense of the adaptability to geometrical complexity to some extent. Alternatively, Ollivier–Gooch [39] proposed a data-dependent reconstruction procedure (DD-ENO) that uses a fixed stencil and a data-dependent weighting to satisfy the ENO condition based on the k−exact least-squares reconstruction. Afterwards, a high-order limiting strategy [41,37] was proposed to improve the monotonicity and convergence properties. More recently, Ivan and Groth [24] developed a central ENO (CENO) scheme that combines high-order k−exact least-square reconstruction based on a fixed central stencil with a monotonicity-preserving limited reconstruction algorithm. Christlieb [6] constructed a WENO scheme to enhance the positivity-preserving property by applying the parameterized limiting technique. Liu [35] designed a WLS-ENO scheme based on Taylor series expansions and weighted least square formulation. In spite of the efforts devoted, some technical barriers still remain for constructing efficient and robust non-oscillatory high-order finite volume schemes on the unstructured grids. High-order reconstructions with limiting projections using the WENO concept usually require a wide stencil which gives rise to substantially increased algorithmic complexity and difficulty for unstructured grids, where choosing admissible stencil cells and applying boundary conditions are not trivial tasks in 3D cases. Being an alternative to the conventional finite volume method, high-order methods using additional degrees of freedom (DOFs) defined within a compact stencil have received particular attention in recent years because of not only their superior convergence property but also the flexibility and easier implementation on unstructured grids. Among others, the discontinuous Galerkin (DG) method [7–10] and the spectral volume/difference method [51–54], as well as their advanced variants, flux reconstruction (FR) [23] and lifting collocation penalty (LCP) [55] methods, can be viewed as the representative schemes of this sort. In these schemes, high-order reconstructions are realized by using locally defined DOFs with structured connections, which greatly improves the efficiency of data management in local and compact stencils for reconstruction, thus well-suited for building high-order schemes on unstructured grids. Despite remarkable success in solving convection-dominant flows with smooth solutions, the compact-stencil high-order methods do have distinct disadvantage against conventional finite volume methods (FVMs) in particularly the following aspects: (i) low productivity per DOF which results in huge computational cost and high storage demand; (ii) more restrictive stability which is generally associated with exacerbated Courant–Friedrichs–Lewy (CFL) condition for explicit time integration; and more profoundly (iii) lack of robustness and reliable limiting projection to suppress spurious oscillations associated with discontinuous solutions. To alleviate the intrinsic weakness of the compact-stencil high-order methods aforementioned, some hybrid schemes have been formulated to accommodate both finite volume method and DG scheme in a unified framework, like PN PM [14,15,30] method and DG/FV [66,67] method. Seeking an approach that can provide accurate and robust solutions without loss of algorithmic simplicity and computational efficiency, we have developed a new class of numerical formulations on unstructured grids that make use of point values (PVs) at cell vertices as computational variables in addition to the volume integrated average (VIA) in the conventional FVM. The so-called VPM (Volume integrated average and Point value based Multi-moment) method [60,59] constructs higher-order reconstruction polynomial with DOFs of VIAs and PVs collected from the target cell and its immediate neighbors. Being conservative quantity, the VIA is updated by finite volume scheme for the governing equations in integral form, while PV is computed by point-wise numerical formulation for differential-form equations. More recently, the multi-moment constrained finite volume method with solution points at cell center and vertices (MCV-SPCV) [61] was

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

561

Fig. 1. The computational grid elements in 2 and 3 dimensions. From left to right is triangular, quadrilateral, tetrahedral, hexahedral, prismatic and pyramidal element.

proposed to further improve the algorithmic simplicity which replaces the VIA with the PV at the mass center of each cell and predict it by time evolution equations derived from the constraint conditions. The VPM/MCV schemes have been used to develop numerical models for incompressible flows [61,63,64] which demonstrate excellent solution quality with greatly improved numerical accuracy, robustness and computational efficiency in comparison with existing finite volume schemes and high-order schemes using local DOFs. Successive efforts have been devoted to implement VPM/MCV schemes to solve compressible flows where particular attention has to be paid for numerical oscillation and dissipation associated with discontinuous solutions. In the previous works [62,17,18], a limiting projection is applied to the trouble cells detected by a smooth indicator, where monotonicity is enforced by replacing high-order reconstruction to a limited linear polynomial and the boundary variation diminishing (BVD) formulations are devised to suppress the excessive numerical dissipation generated in the limiting projection process. As a result, the VPM/MCV schemes implemented with the limiting projection and BVD formulation have shown promising capability in accurately resolving both smooth and discontinuous solutions for various flow conditions. In this paper, we proposed another improved high-order MMFVM formulation, called VPM-Ex (VPM with extendedstencil reconstruction) scheme to solve compressible viscous flows of practical interest on unstructured grids. Different from the previous VPM method that reconstructs solutions of VIA and PV uniquely by multi-moment formulation on local coordinate system, the VPM-Ex scheme firstly determines the reconstruction polynomial of VIA by least-square method directly on the global coordinate, then retrieves its derivative values to construct the multi-moment formulations for updating PV moment. The VIA is predicted by finite volume method with numerical fluxes computed from the global reconstruction, while the PV, which is not required to be numerically conservative, is updated by a point-wise numerical formulation with high efficiency. To improve the solution quality for discontinuous solutions, we have also proposed a smoothness adaptive fitting (SAF) scheme based on two reconstructions obtained separately from the THINC/QQ (tangent of hyperbola interface capturing (THINC) with quadratic surface representation and Gaussian quadrature) method [65] and the limiting projection that adopts either the multi-dimensional limiting process (MLP) [42] or the WENO reconstruction. Compared with WENO scheme based on conventional FVM, the proposed VPM-WENO scheme prepares candidate reconstruction polynomials within more compact stencils. The resultant VPM-Ex formulation with SAF-based limiting projection scheme can resolve smooth solution with 3rd-order accuracy on arbitrary unstructured grids and obtain non-oscillatory discontinuous solution with substantially reduced numerical dissipation error. The rest of this paper is organized as follows. In Section 2, we present two different reconstruction formulations of VPMEx scheme and numerical procedure of THINC/QQ scheme. Then numerical model that uses VPM-Ex to solve Navier–Stokes equations for compressible viscous flows are described in details in Section 3 followed by the newly developed techniques for limiting projection that substantially reduces numerical dissipation. In Section 4, some benchmark tests are carried out to evaluate the performance of the proposed schemes and numerical models. We end this paper with some concluding remarks in Section 5. 2. Preliminary definition and numerical reconstructions on unstructured grids 2.1. The computational grids and variables The computational domain is partitioned into non-overlapping control volume i (i = 1, ..., I ) by single-type or mixed grid elements, including triangular and quadrilateral elements for 2D, and tetrahedral, hexahedral, prismatic and pyramidal elements for 3D as shown in Fig. 1. The vertices of cell i are denoted by θik (k = 1, . . . , K ). The encompassed edges/surfaces are similarly denoted by i j ( j = 1, . . . , J ) which are straight lines in 2D and planar surfaces of triangle or quadrangle in 3D, where K and J stands for the total numbers of vertices and surface segments of the cell. Since we only consider linear element in this work, the outward normal unit vector ni j = (nxi j , n yi j , n zi j ) is assumed to be constant on each surface i j and directed from left (+ ) to right (− ) side. We define the coordinate of vertices θik , middle points θi j of i j and mass   center θic by (xik , y ik , zik ), (xi j , y i j , zi j ) and (xic , y ic , zic ) respectively. We denote the magnitude of boundary area by i j  and the volume of cell i by |i |. We use the following convention to connect the target cell i and its neighbors: i j denotes the neighboring cell adjacent to cell i sharing boundary segment i j ; ikl (l = 1, . . . , L ) denotes the lth cell of the cell union that shares vertex θik . According to the concept of multi-moment finite volume method, we define both VIA and PVs at the vertices of cell i as computational variables

562

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 2. The reconstruction stencil for triangular (left) and quadrilateral (right) elements respectively. The cells colored with gray, green and blue mean the target cell i , neighboring cells ikl around the joint vertex θik and neighboring cell i j sharing boundary segment i j . (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

φ i (t ) ≡

1

| i |

 φ(x, y , z, t )d, (1)

i

φik (t ) ≡ φ(xik , y ik , zik , t ), k = 1, . . . , K where φ stands for any physical variable. 2.2. The reconstruction stencil The stencil selection is not trivial for conventional finite volume method on unstructured grids, particularly for highorder reconstructions, where interpolation functions are constructed by recursively adding the nearest neighbors until the minimum number of DOFs is reached for the pre-determined reconstruction function. Additional attention is also necessary to keep the stencil as symmetric as possible so as to cancel out the truncation error. In the present work, we define the stencil for the spatial reconstruction on i as the union of cells that share the vertices of i , K

L

{is | : s = 0, 1, · · · , S } ≡ ∪ ∪ ikl . k =1 l =1

(2)

Here we re-label the cells with a new set of index “is” among which the target cell i is indexed by i0 . As illustrated in Fig. 2 for two-dimensional case, for any given unstructured grid stencil (2) is uniquely determined. It is compact and involves only the immediate neighbors of the target cell. Stencil (2) circumvents the difficulties for the wide-stencil reconstructions in conventional finite volume method [22,41,33], associated with the admissible cell selection, boundary condition formulation and parallel processing. Meanwhile, it also provides the data to formulate the limiting projection needed to suppress spurious oscillations in presence of discontinuity. Different from our past works [60,59,63] that only include adjacent neighbors that share the surface i j , the present stencil contains more DOFs so as to design more sophisticated non-oscillatory high-order schemes for all-type of unstructured grids. The spatial reconstruction based on this extended stencil effectively enhances numerical accuracy and robustness without loss of computational efficiency, and is thus highly appealing in the practical simulation. To establish multi-moment reconstruction on this extended stencil, we need to reformulate the numerical procedure for reconstruction as described in the following subsections. 2.3. The third-order least-square reconstruction with multiple DOFs The third-order reconstruction can be formulated by piecewise polynomial using k−exact least-square method [3]. We use both PVs of target cell φik (k = 1, . . . , K ) and VIA of all cells φ ikl (k = 1, . . . , K ; l = 1, . . . , L ) as independent DOFs in the supporting stencil. We denote φ is (s = 0, . . . , S ) with φ i0 = φ i in accordance with the index convention of is described above. On regular grid (e.g. Fig. 2), each stencil includes 16 and 13 DOFs in total for triangular and quadrilateral elements which adequately suffice to determine the third-order reconstruction polynomial even for cells at boundaries of computational domain. Here, we present two algorithms based on different numerical procedures as below. It is noted that the number of DOFs, including both VIA and PV moments, in the compact stencil of the target cell and its immediate neighboring cells enables to construct polynomials higher than 3rd order. As we are focusing on a uniform 3rd-order finite volume method in this paper, we use the quadratic polynomial instead as follows. More profoundly, the multiple moments available in the compact stencil support the WENO-type limiting projection shown later.

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

563

2.3.1. Algorithm I: direct reconstruction The third-order accuracy can be gained straightforwardly through a piecewise quadratic reconstruction written in global coordinate system as (1 )

i (x, y , z) = c 000 + c 100 (x − xic ) + c 010 ( y − y ic ) + c 001 ( z − zic ) + c 200 (x − xic )2 + c 020 ( y − y ic )2

(3)

2

+ c 002 ( z − zic ) + c 110 (x − xic )( y − y ic ) + c 011 ( y − y ic )( z − zic ) + c 101 (x − xic )( z − zic ). The barycentric coordinates (xic , y ic , zic ) are used in order to circumvent the scaling efforts [19]. The coefficients crst (0  r , s, t  2) of (3) are derived by enforcing the following constraint conditions, (1 )



1

|is |

i (xik , y ik , zik ) = φik , (k = 1, . . . , K )

(4)

(i 1) (x, y , z)d = φ is , (s = 0, . . . , S ).

is

The volume integration of (3) is performed by numerical quadrature which is elaborated in Appendix A. For simplicity, we further re-write constraint conditions (4) into



1

|is |

|i0 |

is

|i0 |

(1 )

i (x, y , z)d = φik − φ i0 , (k = 1, . . . , K )

(5)

i0



1

(1 )

i (x, y , z)d −



1

(1 )

i (xik , y ik , zik ) −

(1 )

i (x, y , z)d = φ is − φ i0 , (s = 1, . . . , S ), i0

which is cast into a matrix form of overdetermined linear system as

M · A = B,

(6)

where M denotes the coefficient matrix depending only on the mesh geometry, A the vector containing the unknown coefficients crst (r + s + t > 0) for the reconstruction, and B contains solution values as the right-hand-side of (5). Then the unknowns are uniquely determined from



 −1

A = MT M

M T B or

A = (M)−1 B

(7)

directly following least-square method. In practice, a geometric weighting technique is introduced to improve the locality of reconstruction which is ignored here for the sake of simplicity. We prefer singular value orthogonal decomposition (SVD) method to compute the inverse matrix which can be stored for reuse in reconstruction of all variables. It is noted that c 000 is eliminated from the unknown coefficients in (4) thus leads to a one-rank-lower matrix in (6). Once A is solved, we can explicitly compute c 000 by

c 000 = φ i0 −

 

1

|i0 |

c 100 (x − xic ) + c 010 ( y − y ic ) + c 001 ( z − zic ) + c 200 (x − xic )2 + c 020 ( y − y ic )2

(8)

i0

 + c 002 ( z − zic )2 + c 110 (x − xic )( y − y ic ) + c 011 ( y − y ic )( z − zic ) + c 101 (x − xic )( z − zic ) d. The reconstruction function (3) can be also equivalently expressed into a basis function form (1 )

i (x, y , z) =

K 

ψik φik +

k =1

S 

(9)

ψ is φ is .

s =0

The ψik and ψ is are the basis functions of their corresponding DOFs with explicit dependence limited to the coordinate



−1

system, i.e. {ψ} = {x, y , z, x2 , y 2 , z2 , xy , yz, xz} which can be extracted from M T M M T in (7) and (8) straightforwardly. In (9), both PVs and VIAs are used through the least square formulation. It is obvious that we can recover the conventional k-exact finite volume method by removing the PVs related terms in (9). Our numerical experiments show that both numerical errors and elapsed times of the k-exact finite volume method are comparable to those of VPM-Ex scheme with algorithm I, which includes the 3rd-order k-exact finite volume method as a particular case. As aforementioned, we have sufficient DOFs on the same stencil to construct a higher-order, e.g. cubic polynomial reconstruction in the following form, (1 )

i (φ|; x) =

(α +β+γ )3  

α =0 β=0 γ =0

(x − xic )α ( y − y ic )β ( z − zic )γ c α β γ ,

(10)

564

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

which contains additional cubic terms, i.e. x3 , y 3 , z3 , x2 y , xy 2 , x2 z, xz2 , y 2 z, yz2 , xyz in the basis function. The PVs at cell vertices provide the DOFs needed to raise the order of reconstruction polynomial in comparison with the conventional k−exact least-square method which at most allows to build a quadratic polynomial with only the VIAs on the same stencil. In this work, we make use of the quadratic reconstruction (3) rather than (10) so as to be consistent with the following algorithm II, and integrate with the limiting projection schemes described afterwards in section 3.3. 2.3.2. Algorithm II: hierarchical reconstruction The quadratic polynomial can be alternatively constructed by computing first and second order derivatives at cell center using the available DOFs. The reconstruction procedure is arranged to hierarchically implement the following two steps:

• Step 1: compute first order derivatives ∇φik = (φxik , φ yik , φzik ) at the vertices θik . The linear polynomial can be constructed within the sub-stencil consisting of the surrounding cells ikl sharing vertex θik (the green and gray colored cells in Fig. 2)

ik (x, y , z) = φik + φxik (x − xik ) + φ yik ( y − y ik ) + φzik ( z − zik ).

(11)

Following the constraint conditions in terms of the volume integrated average of (11) for the surrounding cells,

1

|ikl |



ik (x, y , z)d = φ ikl , (l = 1, . . . , L ),

(12)

ikl

the unknown (φxik , φ yik , φzik ) can be uniquely determined by least-square method. • Step 2: approximate first and second order derivatives at cell center θic of target cell i . We can construct linear polynomial in the target cell i

ic (x, y , z) = ψic + ψxic (x − xic ) + ψ yic ( y − y ic ) + ψzic ( z − zic ), where ψ = φγ (γ ∈ (x, y , z)) stands for the partial derivative respect to coordinate ditions,

(13)

γ . According to the constraint con-

ic (xik , y ik , zik ) = ψik , (k = 1, . . . , K ),

(14)

as well as (φxik , φ yik , φzik ) at cell vertices obtained in step 1, the unknowns of (14) can be determined by least-square method in the same manner. Then we can get the first-order derivatives (φxic , φ yic , φzic ) and second-order derivatives (φx2 ic , φ y2 ic , φz2 ic , φxyic , φ yzic , φxzic ) at cell center. Here, the second-order derivatives are computed by differentiating the first-order derivatives as

φα β ic =

1



2

 ∂φβ  ∂φα , (α , β = x, y , z). + ∂β ∂ α θic

Once the first and second order derivatives are obtained, we can construct the quadratic polynomial from the Taylor expansion as (1 )

1

1

2

2

i (x, y , z) = φic + φxic (x − xic ) + φ yic ( y − y ic ) + φzic ( z − zic ) + φx2 ic (x − xic )2 + φ y 2 ic ( y − y ic )2

(15)

1

+ φz2 ic ( z − zic )2 + φxyic (x − xic )( y − y ic ) + φ yzic ( y − y ic )( z − zic ) + φxzic (x − xic )( z − zic ) 2

with φic determined to coincide with c 000 in (8). In spirit of [35] where the accuracy of arbitrarily qth-order least-square reconstruction is analyzed, we can prove that both algorithm I and algorithm II described above for spatial reconstruction have 3rd order accuracy. 2.3.3. Implementation details and comparisons between algorithm I and algorithm II In general, the CFD General Notation System (CGNS) convention [49] is employed to represent unstructured mesh using connectivity tables to link topological entities, such as vertices, edges/surfaces and elements. Without pre-existing connectivity information among the cells which are not directly adjacent, implementation of the CGNS convention to generate stencil is (as Fig. 2) for target cell i involves global search via vertex-element connectivity (ikl ) and renumbering to remove the repeated index and communication in parallel implementation. In algorithm I, we use a pre-processing step before the computation to build the connectivity between the target cell i and its neighboring elements as well as their vertices within stencil is (as Fig. 2). The reconstruction in algorithm II is more straightforward and only relies on cellvertex connectivity. The large condition number relevant to high order polynomial is not experienced in the present work which can be further controlled by a simple column-scaling procedure [24]. The curved boundary is not considered here but high order boundary representation can be naturally treated within the present framework according to existing studies [40,24,56].

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

565

Fig. 3. The computational cost and parallel scalability for VPM-Ex scheme with algorithm I and algorithm II for reconstruction measured on CPU cores of Intel Xeon(R) E5-2630, 2.20 GHZ. The solid lines show the elapse time per iteration and MPI speedup obtained in a test problem with fixed DOFs. The dashed lines indicate the ideal parallel speedup.

We compare the major properties of two algorithms as follows. Algorithm I is straightforwardly derived from quadratic polynomial which requires Gaussian quadrature for numerical integration during the reconstruction. In contrast, algorithm II involves only linear reconstruction and enables to compute the volume integration (12) simply by the PV at mass center which is thus more convenient and computationally efficient. The reconstruction of algorithm I using higher order polynomial deals with more unknowns and results in a larger matrix (6), consisting 9 × 9 coefficient entries, while only 4 × 4 matrix is involved in algorithm II. Inversion of the matrix in algorithm I requires the SVD approach which demands larger memory storage and more operation counts, while the matrix in algorithm II can be inverted analytically. Considering the dynamic mesh problem, algorithm II is preferable since the matrix inversion has to be recomputed at each time step. In spite of larger computational cost, the basis function of each DOF in algorithm I gathered in (7) is computed as a preprocessing step and memorized before the time iterations, which facilitates the incorporation of implicit formulation for spatial discretization. Since the present work is restricted to the explicit time integration, we prefer algorithm II in the computations for simplicity except for some verification tests. The computational cost and parallel scalability are the key aspects for a practical numerical method, which are of particular importance for massively parallel processing. We measured the computational cost and MPI speedup of the reconstruction algorithms for a problem of a fixed number of DOFs on a parallel computer with 32 cores. In order to obtain reliable measurement of execution time, we repeated the computation of each algorithm for 106 iterations, and estimated the elapse time per iteration by simple averaging. To assess the parallel performance, we ran the same problem with gradually increased CPU cores from 2 to 32. We plotted the elapse time and parallel speedup in Fig. 3 for both reconstruction algorithms. It is observed that the elapse time of algorithm I on single core is slightly larger than that of algorithm II. On the other hand, both schemes realize satisfactory parallel scalability for this strong scaling test, which indicates a less ratio between data communication overhead and computing intensity. Similar to other high-order methods using additional local degrees of freedom (DOFs), the reconstruction algorithms proposed in this paper are more computing intensive and thus superior in parallel processing. It is noted that in this test a simple MPI implementation is used without any special technique to hide the communication overhead for MPI data transfer. 2.4. The global reconstruction of THINC/QQ method In addition to the polynomial reconstruction in section 2.3, we can also piecewisely construct the hyperbolic tangent function featured with step-like distribution which is preferable for the discontinuous solution. Considering a hyperbolic (2) tangent function i (x, y , z) on Cartesian coordinate in the target cell i ,

(i 2) (x, y , z) = φmin +

φmax − φmin 2

(1 + tanh (β (Pi (x, y , z) + di ))) ,

(16)

where φmin and φmax are minimum and maximum values of VIAs in the stencil is (as depicted in Fig. 2). β is the steepness parameter to control the thickness of the transition layer. The Pi (x, y , z) + di = 0 is the equation of jump interface and approximated by a quadratic polynomial,

566

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Table 1 Switching parameters in (22) for different elements. Element type

Switching parameters

Triangular element Quadrilateral element Tetrahedral element Hexahedral element Prismatic element Pyramidal element

A1 A1 A1 A1 A1 A1

= 1, = 1, = 1, = 1, = 1, = 1,

B1 B1 B1 B1 B1 B1

= 0, = 0, = 1, = 1, = 1, = 1,

A2 A2 A2 A2 A2 A2

= 0, = 1, = 0, = 1, = 0, = 1,

B2 B2 B2 B2 B2 B2

=0 =0 =0 =1 =1 =0

Pi (x, y , z) = c 200 (x − xic )2 + c 020 ( y − y ic )2 + c 002 ( z − zic )2 + c 110 (x − xic )( y − y ic ) + c 011 ( y − y ic )( z − zic ) (17)

+ c 101 (x − xic )( z − zic ) + c 100 (x − xic ) + c 010 ( y − y ic ) + c 001 ( z − zic ), (1)

with coefficients c str (s, r , t = 0, 1, 2 and s + r + t  2) retrieved from i (x, y , z) in (3) or (15). The only unknown di in (16) is solely determined from the VIA φ i in the target cell by,



1

| i |

(i 2) (x, y , z)dxdydz = φ i .

(18)

i

Following [65], we solve (18) by using numerical Gaussian quadrature to approximate the integration of hyperbolic tangent function (16). Denoting the coordinates and weight of Gaussian points by xig = (xig , y ig , zig ) and ωig ( g = 1, 2, . . . , G ) satisfying G 



G g =1



ωig

g =1



ωig = 1 in element i , we approximate (18) by Gaussian quadrature as follows,

   φ i − φmin = 1 + tanh β Pi (xig ) + di . 2 φmax − φmin 1

 

Re-writing tanh β Pi (xig ) + di







tanh β Pi (xig ) + β di =

(19)

in (19) as

tanh(β Pi (xig )) + tanh(β di ) 1 + tanh(β Pi (xig )) · tanh(β di )

,

then we recast (19) into G  g =1

φ i − φmin 1 ωig − =2 , 1 + tanh(β Pi (xig )) · tanh(β di ) φmax − φmin 2 tanh(β Pi (xig )) + tanh(β di )



(20)

and solve di as the only unknown from (20). For brevity, we denote A g = tanh(β Pi (xig )), D = tanh(β di ) and Q = 2(φ i − φmin ) /(φmax − φmin ) − 1, and rewrite (20) as G 

ωig

g =1

Ag + D 1 + Ag D

= Q,

(21)

yielding a high order rational equation regarding the unknown D which can be solved by Newton–Raphson method within a few iterations. See [65] for more details. 2.5. The local reconstruction of VPM method





Provided φξ ic , φηic , φζ ic and



 φξ 2 ic , φη2 ic , φζ 2 ic

at mass center θic , we can build the cell-wise reconstruction of VPM

method. The form of the reconstruction function varies according to the type of mesh cell. We cast it in an unified formulation (22) with switching parameters A 1 , A 2 , B 1 , B 2 summarized in Table 1 for different cell types of unstructured grids.

(i 3) (ξ, η, ζ ) =

K 

ψik φik + ψ φ i + A 1 (ψξ φξ ic + ψη φηic ) + B 1 ψζ φζ ic + A 2 (ψξ 2 φξ 2 ic + ψη2 φη2 ic ) + B 2 ψζ 2 φζ 2 ic , (22)

k =1

The basis functions {ψ} are merely related to the  element shape and  uniquely determined by the constrained conditions   (1) (23). The derivative terms of φξ ic , φηic , φζ ic and φξ 2 ic , φη2 ic , φζ 2 ic are retrieved from i (x, y , z) on global coordinate and then transformed into local coordinate regarding to each grid element. For more details, see [60,59,63,61].

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

(3 )



1

i (ξik , ηik , ζik ) = φik

567

(23)

(3 )

i (ξ, η, ζ )d = φ i

| i |

 (3 )

∂ i (ξ, η, ζ )   = φξ ic (ξic ,ηic ,ζic ) ∂ξ ∂ (i 3) (ξ, η, ζ )   = φηic (ξic ,ηic ,ζic ) ∂η (3 )

∂ i (ξ, η, ζ )   = φζ ic (ξic ,ηic ,ζic ) ∂ζ ∂ 2 (i 3) (ξ, η, ζ )   = φξ 2 ic (ξic ,ηic ,ζic ) ∂ξ 2 (3 )

∂ 2 i (ξ, η, ζ )   = φη2 ic (ξic ,ηic ,ζic ) ∂ η2 (3 ) ∂ 2 i (ξ, η, ζ )   = φζ 2 ic (ξic ,ηic ,ζic ) ∂ζ 2

Remark.

• Within the supporting stencil is (Fig. 2), we can independently construct quadratic reconstruction (i 1) , (i 2) and (i 3) for the target cell i . (1)∼(3) • i are well linked by sharing the first-order derivatives (φxic , φ yic , φzic ) and second-order derivatives (φx2 ic , φ y2 ic , φz2 ic , φxyic , φ yzic , φxzic ) at mass center θic . (1)

• i

(3)

and i , denoted by VPM-Ex reconstruction hereafter, are tightly coupled to compute numerical fluxes for updating VIA and PV.

3. Numerical formulation for Navier–Stokes equations 3.1. The governing equations and spatial discretization We describe unsteady viscous compressible flow by the following Navier–Stokes equations,

∂U + ∇ · F (U) − ∇ · G (U) = 0, ∂t

(24)

in which the conservative variables U, the inviscid fluxes F = (F x , F y , F z ) and the viscous fluxes G = (G x , G y , G z ) are defined as

U = [ρ ,



ρ u , ρ v , ρ w , E]T ,

(25)

T F x (U) = ρ u , ρ u 2 + p , ρ uv , ρ u w , u ( E + p ) , (26)  T F y (U) = ρ v , ρ uv , ρ v 2 + p , ρ v w , v ( E + p ) , (27)  T F z (U) = ρ w , ρ u w , ρ v w , ρ w 2 + p , w ( E + p ) , (28)   ∂T T G x (U) = 0, τxx , τ yx , τzx , u τxx + v τxy + w τxz + κ , (29) ∂x   ∂T T G y (U) = 0, τxy , τ y y , τ yz , u τ yx + v τ y y + w τ yz + κ , (30) ∂y   ∂T T G z (U) = 0, τxz , τ yz , τzz , u τzx + v τzy + w τzz + κ , (31) ∂z where m = (mx , m y , m z ), u = (u x , u y , u z ) = (u , v , w ), ρ , E, T and p denote the momentum, velocity, density, total energy, temperature and pressure respectively. It is noted that G (U) are written in primitive variable form which are computed from conservative variables by related formulations. Besides, the thermal conductivity coefficient κ is computed by

568

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

κ=

μγ cν

with c ν =

Pr

1

γ −1

R

(32)

and the viscous stress tensor is given by



τα β = μ

∂ uβ ∂ uα + ∂β ∂α

−



2

(∇ · u) δα β , (α , β = x, y , z),

3

(33)

where μ is the dynamic viscosity of the fluid, R the gas constant, P r the Prandtl number and δα β the Kronecker delta. For   perfect gas, the pressure is determined from the equation of state as p = E − ρ u2 2 (γ − 1) with γ = 1.4 being the ratio of the specific heats. We define both VIA and PV moments for all computational variables and update simultaneously in time with the third order TVD Runge–Kutta scheme [47]. The approximate Roe solver [44] with entropy fix [31] is used to evaluate inviscid flux in both integral and differential form. The spatial discretization formulations regarding to the solutions of VIA and PV moments are separately described in semi-discrete form as follows. 3.1.1. Solution of Navier–Stokes equations for VIA moment The VIA moment is updated by the integral form of (24) which can be recast as follows according to the divergence theorem,

∂ ∂t







Ud + i

Fi j (U) d =

i j

Gi j (U) d.

(34)

i j

Approximating the surface integration with Gaussian quadrature, we arrive at G     ∂ Ui 1  =− ωig Fi j (U) − Gi j (U) · ni j i j  . | i | ∂t J

(35)

j =1 g =1

The numerical fluxes can be expressed as Fi j (U) = AUi j where A = R Λ L is the Jacobian matrix with R / L and Λ being the matrices of right/left eigenvectors and corresponding eigenvalues respectively. According to Riemann solver of Roe algorithm [44], we compute inviscid flux in (35) by

    F i j ( U) = A U+ + U− − U− + R |Λ| L U+ . ij ij ij ij

(36)

It is noteworthy that all the values in the eigensystem are computed by Roe-averaging approximation with two adjacent VIAs. To evaluate the viscous flux at boundary surface i j , we first compute values (Ui j ) and derivatives (∇ U)i j of conservative variables by arithmetic mean of left and right reconstructed solutions and then compute the required gradients of velocity and temperature. For brevity, the viscous flux in (35) is expressed by

  Gi j (U) = G i j Ui j , (∇ U)i j = G i j

 1 2

 1 

− . + (∇ U ) (∇ U)+ ij ij

U+ + U− , ij ij

2

(37)

Direct substitution of (36) and (37) into (35) gives the final form of finite volume scheme J G      ∂ Ui 1  =− ωig A U+ + U− − U− + R |Λ| L U+ ij ij ij ij | i | ∂t j =1 g =1   1 

  1 − + − −G i j U+ + U U ) + (∇ U ) n i j  i j  . , (∇ ij ij ij ij

2

(38)

2

3.1.2. Solution of Navier–Stokes equations for PV moment Given Jacobian matrices A = F x (U)ik / ∂ Uik , B = F y (U)ik / ∂ Uik and C = F z (U)ik / ∂ Uik , we rewrite (24) into

∂ Uik ∂ Uik ∂ Uik ∂ Uik +A +B +C = (∇ · G (U))ik , ∂t ∂x ∂y ∂z

(39)

where A = R A Λ A L A , B = R B Λ B L B and C = R C ΛC L C are Jacobian matrices. The corresponding left/right eigenvectors and the diagonal matrix of eigenvalues are denoted by L / R and Λ respectively. We evaluate the matrices by the Roe-averaged values φ˜ ik from the surrounding VIAs φ ikl (l = 1, 2, . . . , L ) encompassing vertex θik

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

569

L  

ρ ikl · φ ikl

φ˜ ik =

l =1

(40)

,

L  

ρ ikl

l =1





˜ 2ik )) . ˜ ik and H˜ ik , and the sound speed is computed by a˜ ik = (γ − 1)( H˜ik − 12 (u˜ 2ik + v˜ 2ik + w where φ˜ ik , standing for u˜ ik , v˜ ik , w We solve (39) with inviscid flux computed point-wisely by Roe’s Riemann solver and viscous flux simply approximated by the TEC (Time-evolution Converting) formulation [58]

(∇ · G (U))ik = T EC ((∇ · G (U))ikl ) , l = 1, . . . , L ,

(41)

which yields the semi-discrete scheme for updating PV moment





  ∂ U− ∂ U− ∂ U+ ∂ U+  ˜ ˜ ik ik ik ik ˜A ˜ + − + R A Λ A  L A ∂x ∂x ∂x ∂x



  ∂ U− ∂ U− ∂ U+ ∂ U+ 1  ˜ ˜ ik ik ik ik ˜ ˜ − B + − + R B Λ B  L B 2 ∂y ∂y ∂y ∂y



  ∂ U− ∂ U− ∂ U+ ∂ U+ 1  ˜ ˜ ik ik ik ik ˜ ˜ C + − + R C ΛC  L C − 2 ∂z ∂z ∂z ∂z

∂ Uik 1 =− ∂t 2

where

(42)

+T EC ((∇ · G (U))ikl ) , +

∂ Uik ∂ U− ik ∂α , ∂α



denote the derivatives of variables Uik regarding

α ∈ (x, y , z) direction reconstructed on left-side and

right-side stencil of point θik . The algorithmic details can be found in [44,18]. 3.2. VPM-Ex multi-moment reconstruction formulation For Euler equations, the reconstructions can be applied to any set of dependent quantities of either conservative, primitive or characteristic variables. Nevertheless, the primitive variables are preferable regarding to preserving monotonicity, positivity and efficiency, as commented in [37,24,5,27]. In this work, we employ the primitive variables W = [ρ , u , v , w , T ]T derived from conservatives variables by W = Q−1 U with Q being the transformation matrix [4,46]. We de(1)∼(3) (1)∼(3) (W|; x|ξ ) to represent the reconstructions i in terms of the primitive variables on global coordinate(x) note i (1)

(1)

or local coordinate(ξ ) as described in sections 2.3-2.5. Provided i (W|; x) and i j (W|; x) in two adjacent cells which may be discontinuous at their boundary surface Γi j , the integrated averages as well as derivatives of conservative variables in (36) are computed from left-side



Q+ ij   U+ = ij

Γi j 

(1 )

i (W|; x) dΓi j ,

(43)

Γi j

Q+

ij (∇ U)+ =  ij Γ i j 



∇ (i 1) (W|; x) dΓi j ,

(44)

Γi j

or from right-side

Q− ij   U− = ij

Γi j 



(1 )

i j (W|; x) dΓi j .

(45)

Γi j

Q−

ij (∇ U)− =  ij Γ i j 



∇ (i j1) (W|; x) dΓi j .

(46)

Γi j

The individual entries in transformation matrix Q± is computed from reconstructed primitive variables in left or right side ij stencil. The VIA Ui then is updated by (38) to proceed into the next time step.

570

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Table 2 The cutoff number Sc of indicator function for different type of grid elements. Elements

Triangular

Quadrilateral

Tetrahedral

Hexahedral

Prismatic

Pyramidal

Sc

5 × 104

5 × 105

2 × 105

5 × 105

1 × 106

1 × 105

(3)

For PV at vertices θik , we firstly compute the derivatives of primitive variables from ikl (W|; ξ ) in the local coordinate ξ , which are constructed cells sharing vertices θik . Then we compute derivatives of conservative on global  on all surrounding  ∂U

coordinate ∂ xikl =

with transformation matrices Qik and Mξ →x , i.e.

∂ Uikl ∂ Uikl ∂ Uikl ∂x , ∂ y , ∂z

∂ Uikl = Qik · Mξ →x ∂x



(3 ) ∂ ikl (W|; ξ ) , ∂ξ

(47)

where Mξ →x denotes the Jacobian matrix mapping from local coordinate system ξ = (ξ, η, ζ ) to physical coordinate x = (x, y , z) as given in [61]. It is noted that Qik is simply computed by PVs at vertices θik which are continuous for all adjacent cells, whereas the values on the boundary surface Γi j are not continuous. In order to get the left- and right-side derivative values required in (42), we devised a biased weighting formula L  ∂ U± τ ± ∂ Uikl ik = ωikl , τ = x, y or z ∂τ ∂τ

(48)

l =1

where

τ± ωikl is a weight that adaptively computes the appropriate derivative at θik in left- or right-direction τ ,  −−−→ max 0, nτ ± · θiklc θik τ±  ωikl = −−→ . L τ± · − max 0 , n θ θik iklc l =1

(49)

−−−→

The θiklc θik denotes the vector from the mass center of cell ikl to vertex θik and nτ ± the unit normal vectors, nx± = (±1, 0, 0), n y± = (0, ±1, 0) and nz± = (0, 0, ±1) indicating the “+” or “−” directions of x, y and z components respectively. Once derivatives (48) are determined, PVs can be updated by (42) straightforwardly. 3.3. The smoothness evaluation and limiting projection As suggested in [62,18], the smoothness of high-order reconstruction can be effectively evaluated by examining the total boundary variation T B V ()i of the target cell

J T B V ()i = 1 −

j =1





1 i j  i j

4  (1 ) (1 ) i (x, y , z) d − 1  i j i j (x, y , z) d ij , 4 J  φi − φi j j =1

(50)

which indicates the total differences of the integrated averages reconstructed at the cell boundaries. The smoothness indicator then is structured in form as

S= with

T B V ()i max((1 − T B V ()i ),  )

,

(51)

 = 10−16 being a small value for avoiding zero-division. The total boundary variation T B V ()i ranges between

(−∞, 1] and it rapidly approaches unity for the sufficient smooth solution. Conversely for reconstruction containing discontinuous solution, the T B V ()i will be far from unity which yields a smaller or even negative value for the smoothness indicator S . We identify the oscillatory solution as S < Sc by defining a cutoff number Sc which is case-independent but has optimal values for different grids. We give preferred Sc for different grid elements in Table 2 which are used in the numerical tests throughout this paper. 3.3.1. Multi-dimensional limiters for VIA For the trouble cell in presence of discontinuous solution which is identified by S < Sc , we impose monotonicity by degrading the high-order reconstruction to a limited linear interpolation function (1)

  (x, y , z) = φ i +  φxic (x − xic ) + φ yic ( y − y ic ) + φzic ( z − zic )

(52)   (1) where φxic , φ yic , φzic is the first-order derivative at mass center θic and directly obtained from i (x, y , z) in (3) or (15). The limiter function  ∈ [0, 1] is determined by the multi-dimensional limiting process (MLP) [42] scheme with MLP-u2

i

limiter. The VPM-Ex reconstruction combined with the MLP limiter is referred to as VPM-MLP scheme.

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

571

Fig. 4. The configuration of biased sub-stencils for 3rd-order VPM-WENO scheme using multiple reconstructions.

Besides MLP, we also combine the WENO method to improve the accuracy of numerical solutions in this study. In addition to the primary reconstruction (section 2.3), several biased stencils are used to construct candidate polynomials which are then used to build a polynomial of higher order through a convex combination in the spirit of WENO. (1) Concretely, we build the K quadratic polynomials denoted by ik (x, y , z), 0  k  K

1

1

2

2

k k k (ik1) (x, y , z) = φic + φxic (x − xic ) + φ kyic ( y − y ic ) + φzic ( z − zic ) + φxk2 ic (x − xic )2 + φ ky 2 ic ( y − y ic )2

(53)

1

k k + φzk2 ic ( z − zic )2 + φxyic (x − xic )( y − y ic ) + φ kyzic ( y − y ic )( z − zic ) + φxzic (x − xic )( z − zic )

2

with k = 0 being the primary reconstruction on the biased stencils as shown in Fig. 4. The biased stencils are selected so as to include DOFs defined over a compact support for spatial reconstruction which includes all the VIAs (φikl ) of adjoint vertex and PVs (φik ) of the target cell. The number of DOFs in the sub-stencils are sufficient for quadratic polynomial such that we can readily determine the unknowns at mass center (θic ) by the least-square approximation. Compared with the conventional finite volume WENO schemes (e.g. [12] (Fig. 2)), the stencil used in the proposed scheme is more compact, which significantly improves the algorithmic complexity and numerical efficiency, particularly in practical applications and parallel computations. We find the combination coefficients by measuring the smoothness of reconstruction candidates as in [22,45]

σik =

3 

|l|=1

 | i |ϑ

∂ |l| (1 )  (x, y , z)dxd ydz , ∂ xl1 ∂ yl2 ∂ zl3 ik

i

where l = (l1 , l2 , l3 ) is a multi-index with |l| = l1 + l2 + l3 . The superscript ϑ denotes the normalization factor for dimensionality which is ϑ = |l| − 1 for 2D and ϑ = 23 |l| − 1 for 3D respectively. Following [33], we compute the non-linear weights by

ω˜ ωik =  K ik , ˜ ik k =1 ω

ω˜ ik =

λik ( + σik )4

.

Here,  = 10−8 is a small positive number used for avoiding zero-division. The λik is a pre-specified linear weights chosen as 100 for primary reconstruction (k = 0) and 1 for the rest ones. We come to the final formulation of VPM-WENO reconstruction given by a convex combination of all candidates polynomials prepared in the stencil 1) ( (x, y , z) = i

K 

ωik (ik1) (x, y , z).

k =0

The present scheme provides excellent capability to capture shock wave and contact discontinuities essentially free from numerical oscillation, and more importantly retains uniformly high order accuracy for both smooth and non-smooth solutions, which will be demonstrated in section 4. (1) For conciseness, we denote the limited VPM-MLP|WENO reconstruction by the same notation i (x, y , z). 3.3.2. Limiting constraints for PV Since it is not necessary to ensure conservativeness of PV moment, the overshoots/undershoots of PV can be simply removed by

572

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

   φik = min φik , max φ ikl ,    φik = max φik , min φ ikl , 





(54)



where max φ ikl / min φ ikl indicate the maximum/minimum VIAs (φ ikl ) on the surrounding cells sharing vertex θik . 3.4. The smoothness adaptive fitting reconstruction for discontinuous solution The limiting projection process generally introduces excessive numerical dissipation for discontinuous solution. As illustrated in [62,18], the THINC scheme is desirable for the reconstruction of discontinuities inside cell so as to reduce the numerical dissipation caused by the boundary variations of reconstructed values. To further improve the solution fidelity, we (1) propose a novel concept to combine limited VPM-MLP|WENO reconstruction i (x, y , z) and the THINC/QQ reconstruction (3)

i (x, y , z) via a weight function, which is adaptively determined by minimizing the variations between integrated aver-

age of interpolation function and surrounding VIAs of numerical solutions. As shown later, the present scheme, so-called smoothness adaptive fitting (SAF) reconstruction significantly enhances the solution faithfulness that outperforms the results of previous BVD scheme [62,18]. We firstly describe the solution procedure of SAF for linear advection equation then the Euler equations. (1)

1. Recall the limited VPM-MLP|WENO reconstruction i

(x, y , z) and the unstructured multi-dimensional THINC recon-

(2)

struction as i (x, y , z) in the target cell i . To measure the matching degree of THINC reconstruction, we define the indicator function as

S A F ()i

R=

(55)

,

max((1 − S A F ()i ,  )

with S A F ()i being determined by

S



s =1

S A F ()i = 1 −

1



(2 )

i (x, y , z) d − φ is 4 S  s=1 φ i − φ is

|is | is

4 (56)

,

(2)

which indicates the variations between the i (x, y , z) and surrounding VIAs within the stencil is (Fig. 2). A cutoff number Rc = 100 is given to examine whether the target cell is suitable to be reconstructed by THINC function. We proceed to the next step only if R < Rc holds. (4) (2) (1) 2. Re-formulate the reconstruction function by i (x, y , z) = ωi i (x, y , z) + (1 − ωi )i (x, y , z) where ωi ∈ [0, 1] is (4)

the weighting function, which is determined by minimizing the variation of i the supporting stencil,

⎛

i =

S  s =1

⎜ 1 ⎝ |is |

(x, y , z) with its neighboring VIAs in

⎞2



⎟ (i 4) (x, y , z) d − φ is ⎠ .

(57)

is

Consequently, the weight function

ωi can be determined from

∂ i = 0, ∂ ωi which yields

S

ωi =

s =1



1



|is | is

    (2 ) (1) (1) i (x, y , z) − i (x, y , z) d φ is − |1is | is i (x, y , z) d .  2  S  1   (2 ) (1) (x, y , z) d s=1 |is | is i (x, y , z ) − i



(4)

Once ωi is obtained, we make use of i (x, y , z) to compute integrated averages and derivatives as (43)–(46) in target cell i . It is noted that the above SAF algorithm is only applied for cells containing discontinuities identified by the criterion S < Sc . We refer to the VPM-MLP|WENO reconstruction implementation in conjunction with the SAF technique as VPM-MLP|WENO-SAF reconstruction. For Euler equations, we extend the SAF algorithm to each of the primitive variables directly with a strategy that the velocity components and the pressure variable are only applied to discontinuous vicinity identified with ωi > 0 for density variable. This simple practice gives physically admissible solutions with high-fidelity reconstruction where the numerical oscillation and dissipation are effectively prevented from appearing close to the discontinuities.

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

573

3.5. The summary of solution procedure To facilitate the reader to follow the presented numerical methods, we summarize the solution procedure as follows: n

n

1. Given conservative variables of VIA Ui and PV Unik at step n, we transform them into primitive variables Wi and Wnik (1)−(3)

and construct the piecewise reconstructions i

(W, x) following sections 2.3–2.5.

(1)

2. Project the high order reconstruction i (W, x) to the limited linear interpolation (52) or non-oscillatory quadratic polynomial (53) for trouble cells in presence of discontinuous solutions which are detected by indicator function (51) under the condition of S < Sc . (2) 3. Evaluate the matching degree of THINC reconstruction i (W, x) by (55) and apply the SAF algorithm for cells only if R < Rc and S < Sc hold, as described in section 3.4. 4. Compute the piecewisely reconstructed integrated averages and derivatives on cell boundaries (43)–(46) and derivatives at cell vertices by (48). n 5. Update VIA Ui and PVs Unik by Roe’s approximate Riemann solver (38) and (42), respectively. Remark.

• The proposed reconstruction procedures described in section 2.3 reduces to a 3rd order finite volume scheme if we remove the PVs and related constrained conditions in the linear equation system (6). In this case, the resulted scheme can still be combined with the MLP limiter and SAF algorithm to formulate a conventional finite volume method. However, without the PVs the proposed WENO limiter can not be used along with reconstruction because there is not sufficient DOFs on the biased sub-stencils. 4. Numerical results In this section, we perform presented model with some widely used benchmark tests to verify its performance of numerical accuracy and computational efficiency. In order to quantify the convergence behavior, we define numerical errors in L 1 and L ∞ norms as follows,

Ne

E (L1 ) =

(|φni − φei ||i |) , Ne (|φei ||i |) i =1

i =1

E (L∞ ) =

max |φni − φei | max |φei |

,

(58)

where φni and φei denote numerical and exact solutions respectively. 4.1. Accuracy of advection equation on a sine wave test We firstly evaluate the numerical accuracy and efficiency of scheme for the linear advection equation where a sine function [42] is initially given as follows

φ0 (x, y ) = sin(2π x) · sin(2π y ). The computational domain is [0, 1] × [0, 1] with periodic boundary condition applied in x − y directions. In order to check the grid dependency, we divide the domain by both triangular and quadrilateral elements with gradually increasing resolution. We choose the uniform advection velocity (u , v ) = (1, 2) to avoid mesh alignment. The time step is adjusted with fixing CFL number of 0.2 until t = 1.0. To compare the numerical accuracy of VPM(-Ex) scheme obtained with different reconstruction manner, we present the L 1 and L ∞ errors and convergence rates for VPM schemes implemented with algorithm I and algorithm II as shown in Table 3 which are denoted by VPM-I and VPM-II respectively. It is clearly seen that both schemes have realized desired 3rd order accuracy on all grids and marginal variations in numerical errors are observed between VPM-I and VPM-II schemes, which exhibit an almost theoretically optimal performance as expected. Then we combine VPM-I and VPM-II schemes with MLP|WENO limiters and SAF reconstruction to show their influence on numerical accuracy for smooth solutions. As displayed in Table 4, all the schemes effectively preserve the 3rd order accuracy and give identical errors to the original schemes on the relatively refined grids, which indicates that the smooth indicator effectively detect solution smoothness and prevented limiting projection from degrading the numerical accuracy. On the other hand, a slightly lower level of solution errors are observed for the scheme with WENO limiter on the coarser grids, which justifies its excellent capability to recover the potentially highest accuracy for smooth and non-smooth solution. In order to quantify the computational cost, we perform the above tests on a PC with a single CPU of Intel Xeon(R) E5-2630, 2.20 GHZ and include the elapse time in Table 5. It is seen that algorithm I arrives at nearly same level of accuracy with slightly larger elapsed time which verifies the excellent computational efficiency of both schemes. The treatment of limiting projection and SAF reconstruction only consumes the mildest extra overhead which occupies a small portion of the whole advection computation. The CPU consumption increased by WENO limiter is comparable to the MLP limiter where the smooth indicator (51) effectively avoids the implementation of costly multiple reconstructions on smooth solution. We

574

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Table 3 Numerical errors and convergence orders of VPM-Ex scheme with algorithm I and algorithm II for advection test on quadrilateral and triangular grids. DOFs

VPM-I

VPM-II L 1 error

Order

L ∞ error

Order

Quadrilateral

Grids

1/20 1/40 1/80 1/160 1/320

841 3281 12961 51521 205441

1.179 × 10−1 1.587 × 10−2 2.008 × 10−3 2.517 × 10−4 3.178 × 10−5

– 2.89 2.98 3.00 2.99

1.200 × 10−1 1.594 × 10−2 2.010 × 10−3 2.518 × 10−4 3.177 × 10−5

– 2.91 2.99 3.00 2.99

1.119 × 10−1 1.501 × 10−2 1.898 × 10−3 2.380 × 10−4 3.007 × 10−5

– 2.90 2.98 3.00 2.99

1.138 × 10−1 1.508 × 10−2 1.900 × 10−3 2.381 × 10−4 3.006 × 10−5

– 2.92 2.99 3.00 2.99

Triangular

Type

L 1 error

Order

L ∞ error

Order

1/20 1/40 1/80 1/160 1/320

1241 4881 19361 77121 307841

9.314 × 10−2 1.251 × 10−2 1.583 × 10−3 1.983 × 10−4 2.477 × 10−5

– 2.90 2.98 3.00 3.00

7.104 × 10−2 9.361 × 10−3 1.178 × 10−3 1.474 × 10−4 1.841 × 10−5

– 2.92 2.99 3.00 3.00

9.990 × 10−2 1.345 × 10−2 1.702 × 10−3 2.132 × 10−4 2.667 × 10−5

– 2.89 2.98 3.00 3.00

7.622 × 10−2 1.007 × 10−2 1.268 × 10−3 1.586 × 10−4 1.983 × 10−5

– 2.92 2.99 3.00 3.00

Table 4 The same as Table 3 but for the schemes implemented with MLP|WENO limiters and SAF reconstruction. We omit the notation of SAF in the table for brevity. Type

Grids

Quadrilateral

Triangular

L 1 error

Order

L ∞ error

Order

L 1 error

Order

L ∞ error

Order

VPM-I MLP limiter

1/20 1/40 1/80 1/160 1/320

1.788 × 10−1 1.587 × 10−2 2.008 × 10−3 2.517 × 10−4 3.178 × 10−5

– 3.49 2.98 3.00 2.99

2.418 × 10−1 1.594 × 10−2 2.010 × 10−3 2.518 × 10−4 3.177 × 10−5

– 3.92 2.99 3.00 2.99

9.316 × 10−2 1.251 × 10−2 1.583 × 10−3 1.983 × 10−4 2.477 × 10−5

– 2.90 2.98 3.00 3.00

7.105 × 10−2 9.361 × 10−3 1.178 × 10−3 1.474 × 10−4 1.841 × 10−5

– 2.92 2.99 3.00 3.00

VPM-II MLP limiter

1/20 1/40 1/80 1/160 1/320

1.753 × 10−1 1.501 × 10−2 1.898 × 10−3 2.380 × 10−4 3.007 × 10−5

– 3.55 2.98 3.00 2.99

2.399 × 10−1 1.508 × 10−2 1.900 × 10−3 2.381 × 10−4 3.006 × 10−5

– 3.99 2.99 3.00 2.99

9.994 × 10−2 1.345 × 10−2 1.702 × 10−3 2.132 × 10−4 2.667 × 10−5

– 2.89 2.98 3.00 3.00

7.623 × 10−2 1.007 × 10−2 1.268 × 10−3 1.586 × 10−4 1.983 × 10−5

– 2.92 2.99 3.00 3.00

VPM-I WENO limiter

1/20 1/40 1/80 1/160 1/320

1.271 × 10−1 1.587 × 10−2 2.008 × 10−3 2.517 × 10−4 3.178 × 10−5

– 3.00 2.98 3.00 2.99

1.431 × 10−1 1.594 × 10−2 2.010 × 10−3 2.518 × 10−4 3.177 × 10−5

– 3.17 2.99 3.00 2.99

9.314 × 10−2 1.251 × 10−2 1.583 × 10−3 1.983 × 10−4 2.477 × 10−5

– 2.90 2.98 3.00 3.00

7.104 × 10−2 9.361 × 10−3 1.178 × 10−3 1.474 × 10−4 1.841 × 10−5

– 2.92 2.99 3.00 3.00

VPM-II WENO limiter

1/20 1/40 1/80 1/160 1/320

1.203 × 10−1 1.501 × 10−2 1.898 × 10−3 2.380 × 10−4 3.007 × 10−5

– 3.00 2.98 3.00 2.99

1.318 × 10−1 1.508 × 10−2 1.900 × 10−3 2.381 × 10−4 3.006 × 10−5

– 3.13 2.99 3.00 2.99

9.990 × 10−2 1.345 × 10−2 1.702 × 10−3 2.132 × 10−4 2.667 × 10−5

– 2.89 2.98 3.00 3.00

7.622 × 10−2 1.007 × 10−2 1.268 × 10−3 1.586 × 10−4 1.983 × 10−5

– 2.92 2.99 3.00 3.00

have also examined the solution accuracy computed with different CFL numbers and found only negligible differences in numerical errors and convergence rates. 4.2. Advection test on hybrid unstructured grids To further evaluate the performance of proposed model, we investigate the numerical accuracy and convergence behavior of 3D advection test [50] on hybrid unstructured grids. A smooth initial condition φ0 (x, y , z) = sin(2π x) · sin(2π y ) · sin(2π z) is defined in a unit cube and transported by a uniform velocity of u = (1, 1, 1). The boundaries are specified with periodic conditions and the computations are advanced with adaptively chosen time steps to meet a CFL restriction of 0.2. The unstructured grids used in the calculation include four types, i.e. three grids with pure hexahedral, tetrahedral or prismatic elements, and a hybrid grid blending all element types as depicted in Fig. 5. The grids with different resolution are generated by specifying 20, 30, 40, 50, 60 cells along each boundary edge. The statistical information of each element type is tabulated in Table 6. The total number of elements is denoted as N total which varies for different grids. We carried out computation by VPM(-Ex) and VPM(-Ex)-MLP|WENO-SAF schemes for time at t = 1 and present the L 1 errors and the convergence rates in Table 7 and Table 8 for algorithm I and algorithm II respectively. We first compare the results between VPM and VPM-MLP|WENO-SAF schemes. It is seen that all the schemes are verified with approximately 3rd order accuracy while the former version [18] degrades to 2.2th order on the tetrahedral grids. This improvement is partially attributed to the choice of selective stencil which provides sufficient amount of DOFs to realize exactly 3rd order reconstruction. The numerical accuracy of VPM-MLP-SAF scheme is found slightly worse on the coarsest grids and rapidly improved with grid refinement which leads to a relatively larger convergence rate at beginning. This is in accordance with the convergence behavior of smooth indicator function which is analyzed to be utmost eighth order [62]. Besides, the

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

575

Table 5 The elapse time measured on a single core of Intel Xeon(R) E5-2630, 2.20 GHZ for schemes computed in Table 3 and Table 4. Type

Grids

Quadrilateral

Triangular

algorithm I

algorithm II

algorithm I

algorithm II

VPM

1/20 1/40 1/80 1/160 1/320

0.66 3.49 36.79 370.37 3371.03

0.60 3.30 23.80 345.37 3241.46

0.76 4.37 42.57 454.13 4132.01

0.71 4.03 31.85 442.22 3866.67

VPM MLP limiter

1/20 1/40 1/80 1/160 1/320

0.69 3.75 47.96 416.32 3413.96

0.62 3.43 28.87 370.28 3273.89

0.82 4.54 57.46 471.69 4169.96

0.73 4.13 33.17 465.61 3897.25

VPM WENO limiter

1/20 1/40 1/80 1/160 1/320

0.71 3.87 59.84 421.56 3443.76

0.63 3.52 31.84 393.42 3290.12

0.84 4.60 59.52 509.38 4257.54

0.75 4.16 35.32 499.87 3920.91

Fig. 5. Cutaway section of different type of unstructured grid elements for advection test.

numerical errors of VPM-MLP-SAF schemes are slight larger than VPM scheme particularly for the coarser grids, which may be affected by the solution degradation of linear MLP reconstruction. The loss of accuracy is entirely restored by the VPM-WENO-SAF scheme, to an extent that the error introduced by limiting projection is reduced to negligible level. The results demonstrate the major advantage of VPM-WENO-SAF scheme whose accuracy-preserving capability benefits from its subtle design of nonlinear weighting strategy rather than VPM-MLP-SAF scheme that mainly relies on the use of smoothness indicator function. For better understanding, we visualize the data tabulated in Table 6 for numerical errors and computational costs in Fig. 6 measured on a PC with a single CPU of Intel Xeon(R) E5-2630, 2.20 GHZ. It is seen that there is no discernible difference in terms of numerical errors for a given number of DOFs among different mesh types or multi-moment reconstruction algorithms in the VPM-Ex scheme. On the other hand, the elapse time of algorithm I is slightly larger than that of algorithm II, which shows the superiority in numerical efficiency of the latter scheme. These observations apply to VPM-MLP|WENO-SAF as well, where the VPM-Ex scheme forms the major part of the numerical procedures.

576

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Table 6 Element statistics for advection test on different type of unstructured grids. Element

N edge

N hexa

N tetra

N prsim

N pyra

N total

Hexahedral

20 30 40 50 60

8,000 27,000 64,000 125,000 216,000

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

8,000 27,000 64,000 125,000 216,000

Tetrahedral

20 30 40 50 60

0 0 0

0 0 0

0 0 0

0

58,880 190,535 411,104 867,003 1,540,893

0

0

58,880 190,535 411,104 867,003 1,540,893

Prismatic

20 30 40 50 60

0 0 0

0 0 0

0 0 0

0

0

17,960 60,780 143,520 281,200 485,640

0

17,960 60,780 143,520 281,200 485,640

Hybrid

20 30 40 50 60

6,720 22,680 53,760 105,000 181,440

2,630 8,642 19,580 37,886 68,139

1,776 5,832 13,824 26,940 46,800

256 576 1,024 1,600 2,304

11,382 37,730 88,188 171,426 298,683

Table 7 Numerical accuracy and convergence rate of scheme VPM-MLP|WENO-SAF implemented with algorithm I for 3D advection test on hybrid unstructured grids. Type

Grids

Hexahedral

Tetrahedral

Prismatic

Hybrid

VPM

VPM-MLP-SAF

VPM-WENO-SAF

L 1 error

Rate

L 1 error

Rate

L 1 error

Rate

1/20 1/30 1/40 1/50 1/60

1.710 × 10−1 5.524 × 10−2 2.384 × 10−2 1.233 × 10−2 7.171 × 10−3

−

3.028 × 10−1 5.524 × 10−2 2.384 × 10−3 1.233 × 10−2 7.171 × 10−3

−

1.720 × 10−1 5.524 × 10−2 2.384 × 10−2 1.233 × 10−2 7.171 × 10−3

−

1/20 1/30 1/40 1/50 1/60

3.985 × 10−2 1.274 × 10−2 5.955 × 10−3 2.830 × 10−3 1.595 × 10−3

−

3.990 × 10−2 1.274 × 10−2 5.956 × 10−3 2.830 × 10−3 1.595 × 10−3

−

3.986 × 10−2 1.274 × 10−2 5.955 × 10−3 2.830 × 10−3 1.595 × 10−3

−

1/20 1/30 1/40 1/50 1/60

9.345 × 10−2 2.901 × 10−2 1.244 × 10−2 6.399 × 10−3 3.712 × 10−3

−

1.066 × 10−1 2.906 × 10−2 1.244 × 10−2 6.399 × 10−3 3.712 × 10−3

−

9.349 × 10−2 2.901 × 10−2 1.244 × 10−2 6.399 × 10−3 3.712 × 10−3

−

1/20 1/30 1/40 1/50 1/60

1.555 × 10−1 5.013 × 10−2 2.169 × 10−2 1.134 × 10−2 6.544 × 10−3

−

2.709 × 10−1 5.872 × 10−2 2.349 × 10−2 1.185 × 10−2 6.835 × 10−3

−

1.561 × 10−1 5.010 × 10−2 2.168 × 10−2 1.134 × 10−2 6.544 × 10−3

−

2.79 2.92 2.95 2.97

2.91 2.97 2.99 2.99

2.88 2.96 2.97 2.99

2.83 2.96 2.93 2.97

4.20 2.92 2.95 2.97

2.92 2.97 2.99 2.99

3.20 2.96 2.97 2.99

3.83 3.24 3.09 2.97

2.80 2.92 2.95 2.97

2.91 2.97 2.99 2.99

2.88 2.96 2.97 2.99

2.84 2.96 2.92 2.97

4.3. Solid rotation of a complex profile To assess the ability of present schemes to provide high-fidelity reconstruction for both smooth and non-smooth solution, we consider a two-dimensional benchmark problem obtained by a modification of a complex profile introduced in [28] which consists four shapes located within a circle of r0 = 0.15 on a unit square [0, 1]2 . The initial condition is given as

⎧ 1, if |x − 0.5| < 0.25, or y > 0.85 ⎪ ⎪ ⎪ 1 ⎪ ⎪ ⎨ 4 (1 + cos (π min ( r (x, y )/ r0 , 1))) , φ(x, y ) = 1 − r (x, y )/ r0 , ⎪ ⎪ ⎪ 1 ( F (r (x, y ) + δ, α ) + F (r (x, y ) − δ, α ) + 4F (r (x, y ), α )) , ⎪ ⎪ ⎩6 0,

where

|r (x, y )| < r0 |r (x, y )| < r0 |r (x, y )| < r0 |r (x, y )| < r0 otherwise

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

577

Table 8 The same as Table 7 but for the scheme implemented with algorithm II. Type

Grids

Hexahedral

Tetrahedral

Prismatic

Hybrid

VPM

VPM-MLP-SAF

VPM-WENO-SAF

L 1 error

Rate

L 1 error

Rate

L 1 error

Rate

1/20 1/30 1/40 1/50 1/60

1.712 × 10−1 5.526 × 10−2 2.384 × 10−2 1.233 × 10−2 7.170 × 10−3

−

3.010 × 10−1 5.524 × 10−2 2.384 × 10−3 1.233 × 10−2 7.171 × 10−3

−

1.722 × 10−1 5.526 × 10−2 2.384 × 10−2 1.233 × 10−2 7.170 × 10−3

−

1/20 1/30 1/40 1/50 1/60

4.385 × 10−2 1.414 × 10−2 6.671 × 10−3 3.212 × 10−3 1.843 × 10−3

−

4.709 × 10−2 1.540 × 10−2 7.459 × 10−3 3.621 × 10−3 2.083 × 10−3

−

4.383 × 10−2 1.414 × 10−2 6.671 × 10−3 3.212 × 10−3 1.843 × 10−3

−

1/20 1/30 1/40 1/50 1/60

9.502 × 10−2 2.954 × 10−2 1.267 × 10−2 6.521 × 10−3 3.782 × 10−3

−

1.191 × 10−1 2.987 × 10−2 1.319 × 10−2 6.953 × 10−3 4.093 × 10−3

−

9.508 × 10−2 2.954 × 10−2 1.267 × 10−2 6.521 × 10−3 3.782 × 10−3

−

1/20 1/30 1/40 1/50 1/60

1.560 × 10−1 5.022 × 10−2 2.172 × 10−2 1.124 × 10−2 6.550 × 10−3

−

2.721 × 10−1 5.925 × 10−2 2.392 × 10−2 1.218 × 10−2 6.864 × 10−3

−

1.567 × 10−1 5.022 × 10−2 2.172 × 10−2 1.124 × 10−2 6.550 × 10−3

−

2.79 2.92 2.95 2.97

2.89 2.93 2.94 2.90

2.88 2.96 2.96 2.99

2.84 2.96 2.97 2.92

4.20 2.92 2.95 2.97

2.86 2.83 2.91 2.88

3.40 2.96 2.86 2.91

3.83 3.21 3.05 3.10

2.80 2.92 2.95 2.97

2.89 2.97 2.94 2.90

2.88 2.96 2.96 2.99

2.85 2.96 2.97 2.92

Fig. 6. Numerical errors and CPU time of 3D advection test computed by algorithm I and algorithm II of VPM-Ex scheme on different grid resolutions.

r (x, y ) =

 (x − x0 )2 + ( y − y 0 )2 ,

F (r (x, y ), α ) =



max(1 − α 2 r 2 (x, y )),

α = 5 and δ = 0.01.

The center (x0 , y 0 ) of each shape is specified with (0.5, 0.78), (0.22, 0.5), (0.5, 0.22) and (0.78, 0.5) respectively. The initial profiles are depicted in Fig. 7 to include a more complex profile which is challenging to be accurately resolve by the present scheme. We perform the rotation on 57,518 triangles by the velocity field u(x, y ) = (0.5 − y , x − 0.5) and present the solution shapes at time t = 2π in Fig. 8. It is seen that the VPM-MLP-SAF scheme obtains remarkable results for all the shapes which effectively prevent numerical oscillation in presence of discontinuous solution and meanwhile suppress the excessive dissipation introduced by the limiting projection, despite slightly flattens connective vicinity of smooth and non-smooth solution. We plot the weight function ωi at time t = 2π in Fig. 8 for 3D bird-eye’s view. It is shown that the VPM-MLP-SAF scheme precisely identifies the discontinuities with weight values occupying rightly in the cells containing discontinuities. Different from the previous BVD scheme [62] based on dissipation-minimizing approach, the SAF process is designed in a thoroughly different spirit that adaptively chooses optimal weight function to matching the VIAs of neighboring cells which yields straightforward and faithful reconstruction even for rather complicated solution distribution. As shown in

578

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 7. The exact or initial solution of a complex profile on triangular grids with 57,518 elements.

Fig. 8. The 3D bird’s view of solution body and weight function at t = 2π computed by VPM-MLP-SAF scheme.

Fig. 9, the 2D contour of weight function computed by VPM-MLP-BVD scheme is included for comparison. It is observed that the previous BVD scheme [62] probably incur solution degradation with artificial values on unstructured grids, which demonstrate the superiority of present SAF process to reconstruct the complex profiles mixing with smooth and non-smooth solutions. Moreover, the matching degree (55) is defined to further control the vicinity of THINC reconstruction. As an instance, we specify the threshold values Rc = 0 and Rc = 200 in the computation, and plot the weight functions in Fig. 10 respectively. It is shown that the choice of larger value corresponds to the narrower region of identified discontinuities which leads to a more robust calculation. We have also performed the computation with VPM-WENO-SAF scheme and displayed the solutions and weights in Fig. 11 and Fig. 12 for 3D/2D visualization. It is realized that the non-oscillatory solution is reproduced with sufficiently high geometrical faithfulness. Compared with the results by VPM-MLP-SAF scheme, the flattened deformation is largely improved which is attributed to the sophisticated quadratic reconstruction of WENO limiting process. Moreover, we have checked the numerical results under different time step size and do not see noticeable dependence of solution quality on the CFL number.

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

579

Fig. 9. The 2D contour maps of weight function within range of [0.001, 0.999] by BVD and SAF scheme. It is noted that the value below 0.001 is cutoff for visual clarity.

Fig. 10. The comparison of weight function computed by VPM-MLP-SAF scheme with different cutoff numbers Rc .

4.4. Accuracy test for 2D Euler equations To investigate the accuracy of present scheme in solving the non-linear system of Euler equations, we use the following initial condition [43,69].

ρ (x, y , 0) = 1 + 0.2 sin(π (x + y )), u (x, y , 0) = 0.7, v (x, y , 0) = 0.3, p (x, y , 0) = 1.0. The density is disturbed by sinusoidal perturbation and transported by uniform velocity and pressure fields. The exact solution is given as ρ (x, y , t ) = 1 + 0.2 sin(π (x + y − t )). The test is carried out on a square domain of [−1, 1]2 divided by regular quadrilateral and triangular elements with different grid resolution. We use five set of meshes ranging in size of 2/10, 2/20, 2/40, 2/80 and 2/160. The periodic conditions are imposed for all the boundaries. We perform the computation

580

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 11. The 3D bird’s view of solution body and weight function at t = 2π computed by VPM-WENO-SAF scheme.

Fig. 12. The same as Fig. 11 but for 2D contour map.

by VPM-MLP|WENO-SAF schemes implemented with algorithm II upto time t = 2. The numerical errors in terms of L 1 and L ∞ norms of density are given in Table 9 with convergence rates included for examining the orders of accuracy. We can clearly see that the present numerical model with both schemes can exactly ensure 3rd-order accuracy for 2D Euler equation on refined mesh while return different level of errors on the coarser mesh. Particularly on the grid of 2/10, the smoothness indicator function is triggered in a few cells to switch to MLP or WENO limiting process and the errors obtained with WENO limiter are almost 1/2 ∼ 1/3 fewer than those with MLP limiter. The improved numerical accuracy indicates that the WENO limiting process possesses superior accuracy-preserving capability to achieve the utmost accuracy for both smooth and non-smooth solution. In [69], the same test case is computed by Runge–Kutta DG method with or without WENO limiter. On triangular grids with N cells, DG(P1 ) and DG(P2 ) method requires 3N and 6N DOFs in total respectively whereas VPM-MLP|WENO-SAF methods only use 1.5N DOFs that include the VIAs and PVs at cell vertices. In addition, the MLP and WENO limiting projections are naturally implemented in our framework and don’t entail moment recovering process in RKDG-WENO method. The accuracy of present schemes overwhelmingly outperform DG(P1 ) method with only 1/8 solution errors. Compared with DG(P2 ) method, the present errors are even smaller in L ∞ norm while fall behind in L 1 norm, which owes to its nearly 4 times of more DOFs. To further test the accuracy of inviscid Euler equations, we present an alternative example [32] by solving the manufactured solution using the present schemes. We choose the exact solution as

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

581

Table 9 Numerical errors and convergence orders of the Euler equation for initial conditions of sinusoidal perturbation on quadrilateral and triangular grids. DOFs

VPM-MLP-SAF

VPM-WENO-SAF L 1 error

Order

L ∞ error

Order

Quadrilateral

Grids

2/10 2/20 2/40 2/80 2/160

221 841 3281 12961 51521

6.574 × 10−2 1.867 × 10−2 6.417 × 10−4 8.064 × 10−5 1.009 × 10−5

– 1.82 4.86 2.99 3.00

8.840 × 10−1 3.881 × 10−2 8.380 × 10−4 1.055 × 10−4 1.321 × 10−5

– 1.19 5.53 2.99 3.00

3.297 × 10−2 5.101 × 10−3 6.417 × 10−4 8.064 × 10−5 1.009 × 10−5

– 2.69 2.99 2.99 3.00

4.260 × 10−2 6.942 × 10−3 8.380 × 10−4 1.055 × 10−4 1.321 × 10−5

– 2.62 3.05 2.99 3.00

Triangular

Type

L 1 error

Order

L ∞ error

Order

2/10 2/20 2/40 2/80 2/160

321 1241 4881 19361 77121

2.050 × 10−2 1.028 × 10−3 1.293 × 10−4 1.617 × 10−5 2.022 × 10−6

– 4.32 2.99 3.00 3.00

3.772 × 10−2 1.341 × 10−3 1.691 × 10−4 2.117 × 10−5 2.647 × 10−6

– 4.81 2.99 3.00 3.00

7.875 × 10−3 1.028 × 10−3 1.293 × 10−4 1.617 × 10−5 2.022 × 10−6

– 2.94 2.99 3.00 3.00

1.016 × 10−2 1.341 × 10−3 1.691 × 10−4 2.117 × 10−5 2.647 × 10−6

– 2.92 2.99 3.00 3.00

Table 10 Numerical errors and convergence orders of the Euler equation for manufactured solution on quadrilateral and triangular grids. DOFs

VPM-MLP-SAF

VPM-WENO-SAF L 1 error

Order

L ∞ error

Order

Quadrilateral

Grids

1/10 1/20 1/40 1/80 1/160

221 841 3281 12961 51521

3.906 × 10−2 8.206 × 10−3 5.210 × 10−4 6.665 × 10−5 8.717 × 10−6

– 2.25 3.98 2.97 2.93

5.195 × 10−1 1.351 × 10−2 5.498 × 10−4 6.982 × 10−5 9.129 × 10−6

– 1.94 4.62 2.98 2.94

2.408 × 10−2 4.013 × 10−3 5.207 × 10−4 6.665 × 10−5 8.717 × 10−6

– 2.59 2.95 2.97 2.93

2.652 × 10−2 4.793 × 10−3 5.490 × 10−4 6.982 × 10−5 9.129 × 10−6

– 2.47 3.13 2.98 2.94

Triangular

Type

L 1 error

Order

L ∞ error

Order

1/10 1/20 1/40 1/80 1/160

321 1241 4881 19361 77121

4.760 × 10−3 6.758 × 10−4 8.728 × 10−5 1.129 × 10−5 1.515 × 10−6

– 2.82 2.95 2.95 2.90

5.671 × 10−3 7.387 × 10−4 9.501 × 10−5 1.227 × 10−5 1.645 × 10−6

– 2.94 2.96 2.95 2.90

5.020 × 10−3 6.859 × 10−4 8.728 × 10−5 1.129 × 10−5 1.515 × 10−6

– 2.87 2.97 2.95 2.90

5.627 × 10−3 7.496 × 10−4 9.501 × 10−4 1.227 × 10−5 1.645 × 10−6

– 2.91 2.98 2.95 2.90

ρ (x, y , t ) = 2 + sin (2π (t − (x + y ))) , u (x, y , t ) = 0, v (x, y , t ) = 0, p (x, y , t ) = 2 + sin (2π (t − (x + y ))) , with the corresponding source terms

S = [S 1 , S 2 , S 3 , S 4 ]T

⎡

⎤

2π cos (2π (t − (x + y ))) ⎢ −2π cos (2π (t − (x + y ))) ⎥ ⎥ =⎢ ⎣ −2π cos (2π (t − (x + y ))) ⎦ 2π γ −1 cos (2π (t − (x + y ))) added to the right hand side of equations. This test problem is computed over the domain of [0, 1]2 with periodical boundary condition upto time t = 1. Likewise, we use two types of grids, namely regular quadrilateral and triangular grids with different resolution. The VPM-MLP|WENOSAF schemes are implemented with algorithm II and the results are shown in Table 10. We can see that both schemes can achieve the expected 3rd order convergence rate regarding to the errors of density in terms of both L 1 and L ∞ norms. 4.5. A Mach 3 wind tunnel flow with a step This is a classical and well-known test problem, originally proposed by Woodward and Colella in [57] for verifying the capability of high-resolution schemes in computing unsteady shock waves. The problem under consideration is a Mach 3 flow in a wind tunnel with a step where a uniform Mach 3 flow is blown into a wind tunnel of [0, 3] × [0, 1] and the step of 0.2 units length is located at 0.6 units away from the left end of the tunnel. The inflow and outflow conditions are prescribed at the entrance and exit while the reflective conditions are applied along the walls of tunnel. The numerical experiments are conducted with VPM-WENO-SAF scheme on triangular grids with different resolutions. It is noted that we did not modify our schemes nor locally refine the mesh around the corner. The density contours at time t = 4 are plotted in Fig. 13 and Fig. 14 for grids of 1/160 and 1/320 respectively. The shock waves and contact discontinuities are well resolved without numerical oscillations, and meanwhile the numerical dissipation is effectively suppressed to capture

582

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 13. Density contours for the Mach 3 wind tunnel with a forward step at time of t = 4 on a unstructured grid with triangular elements that have 1/160 edge length along the lower domain boundary.

Fig. 14. Same as Fig. 13, but with 1/320 triangular grid elements.

Fig. 15. Same as Fig. 13, but for weight function with 1/320 triangular grid elements.

the vortex structure waves, which exhibit very robust shock capturing performance. Compared with the results in [31,36], it is clearly seen that present third order scheme provides sharper resolution than that of second order scheme. Moreover, the numerical results converge rapidly on a finer mesh and show competitive solution quality compared with other existing high-order schemes in [22,33,36]. At the last time step, the trouble cells identified by SAF reconstruction technique are displayed in Fig. 15 which seems to act properly. 4.6. Double Mach Reflection The double Mach problem is another popular case from [57] which contains strong shock-vortex interaction and rich small-scale structures that are more difficult to resolve. This example is usually used as a test bed to evaluate the capability of numerical schemes to deal with strong shock waves and to reproduce complex vortical structures due to the contact discontinuities in the right part of the shock. A right-moving Mach 10 shock is imposed on the computational domain D = [0, 4] × [0, 1] with a 60◦ angles inclined with respect to the x-axis. The reflection condition is used at the wall lying at

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

583

Fig. 16. Density contours of the double Mach reflection test at t = 0.2 on the grid of triangular elements of 1/240 resolutions.

Fig. 17. Same as Fig. 16, but for results on the grid of 1/480 resolution.

the bottom starting from x = 1/6, while the exact post-shock condition is prescribed for the rest of the boundaries. The top boundary is assigned values for the exact motion of Mach 10 shock. The domain is divided by triangular mesh with two set of grid resolution, namely 1/240 and 1/480 respectively. The computation is carried out with VPM-WENO-SAF scheme until t = 0.2. We first present the density contours at t = 0.2 in Fig. 16 and Fig. 17 for grid of 1/240 and 1/480 respectively. The vortex structures resulted from Kelvin–Helmholtz instabilities are prominently captured along the slip line which are the key features to evaluate the intrinsic numerical dissipation of the tested scheme. These instabilities further develop to roll-up on a finer mesh where the amount of rolling is considered as qualitative indicator to evaluate the numerical dissipation introduced by the scheme. Our numerical results show excellent solution quality which is competitive to the reference solutions computed by other high-order schemes in [22,13,32,69,16]. To examine which cells are reconstructed with THINC/QQ function, we display the weight function of density variable in 3D bird’s view supplemented with density color contours in Fig. 18 and Fig. 19 for coarse and finer grid resolution. It is clearly seen that the SAF reconstruction technique precisely distinguishes the trouble cells in presence of discontinuous solution and adaptively yields appropriate weight function to reproduce a reconstruction function that minimizes the variation errors with neighboring solution distribution. The vicinity detected by SAF technique almost perfectly match the region of shock waves, contact discontinuities and even the smallscale vortex structures. The results significantly outperform the MOOD concept used in [16] where erroneously detected troubled cells are found behind the main shock wave particularly for lower-order scheme. We can realize the extraordinary performance of SAF technique in recognition accuracy and robustness that can resolve both smooth and non-smooth solution with high fidelity. To further confirm the observation, we present two zooms on the interaction zone containing rolling-up vortices in Fig. 20 and Fig. 21 for grid of 1/240 and 1/480 respectively. The color maps of weight function are plotted in addition to density contours for detailed comparison in close-up view. We can clearly see the compact thickness of the shock wave and pronounced resolution of small-scale vortical structure which are often smeared out by methods with larger numerical dissipation. The vortical structures are resolved with richer small-scale details on the finer grid which is comparable to the results in [16] computed by schemes DG(P5 ) and DG(P9 ) with significantly increasing DOFs. It is also found that the weight functions are credibly identified by SAF technique with desirable solution distribution that almost overlaps with detailed structure of density contours. The results computed by VPM-WENO scheme without SAF technique are also included in Fig. 22 which convincingly illustrate the remarkable effort of SAF technique on suppressing numerical dissipation. In a nutshell, the present VPM-WENO-SAF scheme provides a promising numerical model for Euler equations that is able to effectively preserve accuracy of smooth solution and meanwhile resolve discontinuities with significantly improved solution quality.

584

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 18. Double Mach reflection problem at t = 0.2 on the triangular grid of 1/240. The weight function of density variable identified by SAF technique is supplemented with the density color contours in 3D bird’s view.

Fig. 19. Same as Fig. 19, but for results on the grid of 1/480 resolution.

Fig. 20. Closed-up view around the double Mach stem for density contours and color map of weight function on triangular grid of 1/240.

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

585

Fig. 21. Same as Fig. 21, but for results on the grid of 1/480 resolution.

Fig. 22. Closed-up view of density contours computed by VPM-WENO scheme without SAF technique on the triangular grids with element size of 1/240 and 1/480, respectively.

4.7. Viscous flow past a circular cylinder In this section, we consider the unsteady viscous flow past a circular cylinder to test the capability of present scheme to simulate compressible Navier–Stokes equations. The cylinder obstacle is generated with a diameter of d = 1 and surrounded by an uniform flow incoming with the Reynolds number Re = 150 and the Prandtl number P r = 1. The initial condition in the far field is given as (ρ∞ , p ∞ , u ∞ , v ∞ ) = (1, 1 /γ , 0.2, 0) with γ being 1.4 and the Mach number M ∞ = 0.2. As shown in Fig. 23, the computational domain is a circle with diameter D = 400 which is partitioned by hybrid unstructured grid consisting of 42,632 triangular elements and 14,200 quadrilateral elements. The quadrilateral cells are distributed uniformly by 200 boundary points around the circular cylinder with gradually increasing grid density to accurately resolve the effort of viscosity in the boundary layer. The simulation is performed with Sutherland’s law for the viscosity and carried out in parallel by using six cores. It takes nearly 18 hours to arrive at time t = 1200. The time histories of drag and lift coefficients are plotted in Fig. 24 with the enlarged views being displayed in Fig. 25 for further comparison. The results of drag coefficient, lift coefficient and Strouhal number computed by present model are tabulated in Table 11 which compare favorably with reference solutions found in [15,38,68]. Particularly, the Strouhal number in our computation is 0.182 which is exactly the same with values obtained by Müler [38] and agrees well with the ones obtained by DG/FV scheme [68]. The contour plots of the Mach number and viscosity at instant t = 800 are depicted in Fig. 26 and Fig. 27 respectively which clearly reproduce the von Kárman vortex street behind the circular cylinder. The numerical results confirm the significant capability of present model to solve compressible viscous flows with sufficiently high accuracy.

586

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 23. The computational configuration and zoomed view of unstructured grids used for the flow problem around the circular cylinder.

Fig. 24. Time history of drag and lift coefficients at Reynolds number Re = 150.

4.8. Viscous flow around a sphere The test we choose as the last numerical example is to demonstrate the potential capability of present model to simulate three-dimensional problems with complex geometry. The unsteady flow is considered to pass through a unit sphere at Reynolds number Re = 300 and Mach number M ∞ = 0.3. The initial condition of freestream is specified with density ρ∞ = 1 and pressure p ∞ = 1 /γ with γ again being given as 1.4. The computational domain is defined as the union of a half sphere with radius R = 10 and the cylinder of the same radius with extension of 40 in its perpendicular direction. The unit sphere is located in the origin which is also the center of the half sphere. We generate hybrid unstructured grid to discretize the domain where nearly 992,610 hexahedral elements are used in the cylindrical region, 1,807,030 tetrahedral elements in the vicinity of spherical region and 20,310 pyramidal elements in the rest to build the connections. As shown in Fig. 28, the mesh is significantly refined upto a characteristic size of h = 0.014 for the domain composed of sphere and cylinder with a smaller radius, for the interest of vortical structure behind the sphere. We can realize the huge advantage of present model with hybrid unstructured grids to resolve the flow problems with complex geometry, which substantially cut off the grid numbers to improve the computational efficiency. The simulation is performed by present model with Sutherland’s law for the viscosity until t = 1000. It took nearly 180 hours to compute in parallel by using 120 CPUs at the National Supercomputer TH-A1 (in Tianjin, China). The instantaneous contour of Mach number on the xy-plane is first shown in Fig. 29. The magnitude of the Mach number is approximately 0.3 around the sphere which is coincided with our computational condition. One can clearly see the vortex shedding behind the sphere which has a similar pattern compared with the one observed in the experiments [26]. We draw the time histories

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

587

Fig. 25. Zoomed view of time history of drag coefficient (left) and lift coefficient (right).

Table 11 The computed drag coefficient, peak to peak amplitudes of drag and lift coefficients, Strouhal number in comparison with existing results. Scheme

Time-averaged C d

C d peak to peak

C l peak to peak

Strouhal number

Present DGM(4th order) [68] DG/FV(4th order) [68] Müler [38] Dumbser [15]

1.34 1.348 1.349 1.34 1.33

0.052 0.0519 0.0521 0.05228 0.0484

1.058 1.048 1.049 1.0406 1.01

0.182 0.1834 0.183 0.183 0.182

Fig. 26. Mach number contour generated by the von Kárman vortex street at time t = 800.

Fig. 27. Vorticity contour generated by the von Kárman vortex street at time t = 800.

588

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

Fig. 28. The computational mesh of hybrid unstructured grid elements for the viscous flow around the sphere.

of drag and lift coefficients in Fig. 30 with quantitative results given in Table 12. It is seen that the computed drag and lift coefficients are in close agreement with results in reference. The Strouhal number obtained in this simulation is closer to the one reported in [34] and slightly smaller than the one reported in [15] (St = 0.137) which is directly computed from the wave frequency. To provide a clear picture of the vortical structure of the wave shedding, we plot the Q-criterion in Fig. 31 to reveal the regions of vortical motion in the flow field. It is apparent that the hairpin vortex is attached to the surface of the sphere and continues convecting downstream to form a new vortical structure around the legs of the hairpin, which agrees well with the wave pattern given in [26]. The results demonstrate that our numerical model can accurately resolve complicated shedding process of vortical structures within complex geometries. 5. Concluding remarks We have presented a novel numerical model on unstructured grids for compressible viscous flows which is mainly based on newly developed multi-moment reconstruction method, and numerical approaches for limiting projection and diminishing dissipation. Distinguished from solution procedure used in previous VPM [60] scheme, this new method makes reconstructions by least-square method for VIA on global coordinate straightforwardly. The VIA is predicted by finite volume scheme to ensure rigorous conservation and PV is updated efficient in differential form by multi-moment formulation that can be constructed directly with derivatives obtained in the global reconstruction. The so-called VPM-Ex scheme enhances the numerical accuracy by choosing a relatively broader reconstruction stencil to include more DOFs which facilitates us to establish two different reconstruction methods that are suitable for explicit and implicit schemes. Seeking approaches to resolve discontinuous solution for compressible flows, we developed two different limiting projection methods with choices of MLP and WENO schemes to suppress the numerical oscillation in presence of discontinuous solution. For the trouble cells which are detected by an indicator function, MLP scheme reverts the high-order reconstruction to a limited linear polynomial while WENO scheme makes a convex combination of multiple reconstructions constructed with same-order accuracy on primary and biased stencils. Given both VIA and PV as computational variables, the VPM-Ex scheme integrated with WENO limiter can be implemented on a more compact stencils than the conventional WENO schemes based on FVM which largely improves the algorithmic simplicity and computational efficiency. Moreover, SAF method is developed to diminish dissipation by adaptively choosing hyperbolic tangent function for reconstruction of discontinuous solutions. As substantiated in the numerical tests, the proposed reconstruction procedure, called VPM-MLP|WENO-SAF scheme is shown to attain 3rd-order accuracy and offer excellent quality for both smooth and discontinuous solutions on all unstructured grids tested. We have developed various numerical models based on the VPM-MLP|WENO-SAF reconstruction scheme. To evaluate the accuracy of linear equation, we first computed the scalar advection equation with both algorithm I and algorithm II of VPM-Ex schemes. The results show that both schemes can sufficiently maintain the theoretically 3rd-order accuracy. In regard to computational efficiency, algorithm I can be implemented in both explicit and implicit computations, while algorithm II is computationally cheaper and preferable for the explicit computations. We also solved the solid rotation of a complex profile with the proposed VPM-MLP|WENO-SAF schemes which resolve both smooth and discontinuous profiles with high-fidelity solution quality. As verified in the tests of Euler equations, the present model can obtain theoretically

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

589

Fig. 29. The contour of Mach number at Reynolds number Re = 300 on the xy-plane.

Fig. 30. Zoomed view of time history of drag coefficient (left) and lift coefficient (right).

Table 12 The comparisons of drag coefficient, lift coefficient and Strouhal number. Present Johnson and Patel [26] Gassner [20] Li and Ren [34]

Cd

C d

Cl

C l

Strouhal number

0.668 0.656 0.673 0.674

0.0033 0.0035 0.0031 0.003

−0.066 −0.069 −0.065 −0.055

0.016 0.016 0.015 0.013

0.131 0.137 0.135 0.133

Fig. 31. The vortical structure defined by using Q -criteria with value 1 × 10−4 colored by Mach number.

590

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

3rd-order accuracy which substantiates its sufficient accuracy-preserving capability for smooth solution. For cases containing discontinuous solution, it is observed that the present model effectively captures strong shocks with compact thickness and adequately resolves the small-scale vortical structures with sharper resolution, which demonstrates the excellent shockcapturing capability to suppress numerical oscillation and to diminish numerical dissipation. Being a reliable CFD solver for compressible viscous flows, the present model is also used to solve Navier–Stokes equations for two viscous flow cases in 2 and 3 dimensions. It is seen that the vortex shedding regions are clearly reproduced to appear behind the circular or spherical obstacles and the characteristic quantities obtained in the computation show good agreement with existing results. It is convinced that the present numerical model and reconstruction scheme can be further extended to solve more complicated flows problems which are of great practical significance in the industrial applications. Acknowledgements This work was supported in part by National Natural Science Foundation of China (Grant No. 11802178), as well as the funds from JSPS (Japan Society for the Promotion of Science) under Grant Nos. 17K18838 and 18H01366. Appendix A. The Gaussian quadrature used for arbitrary unstructured grids The volume integration is extensively used in the solution reconstruction of VPM-Ex and THINC/QQ scheme and in the computation of smooth indicator function of WENO limiting process which is performed by numerical quadrature. The choice of quadrature points are properly selected to respect the requirement of both solution accuracy and computational efficiency. For the purpose of simplicity, the arbitrary unstructured grid element i is transformed from global coordinate system x = (x, y , z) into standard reference element denoted by local coordinate system ξ = (ξ, η, ζ ). We tabulate cell vertices θik (k = 1, 2, . . . , K ) in Table A.1 whose locations are denoted by (ξik , ηik , ζik ). For 2D, component in third dimension ζ is discarded by default.   We use the Gaussian points ( p g ) in this work whose coordinate and weights are denoted by ξ g = ξ g , η g , ζ g and ω g as summarized in Table A.2. According to the transformation that map from the local coordinate to the global coordinate system, we can determine the Gaussian points in arbitrary unstructured grids by

xg =





xik Nik ξ g



 

where Nik ξ g denote the basis functions regarding to standard reference elements as elaborated in [61]. Table A.1 Local coordinate for each type of standard element. Triangular

Quadrilateral

θi1 : θi2 : θi3 : θi4 :

θi1 : ( 0, 0 ) θi2 : ( 1, 0 ) θi3 : ( 0, 1 )

Tetrahedral

θi1 : θi2 : θi3 : θi4 :

( 0, 0, 0 ) ( 1, 0, 0 ) ( 0, 1, 0 ) ( 0, 0, 1 )

Hexahedral

θi1 : θi2 : θi3 : θi4 : θi5 : θi6 : θi7 : θi8 :

( −1, −1, −1 ) ( 1, −1, −1 ) ( 1, 1, −1 ) ( −1, 1, −1 ) ( −1, −1, 1 ) ( 1, −1, 1 ) ( 1, 1, 1 ) ( −1, 1, 1 )

Prismatic

θi1 : θi2 : θi3 : θi4 : θi5 : θi6 :

( 1, 0, −1 ) ( 0, 0, − 1 ) ( 0, 1, −1 ) ( 1, 0, 1 ) ( 0, 0, 1 ) ( 0, 1, 1 )

( −1, −1 ) ( 1, −1 ) ( 1, 1 ) ( −1, 1 )

Pyramidal

θi1 : θi2 : θi3 : θi4 : θi5 :

( −1, −1, −1 ) ( 1, −1, −1 ) ( 1, 1, −1 ) ( −1, 1, −1 ) ( 0, 0, 1 )

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

591

Table A.2 Configuration of Gaussian points for different type of grid elements. Element Triangular

Locations

Weights

p 1 = (a1 , a1 ), p 4 = (a2 , a2 ),

a1 = a2 =



8− 8−

√ √

p 2 = (1 − 2a1 , a1 ), p 5 = (1 − 2a2 , a2& ),

10 + 10 −

 

√

38 − 44 2/5

√

38 − 44 2/5

ω1 = ω2 = ω3 = b1 , ω4 = ω5 =ω6 = b2 .

18,

b1 =

18.

b2 =

p 1 = (a1 , a1 ), p 4 = (a2 , a1 ), p 7 = (a3 , a1 ), √ a1 = − 3/5,

Tetrahedral

p 1 = (a1 , a1 , a1 ), p 2 = (a3 , a2 , a2 ), p 3 = (a2 , a3 , a2 ), p 4 = (a2 , a2 , a3 ), p 5 = (a2 , a2 , a2 ), p 6 = (a4 , a4 , a5 ), p 7 = (a4 , a5 , a4 ), p 8 = (a4 , a5 , a5 ), p 9 = (a5 , a4 , a4 ), p 10 = (a5 , a4 , a5 ), p 11 = (a5 , a5 , a4 ). a1 = 1/4, a2 = 1/ 14, a3 = 11/14,  √ √ a4 = 1 + 5/14 4, a5 = 1 − 5/14 4.

ω1 = −148/1875, ω2 = ω3 = ω4 = ω5 = 343/7500 ω6 = ω7 = ω8 = ω9 = ω10 = ω11 = 56/375.

p 1 = (a1 , a1 , a1 ), p 4 = (a2 , a1 , a1 ), p 7 = (a3 , a1 , a1 ), p 10 = (a1 , a1 , a2 ), p 13 = (a2 , a1 , a2 ), p 16 = (a3 , a1 , a2 ), p 19 = (a1 , a1 , a3 ), p 22 = (a2 , a1 , a3 ), p 25 = (a3 , a1 , a3 ), √ a1 = − 3/5,

p 2 = (a1 , a2 , a1 ), p 5 = (a2 , a2 , a1 ), p 8 = (a3 , a2 , a1 ), p 11 = (a1 , a2 , a2 ), p 14 = (a2 , a2 , a2 ), p 17 = (a3 , a2 , a2 ), p 20 = (a1 , a2 , a3 ), p 23 = (a2 , a2 , a3 ), p 26 = (a3 , a2 , a3 ), a2 = 0,

p 3 = (a1 , a3 , a1 ), p 6 = (a2 , a3 , a1 ), p 9 = (a3 , a3 , a1 ). p 12 = (a1 , a3 , a2 ), p 15 = (a2 , a3 , a2 ), p 18 = (a3 , a3 , a2 ). p 21 = (a1 , a3 , a3 ), p 24 = (a2 , a3 , a3 ), p 27 = (a3 , a3 , a3 ). √ a3 = 3/5.

ω1 = ω3 = ω7 = ω9 = ω19 = ω21 = ω25 = ω27 = 125/5832, ω2 = ω4 = ω6 = ω8 = ω10 = ω12 = ω16 = ω18 = ω20 = ω22 = ω24 = ω26 = 200/5832, ω11 = ω13 = ω15 = ω17 = ω5 = ω14 = ω23 = 512/5832,

p 1 = (a1 , a1 , c 1 ), p 4 = (a2 , a2 , c 1 ), p 7 = (a1 , a1 , c 2 ), p 10 = (a2 , a2 , c 2 ), p 13 = (a1 , a1 , c 3 ), p 16 =(a2 , a2 , c 3 ),

p 2 = (b1 , a1 , c 1 ), p 5 = (b2 , a2 , c 1 ), p 8 = (b1 , a1 , c 2 ), p 11 = (b2 , a2 , c 2 ), p 14 = (b1 , a1 , c 3 ), p 17 = (b2 , a2& , c 3 ),

p 3 = (a1 , b1 , c 1 ), p 6 = (a2 , b2 , c 1 ), p 9 = (a1 , b1 , c 2 ), p 12 = (a2 , b2 , c 2 ), p 15 = (a1 , b1 , c 3 ), p 18 = (a2 , b2 , c 3 ),

Prismatic

a1 = a2 =



8− 8−

√

√ √

10 + 10 −

c 1 = − 3/5, Pyramidal

 

p 3 = (a1 , a3 ), p 6 = (a2 , a3 ), p 9 = (a3 , a3 ). √ a3 = 3/5.

√

620+ 213125−53320 10 , 3720  √ 620− 213125−53320 10 . 3720

Quadrilateral

Hexahedral

p 2 = (a1 , a2 ), p 5 = (a2 , a2 ), p 8 = (a3 , a2 ), a2 = 0,

&

p 3 = (a1 , 1 − 2a1 ), p 6 = (a2 , 1 − 2a2 ).

√

38 − 44 2/5

√

38 − 44 2/5 c 2 = 0,

&

ω1 = ω3 = ω7 = ω9 = 25/324, ω2 = ω4 = ω6 = ω8 = 40/324, ω5 = 64/324.

18, b1 = 1 − 2a1 , 18, b2 = 1 − 2a2 , c3 =

√

3/5.

p 2 = (a2 , a1 , a7 ), p 3 = (a2 , a2 , a7 ), p 1 = (a1 , a1 , a7 ), p 4 = (a1 , a2 , a7 ), p 5 = (a3 , a3 , a8 ), p 6 = (a4 , a3 , a8 ), p 7 = (a4 , a4 , a8 ), p 8 = (a3 , a4 , a8 ), p 9 = (a5 , a5 , a9 ), p 10 = (a6 , a5 , a9 ), p 11 = (a6 , a6 , a9 ), p 12 = (a5 , a6 , a9 ). a1 = −0.5352071497417589, a2 = 0.5352071497417589, a3 = −0.3770075514575694, a4 = 0.3770075514575694, a5 = −0.1703170535312202, a6 = 0.1703170535312202, a7 = −0.8540119518537005, a9 = 0.4100044197769968, a8 = −0.3059924679232962.

ω1 = ω2 = ω3 = ω13 = ω14 = ω15 = 5b1 / 18, ω4 = ω5 = ω6 = ω16 = ω17 = ω18 = 5b2 / 18. ω7 = ω8 = ω9 = 4b1 / 9, ω10 = ω11= ω12 = 4b2 / 9.

b1 = b2 =

√

620+ 213125−53320 10 , 3720  √ 620− 213125−53320 10 . 3720

ω1 = ω2 = ω3 = ω4 = b1 , ω5 = ω6 = ω7 = ω8 = b2 , ω9 = ω10 = ω11 = ω12 = b3 .

b1 = 0.1178522707986649, b2 = 0.1096847019448994, b3 = 0.0224630272564355.

References [1] R. Abgrall, On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation, J. Comput. Phys. 114 (1994) 45–58. [2] T.J. Barth, P.O. Frederickson, Higher Order Solution of the Euler Equations on Unstructured Grids Using Quadratic Reconstruction, AIAA Paper, 1990, 90-0013. [3] T.J. Barth, Recent Developments in High Order K-Exact Reconstruction on Unstructured Meshes, AIAA Paper, 1993, 93-0668. [4] J. Blazek, Computational Fluid Dynamics: Principles and Applications, Butterworth-Heinemann, 2005. [5] J.C. Chassaing, S. Khelladi, X. Nogueira, Accuracy assessment of a high-order moving least squares finite volume method for compressible flows, Comput. Fluids 71 (2013) 41–53. [6] A.J. Christlieb, Y. Liu, Q. Tang, Z. Xu, High order parametrized maximum-principle-preserving and positivity-preserving WENO schemes on unstructured meshes, J. Comput. Phys. 281 (2015) 334–351. [7] B. Cockburn, 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, C.W. Shu, TVB Runge–Kutta local projection discontinuous Galerkin finite element method for conservation laws III: onedimensional systems, J. Comput. Phys. 84 (1989) 90–113. [9] B. Cockburn, S. Hou, 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, C.W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws V: the multidimensional systems, J. Comput. Phys. 141 (1998) 199–224. [11] M. Delanaye, Y. Liu, Quadratic Reconstruction Finite Volume Schemes on 3D Arbitrary Unstructured Polyhedral Grids, AIAA Paper, vol. 9, 1999, p. 3259. [12] M. Dumbser, M. Käser, Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems, J. Comput. Phys. 221 (2007) 693–723.

592

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

[13] M. Dumbser, M. Käser, V.A. Titarev, E.F. Toro, Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinear hyperbolic systems, J. Comput. Phys. 226 (2007) 204–243. [14] M. Dumbser, D.S. Balsara, E.F. Toro, C.D. Munz, A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes, J. Comput. Phys. 141 (2008) 8209–8253. [15] M. Dumbser, Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier–Stokes equations, Comput. Fluids 39 (2010) 60–76. [16] M. Dumbser, O. Zanotti, R. Loubère, S. Diot, A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws, J. Comput. Phys. 278 (2014) 47–75. [17] X. Deng, B. Xie, F. Xiao, Multi-moment finite volume solver for Euler equations on unstructured grids, AIAA J. 55 (2017) 2617–2629. [18] X. Deng, B. Xie, F. Xiao, A finite volume multi-moment method with boundary variation diminishing principle for Euler equation on three-dimensional hybrid unstructured grids, Comput. Fluids 153 (2017) 85–101. [19] O. Friedrich, Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids, J. Comput. Phys. 144 (1998) 194–212. [20] G. Gassner, Discontinuous Galerkin Method for the Unsteady Compressible Navier–Stokes Equations, PhD thesis, University Stuttgart, German, 2009. [21] A. Harten, B. Enquist, S. Osher, S.R. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes, III, J. Comput. Phys. 71 (1987) 231–303. [22] C. Hu, C.W. Shu, Weighted essentially non-oscillatory schemes on triangular meshes, J. Comput. Phys. 150 (1999) 97–127. [23] H.T. Huynh, A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods, AIAA Paper, 2007, 4079–2007. [24] L. Ivan, C.P.T. Groth, High-order solution-adaptive central essentially non-oscillatory (CENO) method for viscous flows, J. Comput. Phys. 257 (2014) 830–862. [25] G. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [26] T.A. Johnson, V.C. Patel, Flow past a sphere up to a Reynolds number of 300, J. Fluid Mech. 378 (1999) 19–70. [27] K. Kitamura, E. Shima, Towards shock-stable and accurate hypersonic heating computations: a new pressure flux for AUSM-family schemes, J. Comput. Phys. 245 (2013) 62–83. [28] R.J. LeVeque, High-resolution conservative algorithms for advection in incompressible flow, SIAM J. Numer. Anal. 33 (1996) 627–665. [29] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [30] H. Luo, L. Luo, R. Nourgaliev, V.A. Mousseau, N. Dinh, A reconstructed discontinuous Galerkin method for the compressible Navier–Stokes equations on arbitrary grids, J. Comput. Phys. 229 (2010) 6961–6978. [31] W. Li, Y.X. Ren, G. Lei, H. Luo, The multi-dimensional limiters for solving hyperbolic conservation laws on unstructured grids, J. Comput. Phys. 230 (2011) 7775–7795. [32] W. Li, Y.X. Ren, The multi-dimensional limiters for solving hyperbolic conservation laws on unstructured grids II: extension to high order finite volume schemes, J. Comput. Phys. 231 (2012) 4053–4077. [33] W. Li, Y.X. Ren, High-order k-exact WENO finite volume schemes for solving gas dynamic Euler equations on unstructured grids, Int. J. Numer. Methods Fluids 70 (2012) 742–763. [34] W. Li, Efficient Implementation of High-Order Accurate Numerical Methods on Unstructured Grids, Springer, Berlin Heidelberg, 2014. [35] H.X. Liu, X.M. Jiao WLS-ENO, Weighted-least-squares based essentially non-oscillatory schemes for finite volume methods on unstructured meshes, J. Comput. Phys. 314 (2016) 749–773. [36] Y. Liu, W. Zhang, Accuracy preserving limiter for the high-order finite volume method on unstructured grids, Comput. Fluids 149 (2017) 88–99. [37] C. Michalak, C.F. Ollivier-Gooch, Accuracy preserving limiter for the high-order accurate solution of the Euler equations, J. Comput. Phys. 228 (2009) 8693–8711. [38] B. Müler, High order numerical simulation of Aeolian tones, Comput. Fluids 37 (2008) 450–462. [39] C.F. Ollivier-Gooch, Quasi-ENO schemes for unstructured meshes based on unlimited data-dependent least-squares reconstruction, J. Comput. Phys. 133 (1997) 6–17. [40] C. Ollivier-Gooch, M.V. Altena, A high-order-accurate unstructured mesh finite-volume scheme for the advection–diffusion equation, J. Comput. Phys. 181 (2002) 729–752. [41] A. Nejat, C. Ollivier-Gooch, A high-order accurate unstructured finite volume Newton–Krylov algorithm for inviscid compressible flows, J. Comput. Phys. 227 (2008) 2582–2609. [42] J.S. Park, S.H. Yoon, C. Kim, Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids, J. Comput. Phys. 229 (2010) 788–812. [43] J.X. Qiu, C.W. Shu, Hermite WENO schemes and their application as limiters for Runge–Kutta discontinuous Galerkin method II: two dimensional case, Comput. Fluids 34 (2005) 642–663. [44] A. Rohde, Eigenvalues and Eigenvectors of the Euler Equations in General Geometries, AIAA paper, 2001, 2609. [45] J. Shi, C. Hu, C.W. Shu, A technique of treating negative weights in WENO schemes, J. Comput. Phys. 175 (2002) 108–127. [46] C. Shen, F. Sun, X.L. Xia, Implementation of density-based solver for all speeds in the framework of OpenFOAM, Comput. Phys. Commun. 185 (2014) 2730–2741. [47] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [48] T. Sonar, On the construction of essentially non-oscillatory finite volume approximations to hyperbolic conservation laws on general triangulations: polynomial recovery, accuracy and stencil selection, Comput. Methods Appl. Mech. Eng. 140 (1997) 157–181. [49] C.G.N.S. The, Steering Sub-committee, The CFD General Notation System Standard Interface Data Structures, AIAA, 2002. [50] P. Tsoutsanis, V.A. Titarev, D. Drikakis, WENO schemes on arbitrary mixed-element unstructured meshes in three space dimensions, J. Comput. Phys. 230 (2011) 1585–1601. [51] Z.J. Wang, Spectral (finite) volume method for conservation laws on unstructured grids: basic formulation, J. Comput. Phys. 178 (2002) 210–251. [52] 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. [53] Z.J. Wang, L.P. Zhang, Y. Liu, Spectral (finite) volume method for conservation laws on unstructured grids IV: extension to two-dimensional Euler equations, J. Comput. Phys. 194 (2004) 716–741. [54] Z.J. Wang, Y. Liu, G. May, A. Jameson, Spectral difference method for unstructured grids II: extension to the Euler equations, J. Sci. Comput. 32 (2007) 45–71. [55] Z.J. Wang, H. Gao, A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids, J. Comput. Phys. 228 (2009) 8161–8186. [56] Q. Wang, Y.X. Ren, W. Li, Compact high order finite volume method on unstructured grids II: extension to two-dimensional Euler equations, J. Comput. Phys. 314 (2016) 883–908. [57] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1) (1984) 115–173. [58] F. Xiao, R. Akoh, S. Ii, Unified formulation for compressible and incompressible flows by using multi-integrated moments II: multi-dimensional version for compressible and incompressible flows, J. Comput. Phys. 213 (2006) 31–56.

B. Xie et al. / Journal of Computational Physics 394 (2019) 559–593

593

[59] B. Xie, F. Xiao, Two and three dimensional multi-moment finite volume solver for incompressible Navier–Stokes equations on unstructured grids with arbitrary quadrilateral and hexahedral elements, Comput. Fluids 227 (2014) 40–54. [60] B. Xie, S. Ii, A. Ikebata, F. Xiao, A multi-moment finite volume method for incompressible Navier–Stokes equations on unstructured grids: volumeaverage/point-value formulation, J. Comput. Phys. 227 (2014) 138–162. [61] B. Xie, F. Xiao, A multi-moment constrained finite volume method on arbitrary unstructured grids for incompressible flows, J. Comput. Phys. 327 (2016) 747–778. [62] B. Xie, X. Deng, Z.Y. Sun, F. Xiao, A hybrid pressure-density-based Mach uniform algorithm for Euler equations on unstructured grids by using multimoment finite volume method, J. Comput. Phys. 335 (2017) 637–663. [63] B. Xie, F. Xiao, Accurate and robust PISO algorithm on hybrid unstructured grids using multi-moment finite volume method, Numer. Heat Transf., Part B 71 (2017) 146–172. [64] B. Xie, P. Jin, F. Xiao, An unstructured-grid numerical model for interfacial multiphase fluids based on multi-moment finite volume method formulation and THINC method, Int. J. Multiph. Flow 89 (2017) 375–398. [65] B. Xie, F. Xiao, Toward efficient and accurate interface capturing on arbitrary hybrid unstructured grids: the THINC method with quadratic surface representation and Gaussian quadrature, J. Comput. Phys. 349 (2017) 415–440. [66] L.P. Zhang, W. Liu, L. He, et al., A class of hybrid DG/FV methods for conservation laws I: basic formulation and one-dimensional systems, J. Comput. Phys. 231 (2012) 1081–1103. [67] L.P. Zhang, W. Liu, L. He, et al., A class of hybrid DG/FV methods for conservation laws II: two-dimensional cases, J. Comput. Phys. 231 (2012) 1104–1120. [68] L.P. Zhang, W. Liu, M. Li, X. He, H.X. Zhang, A class of DG/FV hybrid schemes for conservation law IV: 2D viscous flows and implicit algorithm for steady cases, Comput. Fluids 97 (2014) 110–125. [69] J. Zhu, X.H. Zhong, C.W. Shu, J.X. Qiu, Runge–Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes, J. Comput. Phys. 248 (2013) 200–220.