An improved third-order HWCNS for compressible flow simulation on curvilinear grids

Due to the very high requirements on the quality of computational grids, stability property and computational efficiency, the application of high-order schemes to complex flow simulation is greatly constrained. In order to solve these problems, the third-order hybrid cell-edge and cell-node weighted compact nonlinear scheme (HWCNS3) is improved by introducing a new nonlinear weighting mechanism. The new scheme uses only the central stencil to reconstruct the cell boundary value, which makes the convergence of the scheme more stable. The application of the scheme to Euler equations on curvilinear grids is also discussed. Numerical results show that the new HWCNS3 achieves the expected order in smooth regions, captures discontinuities sharply without obvious oscillation, has higher resolution than the original one and preserves freestream and vortex on curvilinear grids.


Introduction
In the past three decades, a series of nonlinear schemes with high accuracy and high resolution were developed for the numerical simulation of complex flow fluid with shock wave and multiscale structures, such as turbulence and aeroacoustics problems. Among them, the weighted non-oscillatory (WENO) scheme and the weighted compact nonlinear scheme (WCNS) are both applied more extensively. The WCNS developed by Deng et al. [10,11,12], which is derived from the classical compact scheme on the base of the weighted considerations, is a finite difference method with high accuracy, including the explicit form and the implicit form. The classical WCNS procedure consists of three steps: (1) node-to-midpoint weighted nonlinear interpolation of flow variables, (2) fluxes evaluation at midpoints, and (3) midpoint-to-node centered flux differencing. It has not only the high precision/high resolution properties of the spectral-type compact difference scheme, but also the capability of capturing discontinuities automatically of the WENO scheme. Therefore, the WCNS gets lots of attention in the field of the computational fluid dynamics. Nonomura et al.[1] compared the freestream and vortex preservation properties of the classic WENO and WCNS on curvilinear grids, and pointed out that the WCNS works better than the WENO in the generalized coordinate system. In 2008, two group [14,15] explored higher order WCNS and expanded the scheme to the ninth order, respectively. To further improve the resolution of the WCNS, Yan and Liu [17] developed the Y-type nonlinear weights and proposed a new seventh order compact nonlinear interpolation method based on the same stencil as the fifth order WENO scheme. To recover the optimal numerical resolution in smooth flow field, Zhang etc. [18] recently introduced an ENO-like stencil selection procedure to the WCNS, which fully abandons the oscillatory stencils crossing discontinuities and directly applies optimal linear weights when the flow field is smooth.
In 2011, Deng et al. developed the hybrid cell-edge and cell-node weighted compact nonlinear schemes (HWCNS) [13] based on the original WCNS, and put forward the conservation metric method (CMM) to make the scheme follow the geometrical conservation law. The HWCNS is a hybrid scheme of cell-edge and cell-node, with advantages of narrower difference stencil, higher efficiency and better stability. Furthermore, the HWCNS satisfies the geometric conservation law, and it is more suitable for the simulation of flow field with complex geometric shape. Considering that the non-uniqueness form of CMM may cause serious computational instability and inaccuracy for high-order schemes, a symmetrical conservative metric method (SCMM) is proposed based on the discussions of the metrics and Jacobian in [16]. It makes the calculations based on the high-order WCNS schemes around complex geometry flows possible and somewhat easy. In [2], Nonomura et al. modified the flux differencing step of WCNS and developed the midpointand-node-to-node difference(MND) scheme. Their modification makes the scheme more robust, but at the same time, more dissipative. And so far, the WCNS has been successfully applied to the numerical simulation of the flow field in the low speed, the supersonic speed and the hypersonic speed, and show the good performance in the computation of the complicated flow field.
Compared with the fifth-order and higher-order WCNS, the advantages of the third-order scheme are robust for shock problem, such as fewer grid points, less computational cost. However, the computational results show that the classical WCNS using the JS weight of Jiang and Shu [4] is too dissipative. To overcome the shortcomings, a stencil-selection procedure is introduced to improve the nonlinear weight of the third-order WCNS in [19]. Since the proposed scheme has reduced the dissipation and dispersion errors, compared with the WCNSs using typical nonlinear weights, there is a parameter to choose during the stencil-selection procedure. In [3], Zhu first design a new type of finite difference WENO schemes basing on the original idea of the multiresolution methods [20,21]. For the third-order nonlinear interpolation, only the information defined on one one-point central stencil and one three-point central stencil are used. If the solution is smooth enough, the information defined on such two central spatial stencils are combined together to obtain a third-order approximation. When there is a discontinuity in the three-point central stencil, the information of the three-point stencil is effectively abandoned and the approximation order degrades to one. In this work we incorporate Zhu's nonlinear interpolation procedure to improve the HWCNS and propose a third-order hybrid cell-edge and cell-node weighted compact nonlinear schemes, which is simple and efficient in practical applications. In addition, the new method is applied to the governing equations under curvilinear coordinate system and its freestream and vortex preservation properties are verified. Our numerical experiments show that the improved method has a stronger robustness and its dissipation at shock wave is slightly smaller than that of the original method.
The other parts of the paper are organized as follows. In Section 2, the conventional HWCNS3 on the curvilinear coordinates for hyperbolic conser-vation laws are reviewed. Section 3 presents the improved HWCNS3. Series of numerical examples are given in Section 4. Concluding remarks are provided in Section 5.

