High-Order Hybrid WCNS-CPR Scheme for Shock Capturing of Conservation Laws

By introducing hybrid technique into high-order CPR (correction procedure via reconstruction) scheme, a novel hybrid WCNS-CPR scheme is developed for e ﬃ cient supersonic simulations. Firstly, a shock detector based on nonlinear weights is used to identify grid cells with high gradients or discontinuities throughout the whole ﬂ ow ﬁ eld. Then, WCNS (weighted compact nonlinear scheme) is adopted to capture shocks in these areas, while the smooth area is calculated by CPR. A strategy to treat the interfaces of the two schemes is developed, which maintains high-order accuracy. Convergent order of accuracy and shock-capturing ability are tested in several numerical experiments; the results of which show that this hybrid scheme achieves expected high-order accuracy and high resolution, is robust in shock capturing, and has less computational cost compared to the WCNS.


Introduction
In computational fluid dynamics (CFD), numerical schemes greatly influence the computational accuracy, efficiency, and robustness. Despite the extra complexity of high-order schemes [1], they can provide more accurate results in delicate simulations compared with traditional second-order schemes. With the increasing demand for delicate simulations, such as large-eddy simulation (LES) of turbulent flows, computational aeroacoustics (CAA), and shock-induced separation flows, high-order schemes have attracted remarkable research attentions in recent years [1][2][3][4][5][6].
Among these high-order schemes, the CPR scheme has several good properties. It is compact, efficient, and applicable to complex unstructured meshes [1,6,7]. It was firstly proposed by Huynh [8] and was generalized to unstructured meshes by Wang and Gao [9], which is now widely applied in scientific researches and engineering predictions [6,[10][11][12]. With specific correction functions, the CPR methods are equivalent to specific discontinuous Galerkin (DG) methods, while their calculation cost is relatively lower because they are based on a differential formulation and avoid expensive integration calculations [8].
However, the shock-capturing techniques of high-order finite element methods, including CPR and DG methods, still can not meet the need of many simulations [1,13]. Shu et al. [14] employed a WENO nonlinear limiter on the CPR scheme, which can maintain high-order accuracy in smooth regions and control spurious numerical oscillations near discontinuities. A parameter-free gradient-based limiter for the CPR scheme on mixed unstructured meshes was developed by Lu et al. [15]. By proposing a p-weighted procedure, new smoothness indicators are developed in the DG methods, which obviously improve the shock-capturing ability [16]. However, obvious oscillations can still be observed near discontinuities. Artificial viscosity methods with subcell shock capturing and cross cell smoothness have been proposed for shock capturing [17,18]. However, they are dependent on some user-defined parameters. Great progress has been made; however, the compact polynomial approximations of solutions in a high-order finite element cell are by design not a good approximation of discontinuities, which easily leads to spurious numerical oscillations near discontinuities.
A way around this problem is to develop a hybrid scheme of high-order finite element methods (CPR or DG) and finite difference (FD) or finite volume (FV) methods, which provide better shock-capturing abilities locally near shocks. In this way, high-order finite element methods are adopted in smooth regions, which maintain compactness and high resolution. Meanwhile, FD or FV schemes are adopted to provide robust shock-capturing abilities. Dumbser et al. [19] developed hybrid DG-FV methods and utilized second-order FV to capture shocks. Cheng et al. developed multidomain hybrid RKDG and WENO methods [20,21] with good geometry flexibility and low computational costs. These methods show good shock-capturing ability because they avoid using high-order finite element methods to capture shocks directly, but utilize schemes with better shockcapturing abilities instead.
A hybrid shock-capturing scheme is developed based on the zonal hybrid WCNS-CPR scheme proposed in [7]. In [7], some fundamental problems of the WCNS-CPR schemes, such as spatial accuracy and geometric conservation laws, have been studied, which makes the scheme efficient and applicable to complex curved meshes. In this paper, a shock detector and a nonlinear weighted interpolation are introduced to the hybrid WCNS-CPR scheme for simulating problems with discontinuities. The nonlinear weighted approach of WCNS is introduced to the hybrid scheme for both the shock detection and the shock capturing. Considering the robust shock-capturing ability of WCNS, we intend to develop a hybrid shock-capturing scheme to combine merits of two schemes. Different from the zonal hybrid strategy in [7], the two schemes are switched dynamically based on the smoothness of the local flow field. The main contributions of this paper are as follows: (1) Nonlinear mechanisms, including shock detection and shock capturing, are introduced into the hybrid scheme and a hybrid scheme with good shockcapturing ability is a developed (2) A hybrid WCNS-CPR scheme with dynamical switching strategy between the two schemes is proposed, in which the strategy of dynamical scheme switching and interface treatment methods are addressed (3) Numerical experiments are conducted to show the good properties in accuracy, efficiency, and shockcapturing ability of the hybrid WCNS-CPR scheme.
The rest of this paper is organized as follows. In Section 2, high-order WCNS and CPR schemes are briefly reviewed. The hybrid WCNS-CPR scheme is constructed in Section 3, where the interface treatments are discussed in detail. Section 4 presents several numerical tests to show the performance of the proposed hybrid method. Finally, concluding remarks are given in Section 5.

Brief Review of WCNS and CPR
In this section, a brief review of WCNS and CPR is presented based on two-dimensional hyperbolic conservation laws where U is a conservative variable vector and F and G are the flux vectors.
2.1. Brief Review of Third-Order WCNS. Equation (1) is discretized by a finite difference scheme at every node ðx i , y j Þ where E i,j ′ is the approximation of the spatial derivatives.
The following fourth-order flux difference operators are used where h is the grid spacing. Here, Þ are the Riemann fluxes at cell edges, as shown in Figure 1. Riemann solvers can be used to compute numerical fluxes, such as Lax Friedrichs, Roe, Osher, AUSM, HLL, and their modifications. We refer to papers [3,4,22] and references therein.
To capture discontinuities, the conservative variables at cell edges U R i−ð1/2Þ,j are obtained by the weighted nonlinear interpolation proposed in [23]. The third-order WCNS interpolation has the following steps: Step 1. Divide the third-order three-point stencil Step 2. Construct an interpolation polynomial p ðmÞ ðxÞ in each substencil S ðmÞ i , m = 1, 2, which satisfies p ðmÞ ðx i Þ = U i,j . We have International Journal of Aerospace Engineering Step 3. Calculate the nonlinear weights ω k based on linear weights d m and a smoothness indicator IS k for each substencil S ðmÞ i , m = 1, 2. The linear weights are The improved Z-type nonlinear weights [24,25] are defined by where ε = 10 −6 is a small number, IS k is a smoothness indicator, τ 5 = jIS 0 − IS 1 j, and C q and p z are parameters. We use p z = 2 and C q = 4 to improve the performance of the nonlinear weights. An improved smoothness indicator [25] is employed.
The formulations of U L i−ð1/2Þ,j is symmetric with U R i−ð1/2Þ,j by i.

Review of Third-Order CPR.
For the third-order CPR framework, there are 3 2 solution points ði, j, l, mÞ, l = 1, 2, 3, m = 1, 2, 3, in the ði, jÞ-th cell, as shown in Figure 2. The nodal values of the conservative variable U at solution points are updated by the following differential equation The reconstructed flux F̃i ,j ðξ, ηÞ and G̃i ,j ðξ, ηÞ are expressed as where σ 1 i,j ðξ, ηÞ and σ 2 i,j ðξ, ηÞ are flux correction polynomials Here, correction functions g L ðξÞ, g R ðξÞ, g L ðηÞ, and g R ðηÞ are all cubic polynomials, andF andĜ are Riemann fluxes corresponding to F and G, respectively.
Moreover, the solution polynomial is where L l ðξÞ and L m ðηÞ are the 1D Lagrange polynomials in the ξ and η directions. The flux polynomial is Remark 1. In this paper, solution points are the Gauss points, which are positioned at ξ 1 = − ffiffiffiffiffi 15 p /5, ξ 2 = 0, ξ 3 = ffiffiffiffiffi 15 p /5 for K = 3: Correction function for g DG scheme is used, which is expressed as

Hybrid Schemes Based on WCNS and CPR
3.1. Strategy of Scheme Switching. The hybrid scheme adopts the nonlinear WCNS with shock-capturing ability in nonsmooth regions and the CPR scheme with high computational efficiency in smooth regions. Thus, a shock detector is needed to judge the location of shocks. Recently, many shock detectors with good properties have been developed [26][27][28]. In this paper, we use the NI (nonlinear indicator) [28] to detect which cells probably contain shocks and then mark them as troubled cells. The formula of the NI shock detector is as follows: where r + 1 represents the number of candidate stencils in WCNS with r = 1 for the third-order WCNS. To maintain that discontinuities do not move out of the troubled cells during the time integration, buffer cells are introduced around the troubled cells, which are also calculated using WCNS. As shown in Figure 3, the direct neighbors of troubled cells are defined as buffer cells, and other cells are called smooth cells. Above all, the strategy of scheme switching is as follows. When the NI shock detector detects shocks, in the cells near shocks (namely, in troubled cells and buffer cells), CPR scheme is switched to WCNS. And when shocks leave the cells, these cells are calculated by the CPR scheme again. That is to say, the shock detector judges three types of cells, and we use WCNS in troubled cells and buffer cells and the CPR scheme in other cells.

Interface Treatment.
The interface treatment influences the accuracy and stability of the overall scheme and should be discussed carefully. The main difficulty is at the WCNS-CPR interface. We denote the solution points of the CPR scheme by c and the solution points of WCNS by w.
3.2.1. 1D Interface Treatment. First, WCNS in the i + 1-th cell needs information of ghost points fg i,k , k = 1, 2, 3g from the neighboring i-th CPR cell, as shown in Figure 4. The solution values of ghost points of WCNS cell are calculated via CPR interpolation using data at fc i,k , k = 1, 2, 3g. That is, Second, we apply u −ðCPRÞ i+ð1/2Þ,interface provided by i-th CPR cell and u +ðWCNSÞ i+ð1/2Þ,interface computed in i + 1-th WCNS cell with the values of ghost points to calculate the numerical fluxes at interface: In addition, information exchange of switching strategy at different moments should also be considered. When i-th cell is a smooth cell at the n-th time step and it is judged as a buffer cell or a troubled cell in the next time step, information of solution points of WCNS can be obtained by Lagrange interpolation via information of CPR. Now, we consider the opposite situation: the cell is a buffer cell or a troubled cell at n-th time step and it is switched to a smooth cell at the next time step, information should be transmitted from WCNS to CPR. The procedure involves nonlinear interpolation cross cells. We take the ði, 1Þ-th solution point as an example.
Step 1. Divide the stencil S = fw i−1,3 , w i,1 , w i,2 g into two sub-  i-th cell i+1-th cell  International Journal of Aerospace Engineering Step 2. Construct interpolation polynomials u 1 ðc i,1 Þ and u 2 ð c i,1 Þ in S 1 and S 2 , respectively, we have where l 1 and l 2 are linear Lagrange polynomials.
Step 3. Calculate the nonlinear weights based on linear weights d 1 and d 2 . The linear weights can be obtained by solving the following equations: Then, nonlinear weights can be calculated by equation (6).

2D Interface Treatment.
In this subsection, we consider the interface treatment between ði, jÞ-th CPR cell and ði + 1 , jÞ-th WCNS cell in 2D computational domain in Figure 5(a).
The solution values at ghost points of WCNS in the CPR cell are calculated via CPR interpolation using data at fðc i,j,i 0 ,j 0 Þ , i 0 = 1, 2, 3 ; j 0 = 1, 2, 3g. Using dimension by dimension interpolation strategy, we do the Lagrange interpolation in the x direction first, as shown in Figure 5(b), and then in the y direction, as shown in Figure 5(c). Using solution polynomial equation (15) of CPR, we have at the ghost points fg i,j,i 1 ,j 1 , i 1 = 1, 2, 3 ; j 1 = 1, 2, 3g. Then, we calculate the numerical fluxes at cell interface. As shown in Figure 6,   Therefore, for the CPR scheme, the numerical flux at f p CPR Meanwhile, the numerical flux at f p WCNS k for WCNS iŝ The 2D information exchange of switching strategy at different moments is similar with the 1D situation. The Nonlinear interpolation is used from information in WCNS to CPR, while the linear interpolation from CPR to WCNS. They are not presented here for clarity.

Numerical Investigation
In this section, we will test spatial accuracy, shockcapturing performance, and computational efficiency of the hybrid WCNS-CPR scheme. Comparisons with the WCNS are also investigated. We firstly take 1D linear wave equation to compare numerical errors of single and the hybrid scheme. 1D shock-tube problems are solved to demonstrate the shock-capturing ability of the hybrid scheme. Then, 2D isentropic vortex problem in Euler equation systems is solved for testing spatial accuracy and computational efficiency of WCNS and CPR. At last, 2D Riemann problem and double Mach reflection problem are simulated to test the shock-capturing performance of the hybrid scheme and comparison of computational efficiency is conducted in the double Mach reflection problem. In this paper, Riemann fluxes are computed by the Lax Friedrichs solver. The explicit third-order TVD Runge-Kutta scheme is used for time integration.
The L ∞ error and L 2 error are computed to measure local solution quality and global solution quality, respectively. The L 2 error is estimated by where U h i and Uðr i Þ are the numerical solution and the exact solution at the i-th solution point. Here, N is the total number of solution points.

1D Linear Wave Equation.
Consider the 1D linear wave equation with the initial condition and periodic boundary conditions. L = 6 and time step ΔT = 0:0001 is used to solve the problem till T = 3:0 by different schemes. Firstly, we compare numerical errors of WCNS and CPR on the uniform grid. It can be seen from Table 1 that both WCNS and CPR achieve desired third-order accuracy. In addition, CPR has smaller error than WCNS under the same degrees of freedom (DOFs).
Secondly, the performance of the hybrid WCNS-CPR scheme is studied. To study the convergence accuracy of the hybrid scheme, a very small δ NI = 2:0 × 10 −8 is used to avoid that purely CPR scheme is used in this smooth problem. As shown in Figures 7(c)-7(e), the gray line with light blue dots represents values of NI at each solution point. The value of cell indicator (in dark blue line) represents the scheme used in every cell. When the indicator is 1, WCNS is employed. Indicator of value 0 means the cell is smooth, and it is calculated by CPR. As the number of DOFs gets larger, the percentage of WCNS cells decreases. When the number of DOFs becomes 480, the solution is smooth enough and CPR is used in all the cells. From numerical errors in Table 1, we can see that the hybrid scheme achieves expected third-order accuracy.
where ρ is the density, v is the velocity, E is the total energy, and p is the pressure, which is related to the total energy by E = ðp/ðγ − 1ÞÞ + ð1/2Þρv 2 with γ = 1:4. Three kinds of shock-tube problems are considered in this subsection.

Shock-Tube Problem of Sod.
Consider the Sod problem with initial conditions [29] ρ, u, p ð Þ= 0:125,0, 0:1 ð Þ , −5 ≤ x < 0, and Dirichlet boundary conditions. Time step ΔT = 0:0001, DOFs = 210, and δ NI = 0:2 are used to solve the problem till T = 3:0 by the hybrid WCNS-CPR scheme. The reference solution is obtained by the WCNS3 scheme with DOFs = 2100. As shown in Figure 8, shocks and contact discontinuities can be well captured. Near the contact discontinuity (x ≈ 3) and the shock (x ≈ 5), no obvious oscillation is observed, which shows the good shock-capturing ability of the hybrid scheme. In addition, the CPR scheme is adopted in the main computational domain. The use of CPR scheme helps save computational cost since the linear CPR scheme is more efficient than the nonlinear WCNS scheme as will be shown in Section 4.3.1.    The reference solution is obtained by the WCNS3 scheme with DOFs = 9000. As shown in Figure 9, CPR is employed in the smooth area. Moreover, the hybrid scheme can capture shocks and discontinuities well. The numerical result of density obtained by the hybrid scheme is in good agreement to the reference, which illustrates the high-resolution property of the hybrid scheme.
For an ideal gas, γ = 1:4. Here, ρ,u,v,p, and E are the density, the x direction velocity, the y direction velocity, the pressure, and the total energy, respectively.

2D Isentropic Vortex Problem.
In this subsection, the hybrid WCNS-CPR scheme is adopted to solve the isentropic vortex problem [32]. The initial condition is a mean flow, fρ, u, v, pg = f1, 1, 1, 1g, with the isotropic vortex perturbations added to the mean flow. There are

Distribution of WCNS cells and CPR cells
where r = ffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 + y 2 p and the vortex strength ε = 5. The computational domain is ½−10, 10 × ½−10, 10 with periodic boundary conditions. The problem is solved till T = 1 with time step ΔT = 0:0001 for accuracy test. The parameter is δ NI = 0:001 in this test. Figure 11 shows distribution of different cells, NI in the x direction (NI in the y direction is central symmetric with it), and error distribution. When DOFs are 120 and 240, the central parts of the computational domain (which are in red) are calculated by WCNS. However, when DOFs = 480, all cells are calculated by CPR. This is because that, as the number of DOFs gets larger, NI gets smaller in this smooth problem. When parameter δ NI is lower than the smallest NI in the computational domain, the whole domain is judged as smooth area and only CPR is employed.
From Table 2, we can see that the numerical errors of CPR alone are smaller than those of WCNS under the same DOFs. In addition, the WCNS scheme achieves third-order accuracy, while the order of accuracy of CPR is about 2.5. The hybrid WCNS-CPR scheme achieves third-order accuracy.
From Table 3, we can see that the computational cost of WCNS with characteristic projection is more than twice than that of CPR under the same DOFs. Considering that the errors of CPR are smaller than those of WCNS as shown in Table 2, the efficiency of CPR is even higher than that of WCNS under the same error.    Figure 12, the results of the hybrid scheme in Figure 13(a) capture the flow field accurately with no obvious oscillation. In addition, small-scale structures are better captured by the hybrid scheme. We can see in Figure 13(b) that most of the computational domain is calculated by CPR, which contributes to the high efficiency of the hybrid scheme.

Double Mach Reflection
Problem. This test is first described by Woodward and Colella [33]. It is aimed at testing the robustness of the high-resolution schemes. A Mach 10 oblique shock is initially set up at x = 1/6 on the lower boundary, and its direction is 60°from the x-axis. Inflow and slip wall boundary conditions are specified at the bottom boundary for x = ½0, 1/6 and x > 1/6, respectively. For the upper boundary (y = 1), a time-dependent boundary based on the analytical propagation speed of the oblique shock is imposed.
The simulation is run until T = 0:2 using time step ΔT = 0:0001.   Figure 14. The remaining area is accurately calculated by the CPR method as well. Moreover, the total number of cells using CPR is much more than that using WCNS. As a result, the computational efficiency of the hybrid scheme is higher than that of WCNS as shown in Table 4. For DOFs = 600 × 150 and δ NI = 0:9, the hybrid scheme is slightly more efficient than WCNS. However, for DOFs = 1200 × 300 and δ NI = 0:9, the hybrid time is about 1.6 times more efficient than the WCNS scheme.

Concluding Remarks
In this paper, a new hybrid WCNS-CPR scheme for simulating conservation laws with discontinuities is proposed by introducing nonlinear weighted approach of WCNS to the linear hybrid scheme. Shock detector based on nonlinear weights is used to detect nonsmooth troubled cells and buffer cells. In these cells, WCNS is employed in order to capture shocks, while in smooth cells CPR is used to maintain high computational efficiency. By introducing buffer cells, the scheme interfaces are by design placed in relatively smooth regions and a linear interpolation method is used for the CPR to provide information for the WCNS scheme. Accuracy tests, shock-capturing tests, and computational efficiency tests are conducted. Accuracy tests in 1D and 2D show that the hybrid scheme achieves designed high-order accuracy. Shock-capturing tests illustrate that the hybrid scheme captures discontinuities without obvious oscillations. Computational efficiency tests show that the CPU cost of the hybrid WCNS-CPR scheme is relatively lower than that of WCNS.
In the future, the hybrid scheme will be generalized to three-dimensional cases and to Navier-Stokes equations.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.

12
International Journal of Aerospace Engineering