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 essentially non-oscillatory (WENO) scheme and the weighted compact nonlinear scheme (WCNS) are both applied more extensively. The WCNS developed by Deng et al. [1][2][3], which is derived from the classical compact scheme on the basis 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 computational fluid dynamics. Nonomura et al. [4] 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 groups [5,6] explored higher order WCNS and expanded the scheme to the ninth order, respectively. To further improve the resolution of the WCNS, Yan and Liu [7] 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 et al. [8] 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) [9] 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 fields with complex geometric shapes. 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 [10]. It makes the calculations based on the high-order WCNS schemes around complex geometry flows possible and somewhat easy. In [11], Nonomura et al. modified the flux differencing step of WCNS and developed the midpoint-and-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 shown good performance in the computation of the complicated flow field.
Compared with the fifth-order and higher-order WCNS, the advantages of the thirdorder scheme are robust for shock problems, 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 [12] is too dissipative. To overcome the shortcomings, a stencilselection procedure was introduced to improve the nonlinear weight of the third-order WCNS in [13]. 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 [14], Zhu first designed a new type of finite difference WENO schemes based on the original idea of the multiresolution methods [15,16]. For the third-order nonlinear interpolation, only the information defined on one one-point central stencil and one three-point central stencil is used. If the solution is smooth enough, the information defined on such two central spatial stencils is 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 scheme, which is simple and efficient in practical applications. In addition, the new method is applied to the governing equations under the curvilinear coordinate system and its freestream and vortex preservation properties are verified. Our numerical experiments show that the improved method has stronger robustness and its dissipation at nonlinear waves 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 conservation laws are reviewed. Section 3 presents the improved HWCNS3. A 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 Eq. (1) is first transformed into a generalized coordinates system (ξ , η, ζ ) with the following relationships The transformed equation can be written as

Original HWCNS
The semi-discrete approximation of the governing Eq.
(2) at a grid point indexed as (i, j, k) is as follows: where F ξ i,j,k , G η i,j,k and H ζ i,j,k 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 [11]: whereF i± 1 2 ,j,k are the numerical fluxes at half points, the local Lax-Friedrichs flux is used here:F where α is the maximum eigenvalue of the Jacobian matrix, the cell-boundary recon-structionsŨ L i± 1 2 ,j,k andŨ R i± 1 2 ,j,k are computed by the third-order nonlinear interpolation in the characteristic space. Here, only the construction ofŨ L i+ 1 2 ,j,k is presented, and the subscripts "j, k" are dropped for simplicity. Denote I m and r m the mth left and right eigenvectors 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: and the nonlinear weights ω 1,m and ω 2,m are determined as follows 4 and d 2 = 3 4 are the linear weights, ε = 10 −6 is a small positive number to avoid the denominator becoming 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 results show that the original HWCNS scheme is too dissipative. To overcome this shortcoming, we propose a new nonlinear reconstruction procedure with similar ideas that was first presented in [14].
We still use the construction of U L i+ 1 2 ,m as an example. Two central stencils are chosen: , the corresponding reconstruction polynomials are denoted as q 1,m (x) and q 2,m (x), For consistency, the reconstruction 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 [12]: 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 Eqs. (8) and (9).
The nonlinear weights are defined by where τ m is defined as the absolute difference 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, and Eq. (2) is simplified as where By substituting Eq. (3) into Eq. (17), we can see that which means the uniform flow conditions hold: The condition Eq. (18) is called the freestream condition. To preserve it numerically, we first rewrite the metrics in an equivalent form as suggested in [17]: Then, Eq. (2) at grid node (i, j, k) is discretized as where For l = 1, 2, 3, D ξ l , D η l and D ζ l denote the finite difference operators in ξ , η and ζ coordinate directions respectively.
As proven in [18], the freestream conditions 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 all use the fourth-order compact differential operator Eq. (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 a smaller time step is taken to ensure that spatial errors dominate.

The order of numerical accuracy
Example 1: We consider the one-dimensional Euler equations on [ 0, 2π] with periodic boundary conditions and initial conditions: 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 errors and orders 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 errors and orders 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 is presented in Figs. 1 and 2 respectively. We can observe that the new method has better resolution than the original HWCNS3 in the neighborhood of shocks for these two problems.

Example 4:
We now consider the interaction of two blast waves. The initial conditions are given as 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. Figure 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.
This problem is characterized by rich small scales in the flow fields, and thus is specifically chosen to evaluate the ability of the proposed numerical schemes and resolve the fine structures. We evolve the solution until time t = 0.8. The results with grid points 400 × 400 and 1200 × 1200 are plotted in Figs. 4 and 5. We can see that the new HWCNS3 can capture the rich small scales for this specific example.
Example 6: We consider the double Mach reflection problem. This problem describes the phenomenon that a Mach 10 shock reflected from the wall with an incidence angle of 60 • . The computational domain is chosen to be [ 0, 4] ×[ 0, 1] and the initial conditions are given as (27) γ = 1.4. The left state value and right state value are specified at the left and right boundary respectively. At the bottom boundary y = 0, the reflective boundary conditions are applied when x ≥ 1 6 and the inflow boundary conditions are applied when x < 1 6 . At the top boundary y = 1, the left state value is specified when x < g(t) and the right state value is specified when x ≥ g(t), where g(t) = 1 6 + √ 3 3 (1 + 20t). Using grid points

Example 7:
The freestream is tested on a three-dimensional wavy grid, as shown in Fig. 7. This grid is generated 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.   Table 3, where T 1 denotes the CPU time for computing the grid metrics and Jacobian and T 2 denotes the CPU time for time iteration. We can clearly observe that the L 2 errors are all at the level of round-off errors for different precisions, and verify the expected freestream preservation property accordingly. Since the grid metrics and Jacobian only need to be calculated once and stored in the whole computational process, its CPU time is negligible compared with the CPU time of time iteration. Fig. 7 The three-dimensional wavy grids for the freestream preservation test Example 8: We consider a two-dimensional moving vortex problem on the following wavy grid to test the vortex preservation property.
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/ρ γ , 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.
with θ = 5π 12 , R x = 3, R y = 6, i max = 81 and j max = 61. 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 improved HWCNS3 is more robust than the original one. First, we consider a Mach 2 supersonic flow moving toward the cylinder from left and taking the time-step as t = 0.01. The new scheme works well while the original scheme will diverge unless the time-step is reduced to 0.001. Figure 10 shows the density and pressure  . 9 The computational grids of supersonic flow past a cylinder contours of the two schemes. We can see that the results have good agreement with that of [20]. We then increase the Mach number of the inflow to 10 and take the time-step as t = 0.001. Although both schemes can work in this situation, the density obtained by the original scheme oscillates slightly as shown in Fig. 11. Finally, we increase the Mach number of the supersonic inflow to 20 and take the time-step as t = 0.001, again the new scheme is still robust while the original scheme is divergent under the same condition. Although the original scheme can converge if we further reduce the time step to 0.0001, the density still has small oscillation as shown in Fig. 12.

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