Governing equations in curvilinear coordinates
In this subsection, the three-dimensional Euler equations in the Cartesian and generalized coordinate systems and the metrics used in the generalized coordinate system are described. In addition, the numerical technique that can preserve the freestream conditions for the HWCNS is also presented.
The three-dimensional Euler equations in the Cartesian coordinate system are expressed as ∂U ∂t where are conservative variables and convective fluxes, respectively. Here ρ, u, v, w, p and E are density,components of velocity in the x, y and z coordinate directions, pressure and total energy, respectively. By assuming that the fluid obeys the perfect gas state equation, p can be calculated as where γ indicates the specific heat ratio of the gas. When a curvilinear grid is used for numerical computation, the governing equation(1) is first transformed into generalized coordinates system (ξ, η, ζ) with the following relationships

The transformed equation can be written as
whereŨ (3)

Original HWCNS
The semi-discrete approximation of the governing equation(2) at a grid point indexed as (i, j, k) is as follows: are the approximation to the first spatial derivatives of the fluxesF,G andH with respect to ξ, η and ζ, respectively. For HWCNS3, (F ′ ξ ) i,j,k is computed by the fourth-order MND formula[2]: whereF i± 1 2 ,j,k are the numerical fluxes at half points, the local Lax-Friedrich flux is used here: where α is the maximum eigenvalue of the Jacobian matrix, the cell-boundary reconstructionŨ L i± 1 2 ,j,k andŨ R i± 1 2 ,j,k are computed by the third-order nonlinear interpolation in characteristic space. Here, only the construction of U L i+ 1 2 ,j,k is presented, and the subscripts "j, k"are dropped for simplicity. Denotes I m and r m the mth left and right eigenvector of the Jacobi matrix A = ∂F ∂U . At first, conservative variables are transformed into characteristic variables, the mth characteristic variable is computed as follows: Then U L i+ 1 2 ,m is constructed by the following nonlinear interpolation: whereÛ and the nonlinear weights ω 1,m and ω 2,m are determined as follows with d 1 = 1 4 and d 2 = 3 4 are the linear weights, ε = 10 −6 is a small positive number to avoid the denominator to become zero, and the smoothness indicators are Finally, the characteristic form of the interpolated value is transformed into the conservative form:Ũ

Node-to-midpoint interpolation
The computational result show that the original HWCNS scheme is too dissipative. To overcome this shortcoming, we propose a new nonlinear reconstruct procedure with similar ideas that firstly presented in [3].
We still use the construction of U L i+ 1 2 ,m as an example. Two central stencils are chosen: T 1 = {i} and T 2 = {i − 1, i, i + 1}, the corresponding reconstruct polynomials are denoted as q 1,m (x) and q 2,m (x), For consistency, the reconstruct polynomials are normalized to obtain the following interpolation polynomials: where γ 1 and γ 2 are linear weights with γ 1 + γ 2 = 1. Based on a balance between the sharp and essentially non-oscillatory shock transitions in nonsmooth regions and accuracy in smooth regions, we set γ 1 = 1 11 , γ 2 = 10 11 . The smoothness indicator β 2,m of p 2,m (x) is computed according to the definition in [4]: While for the zeroth-degree polynomial p 1,m (x), the measurement of smoothness is different: 10 11 , IS 1,m < IS 2,m , with IS 1,m and IS 2,m are defined in (8) and (9). The nonlinear weights are defined by where τ m defined as related to the absolute deference between the smoothness indicators: and ε = 10 −6 is a small number to avoid the denominator becoming zero. The final reconstructed variable at cell edges is given by

Freestream preservation property
In uniform flow regions, U, F, G and H are constants, (2) are simplified as ∂Ũ ∂t where By substituting (3) into (17), we can see that which means the uniform flow conditions are hold: The condition (18) is called freestream condition. To preserve it numerically, we firstly rewrite the metrics in an equivalent form as suggested in [6]:   Then, Eq.
(2) at grid node (i, j, k) are discretized as where For l = 1, 2, 3, D ξ l , D η l and D ζ l denote the finite difference operators in ξ, η and ζ coordinate direction respectively.
As proven in [7], the freestream conditions are hold if D 1 and D 2 are exchangeable. Furthermore, D 2 = D 3 in practical calculations can prevent large numerical errors. Therefore, D 1 , D 2 and D 3 are all used the fourthorder compact differential operator (5). The values at midpoint are obtained analytically or by the fourth-order Lagrangian interpolation, for example,

Numerical experiments
In this section, a variety of benchmark problems are given to illustrate the behavior of the present method, compared with the original HWCNS3. The CFL number is taken as 0.6, except for the accuracy tests where smaller time step is taken to ensure that spatial errors domainate.

The order of numerical accuracy
Example 1: We consider the one-dimensional Euler equations on [0, 2π] with periodic boundary conditions and initial conditions: ρ(x, 0) = 1 + 0.2 sin x, u(x, 0) = 1, p(x, 0) = 1, γ = 1.4. This problem has an exact solution of ρ(x, t) = 1 + 0.2 sin(x − t). We compute the solution up to t = 1 using the new HWCNS3, the error and order of convergence are listed in Table 1. We can see that the third order of accuracy is achieved.
The exact solution of this problem is ρ(x, y, t) = 1 + 0.2 sin(x + y − 2t). The final time is t = 1. Table 2 lists the error and order of convergence of density. We can see that the third order of accuracy is also achieved for this two-dimensional case.

Shock capturing ability
Example 3: We consider two shock tube problems of the one-dimensional Euler equations. The Sod shock tube problem involves a right-moving shock of Mach number 1.7. Its initial conditions are given as The computed density ρ of the two problems are presented in Fig.1 and Fig.2 respectively. We can observe that the new method has better resolution than the original HWCNS3 at the neigbhood of shocks for these two problems.
The boundaries at x = 0 and x = 1 are solid walls with reflecting boundary conditions. At the final time t = 0.038, the flow field has three contact discontinuities. The middle one, which is generated by the two shocks interacting with each other, is very difficult to resolve. Fig.3 shows the results with N = 600 grid points. Once again, it appears that the new HWCNS3 works better than the original one for this example.

Freestream and vortex preservation
Example 7: The freestream is tested on a three-dimensional wavy grid, as shown in Fig.7. This grid is genenrated by where L x = L y = L z = 4, A x = A y = A z = 1, n xy = n xz = n yz = n yx = n zx = n zy = 4, and x min = −L x /2, y min = −L y /2, z min = −L z /2. The grid resolution is set to 21 × 21 × 21.
The initial condition of an ideal gas is given as where γ = 1.4 is the specific heat ratio. We take ∆t = 0.03 and compute 200 steps. The L 2 errors of the velocity components v and w of the flow field are shown in Table 3. We can clearly observe that the L 2 errors are all at the level of roundoff error for different precisions, and verify the expected freestream preservation property accordingly.   We consider a two-dimensional moving vortex problem on the following wavy grid to test the vortex preservation property where L x = L y = 20, A x × ∆x 0 = A y × ∆y 0 = 0.6, n xy = n yx = 8, and x min = −L x /2, y min = −L y /2. Initially, the mean flow is ρ = 1, u = 0.5, v = 0 and p = 1 γ . Then we add an isentropic vortex to it with perturbation centered at (x 0 , y 0 ) = (0, 0) in (u, v) with T = p/ρ, and no perturbation in entropy S = p/ρ γ , (δu, δv) = ετ e α(1−τ 2 ) (sin θ, − cos θ), δT = − (γ − 1)ε 2 4αγ e 2α(1−τ 2 ) , δS = 0.
where τ = r rc ,r = √ (x − x 0 ) 2 + (y − y 0 ) 2 and r c = 1. The exact solution of this problem is smooth and periodic. We set α = 0.204, ε = 0.02 and use 5 different grids to evaluate the scheme convergence. The time-step sizes ∆t respect to those grids are selected as ∆t = 0.2 min(∆x, ∆y). The numerical results at T = 40 are given in Fig.8. It can be observed that the new HWCNS3 is able to resolve the moving vortex. In Table 4, grid convergence rates with the wavy grid are shown. As expected, the formal third order of accuracy is obtained.

Supersonic flow past a cylinder
Example 9: Finally, the supersonic flow past a cylinder is simulated on a randomized curvilinear grid [9]. The grid is generated as: inner points, 0, boundary points.
with θ = 5π 12 , R x = 3, R y = 6, i max = 81 and j max = 61. The Mach 2 supersonic flow moves toward the cylinder from left. The reflective boundary condition is applied on the cylinder surface, the supersonic inflow condition at the left boundary and the supersonic outflow condition at the other boundaries. For this problem, the new HWCNS3 is more robust than the original one. We can take time-step ∆t = 0.01 for the new scheme, while the original scheme will divergent under the same condition. The grid and pressure contours are given in Fig.9. We can see that the results have a good agreement with that of [9].

Conclusions
In this paper, we construct a third-order HWCNS for the compressible Euler equations on curvilinear grids. The new scheme uses the information defined on two nested central spatial stencils to compute the cell boundary values, which makes it more stable. This nonlinear weighting mechanism can also be extended to the WCNS with arbitrary high order accuracy. Various numerical analyses have shown that the presented scheme is capable of simulating complex flow problems, including curvilinear geometric boundary, strong shock waves and small-scale structures. Further applications of the proposed method are expected.

Availability of data and materials
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Funding
The work is supported by the Basic Research Foundation of the National Numerical Wind Tunnel Project (Grant No. NNW2018-ZT4A08) and the National Key Project (Grant No. GJXM92579) of China.

Authors' contributions
Cheng Mingyang was a major contributor in writing the manuscript. Tang Lingyan carried out the numerical simulation. Liu Yu directed the numerical simulation and processed the calculated data. Zhu Huajun reviewed and revised the manuscript. All authors read and approved the final manuscript.